Numerical Study of Crystal Size Distribution in Polynuclear Growth
aa r X i v : . [ c ond - m a t . m t r l - s c i ] M a y Numerical Study of Crystal Size Distribution in Polynuclear Growth
Hidetsugu Sakaguchi and Takuma Ohishi
Department of Applied Science for Electronics and Materials,Interdisciplinary Graduate School of Engineering Sciences,Kyushu University, Kasuga, Fukuoka 816-8580, Japan
The crystal size distribution in polynuclear growth is numerically studied using a coupled maplattice model. The width of the size distribution depends on c/D , where c is the growth rate atinterface sites and D is the diffusion constant. When c/D is sufficiently small, the width W increaseslinearly with c/D and saturates at large c/D . Monodisperse square and cubic crystals are obtainedrespectively on square and cubic lattices when c/D is sufficiently small for a small kinetic parameter b . The linear dependence of W on c/D in a parameter range of small c/D is explained by theeigenfunction for the first eigenvalue in a two-dimensional model and a mean-field model. For themean-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, chemicalengineering, and mineralogy. It is also studied in statistical physics as a typical nonequilibrium phenomenon. Nucle-ation, surface kinetics, and diffusion are important factors in crystal growth [1]. Crystal growth is one of the Stefanproblems for moving interfaces, which is difficult to analyze mathematically. Various simulation methods have beenproposed to study crystal growth far from equilibrium. The phase-field model is a partial-differential equation inwhich the solid and melt phases are expressed with a continuous variable corresponding to the order parameter [2, 3].We proposed a coupled map lattice model for crystal growth [4–6]. A similar type of order parameter is defined on alattice, and coupled maps are used for the time evolution. The numerical simulation of the coupled map lattice is muchfaster than that using the phase-field model. One reason is that the width of the interface between the solid and meltphases is only one lattice constant in the coupled map lattice model. On the other hand, the order parameter changessmoothly 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) usingcoupled 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 studiedcrystal growth starting from a single crystal seed. In this paper, we study polynuclear growth in which there aremany 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¨oring studied the time evolution of the size distribution in the nucleation process [7]. Randolph andLarson discussed stationary CSDs in the research field of chemical engineering [8]. Marsh applied the CSD theoryin geology and found an exponential size distribution of mineral crystals deposited from magma [9]. In their model,the exponential distribution is derived by the balance of the outflow, the inflow, and the growth of crystals. In thesynthesis of micro- and nanocrystals, the final size distribution after the crystal growth is important. Various methodsof characterization and control of the crystal size distribution have been studied. In particular, the control of thecrystal size distribution have been investigated in colloidal sciences [10, 11]. The size distributions similar to thegamma distribution or the log normal distribution are often observed in experiments. The size distribution P ( S )is not symmetric around the average value but has a longer tail for large S . It is practically important to producealmost 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 microcolloidalparticles, which can be applied to various systems such as photonic crystals [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 monodispersecrystals, it is important that the time scales of nucleation and growth processes are separated, and the size distributionof the crystal seeds is almost mono-disperse at the initial stage of the growth process. In this paper, we consider aproblem 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 thecompetition of the diffusion and surface kinetics. Our model might be applied to a two-dimensional crystal growthon a substrate in that the positions of crystals are fixed in time. However, the purpose of this study is to understandqualitatively a condition wherein the final crystal size distribution becomes almost monodisperse. (a) (b) j i 050100 0 50 100 j i0(cid:13)20(cid:13)40(cid:13)60(cid:13)80(cid:13)100(cid:13)0(cid:13) 20(cid:13) 40(cid:13) 60(cid:13) 80(cid:13) 100(cid:13) j (cid:13) i(cid:13) (c) j (cid:13) i(cid:13) (d) FIG. 1: Crystal sites at (a) D = 0 .
01 and c = 1, and (b) D = 0 . c = 0 .
01 in the case that the initial sizes of crystals areall 1. Crystal sites at (c) D = 0 .
01 and c = 1, and (d) D = 0 . c = 0 .
01 in the case that the initial sizes of crystals are 1for i <
70 and 9 for i >
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 ourcoupled map lattice for the solution growth, namely, the order parameter x n ( i, j ) and the concentration u n ( i, j ), where( i, j ) and n denote respectively a lattice site and a discrete time. The order parameter x n ( i, j ) 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 andsurface kinetics. The diffusion is expressed as u ′ n ( i, j ) = u n ( i, j ) + D { u n ( i + 1 , j ) + u n ( i − , j ) + u n ( i, j + 1) + u n ( i, j − − u n ( i, j ) } , (1)where D is a diffusion constant. The surface kinetics is expressed as x n +1 ( i, j ) = x n ( i, j ) + cbu ′ n ( i, j ) ,u n +1 ( i, j ) = u ′ n ( i, j ) − cbu ′ n ( i, j ) , (2)where c is a parameter for the growth rate. The time evolution Eq. (2) is applied only at interface sites satisfying0 < x n ( i, j ) <
1. If x n ( i, j ) goes over a critical value 1, interface sites change into crystal sites and the neighboringsolute sites are assigned to interface sites. By repeating Eqs. (1) and (2), polynuclear crystal growth occurs. Theparameter b is 1 at nonflat sites and takes a small value b at flat sites. The flatness of an interface site is determinedby counting the number of crystal sites in the four nearest-neighbor sites of the interface site. If the number is 1, b is set to be b , and if the number is larger than 2, b is set to be 1. For sufficiently small b , square crystals appear,because kink sites and concave sites are quickly occupied, but it takes a rather long time for flat sites to grow. Theparameter b represents an effect of the surface tension on the surface kinetics. In the case of b = 1, the surfacetension effect disappears. In Eq. (2), the increase in the order parameter x is equal to the absorption of u at interfacesites, that is, the conservation law of x + u is satisfied. In our model of Eqs. (1) and (2), only the deposition processfrom 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 ofsize 1 are randomly set in a square of L × L = 100 × u ( i, j ) is uniform and set to be u = 0 . b is assumed to be 0.005. Periodic boundary conditions are imposed. After a long-time iterationof Eqs. (1) and (2), u becomes zero owing to the diffusion and the absorption of u at interface sites, and the crystalgrowth process is completed. Figure 1(a) shows crystal sites at D = 0 .
01 and c = 1 at the final time. The crystalsize is widely distributed. The crystal sizes are small in a region where initial seeds are densely distributed. In ourmodel, the randomness originates only from the initial configuration of crystal seeds. Figure 1(b) shows crystal sitesat D = 0 . c = 0 .
01. Crystals take a clear square shape, and the sizes are almost uniform. As D is large and c issmall, the concentration u tends to be uniform and the crystal sizes become homogeneous even for the random initial (a) (b) (c) P ( S ) S 5101520 0 0.1 0.2 0.3 0.4 0.5 W c 1101000.1 1 10 100 W c/D (d) P ( S ) (cid:13) S(cid:13)
FIG. 2: Size distributions at (a) D = 0 .
01 and c = 1, and (b) D = 0 . c = 0 .
01. (c) Relationship between W and c for c/D = 5 (rhombi) and c/D = 2 . W and c/D . (a) (b) (c) P ( S ) S00.0050.010.0150.02 0 50 100 150 200 250 P ( S ) S 0.1(cid:13)1(cid:13)10(cid:13)100(cid:13) 0.01(cid:13) 0.1(cid:13) 1(cid:13) 10(cid:13) W (cid:13) c/D(cid:13) FIG. 3: Size distributions at (a) D = 0 .
025 and c = 0 .
64, and (b) D = 0 . c = 0 . W and c/D . configuration of seeds. To show an effect of the initial seed size distribution, the sizes of four crystal seeds at i > × D = 0 .
01 and c = 1 and (d) D = 0 . c = 0 .
01 at a final time. At D = 0 . c = 0 .
01, the crystal sizes for thesix 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 largerby 11 − = 40 than those in Fig. 1(b). At D = 0 .
01 and c = 1, the crystal sizes change for all crystals from thevalues in Fig. 1(a), but the change in size from Fig. 1(a) depends on the position. The final crystal size depends onthe initial seeds, but the dependence of the final size distribution on the initial size distribution is not so simple. Weinvestigate only the case that the initial sizes are all 1 hereafter.We have calculated the CSD using a larger square lattice of 400 × S is defined as the sum of x n ( i, j )in crystal sites and interface sites around each crystal seed. At crystal sites, x n ( i, j ) takes a value of 1. Then, S represents the area of a two-dimensional crystal. Figures 2(a) and 2(b) show the size distribution at (a) D = 0 . c = 1, and (b) D = 0 . c = 0 .
01. (There are a few crystals that collide with neighboring crystals, but thesizes of those crystals are not considered in the size distribution. The collision probability decreases if the system sizeincreases for a fixed number of crystal seeds. In the production of colloidal particles, the electric double layer aroundcolloidal particles prevents the collision.) The size distribution is rather wide at D = 0 .
01 and c = 1. The distributionis asymmetric and has a longer tail for large S . The size distribution is very narrow at D = 0 . c = 0 .
01, whichimplies that crystals are almost monodisperse.The average of crystal size is estimated as u L /N , because the total sum of S is equal to L u owing to theconservation law of u + x . The inhomogeneity of the crystal size is characterized by the width W or the standarddeviation of the size distribution. Figure 2(c) shows the widths W for four different c ’s at the same c/D = 2 . W depends strongly on c/D but hardly depends on c under the condition of c/D =const. If c and D are small and the one-step increments of x and u in Eqs. (1) and (2) are sufficiently small, Eqs. (1) and (2) can beapproximated by differential equations. It can be shown that the final stationary states depend only on the ratio c/D owing to a scale transformation of time if the differential approximation is good. Figure 2(d) shows a relationshipbetween W and the ratio c/D in a logarithmic plot. The diffusion constant D is changed as 0 . , . , . , . . c = 1, and then c is changed as 1 / n ( n = 1 , , · · · ,
6) for D = 0 .
2. The width increases linearly as W ∼ . c/D ) in a range of small c/D and tends to saturate at large c/D . That is, monodisperse crystals areobtained in the limit of c/D = 0, and the dispersion is proportional to c/D when c/D is sufficiently small.In the coupled map lattice model of Eqs. (1) and (2), there is a possibility that neighboring crystals become closerto each other and collide with each other as they grow. To neglect the effect of the collision, we can propose a simplercoupled map lattice model. Equation (2) is rewritten as x n +1 ( i, j ) = x n ( i, j ) + c ′ u ′ n ( i, j ) /r n ( i, j ) × πr n ( i, j ) ,u n +1 ( i, j ) = u n +1 / ( i, j ) − c ′ u ′ n ( i, j ) /r n ( i, j ) × πr n ( i, j ) , (3)where c ′ = cb/ (2 π ), and x n ( i, j ) and r n ( i, j ) = p x n ( i, j ) /π are interpreted respectively as the area and radius ofcircular crystals. Here, the growth velocity of crystals of radius r n ( i, j ) is assumed to be inversely proportional to theradius r n ( i, j ) as v r = c ′ /r . This type of relation is satisfied under the condition of diffusion-limited growth. Equation(3) is expressed as x n +1 ( i, j ) = x n ( i, j ) + cu ′ n ( i, j ) , (4) u n +1 ( i, j ) = u ′ n ( i, j ) − cu ′ n ( i, j ) , (5)where 2 πc ′ is rewritten as c . The form of Eqs. (4) and (5) is almost equal to that of Eq. (2); however, the timeevolution of Eqs. (4) and (5) is applied only at the positions of crystal seeds in this model. That is, in contrastto Eq. (2), the variable x increases indefinitely even if x becomes larger than 1. Growth patterns such as squarepatterns do not appear, because x grows only at the seed points. The magnitude of x is interpreted as crystal size.In other words, we consider a situation that the radii of crystals are sufficiently small in contrast to the distancesamong different crystals, and the growth patterns are invisible. In the following numerical simulation, the number N of initial seeds is set to be 180, and the system size is L × L = 400 × n , when u n ( i, j ) decreases to 0. Figure 3(a) shows the crystal size distribution at D = 0 .
025 and c = 0 .
64. Figure3(b) shows the size distribution at D = 0 . c = 0 . c is small and D is large, the size distributionbecomes narrow. The rhombi in Fig. 3(c) show the relationship between the width W of the size distribution andthe ratio c/D , where c is changed as 5 . / n ( n = 1 , , · · · ,
12) for D = 0 .
2. The dashed line denotes 74( c/D ). Thewidth W increases from 0 linearly in a parameter region of small c/D and tends to saturate at large c/D also in thissimplified model.Because Eqs. (1) and (5) are linear equations, an initial value problem for u can be solved using the eigenvalues andeigenfunctions as u n ( i, j ) = P k a k e k ( i, j )( λ k ) n , where λ k and e k ( i, j ) are eigenvalues and eigenfunctions, and a k is theexpansion coefficient. The first eigenvalue λ and the corresponding eigenfunction e ( i, j ) can be numerically evaluatedby a long-time iteration of the same equation Eqs. (1) and (5) and the normalization u n +1 ( i, j ) → u n +1 ( i, j ) /λ . When c/D is sufficiently small, the first eigenvalue is close to 1 and the first component of a e ( i, j )( λ ) n decays slowly orthe other components P k =1 a k e k ( i, j )( λ k ) n decay quickly. Under this condition, the growth rate of each crystal isexpected to be proportional to the first eigenfunction e ( i, j ) at the crystal-seed point, because the concentration u decays as e ( i, j ) λ n and x ( i, j ) increases to compensate for the decrease in u . The crystal-size distribution can beevaluated approximately using numerically estimated e ( i, j ) by assuming that the crystal size is proportional to theeigenfunction e ( i, j ) and using the fact that the average crystal size is equal to L u /N . The pluses in Fig. 3(c) showthe width of the size distribution evaluated from e ( i, j ). The width is close to the results denoted by rhombi obtainedby direct numerical simulation of Eqs. (1), (4), and (5) for small values of c/D . At large c/D , the contribution fromthe 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, N centralpoints are randomly distributed, and the whole plane is partitioned into territories of each point. The boundary linebetween two territories is determined by the perpendicular bisector of the two central points. Each Voronoi cell is aset of points surrounded by perpendicular bisectors. The Voronoi tesselation is applied to various areas such as grainsof crystals and territories of animals. The size distribution of Voronoi cells is approximately given by the gammadistribution: [13, 14] P ( S v ) ∝ S γv e − αS v , (6)where S v is the area of each Voronoi cell, α is a parameter, and γ ∼ . c/D is sufficiently large. Theabsorption range or the territory of each crystal seed might be approximated by the Voronoi cell. We interpret thatthe central point in each Voronoi cell corresponds to a seed point of a crystal. First, we consider a case that all P ( S ) S (a) P ( S ) S (b) W c/D (c) (d) P ( S ) S FIG. 4: (a) Size distribution determined by the Voronoi tessellation where α is given by Eq. (9). (b) Size distribution determinedusing Eqs. (1), (10), and (11) at D = 0 .
01 and c = 1. (c) Size distribution determined using the mean-field model at D = 1and c = 0 .
02. (d) Relationship between the width W and c/D . the solutes in each Voronoi cell are absorbed into the central seed of a crystal. In that case, the crystal size S isproportional to the area S v of the Voronoi cell. In that case, the crystal-size distribution is expressed as P ( S ) ∝ S γ e − αS . (7)The parameter α is determined from the average value of S . In our model, the average size of crystals is given by u L /N , because of the conservation of the total mass of solute. At γ = 2 .
5, the average value can be calculated as h S i = R ∞ S γ +1 e − αS dS R ∞ S γ e − αS dS = 1 α Γ( γ + 1)Γ( γ ) = 72 α , (8)from the definition of the gamma function. The parameter α is given by α = 7 N u L ∼ . N = 180 , u = 0 .
1, and L = 400. Figure 4(a) shows the distribution of Eq. (7) for N = 180 , u = 0 .
1, and L = 400.The size distributions in Figs. 2(a) and 3(a) are qualitatively close to the size distribution shown in Fig. 4(a), in thatthere is a longer tail in a region of S > h S i . The average value in the size distribution in Fig. 3(a) is close to thatin Fig. 4(a) because the number of initial crystal seeds is the same; however, the width of the size distribution inFig. 3(a) is smaller than that of the size distribution of the Voronoi tessellation.If c/D is not so large, the effect of diffusion needs to be taken into consideration. To understand qualitativelythe dependence of the width of the size distribution on c/D using an even simpler model, we propose a mean-fieldtype model for crystal growth, keeping the Voronoi tessellation in mind. The model is coupled differential equationsexpressed as dx i dt = cu i , (10) du i dt = 1 S vi {− cu i + D (¯ u − u i ) } , (11)where x i denotes the crystal size in the i th Voronoi cell, u i is the concentration of the solute in the i th Voronoi cell,¯ u = (1 /N ) P u i , and S vi denotes the size of the i th Voronoi cell. The crystal is assumed to grow only at the seedpoint located at the center of the i th Voronoi cell with the growth rate cu i , which is similar to the previous modelusing Eqs. (1), (4), and (5). The concentration inside each Voronoi cell is assumed to be uniform and denoted as u i . The concentrations u i in different Voronoi cells are different but becomes uniform owing to the mean-field-typediffusion term D (¯ u − u i ) in Eq. (11). The total mass P i ( x i + u i S vi ) is conserved in the time evolution of Eqs. (10)and (11). In our numerical simulation, S vi is determined by the Voronoi tessellation of 180 seeds in a 400 ×
400 squarelattice. The crystal size S i is defined as x i at the final stationary state of Eqs. (10) and (11). Figure 4(b) is the sizedistribution of S i at the final time at c = 1 and D = 0 .
01. A rather wide size distribution close to the distributionshown in Fig. 4(a) is obtained at this parameter set satisfying c/D = 100. Figure 4(c) is the size distribution of S i at c = 0 .
02 and D = 1. The size distribution is very narrow, when c/D is sufficiently small. Figure 4(d) shows therelationship between the width W of the size distribution and c/D . Similarly to previous models, we find the lineardependence of W on c/D in the small range of c/D and the saturation at large c/D .The concentration u i can be expanded using eigenvalues λ k ( k = 1 , , · · · , N ) and the corresponding eigenfunctions e k,i of the linear equation (11) as u i ( t ) = P Nk =1 a k e k,i e − λ k t where a k ’s are the expansion parameters. Substitution ofthis expression into Eq. (11) yields e k,i = D P Nj =1 e k,j N ( c − λ k S vi + D ) . (12)By the summation of Eq. (12) from i = 1 to N , the eigenvalue λ satisfies1 N N X j =1 Dc + D − λS vj = 1 . (13)This is the N th-order algebraic equation and there are N solutions corresponding to the N eigenvalues. If c = 0, λ = 0 is a solution to Eq. (13). For sufficiently small c and small λ , Eq. (13) is rewritten as1 − cD + λD N N X j =1 S vj = 1 , using the Taylor expansion. The first eigenvalue λ is therefore expressed as λ = c ¯ S v , (14)where ¯ S v = (1 /N ) P Nj =1 S vj . If c is sufficiently small, the first eigenvalue λ is small, and u i can be approximated at u i = a e ,i e − λ t because the other terms decay quickly. Owing to Eq. (12), e ,i is proportional to 1 / ( c − λ S vi + D ),and u i ( t ) is expressed as u i ( t ) = u i c + D − cS vi / ¯ S v e − λ t . (15)Since u (0) = u , u i = u ( c + D − cS vi / ¯ S v ) ∼ u D for a sufficiently small c . The final stationary value of x i is givenby the time integration of cu i ( t ) as S i = Z ∞ cu i c + D − cS vi / ¯ S v e − λ t ∼ ¯ S v u Dc + D − cS vi / ¯ S v ∼ ¯ S v u (cid:18) cD S vi − ¯ S v ¯ S v (cid:19) . (16)The average value of S i is ¯ S v u and the width W of the size distribution of S i is given by W = cD δS v u , (17)where δS v is the root mean square of the distribution of S vi . The width W is evaluated as W = 47 . c/D ), since δS v u is evaluated as p / L /N ) u = 47 . c/D . We have shownagain that the eigenfunction for the largest eigenvalue determines the width of the size distribution for sufficientlysmall c/D . For sufficiently large c/D , the total solute u S vi in the i th Voronoi cell is absorbed to the i th seed at thecenter. This is the special situation corresponding to Fig. 4(a). The crystal size is evaluated as S vi u . Then, thewidth of the size distribution is expected to be W = δSu = 47 .
5, which is plotted by the dotted line in Fig. 4(d).In this simple model, x i increases in proportion to u i and the dynamics of u i is determined only by { u j } . Therefore,if the initial sizes x i (0) are not zero but randomly distributed, the final crystal sizes x if are calculated as x if = x if + x i (0), where x if is the final crystal size in the case of x i (0) = 0. That is, the final size distribution isdetermined by the convolution of the initial size distribution of x i (0) and the distribution of x if . IV. THREE-DIMENSIONAL MODEL
In this section, we study a three-dimensional coupled map lattice, which is more realistic than the two-dimensionalmodel. The three-dimensional extension of Eq. (1) on a cubic lattice is expressed as
0 20 40 60 0 20 40 60 0 20 40 60 ijk (a) (b) P ( S ) S (c) (d) W c/D0(cid:13)0.05(cid:13)0.1(cid:13)0.15(cid:13)0(cid:13) 100(cid:13) 200(cid:13) 300(cid:13) 400(cid:13) 500(cid:13) 600(cid:13) P ( S ) (cid:13) S(cid:13)
FIG. 5: (a) 3D plot of the crystalized site at D = 0 .
15 and c = 0 . D = 0 . c = 1. (c) Size distribution at D = 0 .
15 and c = 0 . W and c/D for c = 1. u ′ n ( i, j, k ) = u n ( i, j, k ) + D { u n ( i + 1 , j, k ) + u n ( i − , j, k ) + u n ( i, j + 1 , k ) + u n ( i, j − , k )+ u n ( i, j, k + 1) + u n ( i, j, k − − u n ( i, j, k ) } . (18)The three-dimensional extension of Eq. (2) is expressed as x n +1 ( i, j, k ) = x n ( i, j, k ) + cbu ′ n ( i, j, k ) ,u n +1 ( i, j, k ) = u ′ n ( i, j, k ) − cbu ′ n ( i, j, k ) , (19)where ( i, j, k ) denotes a lattice point in a cubic lattice, and the parameter b is 1 for nonflat sites and takes a smallvalue b = 0 .
005 for flat sites. Here, the flat interface site is a site that has only one crystal site and five solution sitesas its nearest-neighbor sites. The system size is set to be 60 × ×
60. Figure 5(a) shows a 3D plot of the crystal siteobtained in the numerical simulation of the three-dimensional coupled map lattice at D = 0 .
15 and c = 0 .
1. Initially,11 crystal seed points are randomly distributed. Clear cubic crystals are created, and the sizes of crystals are almostuniform. The size distributions are calculated for 8 samples of initial 100 seeds of crystals. Figure 5(b) shows the sizedistribution at D = 0 . c = 1. Figure 5(c) shows the size distribution at D = 0 .
15 and c = 0 . c/D increases. Figure 5(d) shows the relationshipbetween the width W of the size distribution and c/D . Here, D is changed as D = 0 . / n ( n = 1 , , · · · ,
5) for c = 1,and then c is changed as 1 / n ( n = 1 , , · · · ,
5) for D = 0 .
15. The width increases linearly with c/D , which is similarto 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 isreduced by the effect of the diffusion. We have found that the width of the size distribution is approximately determinedby the parameter c/D , where c is the growth rate at the interface sites and D is the diffusion constant. When c/D is sufficiently small, the width W increases linearly with c/D and tends to saturate at large c/D . The numericalsimulations were performed on a square lattice and a cubic lattice. Monodisperse square and cubic crystals are obtainedwhen c/D are sufficiently small and the kinetic parameter b is sufficiently small. The research and development ofthe manufacturing process of monodisperse crystals, especially monodisperse nanocrystals, are important for variousapplications. The linear dependence of W on c/D in the parameter range of small c/D is explained by the eigenfunctionfor the first eigenvalue in a simplified two-dimensional model and an even more simplified mean-field-type model. Forthe mean-field model using the Voronoi tessellation, the slope of the linear dependence has been evaluated theoretically. [1] Y. Saito, Statistical Physics of Crystal Growth , World Scientific (Singapore, 1996).[2] A. Karma and W. J. Rappel, Phys. Rev. E , 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. , 96 (1998). [5] H. Sakaguchi and M. Ohtaki, Physica A , 300 (1999).[6] H. Sakaguchi, M. Gondou, and H. Honjo, J. Phys. Soc. Jpn. , 124004 (2012).[7] R. Becker and W. D¨oring, Ann Physik , 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. , 277 (1988).[10] E. Matijevi´c, Langmuir , 8 (1994).[11] T. Sugimoto, Adv. Colloid Interface Sci. , 65 (1987).[12] I. I. Tarhan and G. H. Watson, Phys. Rev. Lett. , 315 (1996).[13] S. B. DiCenzo and G. K. Wertheim, Phys. Rev. B 39, 6792 (1989).[14] M. Tanemura, Forma18