Cell size distribution in a random tessellation of space governed by the Kolmogorov-Johnson-Mehl-Avrami model: Grain size distribution in crystallization

The space subdivision in cells resulting from a process of random nucleation and growth is a subject of interest in many scientific fields. In this paper, we deduce the expected value and variance of these distributions while assuming that the space subdivision process is in accordance with the premises of the Kolmogorov-Johnson-Mehl-Avrami model. We have not imposed restrictions on the time dependency of nucleation and growth rates. We have also developed an approximate analytical cell size probability density function. Finally, we have applied our approach to the distributions resulting from solid phase crystallization under isochronal heating conditions.


I. INTRODUCTION
In this paper we consider the subdivision of a D-dimensional Euclidean space into disjoint regions created after a process of random nucleation and growth. Random subdivisions can be obtained by several different methods, amongst which Poisson-Voronoi and Johnson-Mehl tessellations 1 have been widely studied. The Poisson-Voronoi tessellation is obtained by randomly picking several points, the seeds P i , by a Poisson process. Next, the space is subdivided in cells, C i , by the rule: C i contains all points in space closer to P i than to any other seed. This cellular structure is extensively applied in many diverse scientific fields including biology, 2,3 computer science, 4,5 materials science, 6,7 astrophysics, 8,9,10 medicine, 11 agriculture, 12 quantum field theory, 13 and sociology. 14 The space tessellation can be fully characterized by means of the probability density function (PDF), f (s), which is the probability that a cell has a size between s and s + ds.
The properties of the PDF of the Poisson-Voronoi tessellation have been extensively studied both theoretically 1,15,16 and numerically. 7,16,17,18,19,20,21,22 It is well known that the Poisson-Voronoi tessellation PDF is described by a gamma distribution where Γ is the gamma function, ν is a parameter that is dependent on the dimension D, i.e. ν = 2, 3.584 and 5.586 for D = 1, 2 and 3 respectively, and E is the expected cell size, It is worth mentioning that Eq. (1) has been analytically derived for the one-dimensional case where ν=2 is an exact result. 15 Conversely, for the two and three-dimensional cases, the validity of Eq. (1) is supported by analytical approximations and numerical fits.
Our main interest is the characterization of grain morphology related to crystallization.
In general, the crystallization of most materials takes place by means of a nucleation and growth mechanism: nucleation starts with the formation of small atom clusters of the new stable phase in the metastable phase. Subsequently, clusters with sizes greater than the critical, or nuclei, start to grow by incorporating neighboring atoms of the metastable phase.
During this growth, grains impinge upon each other. Finally, the structure of the new stable phase consists of disjoint regions or crystals separated by grain boundaries. The evolution of crystallization and grain size distributions is entirely determined by the nucleation rate density I and the grain linear growth rate G. When nucleation takes place for a very short time, its rate may vanish before the onset of particle growth (site saturated nucleation). 23,24 In this case, the crystal structure is equivalent to a Poisson-Voronoi tessellation provided that nucleation is Poissonian through the whole space and growth is isotropic.
Conversely, continuous nucleation takes place when nucleation and growth occur at the same time. In general, there is an energy barrier for nucleation and growth to happen.
Thus, I and G depend on temperature. For the particular case of isotropic and isothermal transformations, where I and G are constant, the resulting crystal structure corresponds to the well-known Johnson-Mehl tessellation. 1 For this tessellation, Axe et al. 25 43 and film growth on solid substrates. 44 In Section II we will describe the basic concepts of KJMA theory and will focus our attention on those aspects that are useful to the development of our work. Section III is devoted to the calculation of the expected value and variance of the distributions related to the KJMA tessellations. In Section IV we will derive a simpler approximate relation for the variance and will check its accuracy. As an application of the previous results, in Section V we will derive an approximate grain size PDF which is the superposition of gamma distributions. Finally, at the end this section we will verify that the grain radius PDF can be expressed as well as the superposition of Gaussian distributions.

II. THE KOLMOGOROV-JOHNSON-MEHL-AVRAMI THEORY
The KJMA theory 45,46,47 describes in a very simple form the kinetics of transformations governed by nucleation and growth that satisfy the following assumptions: i nucleation must be Poissonian through the entire space; ii the volume of an arbitrary grain is much smaller that the volume of the system; iii the crystal growth rate is isotropic.
On the basis of these premises, Kolmogorov calculated the evolution of the transformed fraction, X(t), through the probability, p(t), that an arbitrary point O has not crystallized, i.e., the probability that no nuclei able to transform O will be formed during the time interval [0, t], where g D is a geometrical factor related to the shape of the crystal -for a D-dimensional sphere g D = π D/2 /Γ(D/2+1) -and r(t, τ ) is the minimum distance between O and a nucleus created at τ , so that the nucleus would not transform O.
Based on geometrical arguments, Avrami deduced the following relation: where ∂ t v(t, τ ) and ∂ t v ex (t, τ ) are respectively the actual and extended average volumetric growth rate at time t for grains nucleated at time τ . The word extended refers to the volume a grain would attain if nuclei grew through each other and overlapped without mutual interference.
The integration of Eq. (4) leads to 48 Finally, integration of Eq. (5) gives Avrami's well-known formula The calculation of X ex (t) is straightforward and obtained by simply neglecting the impingement between nuclei The combination of Eqs. (6) and (7) gives Eq. (3). As it is well known, Avrami and Kolmogorov deduced the same relation using different approaches.
Note that in Eq. (7) it is assumed that the nucleation rate is not affected by the shrinking of the untransformed phase. In the calculation of X ex (t) the phantom nuclei are taken into account. Avrami designated as phantom nuclei those nuclei that are formed in the transformed fraction and therefore do not contribute to the formation of new grains. Indeed, the actual nucleation rate can be defined as Concerning the limitations of the KJMA theory, it also holds in the case of anisotropic growth provided that the grains have a convex shape and are aligned in parallel. 31 Moreover, the KJMA theory provides a good approximation when the anisotropy is moderate or for soft impingement. 49 However, KJMA theory fails when nucleation is non-random, 50 when growth is anisotropic, 33,51,52 when growth stops before crystallization is complete 49 and when the incubation time is not negligible. 53

III. STATISTICAL PROPERTIES OF THE KJMA CELL SIZE DISTRIBUTION
The cell size distribution is characterized by its PDF, f (s), the probability that a cell has a size between s and s + ds. From its definition it is obvious that f (s) must be normalized ∞ 0 f (s)ds = 1.
To analyze the properties of the cell size distribution, we will consider the contribution of the crystals formed at a time τ over a time interval dτ (τ -crystals). We will call the cell size distribution of the τ -crystals the τ -distribution. Accordingly, we define the PDF of the τ -crystals, f τ (s), as the probability that a τ -crystal has a volume between s and s + ds.
From the definition of f τ (s), it is also apparent that f τ (s) must be normalized: f (s) is simply the addition of the contributions of the τ -crystals over the time interval in which their nucleation takes place Note that the denominator in Eq. (11) ensures that f (s) is normalized if all f τ (s) are normalized.
In the following sections, we will present the expected grain size and the variance of the cell size distribution and their relationship with the equivalent parameters of the τ -distributions.

A. Expected grain size
It is well known that the expected grain size, E, is the inverse of the final grain density: Likewise, the expected value, E τ of a τ −distribution is simply the final average grain size of a τ -crystal normalized to the total volume: Introducing Eq. (4) into Eq. (13) leads to where the extended average growth rate is given by Note than once the evolution of the transformed fraction, X(t), is known -i.e. the solution of Eq. (3) -the calculation of E τ is straightforward.
Besides, the final space fraction occupied by the τ -crystals, X τ , can be calculated from the integration over the entire space of the probability that a point P in the space belongs to a tau crystal nucleated at O. Since the system is homogeneous and isotropic, this probability only depends on the distance b between O and P . Therefore, where P τ (b) is the probability that a point P , separated by a distance b from the nucleus O, belongs to the crystal nucleated at O. To calculate P τ (b), we will use the same approach that Kolmogorov used for the deduction of Eq. (3). Since nucleation is Poissonian, P τ (b) is given by the probability that no nucleus is formed that could transform P before O does so.
Thus, the nuclei formed at z that could transform P before O does so, are located in a D-sphere of radius r(t b , z) around P . Therefore, according to Eq. (3), P τ (b) is given by The previous integral spans the time interval [0, t b ] since no nucleus formed after t b could transform P . Comparison of Eq. (7) with Eq. (18) gives Finally, if we introduce the value of P τ (b) given by Eq. (18) into Eq. (16) and we change the variable b by t b , we obtain Alternatively, the expected value, E τ is the ratio between the space fraction occupied by the τ -crystals and the density of τ -crystals: As expected, substitution of Eqs. (20) and (8) into Eq. (21) delivers Eq. (14). Moreover, the integration of X τ over the whole time interval where nucleation takes place gives the total transformed fraction, 1: We will end this subsection verifying that the value of E evaluated from the τ -distributions coincides with the value given at the beginning of this subsection [Eq. (12)]:

B. Variance of the grain size distribution
To determine the variance we will adapt the development of Gilbert 1 for a Poisson-Voronoi tessellation to our case. First, we define a new PDF, f * (s), as the PDF of the crystals that is proportional to f (s) and to s, because a large crystal has a proportionally greater chance of containing the point O. Therefore, The constant of proportionality, E −1 , has been deduced by imposing normalization: From the definition of f * (s) it can be easily proved that where E * is the expected value of f * (s).
Once E * is known, the calculation of the variance is simple. To obtain E * we first analyze the contribution of the τ -crystals. To do so, we define f * τ (s) as the PDF of the τ -crystals that contain a given arbitrary point O: i.e., if we pick an arbitrary point O, f * τ (s) is the probability that a τ -crystal has a volume between s and s + ds and contains the point O.
Accordingly, f * τ (s) is proportional to f τ (s) and to s: On the other hand, the integration of f * τ (s) over all possible volumes is the probability that an arbitrary point O belongs to a τ -crystal. This probability is the fraction of the space occupied by the τ -crystals X τ : taking into account that, and combining Eqs. (27), (28) and (21) we obtain f * τ (s): Then the expected value of f * τ , E * τ , is It can be easily verified that E * is related with E * τ through Therefore the contribution of the τ -crystal to the expected value E * is X τ E * τ . In addition, this contribution is the integration over the entire space of the probability that a differential volume around a point P in the space belongs to the same τ -crystal as O. Since the system is homogeneous and isotropic, this probability only depends on the distance b between O and P , where P * τ (b) is the probability that two points, O and P , separated by a distance b belong to the same τ -crystal, The integration domain covers the entire space, dV Q is the D-volume differential around a point Q, I(τ )dV Q dτ is the probability that a nucleus is formed at Q at the time τ during the time interval dτ and P * τ (b, Q) is the probability that both O and P belong to the same crystal nucleated at Q (see Fig. 1 is given by the probability that no nucleus is formed that could transform O or P before Q does, then where r(t O , z) and r(t P , z) are the minimum distance between O, P and a nucleus created at the time z, so that the nucleus would not transform O and P respectively (see Fig. 1). V I is the volume intersection between two D-spheres of radius r(t O , z) and r(t P , z) centered at O and P , respectively (gray region in Fig. 1). The subtraction of the term V I is in accordance with the fact that it has been accounted twice in the first and second integrals in Eq. (35a).
It can be easily proved that the variance of the τ -distributions, var τ , is given by Finally, we will check if the variance of the distribution determined from the decomposition of f (s) into τ -PDF, Eq. (11), gives the expected result, Eq. (26): At this point, we would like to point out that the results obtained so far are exact and general, i.e., we have not made any assumption concerning f (s) and f τ (s).

IV. APPROXIMATE VARIANCE
According to our previous analysis, the exact calculation of the variance is reduced to the calculation of the parameters E and E * in Eq. (26). While the calculation of E is straightforward, the evaluation E * is more cumbersome. Indeed, when compared to Monte-Carlo algorithms, 34 its numerical calculation is more complex without representing any significant reduction in computing time. That is because there are several integrals nested and, in particular, the calculation of the intersection volume V I is complex. When r(t O , z) ≫ b or b ≫ r(t O , z), V I tends towards being a D-sphere of radius r(t O , z) and 0, respectively. On the other hand, when r(t O , z) ≈ r(t P , z) the shape of V I roughly approaches a D-sphere.
Since the width of V I (see Fig. 1) is r(t O , z) + r(t P , z) − b, we approximate V I by a D-sphere of diameter r(t O , z) + r(t P , z) − b: It is worth noting that the previous approximation also works for the limiting cases z). Furthermore, for the one-dimensional case it can be easily verified that Eq. (39) is exact (in Appendix A we derive P * τ for D = 1). Finally, the approximate solution (from here on approximation I) is obtained by substitution of Eqs. (7) and (39) into Eq. (36): Therefore, the calculation of P * τ (b, Q) is simple provided that the evolution of the transformed fraction, X(t), is known. Analytical exact solutions for X(t) are restricted to three particular situations under isothermal conditions: time-independent growth and nucleation rates, time-independent growth rate and nucleation rate proportional to a power of time, 54 and site saturated nucleation. A quasi-exact solution of the KJMA model has recently been obtained under continuous heating conditions. 55 Moreover, there are numerical methods which allow a simple and fast calculation of X(t) for an arbitrary time dependence of the nucleation and growth rates. 34 We have analyzed the distribution emerging from solid phase crystallization under isochronal heating conditions, i.e. heating at a constant rate, to check the accuracy of approximation I, Eq. (40), in the case of time dependent nucleation and growth rates. To work with realistic parameters we have taken those of amorphous silicon crystallization, 38,39,56 in which the nucleation and growth rates are described by an Arrhenius temperature dependence where T is the temperature in Kelvin and k B is the Boltzmann constant. In Table I we summarize the corresponding parameters. When the temperature is raised at a constant  (22). Calculations that exhibit discrepancies larger than 10 −6 were rejected.
As is apparent from Fig. 1, the approximation of V I by a D-sphere of diameter equal to its width, Eq. (39), results in an underestimation of V I , which leads to an undervaluation of E * τ and of var τ . The latter conclusion can be verified in Fig. 2, where the evolution of E * τ and var τ with τ is shown. Although approximation I gives an accurate value of E * τ , the approximate value of var τ shows a significant deviation from the exact value. The reason is that in the evaluation of var τ , Eq.  through the relation between the exact and approximate values of var 0 . To cover a wide range of distributions, we will recall the results given in Ref. 35. In this work it was shown that the shape of the grain size PDF was practically insensitive to the heating rate, but it depends mainly on the ratio E N /E G , i.e., the relative evolution of the nucleation and growth rates with time. The limit E N /E G → 0 corresponds to site saturated nucleation while E N /E G = 1 coincides with the isothermal case.
In Fig. 3 we have plotted the exact and approximate values of var 0 as well as their ratio. This result is due to approximation I, which is based on a geometrical approach that is fairly insensitive to the relation between nucleation and growth rates.
Since the ratio between the exact and the approximate var τ is nearly constant, we can obtain a significantly more accurate approximate value for var τ by simply multiplying it by the corresponding proportionality constant. This constant only depends on the growth dimensionality. We have chosen the values of 2.07 and 1.32 for the 3D and 2D cases, respectively. These values correspond to E N /E G = 1, i.e., they correspond to the isothermal case with E G and E N constant in time. With this correction (from now on approximation II ), the relative error in the calculation of the variance diminishes to less than 2% and for E N /E G ≥ 1 the relative error is less than 0.2%. In Fig. 4 we have plotted the exact and the approximate values of the distribution variance and E * with respect to the ratio in general, the inaccuracies in the evaluation of the variance and E * are significantly smaller.
The reason is that both parameters depend exclusively on E τ and E [Eqs. (26) and (32)].
From Fig. 2 it is clear that the error related to E τ is significantly smaller than the error in the evaluation of var τ , while the calculation of E is exact.
From now on, when will always use approximation II in the calculation of the approximate values of E * and E * τ .

V. APPROXIMATE CELL SIZE PROBABILITY DENSITY FUNCTION
One application of the preceding analysis of the statistical properties of grain size distribution is the derivation of a PDF. If we choose a set of f τ (s) such that their expected value coincides with the result of Eq. (14) and their variance is equal to the value given by Eq. (37), then the variance and expected values of the PDF obtained from Eq. (11) will be exact, i.e., the PDF obtained from Eq. (11) will have the same variance and expected values as the actual PDF. Indeed, Pineda et al., 16 apply this approach in the case of tessellations generated by random nucleation processes where the growth rate was assumed to be constant. The agreement between their approximate PDF and Monte-Carlo simulations was remarkable. However, they did not notice that the expected value and the variance of their approximate PDF were exact. Concerning the particular choice of the f τ (s) functions, we will consider two cases: gamma and Gaussian distributions.

A. Gamma distribution
Given that in a τ -distribution the nucleation events are simultaneous and the τ -nuclei are randomly distributed, we can assume that f τ (s) is similar to the PDF resulting from a process of site saturated nucleation, i.e. that the PDF is that of a Poisson-Voronoi tessellation.
As explained in the introduction, the gamma distribution is the exact PDF for the onedimensional case while it provides a very accurate result for the two-and three-dimensional cases: Since E τ is already the expected value, the problem comes down to the determination of the exponent ν τ . Indeed, ν τ can be calculated from the following property of the gamma To check the accuracy of the PDF we compare them to some Monte-Carlo simulations.
The Monte-Carlo algorithm 34 consists in dividing the space into a cubic lattice. Cells are assigned to nuclei randomly. The nucleation time of each nucleus is precisely calculated from the nucleation rate and it is recorded to evaluate the exact evolution of the grain growth.
Finally, each cell is assigned to the nucleus that first reaches this cell. The evolution of the grain growth transformation is checked whenever the grain growth is equal to the size of a cell in order to avoid incorrect cell assignation due to shielding effects. In particular we consider the crystallization of amorphous silicon under isochronal heating at 40 K/min (Table I). The results of the calculations for three-and two-dimensional growth are given in Figs. (5) and (6), respectively. We have calculated the PDF using the exact, Eq. (36), and approximate, Eq. (40), values of E * τ . We will refer to these PDFs as the approximate III and approximate IV PDFs, respectively. The validity of the selection of a gamma distribution for the f τ (s) functions is confirmed by the good agreement between the calculated PDF and the Monte-Carlo simulations. The excellent agreement between the approximate III and the approximate IV PDFs is also noteworthy -the relative difference is less than 0.1% for s/E > 0.05. Therefore, the approximate calculation of E * τ is useful to obtain a simple and accurate PDF for a KJMA tessellation. However, the complexity and the computing time required for their evaluation is significantly different. For instance, the calculation of the approximate III PDF typically takes more than thirty times the time required for the calculation of the approximate IV PDF. The reason for this significant simplification and reduction in computing time is that the approximate IV PDF saves us from having to evaluate the inner and more complex integral of Eq. (36).
To confirm this last conclusion, we have calculated the PDF for the case of site saturated nucleation. This case corresponds to E N /E G = 0 in Figs. (3) and (4) and is the case that exhibits the greatest discrepancy between the exact and approximate values of the variance and E * . Therefore, it should give the worst agreement between approximate III and the approximate IV PDFs. As is apparent from Fig. (7) here again the agreement between the PDFs obtained from the exact calculation of E * τ , from the approximate calculation of E * τ and from Monte-Carlo simulation is excellent. Only small deviations of the approximate III from the approximate IV PDF are distinguishable for s ≈ E. Finally, for the one-dimensional case and for site saturated nucleation, both the approximate III and approximate IV PDFs turn into the exact PDF (see Appendix B).

B. Gaussian distribution
When analyzing the crystallization morphology, it is often better to use the grain radius distribution instead of the grain size distribution. The grain radius, r, of a grain of size s, is defined as the radius that will have a D-sphere of volume s. For instance, for three-dimensional growth Given the grain size PDF, the grain radius PDF, g(r), can be easily derived; for threedimensional growth, On the other hand, for site saturated nucleation and three-dimensional growth it has been shown that g(r) is accurately described by a Gaussian distribution 34 where E g is the expected value of the grain radius PDF and σ is the standard deviation where var g is the variance of the grain radius distribution. The fitted parameters were actually E g = 0.608/n 1/3 0 and σ = 0.0883/n 1/3 0 , where n 0 is the nuclei density. Conversely, g(r) can be determined by means of Eq. (44) where f (s) is a gamma distribution with ν = 5.581. 16 Under this approach, it can easily be proved that where E = 1/n 0 is the expected value of the grain size distribution (see Appendix B).
Therefore, for a KJMA tessellation we can obtain g(r) as the superposition of g τ (r) PDFs. In contrast with the previous subsection, we will now assume that g τ (r) are Gaussian distributions: To evaluate the g τ (r) PDFs we need to know their expected value, E g,τ , and their variance, var g,τ . These parameters can be easily derived from the statistical parameters of the τ -distributions, E τ and var τ , by means of Eq. (47). In Fig. (8) we have plotted the grain  For the one dimensional case, the extended transformed fraction is and the expected value E τ becomes With regard to the calculation of P * τ (b), we have split the entire space into three regions. The first corresponds to Q located at the left side of O, in this case r P = r O + b and The second region corresponds to Q located between O and P , then r P + r O = b and And the third region corresponds to Q located at the right side of P , r O = r P + b and Finally, combining Eqs. (A3), (A4) and (A5) with (34) we obtain Once P * τ (b) is known, Eq. (33) combined with Eq. (37) delivers var τ .

APPENDIX B: SITE SATURATED NUCLEATION
When nucleation is completed prior to crystal growth, the nucleation rate can be approximated to I(t) ≈ n 0 δ(t) where n 0 is the density of nuclei and δ is the Dirac delta function.
In this case, the extended transformed fraction becomes Moreover, the fraction of space occupied by a τ -crystal, X τ , is reduced to δ(τ ), as expected.
With respect to var 0 , its value is determined by E * 0 and E 0 by means of Eq. (37). E * 0 is given by Eq. (33). We therefore need to evaluate P * 0 (b), the value of which depends on which relation for P * 0 (b, Q) we use: the exact one Eq. (36) or the approximate one Eq. (40). We will evaluate P * 0 (b) for the one-dimensional case because in this case the approximate solution coincides with the exact one. Specifically, substitution of Eq.