By coupling statistics and heat transfer, we investigate numerically laser-induced crystal damage by multi-gigawatt nanosecond pulses. Our model is based on the heating of nanometric absorbing defects that may cooperate when sufficiently aggregated. In that configuration, they induce locally a strong increase of temperature that may lead to a subsequent damage. This approach allows to predict cluster size distribution and damage probabilities as a function of the laser fluence. By studying the influence of the pulse duration onto the laser-induced damage threshold, we have established scaling laws that link the critical laser fluence to its pulse duration τ. In particular, this approach provides an explanation to the deviation from the standard τ 1/2 scaling law that has been recently observed in laser-induced damage experiments with KH2PO4 (KDP) crystals [J.J. Adams et al, Proc. of SPIE 5991, 5991R-1 (2005)]. In the present paper, despite the 3D problem is tackled, we focus our attention on a 1D modeling of thermal diffusion that is shown to provide more reliable predictions than the 3D one. These results indicate that absorbers involved in KDP damage may be associated with a collection of planar defects. First general comparisons with some experimental facts have been performed.
© 2007 Optical Society of America
For more than a decade, large aperture laser facilities like NIF in USA or LMJ in France are being developed in order to produce inertial confinement fusion . The final part of the laser chain consists of frequency conversion that is achieved by potassium dihydrogen phosphate (KDP) crystals. Nanosecond multi-gigawatt laser pulses induce localized and probabilistic bulk damages in these crystals whose optical properties are subsequently altered. In order to control the damage creation or to modify the crystalline growth process, there is a need for understanding mechanisms involved in damage processes. Now, it is currently admitted that damage results from a local increase of temperature that can be as high as 10000K . It gives rise to a subsequent shock wave which propagates over a few micrometers then leading to the characteristic size of laser-induced defects that are observed experimentally by optical techniques [3, 4]. The local increase of temperature is mainly associated with defects that enhances the laser absorption. For instance, charged hydrogen or oxygen vacancies, for which light absorption takes place on a distance of a few nanometers, seem to be good candidates as pointed out by Liu et al [5, 6]. It is reasonable to think that such a nanometric defect alone is not able to produce the critical temperature and thus, the fact that observed damages result from a local defect aggregation seems to be a realistic assumption [7, 8, 9].
In order to characterize experimentally the resistance of these crystals to optical damaging, a standard measurement consists in plotting the bulk damage probability as a function of the laser fluence F  that gives rise to the so-called S-curves [10, 11]. In order to explain this behavior, thermal models based on an inclusion heating have been proposed [12, 13, 14]. In these approaches, statistics (Poisson law) and inclusion size distributions are assumed. On the other hand, pure statistical approaches mainly devoted to the onset determination and that do not take into account thermal processes have been considered [15, 16, 17, 18, 19]. On the basis of the above-mentioned assumption of defects aggregation, we propose a model where Absorbing Defects of Nanometric Size, hereafter referred to as ADNS, are distributed randomly and may cooperate to the temperature rise ΔT through heat transfer within a given micrometric region of the bulk that corresponds to an heterogeneity. Since this approach combines statistics and heat transfer, it allows to provide the cluster size distribution, damage probability as a function of fluence, and scaling laws without any supplementary hypothesis. The present paper aims at introducing the general principle of this model and giving first main results that are compared with experimental facts. A particular attention has been payed to scaling laws since they are very instructive in terms of physical mechanisms. More precisely, a deviation from the standard τ 1/2 law has recently been observed within KDP crystals  and here we propose a plausible explanation of this fact based on thermal cooperation effects and statistics. Despite the 3D representation is tackled, we focus on a one dimensional modeling that gives a good insight about physics and seems to provide a nice counterpart to experimental tendencies.
The paper is organized as follows : in order to support discussions about the statistical modeling, we consider in Section 2 the one cluster case, i.e. a well-defined crystal region filled with ADNS that do not feel any thermal contribution of other places, that may be seen as a single inclusion absorbing the laser energy. Also, in that section, different cluster geometries and their associated temperature evolution with respect to cluster size and time are considered. Section 3 deals with the model coupling statistics and heat transfer. In a first part, the principle of the approach is exposed. Secondly, cluster size distributions obtained from numerical simulation based on statistics are analyzed. Thirdly, we provide numerical predictions of the model in terms of damage probability S-curves and temporal scaling laws. Results are discussed and compared with experimental facts in Section 4. Conclusions are drawn in Section 5. For the reader convenience, details of analytical derivations have been reported in Appendix.
2. One cluster case
We address hereafter considerations about the laser heating of only one cluster that may be seen as a single inclusion in the bulk. Thermal properties of this cluster are assumed to be identical to its surrounding. First, let us dealing with the temperature of ball-like absorbing objects. Since we are interested in the highest temperature, we consider it at the center of the cluster. In Appendix A.1, for the 1D case, it is shown that the temperature rise ΔT at position x = 0 and time t reads:
where erfc is the complementary error function, D the termal diffusivity of bulk, l the cluster radius, ρ the cluster density, C the cluster specific heat capacity and A the absorbed power per unit of volume (assumed to be constant in the course of interaction). In the 3D case, the temperature rise can be written as:
By using the following parameters τ = 1ns, A/ρC = 1013 K.s -1 and DKDP = 6.5×10-7 m 2.s -1, the evolution of these temperatures is plotted as a function of the cluster size l on Fig. 1. In the region where the saturation has not been reached, whatever the number of dimensions allowed for spatial diffusion, it appears that the bigger the cluster, the higher the temperature reached. Further, the influence of the number of spatial dimensions lies in the curve slope. Actually, ΔT 1D goes up faster than ΔT 3D. It is noteworthy that the amount of matter needed to form the cluster increases relatively faster in the 3D case (as l 3) than in the 1D one (as l).
Since the present paper focuses notably on temporal scaling law, it is of interest to study how the temperature evolves as a function of time and cluster geometry. From a general point of view, in condition of linear absorption, the temperature may be written as T(t) = ζIf(t) = ζFf(t)/τ where I is the laser intensity, and where ζ and f(t) are a constant and a function of time respectively that depends on the cluster thermal properties, geometry and size. For a given cluster, assuming that the critical temperature is reached at the end of the laser pulse, the scaling law simply reads:
while the analytical expression of f(τ) remains unchanged, i.e. for a particular diffusion regime. The function f(τ) now has to be evaluated. To do so and in order to keep the analogy with a cluster composed of a multitude of tiny absorbing area, we use the Green’s formalism that consists in summing point source contributions over the appropriate cluster geometry. For instance, the temperature expression has been derived for a sphere, a flat cylinder and the surface of a cone in Appendix B. These different temperature evolutions as a function of time are reported on Fig. 2 in the case where the characteristic cluster size l has been set to 25 nm and with A/ρC = 1013 K.s -1. As it is commonly known, these curves exhibit two kinds of behavior corresponding to the limit cases l ≫ √Dt and l ≪ √Dt, the transition corresponding to time close to tT ≃ l 2/D ≃ 1 ns. Moreover, the temperature rise depends on the cluster geometry. These temporal behaviors may be deduced from a simple following picture in the case l ≫ √Dt: when the time elapsed after the beginning of the laser illumination is of the order of r 2/D, a point source situated at r starts to contribute efficiently to the temperature rise at r = 0 (where it is assumed to be maximum for a spherical cluster). After a few r 2/D of time, this contribution becomes almost stationary since the temperature rise θ 3D of a point source evolves asymptotically as (1 - r/√πDt)/r ≃ 1/r as it is shown in Appendix A.2. Despite the contribution of one source rapidly does not depend on time anymore, the temperature at r = 0 always evolves thanks to contributions coming from sources that are more and more distant when time elapses. Actually, the temperature shape evolves accordingly to the spatial number of point sources between r = √Dt and r + δr. The temperature variation δT induced by one source per interval of time δt = 2rδr/D is roughly given by δT/δt ≃ (1/r) × (D/2rδr). Let us introduce N(r) the number of point sources between r and r + δr. The variation of temperature induced by sources located between r and r + δr then is DN(r)/2r 2 δr. For a spherical cluster whose point source density is constant, N(r) = 4πr 2 δr implying δT/δt = cst, i.e. T ∝ t and thus the temperature increase linearly with respect to time. If we consider a non-constant point source density, conclusions differ. In the case of a spherical cluster whose density goes down as 1/r with respect to the radial distance, N(r) becomes N(r) = 4πr 2 δr × (1/r), that leads to a temperature evolving as √t. Now, in the case of a flat cylinder or a cone surface (without its base), N(r) ∝ rδr thus leading to δT/δt ∝ 1/r or equivalently δT/δt ∝ 1/√t. The integration over time gives T ∝ √t. Also, this simple analysis shows clearly the role of the cluster geometry and density onto temperature rise dynamics when l ≫ √Dt. In case where l ≪ √Dt, as soon as the most distant source has contributed (for time longer than a few l 2/D) to the temperature rise into the cluster, the temperature then remains constant whatever the 3D cluster geometry. For the 1D diffusion, since a point source contribution to the temperature rise is never constant with respect to time but varies as √t (see Appendix A.2), the cluster temperature is never time-independent even if l ≪ √Dτ.
From the above discussion, using Eq. (3), we are able to establish scaling laws according to the considered cluster geometry. In case where l ≫ √Dτ, for a density-constant spherical cluster, one has Fc ∝ τ 0. It means that the relevant parameter is the absorbed energy whatever the time needed to do so. In case of a non constant density that follows 1/r, the scaling law becomes Fc ∝ τ 1/2. The scaling law also reads Fc ∝ τ 1/2 for a flat cylinder or a cone where laser energy is absorbed onto its surface. In case where l ≪ √Dτ, all 3D absorbing objects obey to Fc ∝ τ.
3. Random distribution of absorbing defects
Let us consider a set of ADNS that are distributed on a spatial domain. When a crystal cell contains a defect, it is called an alterate cell which may absorb laser energy much more efficiently than a pure cell crystal. Therefore, from heat transfer point of view, the alterate cell may be seen as a very tiny source inducing temperature rising. This source of size a is assumed to be close to the characteristic crystal cell dimension, i.e. one nanometer. Within a 1D framework, the domain of size Na, that is assumed to correspond to an heterogeneity, then is composed of two kinds of cells. The temperature evolution is governed by the Fourier equation:
where xi ∈ [0,Na] stands for the position of ADNS or alterate cells whose number is nADNS and A is the absorbed power per unit of volume. Material physical parameters such as thermal diffusivity D and conductivity λ, density ρ or specific heat capacity C, linked by the relation D = λ/ρC, are assumed to remain constant in the course of interaction. The function ∏ is defined as follows:
The way to distribute defects is addressed in the next paragraph. A general solution of (4) is given by :
where T 0 = 300K is the initial temperature of the crystal and where the temperature rise ΔT (i)(x,t) induced by one ADNS is solution of the following equation:
Now, since we work in conditions where a ≪ √Dt (for DKDP = 6.5 × 10-7 m 2.s -1 and τ = 1ns, we have √Dt ≃ 25nm), in order to deal with simple formula allowing fast numerical calculations, the function ∏(x) is replaced by the Dirac delta function δ(x), i.e. we impose that the energy absorbed by a finite-size defect is in fact absorbed by a point source. The ADNS then may be seen as a heating point source. Within this framework, a good approximated solution of Eq. (7) is given by (see Appendix A.2):
We have checked the reliability of this approximation performing a full numerical resolution of Eq. (4) based on a finite differences scheme. The procedure in a 3D space is very similar. The only difference lies in the expression of the temperature rise induced by a point source. It is given by (see Appendix A.2):
Since θ 3D(r,t) tends to infinity when r goes to zero, we take θ 3D(r ≤ a,t) = θ 3D(a,t). In order to illustrate the main principle of the model, we plot on Fig. 3 the evolution of the temperature (6) as a function of the 1D spatial coordinate in a case for which nADNS = 15, A/ρC = 1013 K.s -1 and τ = 1 ns. This graph shows a characteristic spatial behavior where it appears clearly that cooperative effect leads to a locally higher temperature.
Now, from these calculations, we are able to construct a damage probability law. For a given fluence F, a number ndraw of ADNS distribution are generated, and we check whether or not each ADNS distribution induces a temperature higher than the critical temperature Tc. Let ndam be the number of ADNS configurations leading to T ≥ Tc. The damage probability then is simply given by P(F) = ndam(F)/ndraw. In order to plot the damage probability as a function of fluence, we have to evaluate the source term of Eq. (7). Since the defect nature is unknown, the absorbed power per unit of volume A cannot be evaluated theoretically thus leading us to an empirical evaluation from experimental data. From dimensional analysis, we write A = ξF/lτ where ξ and l may be identified to a dimensionless absorption efficiency and an effective cluster size respectively. In conditions where thermal diffusion is not taken into account, the temperature rise induced by the laser pulse reads ΔT = ΔρE/ρC where ΔρE = ξF/l and C = 0.023 × 900 ≃ 2 J.K -1.cm -3 are the absorbed energy per unit of volume and the heat capacity per unit volume of KDP respectively. From black body measurements , it turns out that the temperature rise associated with damage induced by a 3ns laser pulse with F ≃ 10J.cm -2 is roughly 10000K. That allows to estimate the unknown ratio ξ/l to about 2 × 103 cm -1. In fact, thermal diffusion processes take place, and ξ/l must be higher. Therefore, we choose ξ/l = 104 cm -1 and we will check in Section 3 that it is consistent with experimental data. Since we are working within the Rayleigh regime for which ξ ∝ l (defects clusters are assumed to be well smaller than λ = 351nm), we assume ξ/l to remain almost unchanged when l varies. Also, we use A = 104 × F/τ as source term empirical expression of Eq. (7) in numerical calculations.
In following calculations, we choose Tc = 5000 K  and we define the Laser-Induced Damage Threshold (LIDT) as the value of the fluence Fc such that the damage probability equals 10%. As usually in Physics, this choice indicates actually the departure from perturbative conditions. It is worth noting that taking 5% or 20% will not affect the main conclusions of the present paper.
3.2. Cluster size distribution
Since we associate laser damage with clustering of ADNS, it is of interest to consider the size distribution of clusters. For the sake of simplicity, we only deal with the 1D case. We define a cluster as an ADNS set that may cooperate in order to produce a temperature rise that cannot be reached by a single absorbing defect. Assuming the reasonable criterion for which a defect contributes to the temperature rise up to a distance of about 2√Dτ (10% of the maximum of the temperature at x = xi with τ = 1ns), in this approach, the maximum cluster size is roughly 4√Dτ≃ 100 nm. Since we wish to deal with clusters where important cooperation effects take place, they should be composed of a significant number of ADNS, say at least a ten. Also we introduce the distance d between two adjacent defects below which our criteria are satisfied. For a standard case, we then choose d = 10 nm. In order to define numerically a cluster, we proceed as follows: the domain is scanned from one border to the other one and when an ADNS is detected, the associated cluster size is evaluated. To do so, moving forward, while the distance between two consecutive defects is smaller than d, the cluster survives. The break of this condition then defines the cluster size. Counting the number of clusters of a given size in a domain thus allows to construct the cluster size distribution.
Figure 4 shows the cluster size distribution in case where ADNS are distributed randomly. Calculations have been performed with n ADNS = 50,100,200 and 400 in a domain of size N = 10000 (i.e. of characteristic length 10 μm ) that implies nADNS ≪ N. For each defects density, 10000 drawings have been done. Sub-figure of Fig. 4 displays the average (on drawings) local ADNS density of clusters as a function of their size for n ADNS = 400. Whatever the value of n ADNS, Figure 4 reveals a particular behavior of the cluster size distribution: curves are composed of a plateau (up to a cluster size of 12nm) followed by a decrease. It is worth noting that fluctuations are due to statistical uncertainties. The plateau corresponds to a region where no much more than two defects compose the cluster. Indeed, as long as the cluster is smaller than d+2 = 12nm, the cluster density can be very well approximated by the function 2/l where l is the cluster size and the higher n ADNS, the better the approximation. Since defects are distributed randomly, in conditions where nADNS << N, the probability of having two defects separated by a given distance Δx is the same whatever the Δx value as long as Δx ≪ Na. That explains the existence of the plateau for l ≤ d + 2 and the subsequent break for l > d + 2 where at least 3 defects are needed to make a cluster. In this region, the ADNS number is roughly proportional to the cluster size as the standard thinking. Sub-figure of Fig. 4 strengthens actually this fact since ADNS density variations in clusters are very smaller than in the plateau region. After the plateau, it is clear that the probability of having a given cluster size decreases with the number of involved ADNS. Therefore, the number of cluster falls down when their size increases. In this region, the distribution behavior may be fitted by exp(-αl) where α is a coefficient depending on the defects density. Further, the higher n ADNS, the better the fit.
Now, on the basis of energy minimization considerations, one can expect the ADNS not to be distributed randomly [9, 21]. Indeed, during the crystalline growth or after, ADNS may migrate until reaching the lowest energy level that corresponds to local periodic structure or particular molecular combinations. This behavior may have dramatic consequence on the damage process since we have shown the bigger the cluster, the higher the temperature reached (see Fig. 1). Therefore, it is important to know whether the cluster size distribution is affected by an ADNS aggregation assisted factor. In order to produce this distorted distribution, when each ADNS is arranged successively on the domain, we impose a larger probability (Pd) in a place close to a defect than in a defectless place (probability Pdl). In the previous case for which equiprobability was satisfied, one has Pd/Pdl = nd/ndl where nd and ndl stands for the number of cells close to alterate cells and those that are not respectively. For a biased statistics, we assume Pd/Pdl = bnd/ndl where b is the bias. Since probabilities must satisfy the consistency relation Pd+Pdl = 1, simple algebra leads to the following expressions:
Since we have no information about the defect nature which allows to perform a physical construction of the distorted distribution, we adopt that simple mathematical method to build the distribution. Cluster size distributions are given on Fig. 5 for values of b ranging from 1 to 106, with N = 10000, n ADNS = 100 (i.e. 1% of damaged cells) and for d = 0 (in this case, we only consider adjacent defects). Note that even if the highest value of b seems not to be realistic, it is of interest to observe its influence onto the cluster size distribution. Up to b = 10000, distributions are standard in sense where the smaller the cluster, the higher the number of clusters of this size. The influence of b is to shift the distribution towards larger values of the cluster size. Since the overall ADNS density is constant, the peak distribution decreases. When b becomes greater than 10000, the distribution tends to be uniform except at l = 100nm where all ADNS form a single cluster. One would expect a dirac-like distribution located at l = 100nm when b goes to infinity. Also, we have observed a strong influence of the bias onto the cluster distribution and its influence onto the damage probability curves is addressed hereafter.
3.3. Influences of laser pulse duration and defects arrangement on the damage probability
Let us begin with the 1D case. The damage probability as a function of the laser fluence is shown on Fig. 6 for τ = 250ps, 1ns, 4ns and 16ns. In all cases, we choose nADNS = 100, N = 10000 and statistics is supported by 200 drawings. The general shape of curves exhibits a standard behavior for which damage probabilities increase monotically with the fluence. Further, calculations for τ = 3 ns show that the LIDT is close to 9 J.cm -2 that is in a good agreement with Carr et al’s experiment . Also, a posteriori, this result confirms the reliability of our evaluation of the coefficient ξ/l that has been set to 104 cm -1. For a given pulse duration, probability becomes non-zero when, within at least one drawing, a group of ADNS is sufficiently aggregated to form a cluster whose size and density are enough to reach locally the critical temperature. When the fluence goes up, the critical temperature may be reached by smaller or less dense cluster (see Section 2). Since they are in a larger number, the probability increases itself. A probability close to one corresponds to the lowest cooperative effect, i.e. involving the smallest number of ADNS around the place where T ≥ Tc. This number is determined counting ADNS that contribute significantly to the place xc where T ≥ Tc. To do so, we make the counting in the range [xc - 4√Dτ, xc + 4√Dτ]. The evolution of this number as a function of the laser fluence is shown on Fig. 7 and confirms the latter explanation: for a given pulse duration, the larger the fluence, the lesser the number of ADNS involved in the damage.
Now, we focus on the influence of the pulse duration on the damage probability curves and, more precisely, the scaling law linking the fluence to the pulse duration. Figure 6 shows clearly that the LIDT is shifted towards higher fluence when τ increases. Actually, for long pulse, the thermal diffusion process is more efficient and more energy is needed to reach the critical temperature for a given ADNS configuration. For instance, as shown in Section 2 for 1D diffusion in case where i ≪ √Dτ, the temperature rises as √τ  thus implying the scaling law F c1/F c2 = √τ 1/τ 2. In order to establish the scaling law in our model, we assume that we may write F c1/F c2 = (τ 1/τ 2)x. The exponent x is determined from the knowledge of F c1, F c2, τ 1 and τ 2. In fact, we choose τ 2 = 4τ 1 and we evaluate the x exponent as a function of τ 2 ranging from 1 ns to 16 ns. The values of Fc and x are displayed for the above-mentioned values of τ in Table 1. Note that more drawings than in Fig. 6 have been performed in order to determine Fc with a correct numerical precision. It appears that x differs from the expected 1/2 value and takes values that depends on the pulse duration. For the sake of completeness, the values of x for a large range of pulses ratios is plotted on sub-figure of Fig. 6. The evolution of x can be fitted by the empirical following expression:
with α = 2.92 × 10-2 and β = 0.3592. The fit error is close to 2% and may be due to statistical uncertainties. We have checked that a different ratio τ 2/τ 1 leads to the same general expression of x but where the coefficients α and β differ slightly. For a given ADNS spatial configuration, this behavior of x can be understood by the following consideration: for a single ADNS, since ΔT ∝ τ 1/2, the according scaling law reads Fc ∝ τ 1/2. For a given pulse duration, since several ADNS are involved, whose effective cluster size is non negligible in comparison with the diffusion length, the temperature goes up faster than τ 1/2 and the scaling law exponent is accordingly lower than 1/2. Now, when τ increases, since the contribution length is proportional to √Dτ, more and more defects take part in the temperature rise. It is like if the laser pulses sees different target sizes with respect to its duration. It then turns out that less energy is needed than in a situation where all target components always contribute for any interaction time. Consequently, the scaling law deviates from the standard x = 1/2 value as the pulse duration goes up. This argument is also supported by results of Fig. 7 that shows an increasing number of involved ADNS as the pulse duration becomes longer.
In order to compare with experiments on KDP for which the scaling law exponent is a constant in a narrow range of pulse durations, and thus where, from a theoretical point of view, x does not change significantly, we have evaluated Fc for τ = 1,2,3,…, 10 ns. A power law fit of the evolution of Fc as a function of τ expressed in nanoseconds then leads to Fc = 6.34 τ 0.31 with an error of less than 1%.
The influence of the ADNS density is shown on Fig. 8(a) where n ADNS is set to the following values: 50, 100, 200 and 400. The laser pulse duration is fixed to 1ns and we always work with N = 10000 and ndraw = 200. As it can be expected, the LIDT shift towards lower values when the defect density increases. Indeed, the probability of having larger cluster increases itself thus needing less laser energy to reach the critical temperature. Again, n ADNS = 100 should be retained in order to stick to experimental data .
The same kind of influence is obtained when one considers a distorted distribution as defined in the latter paragraph. Figure 8(b) shows that the higher the bias, the lower the LIDT. Actually, one increases the probability of having large clusters as b increases itself.
Now, let us focus on the 3D case. Figure 9 shows the evolution of the damage probability as a function of fluence for τ ranging from 125 ps to 4 ns. The domain size is 50 × 50 × 50 nm 3 where 1250 ADNS have been distributed (corresponding to 1% of the domain in terms of absorbing volume) and 200 drawings have been performed for each fluence. In order to deal with realistic fluence of a few Joules per square centimeter, the coefficient ξ/l has been set to 106 cm -1. Indeed, a 3D source saturates rapidly with respect to time whereas a 1D source always increases its temperature as shown by expressions (9) and (8) respectively. As regards damage probability curves, Fig. 9 indicates similar main behavior in comparison with 1D results. Nevertheless, two differences appear. The first one lies in the variation of the probability. Within the 1D modeling, in the case of a 1 ns laser pulse for example, an increase of the fluence of 4.5 J.cm -2 is needed to go from a probability of 10-2 to 0.99 resulting in an average slope of almost 0.2 cm 2.J -1. In the 3D case, the average slope is roughly 4 cm 2.J -1 thus 20 times larger than the 1D case. This fact may be explained by simple probabilistic considerations. Despite equivalent ADNS density, the probability of making a cluster in a 3D space is very smaller than in a 1D space. A corollary is that the cluster size distribution is sharper in 3D. Therefore, for a finite number of drawings, after the required fluence for a non-zero probability has been reached, a small amount of additional energy allows to heat a large number of cluster that leads rapidly to a high damage probability. It is worth noting that slopes depend weakly on τ in the pulse range of interest but that does not affect the latter conclusion. The second difference lies in the temporal scaling law. As shown by Fig. 9(f), the value of the exponent x goes up with respect to the pulse duration whereas the opposite behavior has been observed in the 1D case. For a single ADNS, since ΔT does not depend on time because l ≪ √Dτ (see Section 2), the according scaling law reads Fc ∝ τ. In the statistical model framework for the shortest interaction times, since several ADNS are involved, a significant cluster size implies that the temperature evolves with respect to time and the scaling law exponent is well below 1. Now, when τ increases, unlike the 1D case, the contribution of others ADNS is negligible because of the above-mentioned probabilistic considerations and the rapid decrease of θ as a function of time in the case r = √Dτ (θ 3D (r = √Dτ, τ) ∝ Fτ -3/2 whereas θ 1D (x = √Dτ, τ) ∝ Fτ -1/2, see Appendix). Also, in the 3D case, a cluster can really be seen as a well defined inclusion whose temperature saturates when l becomes shorter than √Dτ. It follows that the longer the pulse duration, the larger the scaling law exponent. This is enhanced by the finite domain size. Indeed, for τ≥ 1 ns, the according diffusion length is longer than 25 nm that is the same order of magnitude of the domain size. Therefore, for times longer than 1 ns, no more ADNS contribute to the temperature rise of the cluster.
Finally, numerical studies reveal that the influence of the total ADNS density and of a biased ADNS distribution is similar to the 1D case : for a given domain size, the higher nADNS or b, the lower the LIDT.
In order to analyze our results and to compare them with experimental data, we have to associate the different studied systems (geometry and dimension) with physical objects of the 3D space. As regards the one cluster case, a sphere is the simple model for representing an inclusion embedded in the bulk material. It is a reasonable assumption of considering aggregation of impurities in such a way . Furthermore, they can be produced artificially as gold in glass . Concerning flat cylinders, they may be associated with planar defects as stacking faults, grain boundaries or an array of dislocations [24, 25, 26]. As regards 1D models, i.e. thermal diffusion in only one direction, they may be identified to a slab whose transverse size (with respect to the thermal flow direction) are larger than the maximum value between the parallel size and a few diffusion lengths. More precisely, by using the Green’s functions formalism (see Appendix), we can show that the temperature induced by a disk of radius a is proportional to (1 + (z/a)2)-1/2erf(((1 + (z/a)2)/Dτ)1/2 a/2) at distance z along the revolution axis. For z << a, the temperature is almost constant whereas it decreases as 1/z for z >> a. The transition from a constant to a 1/z behavior indicates the deviation from a 1D thermal diffusion. Regarding the latter mathematical expression, the transition occurs for a value of a close to z. Since cooperative effects take place while the distance between defects is shorter than a few diffusion lengths, the greatest value of z is roughly 100 nm. It follows that the lateral dimension has to be greater than about 100 nm in order to support a 1D thermal diffusion within a nanosecond regime. Experiments in KDP show that set of planes scatter the light . This fact tends to support the planar defect existence. Also, from this point of view, the 1D statistical approach consists in considering the aggregation of planar defects along a particular direction of the 3D space. On the contrary, the 3D statistical approach well corresponds to the intuitive picture for which the ADNS is nothing but a point defect.
A numerical evaluation of the cluster size distribution based on a random repartition of ADNS has shown that the bigger the cluster, the lower the number of it. This result coincides with the common sense as it is used, in works based on thermal transport within a single inclusion embedded in a matrix, to assume analytical expressions containing fit parameters [13, 14, 28]. Also, by using power laws of absorbing inclusion size distribution, these approaches allow to reproduce experimental tendencies . As regards our statistical model, an exponential-like cluster size distribution is predicted. An exponential function being nothing but a power series expansion, it should lead to the similar above-mentioned tendencies. Nevertheless, we emphasize that the distribution obtained from our model may be adapted, introducing a biased statistics, to particular physical conditions for which clustering is made easier.
By coupling statistics with heat transfer, we have been able to predict the damage probability as a function of fluence. In all cases, we predict S-curve shapes that are in good agreement with experimental ones (see for example [3, 10]). Within the one dimensional calculations, the empirical formula A = ξF/lτ with ξ/l = 104 cm -1, that provides the power density absorbed by a single defect, allows to obtain LIDT of a few Joules per square centimeter (for a few ns laser pulse) comparable with experimental data . In the 3D space, reasonable LIDT are obtained setting ξ/l to 106 cm -1, i.e. two orders of magnitude higher than the 1D-model value. This is due to the fact that, firstly, ADNS temperature saturates rapidly with respect to time, secondly, the probability of making cluster is small and, thirdly, the interaction distance is very finite (~ 10nm). Then, from this point of view, it appears that a 1D diffusion mechanism is the best candidate to reproduce experimental data.
From experimental results and for a given laser set up, one observes large variations of the LIDT [11, 30]. This fact can be easily explained in this model framework using different defect density or cluster size distribution, and may be directly linked to crystalline growth. Indeed, one may expect different defects distributions with respect to the quantity of catalysts  or the growth velocity which influences species migrations and subsequently clustering.
Now, average slope of damage probability curves also may inform about physical processes. From experimental data, it appears that the slopes never exceed a few tenth of square centimer per Joule, see e.g. Refs. [10, 13]. For instance, with a 7.5 ns laser pulse, one can state that the average slope is close to 0.05 cm 2.J -1 or so [11, 30]. In our calculations, with an ADNS density of 1%, slopes are roughly 0.1 cm 2.J -1 and 1 cm 2.J -1 for the 1D and 3D cases respectively. Again, it appears clearly that the 1D diffusion mechanism provides the best predictions.
Historically, a large importance has been given to the scaling law Fc ∝ τx with x = 1/2 since it can be established by simple physical considerations about heat transfer [31, 32]. In particular, it is supported by an important experimental work performed by Stuart et al  showing that the scaling law exponent x fits remarkably the 1/2 value for pulse duration longer than a few tens of picoseconds, i.e. when thermal diffusion is allowed. Nevertheless, it is worth noting that surface damages of silica and CaF2 have been considered in this experiment. Dealing with damages in KDP bulk, scaling laws seems not to be establish so clearly and x takes a value close to 0.3 . The 1/2 exponent can be retrieved considering a flat cylinder, a cone surface absorption geometry or a single inclusion of spherical aspect whose density evolves as 1/r. If one considers a set of homogeneous inclusions of different size heated through surface energy absorption, it can be shown that radius of inclusions leading to damage [11, 14] is close to 2√Dτ. Then, it turns out the scaling law reads Fc ∝ τ 1/2. In such approaches, the absorption coefficient may be established using the Mie’s theory that leads it to depend on the inclusion radius as Rm where m ≤ 1. In this case, x becomes (1-m)/2. Depending on the refractive indexes (n,k) values, x can reach the observed experimental value of 0.3 or so . Here we have underlined a new mechanism, based on cooperative effects within a heat transfer framework, that allow to predict non-one-half scaling law exponent. Since the 3D model provide x values rather larger than 1/2, it has to be turned down. As regards the 1D model, it appears for providing x values that are in a good agreement with experimental data when one considers standard values of τ between roughly 1 ns and 10 ns. Finally, even for different defect geometry as for example a non-uniformely ADNS-filled sphere with ρ(r) behaving as 1/r, one may expects that cooperative effects take place between such defects leading to a scaling law exponent below 1/2, as we have shown in the statistical model framework.
We have developed a model coupling statistics and heat transfer within a spatial domain filled partially with absorbing defects. Whereas a single nanometric defect cannot be heated up to the critical temperature leading to damage, when a set of them are allowed to cooperate, the critical temperature can be reached. A statistical spatial distribution of these defects then permits to predict damage probabilities as a function of fluence. We then interpret damage as a result of defect clustering. Since cluster distribution goes down monotically with respect to their size, the higher the damage probability, the smaller the cluster sizes involved in damage. Also, the defect density or their ability to aggregate have been shown to shift the LIDT. By changing the pulse duration τ, we have been able to establish scaling laws that link the critical laser fluence to τ. Results depend on the problem dimensionality and the geometry of clusters that may be formed. By comparing numerical data with experimental tendencies, the 1D model has been shown to provide more reliable results than the 3D. This result indicates that defects may be a collection of planar defects. Within this framework, a temporal scaling law with an exponent value close to 0.3 that deviates from the standard square root law is refound.
A. Analytical resolution of heat equations
We propose hereafter to solve the following set of 1D-equations:
In the Laplace space, ΔT(x,t) transforms into ΔT̄(x,s), and the latter set of equations can be re-written as:
where the upperscript ʺ stands for the second spatial derivative and q = √s/D. Solutions of this set of equations reads:
Since solution has to be symmetric with respect to the operation x → -x, it follows that α = δ and β = γ. The two final constants can be finally determined considering continuity of temperature and flux at x = a. Following this procedure, in the Laplace space, the temperature for x ≥ a is given by:
Using the following result 
the temperature outside the sphere where x ≥ a reads:
The temperature for x ≤ -a may be obtained by changing x into -x. Inside the 1D-sphere, similar calculations leads to:
Now, let us consider a sphere in a 3D space. Goldenberg and Tranter  gave the following expression for temperature at the center of the absorbing sphere:
where subscripts 1 and 2 stands for the sphere and the surrounding medium respectively. Parameters are : . Since we assume λ 1 = λ 2 and D 1 = D 2, it follows that c = 0 and b = 1. By using the formal result 
the temperature rise at the center of the sphere in this particular case reads:
Outside the sphere, the temperature rise may be written as :
In case where a ≪√Dt, some algebra leads to:
A.2. Heating point source
For a 1D one point source, the heat equation reads :
where D, ρ and C are thermal diffusivity, density and specific heat capacity respectively. E has the dimension of an energy per surface unit and F(t) is the temporal profile of the laser pulse. In case where F(t) = δ(t), using Fourier and Laplace transformations, we show that :
Now, instead of a temporal Dirac source, we assume a time-dependent rectangular laser pulse shape such that:
Using convolution properties of Laplace transformation, the new temperature rise reads:
where (1-e -τs)/s is nothing but the Laplace transformation of the rectangular laser pulse. The latter expression may be written in the following form:
We are only interested in the temperature at t = τ thus implying the last term of Eq. (28) can be omitted (since it contributes only for t > τ). Finally, the temperature induced by a 1D point absorber is:
When setting E/τ = Aa, it is worth noting that the latter expression is consistent with Eq. (17) when the sphere radius a goes to zero.
Now, let us consider a point source in a 3D space. The procedure is strictly identical to the previous calculations. Also, the temperature rise induced by a 3D point absorber is:
where erfc denotes the complementary error function and E the absorbed energy. By setting (energy absorbed by an equivalent sphere of radius a during the interaction time τ) and in conditions where r ≪ √Dt, we can refind the expression (23). It is worth noting that θ(r,t) is nothing but the Green’s function that allows to recover the temperature inside an inclusion or a cluster performing the integral over the source of the Green’s function times the source strength .
B. Temperature for different cluster geometries
In this section, we propose to derive the expression of the temperature inside different cluster geometries using the above-mentioned Green’s formalism. Let us begin by the simplest case, a 3D spherical cluster, i.e. evenly distributed sources on spherical shells. This case is a reference since we have to find again the standard result of an absorbing sphere of radius a. Let us focus on the temperature rise at the center of the spherical cluster. This value is obtained by summing all point source contributions in the following way :
where a is the cluster radius and θ(r⃗,t) is the contribution of one point absorber. In order to compare this procedure with the reference, we assume that . Eq. (32) then reads:
where the angular integral has been reduced to 4π since we are dealing with spherical symmetry. Using formula 
Eq. (33) may be written as:
with α = 1/2√Dt. A rearrangement of the latter equation terms leads to the Goldenberg and Tranter  expression of temperature rise (21) at the center of the sphere:
Now, our procedure being validated, we are able to investigate cluster geometries of interest in the text. Let us begin by a very flat cylindric cluster. It is defined by a radius a and a thickness Δz. We emphasize we work in the case Δz << a. Under this condition, the temperature at the barycentric position is roughly given by:
Using Eq. (34), simple calculations leads to:
Finally, let us consider point sources arranged on a cone surface of negligeable thickness ΔR. Let a and h be the radius at the base of the cone and the height cone respectively. We assume that the revolution axis is along the z-axis of the (x,y,z) space. The temperature rise at the cone top induced by the point source distribution may be written as:
where R(z) is the cone radius defined by the intersection of the cone and a plane such that z = cste. So we have R(z) = az/h and r(z) = √z 2+R(z)2 that leads to:
with . It is worth noting that in case where a>> h, the latter expression tends to the flat cylinder cluster case.
The authors are grateful to Franck Enguehard for helpful discussions about heat transfer. Hervé Mathis is aknowledged for his comments about the manuscript. The authors also thank Pierre Grua and Jean-Pierre Morreuw (CEA/CESTA) for fruitful discussions. A. Dyan thanks the Commissariat à l’Energie Atomique for providing a post-doctoral grant.
References and links
1. J. J. De Yoreo, A. K. Burnham, and P. K. Whitman, “Developing KH2PO4 and KD2PO4 crystals for the world’s most powerful laser,” Int. Mater. Rev. 47, 113–152 (2002) [CrossRef]
2. C. W. Carr, H. B. Radousky, A. M. Rubenchik, M. D. Feit, and S. G. Demos, “Localized dynamics during laser-induced damage in optical materials,” Phys. Rev. Lett. 92, 087401 (2004) [CrossRef] [PubMed]
3. J. J. Adamset al, “Wavelength and pulselength dependence of laser conditioning and bulk damage in doubler-cut KH2PO4,” in Laser-induced damage in optical materials, Proc. SPIE 5991, 59911R-1 (2005)
4. C. W. Carr, M. D. Feit, M. A. Johnson, and A. M. Rubenchik, “Complex morphology of laser-induced bulk damage in KH2-xD2xPO4 crystals,” Appl. Phys. Lett. 89, 131901 (2006) [CrossRef]
5. C. S. Liu, C. J. Hou, N. Kiousis, S. G. Demos, and H. B. Radousky, “Electronic structure calculations of an oxygen vacancy in KH2PO4,” Phys. Rev. B 72, 134110 (2005) [CrossRef]
7. S. G. Demos, M. Staggs, J. J. De Yoreo, and H. B. Radousky, “Imaging of laser-induced reactions of individual defect nanoclusters,” Opt. Lett. 26, 1975–1977 (2001) [CrossRef]
8. S. G. Demos, M. Staggs, M. Yan, H. B. Radousky, and J. J. De Yoreo, “Investigation of optically active defect clusters in KH2PO4 under laser photoexcitation,” J. Appl. Phys 85, 3988–3992 (1999) [CrossRef]
9. J. Dijon, T. Poiroux, and C. Desrumaux, “Nano absorbing centers: a key point in the laser damage of thin film,” in Laser-induced damage in optical materials, Proc. SPIE 2966315 (1997)
11. M. D. Feit, A. M. Rubenchik, M. R. Kozlowski, F. Y. Génin, S. Schwartz, and L. M. Sheehan“Extrapolation of damage test data to predict performance of large area NIF optics at 355nm,” in Laser-induced damage in optical materials, Proc. SPIE 3578, 226–234 (1998)
12. R. W. Hopper and D. R. Uhlmann, “Mechanism of inclusion damage in laser glass,” J. Appl. Phys. 41, 4023–4037 (1970) [CrossRef]
13. M. D. Feit and A. M. Rubenchik, “Implications of nanoabsorber initiators for damage probability curves, pulse-length scaling and laser conditioning,” in Laser-induced damage in optical materials, Proc. SPIE 5273, 74–82, (2004)
14. A. Dyan, F. Enguehard, S. Lallich, H. Piombini, H. Mathis, and G. Duchateau, “Laser-induced damage in KH2PO4 and D2xKH2(1-x)PO4: influence of the laser probed volume and optimization of the conditioning through a revisited thermal approach,” submitted to Optics Communications
18. L. Gallais, J. Y. Natoli, and C. Amra, “Statistical study of single and multiple pulse laser-induced damage in glasses,” Opt. Express 10, 1465–1474 (2002) [PubMed]
20. H.S. Carslaw and J.C. Jaeger, Conduction of Heat in Solids, Oxford Science Publications, Second Edition (1959)
21. K. E. Montgomery and F. P. Milanovich, “High-laser-damage-threshold potassium dihydrogen phosphate crystals,” J. Appl. Phys. 683979-82 (1990) [CrossRef]
22. P. K. Bandyopadhyay and L. D. Merkle, “Laser-induced damage in quartz: a study of the influence of impurities and defects,” J. Appl. Phys. 631392–1398 (1988) [CrossRef]
23. F. Bonneauel al, “Study of UV laser interaction with gold nanoparticles embedded in silica,” Appl. Phy. B 75, 803–815 (2002) [CrossRef]
24. N. W. Ashcroft and N. D. Mermin, Solid state physics, Brooks Cole, first edition, 1976
25. C. Kittel, Introduction to solid state physics, Wiley, seventh edition, 1995
26. F. Y. Génin, A. Salleo, T. V. Pistor, and L. L. Chase, “Role of light intensification by cracks in optical breakdown on surfaces,” J. Opt. Soc. Am. A 182607-16 (2001) [CrossRef]
27. I. Drevensek, M. Zgonik, M. Copic, R. Blinc, A. Fuith, W. Schranz, M. Fally, and H. Warhanek, Phys. Rev. B49, 3082–3088 (1994) [CrossRef]
28. Z. L. Xia, Z. X. Fan, and J. D. Shao, “A new theory for evaluating the number density of inclusions in films,” Appl. Surf. Sci. 252 (23), 8235–8238 (2006) [CrossRef]
29. S. Papernov and A. W. Schmid, “Correlations between embedded single gold nanoparticles in SiO2 thin film and nanoscale crater formation induced by pulsed-laser radiation,” J. Appl. Phys. 92, 5720–5728 (2002) [CrossRef]
30. M. Runkel, J. DeYoreo, W. Sell, and D. Milam, “Laser conditioning study of KDP on the optical sciences laser using large area beams,” in Laser-induced damage in optical materials, Proc. SPIE 3244, 51 (1998)
31. E. S. Bliss, “Pulse duration dependence of laser damage mechanism” Opto-electronics 3, 99 (1971) [CrossRef]
32. R. M. Wood, Laser-induced damage of optical materials, Institute Of Physics publishing series in optics ans optoelectronics, Bristol and Philadelphia (2003)
33. B. C. Stuart, M. D. Feit, S. Herman, A.M. Rubenchik, B. W. Shore, and M. D. Perry, “Nanosecond-to-femtosecond laser-induced breakdown in dielectrics,” Phys. Rev. B 53, 1749–1761 (1996) [CrossRef]
34. A. Dyan and G. Duchateau, CEA, Centre d’Etudes du Ripault, BP 16, 37260 Monts, France, are preparing a manuscript to be called “Laser-induced damage by a nanosecond pulse: a method coupling heat transfer, Mie’s theory and microscopic processes”
35. H. Goldenberg and C.J. Tranter, “Heat flow in an infinite medium heated by a sphere,” Br. J. Appl. Phys. 3, 296–298 (1952) [CrossRef]
36. M. Sparks, “Theory of laser heating of solids: metals,” J. Appl. Phys. 47, 837–849 (1976) [CrossRef]
37. I.S. Gradshteyn and I.M. Ryzhik, Table of Integrals, Series, and Products, Alan Jeffrey Editor, Fifth Edition (1994)