Simulation of 2D Nucleation and Crystal Growth BY GEORGE H. GILMER Bell Laboratories Murray Hill New Jersey 07974 U.S.A. Received 1st September 1977 The kinetic Ising model of a crystal-fluid interface is investigated by the Monte Carlo method. The growth rate transient is calculated in the case of a perfect crystal and good agreement with experi- mental data is obtained for the first time. The steady-state growth rates are compared with values calculated using nucleation theory. Large growth rates at high temperatures are related to the small step free energies. The effect of surface mobility is considered. The growth rates obtained for low-index faces of the simple cubic (s.c.) and face centred cubic (f.c.c.) lattice correlate well with the a-factors for these faces.Computer models have been applied to crystal growth and dissolution processes for a number of years. Complex situations that proved difficult to treat analytically could often be modelled by numerical techniques. Early applications included the dynamics of ledge formation by a one-dimensional array of parallel steps.l Two-dimensional (2D) arrays were used to represent the complex surface structures that arise from the nucleation and spreading of 2D clusters on a growing crystal.2 Later lattice models were developed that included the atomic processes of con-densation evaporation and migrati~n.~-~ These models include the long range correlations that are basic to such phenomena as nucleation and second order phase transtions. Here a lattice model is applied to the deposition of atoms of a single component onto a low-index face of a crystal.THE MODEL The Ising (or lattice gas) model provides an atomic-scale representation of the growing crystal yet it is simple enough to permit rapid computation. A fixed array of lattice sites is assumed; each site may either be occupied by a single atom or be vacant. Surface configurations are represented by a square array of integers which specifies the numbers of atoms in columns of sites perpendicular to the plane of the interface. Vacancies and overhanging atomic configurations are not permitted; hence this is a restricted version of the Ising model. Crystal growth kinetics are simulated through the exchange of atoms between the crystal and the adjacent fluid.Atoms impinge at randomly selected sites in the square array at a rate k+ = v exp up) where v is the frequency factor is the reci- procal of Boltzmann’s factor multiplied by the absolute temperature and ,u is the chemical potential. The dissolution (or evaporation) of surface atoms occurs at a rate that depends on the coordination k; = v exp (-Pnp) where n is the number of nearest neighbours and q is the energy of the bond between a pair of such neigh- bours. The frequency factor v includes any retardation due to the formation of ac- tivated complexes during the transfer of atoms between the fluid and sites on the crystal surface. Equilibrium is achieved with ,u = -2q/2 where z is the coordination num- ber of the lattice; this assures that impingement and evaporation at a kink site occur SIMULATION OF 2D NUCLEATION AND CRYSTAL GROWTH at equal rates.The simulation proceeds by the generation of a chain of deposition and evaporation events with the calculated frequencies. Surface migration can also be included; an atom on the surface is then allowed the additional option of a move to an adjacent lattice site. The model is capable of exhibiting a large number of differ- ent growth mechanisms since the fundamental nature of its operation affords con- siderable flexibility. However the model does have some important limitations. For example the lattice structure must be chosen apriori; hence this model can not be used to investigate the formation of dislocations or extended lattice defects.However if a lattice containing defects is chosen initially their effect on the kinetics can be measured. INITIAL TRANSIENT Low index faces of perfect crystals are often quite flat and smooth with only iso- lated adatoms and vacancies. This structure is expected at temperatures well below the melting point if the crystal is in equilibrium with the surrounding fluid. If a crystal bounded by low-index surfaces is subjected to a large driving force for crystal- lization rapid changes in the surface structures occur. The adatom concentration rises very quickly to a new level and on a somewhat slower time scale stable clusters of atoms are nucleated. These clusters expand and eventually entire layers are deposited. This process is illustrated in fig. 1.Fig. l(a) shows a typical structure after the deposition of a quarter of a monolayer. Even at this early stage some stable clusters have nucleated in the second layer. Fig. l(b) is a drawing of the same system after several monolayers have been deposited. Holes are present as a result of the incomplete mergers of large clusters at the lower levels and small clusters predominate at the upper levels. Fig. l(c) illustrates the effect of a larger driving force. The rate of nucleation generally increases more rapidly with driving force than does the edge velocity of a stable cluster and this causes the higher cluster density and rougher appearance of this surface. Fig. 2 is a plot of the average growth rates measured on ninety different 60 x 60 sections of S.C.(100) faces. Each section was subjected to a sudden application of the driving force. The open triangles represent the growth rates simulated without surface migration; the temperature and driving force correspond to those of fig. l(c). The open squares represent growth rates obtained in the presence of surface migration where the migration rate to sites of equal coordination is 7.4 times the evaporation rate. At the instant the driving force is applied there is a very large surge that persists only until the new adatom concentration is established. This surge is reflected in the first data point of fig. 2 which appears somewhat higher than anticipated. (The corre- sponding relaxation time is the adatom evaporation rate ki = 4.1 x 10-3k+ whereas the deposition rates are averaged over the interval At between data points; At = 0.34k+ for the upper curve.) The growth rates exhibit damped oscillations around the asymptotic rate R repre-sented by the dashed lines.The thickness of the deposit in monolayers is indicated above the curves. The minima correspond roughly to the points where a layer is complete and only small clusters occupy the next level. The amplitude of the oscillations is a measure of the correlation in the surface heights at different sites in the array. At the beginning the majority of the sites are at the same level but later they are distributed over a range of levels as a result of statistical variations in the nucleation rates at the various locations. Sites in close proximity remain highly corre- GEORGE H.GILMER FIG.1.-Computer drawings of the surfaces of perfect crystals during growth by the Monte Carlo method. (a)was calculated with v = 4/pand Ap = 2/p. (b) illustrates the same crystal at a later stage of growth. In (c) cp = 4/p and Ap = 2.51b. SIMULATION OF 2D NUCLEATION AND CRYSTAL GROWTH lated in height since a cluster which nucleates near one site quickly spreads to the other. But sites that are widely separated in the lateral dimension may eventually have large differences in height7 In this case there should be little variation in the growth rate since the different portions of the surface are experiencing the fast and slow growth at different times. When the transient decays and the growth rate approaches R,the mean squared height deviations between such sites should be large and probably exceed one layer spacing.k‘t 3.15 -. =% =L 0.10 I I I 20 40 60 k‘t FIG.2.-The transient growth rates normalized by the impingement rate multiplied by the layer spacing. The squares represent data calculated with mobile surface atoms and the figures at the left and above apply to this case. The triangles represent growth with immobile atoms and the figures on the right and below apply here. Note that the growth rate scale does not extend to zero in this case. As in fig. 1 p = 4/p Ay = 2.5IP. 0kJk-= 7.4 A k,lk-= 0. Note the larger amplitude and greater persistence of the oscillations in the presence of mobile surface atoms. This is a consequence of the increase in the capture region near the cluster edges.Atonis that impinge within a distance of -3 atomic diameters of the edge of a cluster have a good chance of migrating to the edge and of being cap- tured. The adatom concentration in this region is depleted. The surface heights of different sites remain correlated even after several layers have been deposited since nucleation is suppressed on top of the smaller clusters. Many of the atonis that impinge on top of a cluster migrate to the edge and are captured at the lower level. The nucleation of clusters in the second layer for example generally occurs at a much higher coverage than was observed without surface migration [fig. 1(a)]. Bertocci calculated transients with oscillations of large amplitude using a simulation model in which nucleation is excluded near the edges of cl~isters.~ Surface migration also causes a more rapid increase in the deposition rate at the start of the process The average rate during the deposition of the first monolayer is (0.69 & 0.02)R without migration and (0.81 & 0.02)R with migration.The transients of fig. 2 are consistent with experimental observations. Darbing-haus and Meyer observed large amplitude oscillations in the flux of atoms evaporat- ing from a KCI crystal when an impinging beam was suddenly removed.8 (Evapora-tion or dissolution should display phenomena similar to those observed during growth since the nucleation of negative clusters or holes is necessary for the removal of suc-cessive layers of the crystal.) Surface mobility is expected to be quite large in this GEORGE H.GILMER case. Rather large amplitudes were observed in electrochemical systems although a slightly smaller mobility would be expected because of the dense fluid pha~e.~.~ The data of ref. (9) were obtained with crystals free of screw dislocations; these should be directly comparable with the simulations. The growth rates in fig. 2 approach the asymptotic values quite rapidly. For example the average growth rate during the deposition of the third monolayer is -0.97R without migration and -0.99R with migration. The maxima located at about seven-tenths of a monolayer are greater than the asymptotic rate. This is in agreement with the experimental measurements mentioned above.The more gradual increase derived froin previous model^^*^^ is probably a consequence of the idealized cluster shapes that were assumed. The clusters of fig. 1 have a much longer periphery than the squares or circles of the earlier models. As a result a uniform coverage of the surface with cluster edges is accomplished with fewer clusters and a shorter transi- ent. Finally it should be mentioned that a rough initial crystal surface does not exhibit the oscillatory transient. Later we consider certain crystal faces that have weak bond networks within the surface layer of atoms. These surfaces may have a disordered multi-level structure even while they are in equilibrium with the fluid. In fact faces with strong bond networks may also disorder in this way provided the surface roughening temperature can be exceeded.ll Non-oscillatory transients were com- puted above the roughening point of the S.C.(100) face using a pair approximation of the Ising rnode1.l2 STEADY-STATE GROWTH RATES The growth rates during the transient period can provide valuable information about the basic crystal growth process but these measurements are often not feasible. More commonly available in the literature are growth rates averaged over a large number of layers and measured at different values of the driving force. The functional relation between these two variables can also yield useful information. Computer simulation measurements of the asymptotic rate R are shown in fig. 3.13 The effect of the transient was minimized by depositing ten to twenty layers at each value of Ap and omitting the first two layers from the average.The numbers ad- jacent to the curves are the values of L/kT where L = zy/2 is the binding energy of the crystal. The normalized growth rates R/k+d were calculated without surface migration. The product k+dis the rate of growth that would occur if every impinging atom remained with the growing crystal and therefore R/k+d is also the condensa- tion coefficient and has a maximum value of unity. According to the data of fig. 3 the normalized growth rate increases rapidly with temperature. The absolute rate R should be even more sensitive to temperature since the parameter k+ includes a temperature dependent exponential in most cases.Two different growth mechanisms are evident in fig. 3. The roughening transition occurs at L/kT = 5.0 and only the data with L/kT= 4.5 depict the normal growth pro- cess. This surface is disordered and is distributed over a large number of levels even in equilibrium. Impinging atoms can easily be incorporated at the edges of existing clusters. Only in this case does the growth rate increase linearly with Ap at the origin. At the lower temperatures a coherent surface is present in equilibrium and new layers must be initiated by a nucleation event. The transition from normal growth to nucleation kinetics can be shown to take place over a very small range of ternperat~res.~~ This transition is a manifestation of surface roughening. Crystal growth kinetics at low temperature should be described by the theory of SIMULATION OF 2D NUCLEATION AND CRYSTAL GROWTH growth by the 2D nucleation and spreading of clusters.In every case a number of different nuclei are formed in each layer and the “large-crystal ” model is ~a1id.l~ Then the growth rate is related to the nucleation rate J (per unit area) and the edge velocity u of stable R = d(~J~~/3)~1~. (1) Here 21 is assumed to be independent of the cluster size. The nucleation rate J may be evaluated over most of the range of Ap by means of an atomistic expression,16 since the numbers of atoms in the critical clusters are small at these relatively high driving forces. For simplicity we include only the lowest energy clusters in each size class.We calculate the rate of formation of adatoms dimers trimers etc. from the rate equations governing the addition and removal of 3.4 P -t \ Q 0.2 4/kT FIG.3.-Normalized growth rates of perfect S.C. (100) faces. The dashed lines represent theatomistic nucleation model. single atoms from these clusters. [See ref. (13) for a list of the clusters and their con- centrations in equilibrium.] Since k+ and k; are exactly known for each simulation all parameters necessary for a calculation of J are available. The dashed lines in fig. 3 are the growth rates predicted by eqn (1) with the values of J obtained from the atomistic calculation outlined above,13 and v calculated from an expression derived by Temkin.17 The simulation results are in good agreement with the theory at LIKT = 12 but at the higher temperatures the theory exaggerates the nucleation depression.This theory does not predict a finite slope at the origin at any temperature. Clusters other than the lowest energy configurations are present in significant numbers at the higher temperatures ; apparently these must be included in the nucleation theory. Instead of attempting to enumerate all of the different cluster configurations we take a different approach. The rapid nucleation and growth at high temperature can be qualitatively explained using classical nucleation theory. The nucleation rate is related to the edge free energy of clusters by an expression of the form J = k+(Ap/kT)’i’exp (-4F2/kTAji) (2) GEORGE H. GILMER <001> step free energy 1.2 '" 021 0 1IIII \\\ /R -I '---.-! -.0 0.4 0.8 1.2 1.6 2.0 2.4 2. k7/j FIG.4.-The excess free energy of a (001 > step on a S.C. (100) face calculated by the mean field and pair approximation methods. Onsager's expression for the free energy of the 2D Ising model inter- face is included. where F is the free energy of a segment of the edge equal to an atomic diameter.18 Normal growth kinetics can only occur when the edge free energy of the clusters is zero. The excess free energy associated with an isolated (001) step on a S.C. (100) face can be calculated using the mean field and pair approximations.19 The step excess F is defined as the difference between the free energy of a surface containing a single step and that of a (100) surface.The results are shown in fig. 4 where Fis plotted as a function of the temperature. Both are normalized by the energy of a broken bond j = q/2. For comparison the free energy of a (001) interface between the phases in a 2D Ising model is also displayed.20 The edge free energy is approximately equal to the <OW step free energy kT/ j FIG.S.-The excess free energy of a 45" step. The free energy of the interface in the 2D Ising model at this angle was obtained in ref. (21). 66 SIMULATION OF 2D NUCLEATION AND CRYSTAL GROWTH zero Kelvin value (F -j)over a finite range of low temperatures. The kinetic data of fig. 3 at L/kT = 12 (kT = 0.5j) correspond to the upper limit of this range and here F -0.9j.At intermediate temperatures kink sites and other defects appear in significant quantities and cause a drastic reduction in the free energy. This accounts for the fast simulated rates at small values of the driving force in fig. 3 for LIkT = 9 and L/kT = 6. In the regime where kT >1.2j the calculated edge free energy is small but finite. The more accurate pair approximation yields much smaller free energies than the mean field method. This suggests that higher order approximations may converge to zero above T, in agreement with the kinetic data of fig. 3. A cluster may be bounded by step segments of different orientations; variations in the edge free energy are expected. Fig. 5 illustrates the free energy of a step at a 45" angle to the close-packed (001) direction.The free energy at low temperatures is appreciably higher than that of the (001) step since kink sites are present even at T = 0. However when the edge begins to disorder the difference between the two orientations diminishes and at high temperatures the free energy is essentially isotropic. Again these results indicate that the exact Ising model has a unique roughening temperature where the free energy of any step on the (100) face vanishes. GROWTH RATES OF AN F.C.C. MODEL The sensitivity of growth kinetics to the free energy of steps also implies that differ- ent faces of the same crystal may have very different growth rates. A weak bond net- work connecting atoms in a surface layer implies low free energies of steps formed on this surface.Jackson 22 has suggested that corresponding states of different faces may scale with the dimensionless temperature parameter a =c(L/kT),where 6 is the fraction of the zbulk nearest neighbours that are in the same layer as the atom. Here we apply this criterion to low-index faces of the f.c.c. lattice. The lattice model of the f.c.c. crystal is similar to the S.C. model described above. Only nearest neighbour interactions are included and the transition probabilities are given by the same equations. The solid-on-solid restriction is applied to columns oriented along the (1 10) direction.23 The calculated growth rates are shown in fig. 6. A very large anisotropy is pre- dicted at small values of Ap/kT,where the nucleation rates are most sensitive to the edge free energy.The faces with the larger a factors grow more slowly as expected. The kinetics on the (110) face are typical of a normal growth law. Nucleation of 2D clusters is not necessary on this face since the atoms in a (1 10) layer are in isolated rows which are equivalent to the edges of steps on the (100) face. The growth rates of different faces of the f.c.c. and S.C. lattice are compared in fig. 7.23924 The open circles and triangles are the S.C. (100)growth rates at the indicated values of LlkT. The closed circles are f.c.c. (111) growth rates and the closed triangles are f.c.c. (100) rates both with L/kT = 12. The circles represent faces with a = 6 in both cases and the triangles a = 4. The agreement is very good; the slightly lower growth rates in the case of the two f.c.c.faces is probably a result of the larger bulk coordination and binding energy of that lattice. The logarithmic plot in fig. 7 also helps to delineate the region where the growth rate is determined by a nucleation process. If we assume that v cc Ap then from eqn (1) and (2) we have R = A(Ap/kT)516exp (-4F2/3kTAp) (3) where A is a constant. A plot of the logarithm of R from eqn (3)against (Ap/kT)-' is approximately linear and the slope is roughly equal to -4F2/3(lcT)2. The S.C. data at GEORGE H. GILMER A,u/kT FIG.6.-Crystal growth rates on three faces of an f.c.c. model. The growth rates are normalized by the product of k+ and the lattice constant a. A (Ill) 0(loo) 0(1 10). L/kT = 20. L/kT = 6 in fig.7 are approximately linear for kT'/A,u > 2 and the slope in this region affords an estimate of F. A value of F is chosen such that the slope of a plot of eqn (3) is in accord with that of the data. This yields the value F/j = 0.40 5 0.05 and the pair approximation gives F/j = 0.47 at this temperature. A similar procedure with L/kT = 9 yields F/j = 0.90 &-0.05 and the pair method gives F/j = 0.78. The small discrepancies here are partly a result of the pair approximation but apparently the values of F calculated by eqn (1 3) are too large. kl/A/ FIG.7.-Logarithm of the normalized growth rates on several S.C. and f.c.c. faces (see text). Here all growth rates are normalized by the product of kf and the corresponding layer spacing. Open triangles LIkT = 6; open circles L/kT = 9.SIMULATION OF 2D NUCLEATION AND CRYSTAL GROWTH SURFACE MIGRATION The inclusion of surface migration may increase the growth rate dramatically as already observed in the transient calculations of fig. 2. The steady-state growth rates with a migration to evaporation ratio of 7.4 and 54.6 are shown in fig. 8 for L/kT = 12. For comparison the data of fig. 3 without surface mobility are also included. The effect of mobility is accentuated at small values of Ap where the nucleation rate is small A) 1kT FIG.8.-Normalized growth rates on the S.C. (100) face at L/kT = 12. The values of the surface mobility are indicated by the ratios k,/k-next to the curves. The curves were calculated using eqn (3) with Flj equal to 1.2,O.g and 0.75 for k,/k-equal to 0 7.4 and 54.6 respectively.clusters are far apart and the competition for the adatom flux is not important. This preferential enhancement at small driving force causes a reduction in the slope of a plot of the logarithm of the rate and hence calculations of F yield erroneous values when the surface atoms are mobile. CONCLUSlONS Realistic models of crystal growth are essential for a detailed understanding of the process. Crystal growth kinetics are directly related to the structure of the crystal surface and accurate representations of this surface are required. The Ising model simulation of the growth rate transient is in much better accord with measurements than were the previous models. Cooperative interactions among large groups of atoms are accurately represented in this model.A qualitative change in the kinetics is observed at a critical roughening temperature. Above this point the normal growth law obtains but at lower temperatures the perfect crystal grows by a nucleation mechanism. A reduction in the free energy of a close packed step is observed at temperatures as low as half of the roughening point and this tends to increase the growth rate. Some low index faces may roughen at very low temperatures if the atoms in the surface layer are not connected by a strong network of bond chains. These temperatures can be calculated from the a factors of the crystal faces. GEORGE H. GILMER The author thanks J. D. Weeks and IS.A. Jackson for helpful discussions and V.Bostanov for suggesting a calculation of the transient. L. D. Hulett and F. W. Young J. Electrochem. Soc. 1966,113,410. ’A. Bewick M. Fleischmann and H. R. Thirsk Trans. Farday Soc. 1962,58,2200. U. Bertocci J. Electrochem. SOC.,1972 119 822. F. F. Abraham and G. H. White J. Appl. Phys. 1970,41 1841. V. V. Solovev and V. T. Borisov Sou. Phys. Crystallography 1973,17 814. G. H. Gilmer and P. Bennema J. Appl. Phys. 1972,43,1347. J. D. Weeks G. H. Gilmer and K. A. Jackson J. Chem. Phys. 1976,65,712. H. Darbinghaus and H. J. Meyer J. Crystal Growth 1972 16 31. V. Bostanov R. Roussinova and E. Budevski J.Electrochem. Soc. 1972,119,1346. lo L. A. Borovinskii and A. N. Tsindergozen,Sou. Phys. Crystallography 1969,13,1191. J. D. Weeks G.H. Gilmer and H. J. Leamy Phys. Rev. Letfers 1973,31 549; also see H. J. Leamy G. H. Gilmer and K. A. Jackson in Surface Physics of Materials I ed. J. M. Blakeley (Academic Press New York 1975) p. 121. J. P. van der Eerden R. L. Kalf and C. van Leeuwen J. Crystal Growth 1976,35,241. l3 G. H. Gilmer in Computer Simulation.for Materials Applications (National Bureau of Standards Gaithersburg 1976) p. 964. l4 S. W. H. de Haan V. J. A. Meeussen B. P. Veltman P. Bennema C. van Leeuwen and G. H. Gilmer J. Crystal Growth 1974,24/25,491. l5 A. E. Nielsen Kinetics of Prec@itation (Pergamon Oxford 1964). l6 D. Walton J. Chem. Phys. 1962,37,2182. l7 D. E. Temkin Sou. Phys. Crystallography 1969 14 179. l8 G. H. Gilmer J. Crystal Growth 1976 35 15; B. Lewis J. Crystal Growth 1974,21,29.l9 J. D. Weeks and G. H. Gilmer J. Crystal Growth 1977 in press; G. H. Gilmer and J. D. Weeks J. Chem. Phys. 1978 in press. 2o L. Onsager Phys. Rev. 1944,65 117. *’ M. E. Fisher and A. E. Ferdinand Phys. Rev. Letters 1967,19 169. ”K. A. Jackson in Liquid Metals and Solidification (American Society for Metals Metals Park Ohio 1958) p. 174. 23 G. H. Gilmer and K. A. Jackson in Crystal Growth and Materials (North-Holland Amsterdam 1976) p. 79. ”U. Bertocci J. Crystal Growth 1974 26 219.