# Numerical Study of Crystal Size Distribution in Polynuclear Growth

###### Abstract

The crystal size distribution in polynuclear growth is numerically studied using a coupled map lattice model. The width of the size distribution depends on , where is the growth rate at interface sites and is the diffusion constant. When is sufficiently small, the width increases linearly with and saturates at large . Monodisperse square and cubic crystals are obtained respectively on square and cubic lattices when is sufficiently small for a small kinetic parameter . The linear dependence of on in a parameter range of small is explained by the eigenfunction for the first eigenvalue in a two-dimensional model and a mean-field model. For the mean-field model, the slope of the linear dependence is evaluated theoretically.

## I Introduction

Crystal growth has been extensively studied in various research fields such as applied physics, metallurgy, chemical engineering, and mineralogy. It is also studied in statistical physics as a typical nonequilibrium phenomenon. Nucleation, surface kinetics, and diffusion are important factors in crystal growth rf:1 (). Crystal growth is one of the Stefan problems for moving interfaces, which is difficult to analyze mathematically. Various simulation methods have been proposed to study crystal growth far from equilibrium. The phase-field model is a partial-differential equation in which the solid and melt phases are expressed with a continuous variable corresponding to the order parameter rf:2 (); rf:3 (). We proposed a coupled map lattice model for crystal growth rf:4 (); rf:5 (); rf:6 (). A similar type of order parameter is defined on a lattice, and coupled maps are used for the time evolution. The numerical simulation of the coupled map lattice is much faster than that using the phase-field model. One reason is that the width of the interface between the solid and melt phases is only one lattice constant in the coupled map lattice model. On the other hand, the order parameter changes smoothly at the interface in the phase-field model owing to the diffusion term in the equation for the order parameter. We performed numerical simulation of crystal growth of dendrites and diffusion-limited aggregations (DLA’s) using coupled map lattices. These complicated fractal patterns can be easily reproduced in the coupled map lattice model, because the surface tension effect can be set to zero in the coupled map lattice model. In previous papers, we studied crystal growth starting from a single crystal seed. In this paper, we study polynuclear growth in which there are many crystal seeds initially, and investigate the crystal size distribution at the final stationary state.

Crystal size distribution (CSD) is one of important topics in crystallography. Various theories of CSD were proposed. Becker and Döring studied the time evolution of the size distribution in the nucleation process rf:7 (). Randolph and Larson discussed stationary CSDs in the research field of chemical engineering rf:8 (). Marsh applied the CSD theory in geology and found an exponential size distribution of mineral crystals deposited from magma rf:9 (). In their model, the exponential distribution is derived by the balance of the outflow, the inflow, and the growth of crystals. In the synthesis of micro- and nanocrystals, the final size distribution after the crystal growth is important. Various methods of characterization and control of the crystal size distribution have been studied. In particular, the control of the crystal size distribution have been investigated in colloidal sciences rf:10 (); rf:11 (). The size distributions similar to the gamma distribution or the log normal distribution are often observed in experiments. The size distribution is not symmetric around the average value but has a longer tail for large . It is practically important to produce almost monodisperse crystals, because the monodisperse crystals are useful for various applications such as catalysts, magnetic fluids, and drug delivery systems. Colloidal crystals can be produced from monodisperse microcolloidal particles, which can be applied to various systems such as photonic crystals rf:12 ().

The crystal size distribution depends on many factors such as nucleation processes, surface conditions of crystals, and various experimental conditions. The CSD problem is still an open problem. For the synthesis of monodisperse crystals, it is important that the time scales of nucleation and growth processes are separated, and the size distribution of the crystal seeds is almost mono-disperse at the initial stage of the growth process. In this paper, we consider a problem that the size distribution of the crystal seeds is almost monodisperse at the initial stage of the crystal growth, but the initial positions of the crystal seeds are randomly distributed. The final size distribution changes owing to the competition of the diffusion and surface kinetics. Our model might be applied to a two-dimensional crystal growth on a substrate in that the positions of crystals are fixed in time. However, the purpose of this study is to understand qualitatively a condition wherein the final crystal size distribution becomes almost monodisperse.

## Ii Two-Dimensional Coupled Map Lattice Models for Polynuclear Growth

We first study a two-dimensional coupled map lattice model on a square lattice. There are two variables in our coupled map lattice for the solution growth, namely, the order parameter and the concentration , where and denote respectively a lattice site and a discrete time. The order parameter is 1 at crystal sites, 0 at solution sites, and between 0 and 1 at interface sites. The time evolution has two stages, namely, diffusion and surface kinetics. The diffusion is expressed as

(1) |

where is a diffusion constant. The surface kinetics is expressed as

(2) |

where is a parameter for the growth rate. The time evolution Eq. (2) is applied only at interface sites satisfying . If goes over a critical value 1, interface sites change into crystal sites and the neighboring solute sites are assigned to interface sites. By repeating Eqs. (1) and (2), polynuclear crystal growth occurs. The parameter is 1 at nonflat sites and takes a small value at flat sites. The flatness of an interface site is determined by counting the number of crystal sites in the four nearest-neighbor sites of the interface site. If the number is 1, is set to be , and if the number is larger than 2, is set to be 1. For sufficiently small , square crystals appear, because kink sites and concave sites are quickly occupied, but it takes a rather long time for flat sites to grow. The parameter represents an effect of the surface tension on the surface kinetics. In the case of , the surface tension effect disappears. In Eq. (2), the increase in the order parameter is equal to the absorption of at interface sites, that is, the conservation law of is satisfied. In our model of Eqs. (1) and (2), only the deposition process from solution to crystal is assumed, that is, the inverse process of dissolution is not taken into consideration.

We have performed numerical simulation of the two-dimensional coupled map lattices. Initially, 10 crystal seeds of size 1 are randomly set in a square of . The initial value of is uniform and set to be , and the parameter is assumed to be 0.005. Periodic boundary conditions are imposed. After a long-time iteration of Eqs. (1) and (2), becomes zero owing to the diffusion and the absorption of at interface sites, and the crystal growth process is completed. Figure 1(a) shows crystal sites at and at the final time. The crystal size is widely distributed. The crystal sizes are small in a region where initial seeds are densely distributed. In our model, the randomness originates only from the initial configuration of crystal seeds. Figure 1(b) shows crystal sites at and . Crystals take a clear square shape, and the sizes are almost uniform. As is large and is small, the concentration tends to be uniform and the crystal sizes become homogeneous even for the random initial configuration of seeds. To show an effect of the initial seed size distribution, the sizes of four crystal seeds at are set to be , and the sizes of the other seeds are set to be 1. Figures 1(c) and 1(d) show crystal sites at (c) and and (d) and at a final time. At and , the crystal sizes for the six seeds of size 1 are the same as those in Fig. 1(b), and the crystal sizes for the four larger seeds of size 9 are larger by than those in Fig. 1(b). At and , the crystal sizes change for all crystals from the values in Fig. 1(a), but the change in size from Fig. 1(a) depends on the position. The final crystal size depends on the initial seeds, but the dependence of the final size distribution on the initial size distribution is not so simple. We investigate only the case that the initial sizes are all 1 hereafter.

We have calculated the CSD using a larger square lattice of . The size is defined as the sum of in crystal sites and interface sites around each crystal seed. At crystal sites, takes a value of 1. Then, represents the area of a two-dimensional crystal. Figures 2(a) and 2(b) show the size distribution at (a) and , and (b) and . (There are a few crystals that collide with neighboring crystals, but the sizes of those crystals are not considered in the size distribution. The collision probability decreases if the system size increases for a fixed number of crystal seeds. In the production of colloidal particles, the electric double layer around colloidal particles prevents the collision.) The size distribution is rather wide at and . The distribution is asymmetric and has a longer tail for large . The size distribution is very narrow at and , which implies that crystals are almost monodisperse.

The average of crystal size is estimated as , because the total sum of is equal to owing to the conservation law of . The inhomogeneity of the crystal size is characterized by the width or the standard deviation of the size distribution. Figure 2(c) shows the widths for four different ’s at the same and . The width depends strongly on but hardly depends on under the condition of const. If and are small and the one-step increments of and in Eqs. (1) and (2) are sufficiently small, Eqs. (1) and (2) can be approximated by differential equations. It can be shown that the final stationary states depend only on the ratio owing to a scale transformation of time if the differential approximation is good. Figure 2(d) shows a relationship between and the ratio in a logarithmic plot. The diffusion constant is changed as and for , and then is changed as () for . The width increases linearly as in a range of small and tends to saturate at large . That is, monodisperse crystals are obtained in the limit of , and the dispersion is proportional to when is sufficiently small.

In the coupled map lattice model of Eqs. (1) and (2), there is a possibility that neighboring crystals become closer to each other and collide with each other as they grow. To neglect the effect of the collision, we can propose a simpler coupled map lattice model. Equation (2) is rewritten as

(3) |

where , and and are interpreted respectively as the area and radius of circular crystals. Here, the growth velocity of crystals of radius is assumed to be inversely proportional to the radius as . This type of relation is satisfied under the condition of diffusion-limited growth. Equation (3) is expressed as

(4) | |||||

(5) |

where is rewritten as . The form of Eqs. (4) and (5) is almost equal to that of Eq. (2); however, the time evolution of Eqs. (4) and (5) is applied only at the positions of crystal seeds in this model. That is, in contrast to Eq. (2), the variable increases indefinitely even if becomes larger than 1. Growth patterns such as square patterns do not appear, because grows only at the seed points. The magnitude of is interpreted as crystal size. In other words, we consider a situation that the radii of crystals are sufficiently small in contrast to the distances among different crystals, and the growth patterns are invisible. In the following numerical simulation, the number of initial seeds is set to be 180, and the system size is . The crystal growth stops for sufficiently large , when decreases to 0. Figure 3(a) shows the crystal size distribution at and . Figure 3(b) shows the size distribution at and . When is small and is large, the size distribution becomes narrow. The rhombi in Fig. 3(c) show the relationship between the width of the size distribution and the ratio , where is changed as () for . The dashed line denotes . The width increases from 0 linearly in a parameter region of small and tends to saturate at large also in this simplified model.

Because Eqs. (1) and (5) are linear equations, an initial value problem for can be solved using the eigenvalues and eigenfunctions as , where and are eigenvalues and eigenfunctions, and is the expansion coefficient. The first eigenvalue and the corresponding eigenfunction can be numerically evaluated by a long-time iteration of the same equation Eqs. (1) and (5) and the normalization . When is sufficiently small, the first eigenvalue is close to 1 and the first component of decays slowly or the other components decay quickly. Under this condition, the growth rate of each crystal is expected to be proportional to the first eigenfunction at the crystal-seed point, because the concentration decays as and increases to compensate for the decrease in . The crystal-size distribution can be evaluated approximately using numerically estimated by assuming that the crystal size is proportional to the eigenfunction and using the fact that the average crystal size is equal to . The pluses in Fig. 3(c) show the width of the size distribution evaluated from . The width is close to the results denoted by rhombi obtained by direct numerical simulation of Eqs. (1), (4), and (5) for small values of . At large , the contribution from the other eigenfunctions cannot be neglected, and some deviation is observed.

## Iii Voronoi Tessellation and a Mean-Field Model of Polynuclear Growth

The Voronoi tesselation is a method of space partitioning. In the construction of the Voronoi tessellation, central points are randomly distributed, and the whole plane is partitioned into territories of each point. The boundary line between two territories is determined by the perpendicular bisector of the two central points. Each Voronoi cell is a set of points surrounded by perpendicular bisectors. The Voronoi tesselation is applied to various areas such as grains of crystals and territories of animals. The size distribution of Voronoi cells is approximately given by the gamma distribution: rf:13 (); rf:14 ()

(6) |

where is the area of each Voronoi cell, is a parameter, and .

In our model, the crystal grows fast by absorbing the surrounding solute when is sufficiently large. The absorption range or the territory of each crystal seed might be approximated by the Voronoi cell. We interpret that the central point in each Voronoi cell corresponds to a seed point of a crystal. First, we consider a case that all the solutes in each Voronoi cell are absorbed into the central seed of a crystal. In that case, the crystal size is proportional to the area of the Voronoi cell. In that case, the crystal-size distribution is expressed as

(7) |

The parameter is determined from the average value of . In our model, the average size of crystals is given by , because of the conservation of the total mass of solute. At , the average value can be calculated as

(8) |

from the definition of the gamma function. The parameter is given by

(9) |

for , and . Figure 4(a) shows the distribution of Eq. (7) for , and . The size distributions in Figs. 2(a) and 3(a) are qualitatively close to the size distribution shown in Fig. 4(a), in that there is a longer tail in a region of . The average value in the size distribution in Fig. 3(a) is close to that in Fig. 4(a) because the number of initial crystal seeds is the same; however, the width of the size distribution in Fig. 3(a) is smaller than that of the size distribution of the Voronoi tessellation.

If is not so large, the effect of diffusion needs to be taken into consideration. To understand qualitatively the dependence of the width of the size distribution on using an even simpler model, we propose a mean-field type model for crystal growth, keeping the Voronoi tessellation in mind. The model is coupled differential equations expressed as

(10) | |||||

(11) |

where denotes the crystal size in the th Voronoi cell, is the concentration of the solute in the th Voronoi cell, , and denotes the size of the th Voronoi cell. The crystal is assumed to grow only at the seed point located at the center of the th Voronoi cell with the growth rate , which is similar to the previous model using Eqs. (1), (4), and (5). The concentration inside each Voronoi cell is assumed to be uniform and denoted as . The concentrations in different Voronoi cells are different but becomes uniform owing to the mean-field-type diffusion term in Eq. (11). The total mass is conserved in the time evolution of Eqs. (10) and (11). In our numerical simulation, is determined by the Voronoi tessellation of 180 seeds in a square lattice. The crystal size is defined as at the final stationary state of Eqs. (10) and (11). Figure 4(b) is the size distribution of at the final time at and . A rather wide size distribution close to the distribution shown in Fig. 4(a) is obtained at this parameter set satisfying . Figure 4(c) is the size distribution of at and . The size distribution is very narrow, when is sufficiently small. Figure 4(d) shows the relationship between the width of the size distribution and . Similarly to previous models, we find the linear dependence of on in the small range of and the saturation at large .

The concentration can be expanded using eigenvalues ( and the corresponding eigenfunctions of the linear equation (11) as where ’s are the expansion parameters. Substitution of this expression into Eq. (11) yields

(12) |

By the summation of Eq. (12) from to , the eigenvalue satisfies

(13) |

This is the th-order algebraic equation and there are solutions corresponding to the eigenvalues. If , is a solution to Eq. (13). For sufficiently small and small , Eq. (13) is rewritten as

using the Taylor expansion. The first eigenvalue is therefore expressed as

(14) |

where . If is sufficiently small, the first eigenvalue is small, and can be approximated at because the other terms decay quickly. Owing to Eq. (12), is proportional to , and is expressed as

(15) |

Since , for a sufficiently small . The final stationary value of is given by the time integration of as

(16) |

The average value of is and the width of the size distribution of is given by

(17) |

where is the root mean square of the distribution of . The width is evaluated as , since is evaluated as using a formula of the gamma function. The theoretical result is shown by a dashed line in Fig. 4(d). Fairly good agreement is seen in a region of sufficiently small . We have shown again that the eigenfunction for the largest eigenvalue determines the width of the size distribution for sufficiently small . For sufficiently large , the total solute in the th Voronoi cell is absorbed to the th seed at the center. This is the special situation corresponding to Fig. 4(a). The crystal size is evaluated as . Then, the width of the size distribution is expected to be , which is plotted by the dotted line in Fig. 4(d).

In this simple model, increases in proportion to and the dynamics of is determined only by . Therefore, if the initial sizes are not zero but randomly distributed, the final crystal sizes are calculated as , where is the final crystal size in the case of . That is, the final size distribution is determined by the convolution of the initial size distribution of and the distribution of .

## Iv Three-dimensional model

In this section, we study a three-dimensional coupled map lattice, which is more realistic than the two-dimensional model. The three-dimensional extension of Eq. (1) on a cubic lattice is expressed as

(18) | |||||

The three-dimensional extension of Eq. (2) is expressed as

(19) |

where denotes a lattice point in a cubic lattice, and the parameter is 1 for nonflat sites and takes a small value for flat sites. Here, the flat interface site is a site that has only one crystal site and five solution sites as its nearest-neighbor sites. The system size is set to be . Figure 5(a) shows a 3D plot of the crystal site obtained in the numerical simulation of the three-dimensional coupled map lattice at and . Initially, 11 crystal seed points are randomly distributed. Clear cubic crystals are created, and the sizes of crystals are almost uniform. The size distributions are calculated for 8 samples of initial 100 seeds of crystals. Figure 5(b) shows the size distribution at and . Figure 5(c) shows the size distribution at and . Similarly to the two-dimensional models, the size distribution becomes wider as increases. Figure 5(d) shows the relationship between the width of the size distribution and . Here, is changed as () for , and then is changed as () for . The width increases linearly with , which is similar to the results found in the two-dimensional models.

## V Summary

We have numerically studied the crystal size distribution in polynuclear growth using coupled map lattice models. The size distribution of crystals originates from the initial random configuration of crystal seeds, and the dispersion is reduced by the effect of the diffusion. We have found that the width of the size distribution is approximately determined by the parameter , where is the growth rate at the interface sites and is the diffusion constant. When is sufficiently small, the width increases linearly with and tends to saturate at large . The numerical simulations were performed on a square lattice and a cubic lattice. Monodisperse square and cubic crystals are obtained when are sufficiently small and the kinetic parameter is sufficiently small. The research and development of the manufacturing process of monodisperse crystals, especially monodisperse nanocrystals, are important for various applications. The linear dependence of on in the parameter range of small is explained by the eigenfunction for the first eigenvalue in a simplified two-dimensional model and an even more simplified mean-field-type model. For the mean-field model using the Voronoi tessellation, the slope of the linear dependence has been evaluated theoretically.

## References

- (1) Y. Saito, Statistical Physics of Crystal Growth, World Scientific (Singapore, 1996).
- (2) A. Karma and W. J. Rappel, Phys. Rev. E 53, R3017 (1996).
- (3) N. Provatas and K. Elder, Phase-Field Methods in Material Science and Engineering (Wiley-VCH, 2010).
- (4) H. Sakaguchi, J. Phys. Soc. Jpn. 67, 96 (1998).
- (5) H. Sakaguchi and M. Ohtaki, Physica A 272, 300 (1999).
- (6) H. Sakaguchi, M. Gondou, and H. Honjo, J. Phys. Soc. Jpn. 81, 124004 (2012).
- (7) R. Becker and W. Döring, Ann Physik 24, 719 (1935).
- (8) A. D. Randolph and M. A, Larson, Theory of Particulate Processes (Academic Press, New York, 1988).
- (9) B. D. Marsh, Contrib. Mineral. Petrol. 99, 277 (1988).
- (10) E. Matijević, Langmuir 10, 8 (1994).
- (11) T. Sugimoto, Adv. Colloid Interface Sci. 28, 65 (1987).
- (12) I. I. Tarhan and G. H. Watson, Phys. Rev. Lett. 76, 315 (1996).
- (13) S. B. DiCenzo and G. K. Wertheim, Phys. Rev. B 39, 6792 (1989).
- (14) M. Tanemura, Forma 18, 221 (2003).