- Email: [email protected]
Kinetic Monte Carlo study of crystal growth from solution
Computer Physics Communications 138 (2001) 250–263 www.elsevier.com/locate/cpc
Kinetic Monte Carlo study of crystal growth from solution Mirosława Rak ∗ , Marek Izdebski, Andrzej Brozi Institute of Physics, Technical University of Łód´z, ul. Wólcza´nska 219, PL-93-005 Łód´z, Poland Received 7 June 1999; received in revised form 26 March 2001
Abstract A study of both surface micromorphology and growth kinetics of the (001) face of a crystal growing from solution is performed on the basis of Monte Carlo simulations taking into account surface diffusion. Special attention is paid to growth simulations performed for layer-by-layer growth mechanism. An influence of the size of two-dimensional array describing crystal surface on simulation results is also discussed. The presented results show that the size of the surface array affects both growth mechanism and growth rate. The simulation results are confronted with predictions of crystal growth theories. The influence of surface diffusion on growth process is compared with some results reported for vapor growth, and it is found that for solution growth, the influence is much weaker. An implementation of fast simulation algorithm with explicit growth time is presented. It is revealed that calculations of the average crystal growth rate using time proportional to the number of events fail in the case of the perfect face growing at low supersaturations. 2001 Elsevier Science B.V. All rights reserved. PACS: 81.10.Aj; 81.10.Dn; 02.70.Lq; 07.05.Tp Keywords: Monte Carlo simulations; Crystal growth models; Solution growth; Layer-by-layer growth mechanism; Surface diffusion; Two-dimensional nucleation
1. Introduction Monte Carlo (MC) simulations are a convenient method for investigation of both surface micromorphology [1,2] and growth kinetics (e.g., Refs. [1–4]) of crystals. The model of Gilmer and Bennema [1,2] takes special position among numerous works in this field. It is based on the kinetic Ising’s model with the (001) face of Kossel crystal and on the solid on solid assumption. The model proposed in Refs. [1,2] is still used as a basis for construction of new ones (e.g., Refs. [5–9]) and its results are treated as a standard other works are compared to [10–12]. * Corresponding author.
E-mail address: [email protected] (M. Rak).
Although the statistical mechanical background of the simulations of crystal growth from solution were presented if Ref. [11], most of the published MC simulation results concern vapor growth, e.g., Refs. [1– 6,10]. The works concerning solution growth (e.g., Refs. [7–9,11,13]) are not numerous and they rather do not contain results obtained for growth conditions close to those used in practice of single crystal growth, i.e. for low supersaturations. The model proposed by Gilmer and Bennema was also designed for simulations of vapor growth, however, after modification of the formulae describing frequencies of elementary events (creation, annihilation and surface diffusion), it may be used for simulations of solution growth [7]. The aim of this paper is to study, using MC simulations, the growth kinetics and micromorphology of
0010-4655/01/$ – see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 1 ) 0 0 2 3 8 - 7
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
the (001) face of a crystal growing from solution at supersaturations varying in wide range. The simulations are performed for both perfect crystal face and crystal face with steps. Special attention is paid to the growth by two-dimensional nucleation at supersaturations as low as possible. In our work, we also analyse the influence of the size of the two-dimensional array describing crystal surface (further on called the surface array) on the simulation results. The influence of surface diffusion on growth kinetics is also studied. Performing growth simulations for perfect crystal face growing at low supersaturations requires generation of enormous number of events. For this reason, we have applied a simulation algorithm which is much more efficient than the one used by Gilmer and Bennema [1,2]. In this algorithm, growth time is defined explicitly, and therefore the average growth rate can be calculated as a quotient of the normal growth of the face and the time of the growth. The measure of elapsed time used in Refs. [1,2] is not applicable in our model.
2. The model 2.1. Statement of the problem In this study, we simulate the growth of the (001) face of a simple cubic crystal. The crystal is a tetragonal Kossel crystal, thus we take into account only first neighbor interactions. The solid on solid (SOS) condition, which excludes solid overhangs, fluid inclusions in crystal and vacancies, is imposed. A two-dimensional (2D) square array represents a crystal surface. Growth units in the crystal-solution interface are represented by solid tetragonal blocks. Three kinds of elementary events are possible to occur on the crystal surface, namely: (1) creation (addition of a solid block to the surface) which represents adsorption of a growth unit at crystal surface; (2) annihilation (subtraction of a solid block from the surface) which represents desorption of a growth unit from crystal surface; (3) surface diffusion (subtraction of a surface solid block from a site of the surface array and immediate addition of it on one of the four nearest neighbor surface positions).
251
In the investigation presented here, we do not take into account the solvation and desolvation phenomena because their influence on growth kinetics and surface morphology was analyzed in our earlier paper [13]. It should be noted that the above presented assumptions, used in this model, are coincident with those used in Refs. [1,11], where they were discussed in detail. 2.2. Creation and annihilation frequencies The exchange of solid blocks between the solid and fluid phases is governed by the frequencies of creation and annihilation. The frequencies are equal to the transition probabilities per unit time for creation or annihilation and are related to the energy change of the system due to creation or annihilation. In Refs. [1, 2], the frequencies of creation ki+ and annihilation ki− were assumed to be independent of the position on the surface and supersaturation, respectively. In the present study concerning solution growth, we assume that the frequencies of creation ki+ and annihilation ki− are given by the formulae proposed by Binsbergen [1, 14]: β α ki+ = ft · exp − (2 − i) + , (1) 4 2 β α , (2) ki− = ft · exp (2 − i) − 4 2 where i is the number of solid lateral neighbors (i = 0, 1, 2, 3, 4). We do not consider vertical bonds because the SOS condition is imposed. The frequency ft denotes the number of trials to occur per unit time. It depends on the free energy Gη of activation of viscous flow [14]: Gη kT exp − ft = , (3) h kT where h is the Planck constant, k is the Boltzmann constant and T is temperature. The parameter α, occurring in Eqs. (1) and (2), describes thermal roughening of crystal surface and is defined [1,2,12,15] as 2(Φss + Φff ) − 4Φsf , (4) kT where Φss is the bond energy of a solid–solid nearest neighbor pair, Φsf is the energy of solid–fluid pair, Φff is the interaction energy between the average α=
252
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
contents of the neighboring blocks of fluid. Expression for the parameter α for a real solid–solution system is derived in Ref. [11]. Dimensionless parameter β, occurring in Eqs. (1) and (2), is sometimes called the kinetic roughening coefficient [12] and represents driving force for crystallization. It is proportional to the difference µ of chemical potentials [16,17] between the solute µfs and the solid phase µs [11]: β=
µ µfs − µs = , kT kT
(5)
where it is assumed that µs = µfs eq (corresponding to equilibrium, i.e. saturated solution) [11]. For perfect [16,17] and ideal [11,16,17] solutions it was derived (e.g., Ref. [11]) that µ = ln(1 + σ ), (6) kT where σ denotes the relative supersaturation. However, Eq. (6) is not valid for non ideal solutions since in such a case, chemical potentials depend not only on concentration of solute molecules but also on activity coefficients [11,16,17] which describe mutual interaction of solute particles. It should be pointed out that in our simulations, µ is assumed to be constant during growth process. To keep this quantity constant requires either an additional reservoir of solid particles in the crystallizer or evaporating solvent from the solution to the gas reservoir. It can be easily seen that the frequencies given by Eqs. (1) and (2) satisfy the condition ki+ α (7) = exp (i − 2) + β . 2 ki− On the other hand, Eq. (7) relating creation and annihilation frequencies was obtained [11] as a result of statistical mechanical arguments and the principle of microscopic reversibility which was supposed to be also valid for non-equilibrium situations. The fact, whether µ is produced by fluctuations or by growth conditions in a crystallizer does not make any difference [11]. It should be mentioned that Eq. (7) was also derived by Gilmer and Bennema [1]. However, the authors did not use statistical mechanical arguments. They considered [1] the principle of microscopic reversibility for equilibrium ( µ = 0) and took into
account the fact that in equilibrium, frequencies of creation and annihilation at kink sites (i = 2) must equal each other. In this way it was found that for equilibrium: ki+ α (i − 2) . (8) = exp 2 ki− For the case of non-equilibrium, according to rate theory it was shown [1] that the ratio ki+ /ki− increases by a factor exp(+β), which leads to Eq. (7). Eq. (7) enables a simple kinetic interpretation of driving force for crystallization. For a kink site (i = 2), we obtain k2+ /k2− = exp(β), which shows how unbalanced is the creation/annihilation of kink sites. Thus parameter β characterizes the growth kinetic of a 2D layer. It should be pointed out that our simulations concern layer-by-layer growth mechanism occurring at very low supersaturations, which means that in this paper, we study near-equilibrium growth. The kinetic roughening coefficient β is varied from the lowest value at which growth is observed to 0.8, and the parameter α is varied from 4 to 8. 2.3. Surface diffusion frequency In our simulations a surface diffusion jump is realized in a way similar to that presented in Ref. [1]. Therefore, the frequency kij of surface diffusion jump from a site with i lateral solid neighbors to a site with j lateral solid neighbors is calculated using the method presented in Ref. [1]. Since surface diffusion is represented by the combination of subtracting a solid block from a site and adding it on one of the four neighbor sites, the frequency kij can be assumed [1] to be equal to: kij = bki− kj+ ,
(9)
where b is a constant with the dimension of time. For i = j = 0 it follows from Eq. (9) that k00 (10) b = − +. k0 k0 Since an isolated (i = 0) adsorbed growth unit jumps to a new isolated (j = 0) site with frequency k00 , the average time τsdiff between surface diffusion jumps is given by: 1 . (11) τsdiff = k00
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
Additionally, the frequency k0− of annihilation from isolated position can be assumed [1] to be equal to the reciprocal of the mean lifetime τdeads of such a growth unit on the surface 1 (12) τdeads = − . k0 On the other hand, the mean diffusion distance Xs is expressed (e.g., Refs. [1,15]) by Xs = Ds τdeads , (13) where the surface diffusion coefficient is given by the Einstein equation (e.g., Refs. [1,15]): Ds =
a2
, (14) τsdiff where a denotes the lattice spacing. Combining Eqs. (9)–(14) we obtain 2 − + Xs ki kj . (15) kij = a k0+ After substituting Eqs. (1) and (2) into Eq. (15), we find [7] 2 β α Xs kij = ft exp (2 + j − i) − , (16) a 4 2 which differs from that used in Ref. [1] for vapor growth. In our simulations we consider values of Xs /a from 0 to 3. It should be pointed out that Eq. (16) satisfies the principle of microscopic reversibility as kij α (17) = exp (j − i) . kj i 2 In our simulations, the number j of solid lateral neighbors of a diffusing block at the target location is calculated after the removal of the block from its starting position. It seems natural, although the reverse possibility is also worth considering. The calculation of the neighbors number at the target position before the removal of the diffusing block enables preference of diffusion jumps along the shortest path, i.e. without change of the layer.
3. Method of simulation In the Monte Carlo method, selection of both the event to perform and its surface coordinates is
253
realized by comparison of random numbers with some number intervals. In the algorithm presented in Refs. [1,2] the lengths of the intervals are not varying during the simulation. The adaptation of the probability distribution of the currently generated events to the current state of the crystal surface is obtained by introducing another random number selection deciding whether to realize or abandon the event just chosen. The efficiency of the Gilmer and Bennema algorithm [1,2] depends strongly on the value of the parameter α. The ratio of the number of generated events to the number of the algorithm loop repetitions depends on the selection of formulae for the frequencies of elementary events, though in general, it decreases quickly with increasing value of parameter α. Apart from this, if we intend to take into account the influence of many factors (e.g., mechanical stresses, solvation, more than one kind of growth units) on the frequencies of elementary events, then the algorithm includes additional random number selections, which may cause abandoning the event. Every subsequent selection lowers the algorithm efficiency in geometrical progression. The algorithm used in our simulations provides both high effectiveness (independent of the parameter α and the solution supersaturation) and ability to take into account new parameters and phenomena, with only linear increase in simulation cost. In order to achieve this we introduced an array containing the sums of frequencies of all possible individual events for every possible surface coordinates (x, y). After every event some elements of the array situated near the site of the last event must be updated. Selection of the coordinates, where some event is to occur, is realized by selecting one of the number intervals (the widths of which are stored in an array) with the use of random number generator with uniform distribution. Then the configuration of the selected location is checked and one of the possible events is chosen with probability proportional to the frequency of this event in the location. Inspection of subsequent frequency array elements as a means for selection of the event coordinates is a very ineffective solution. On the other hand, the use of a binary tree structure is not a good solution either. The binary tree enables fast finding of coordinates, but its updating after every event is too time-consuming. The best results are obtained by dividing the linear search
254
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
into 2 to 4 levels, in which the coordinate ranges are being contracted. Every level of the search has its corresponding array storing sums of frequencies of events in the surface areas which may be chosen in this level. The implementation of the algorithm described above is presented in Appendix A. Detailed discussion of the efficiency of some MC algorithms together with ideas for its improvement may be found in Refs. [3,4]. The algorithm can be extended to take into account other kinds of elementary events with linearly increasing simulation cost. For example, the results of simulations taking into account solvation and desolvation phenomena were presented in Ref. [13]. In the performed simulations, the main loop of the algorithm was repeated until 11 layers grew on the crystal face. As the simulations started from perfectly smooth surface, the growth of the first layer was not taken into account in the calculations of the growth rate. In this way, the influence of the transition state occurring at the beginning of the simulation is limited.
algorithm of Gilmer and Bennema. The number of creations is a good measure of elapsed time only in the case of creation frequency independent of the number of lateral neighbors. In Section 5.1, we will discuss influence of some approximated measures of time on calculated average growth rates. In our simulations, time tn passing between the nth, and the next event is defined as an reciprocal of the sum K of the frequencies of all possible events in all sites, calculated for the surface configuration determined after the n events. The sum K is defined by Eq. (A.3) in Appendix A. Hence we find: tn =
1 . K(n)
(19)
Consequently, the time t passing in the system of the crystal since the beginning of the simulation is a sum of all elementary times tn calculated after occurrence of every event: t= tn . (20) n
4. Normal crystal growth rate In the model of Gilmer and Bennema, the average growth rate of a crystal face is calculated [1,2,10] on the basis of the difference between the averaged fluxes of the added and removed blocks (I + and I − ): d + (I − I − ) N 4 d + F = ki N i − ki− N Si , N
R =
(18)
i=0
where N Fi and N Si are the average numbers of fluid and solid blocks at the interface, respectively, having i solid lateral neighbors, d is the thickness of a one block layer and N is the number of all blocks in one layer. The average values N Fi , N Si should be found as an arithmetic mean of the instantaneous values NiF , NiS calculated every preset constant time period. However, the time measure used in the model presented in Ref. [2] causes the mean values to be approximate. It is possible to use the number of events [18], the number of creations [2,9] or the number of repetitions of the main algorithm loop as a measure of elapsed time. The latter method is suitable only for the original
The algorithm described in Appendix A contains the sum of frequencies K of all the possible events, which is updated after every occurring event. Therefore, the additional cost of calculation of the growth time defined by Eqs. (19) and (20) is negligible. Using the defined time, values of N Fi and N Si can be calculated according to the definition of the mean value. However, it is equivalent but simpler to calculate the average crystal growth rate as a quotient of normal growth of crystal face and the growth time: d N , (21) tN where N is an increment in the solid block number during time t. The increment N is taken from the results of simulation. The growth rate of the face depends directly on the thickness d of one layer of growth units and indirectly on the free energy Gη of viscous flow activation appearing in the definition of the frequency ft (Eq. (3)). As the dependence of the growth rate on the coefficient β at various values of parameters α and Xs is of particular interest to us, we introduce dimensionless average growth rates defined as: r=
R =
R dft
(22a)
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
255
and r =
r , dft
(22b)
where R and r are given by Eqs. (18) and (21), respectively. The dimensionless growth rate r is proportional to the growth rate expressed in physical units, e.g., in cm/s, if the activation energy Gη of viscous flow and temperature (and hence α) are constant, while the parameter β (corresponding to the supersaturation of the solution — Eq. (6)) may vary as it does not affect the proportionality coefficient between the two growth rates. In earlier papers [1,2,12,19], the dimensionless growth rate was calculated using the expression: R =
R . dk0+
(23)
It should be noted, however, that the frequency k0+ depends on the parameter β, and this dependence appears in both the vapor growth model [1,2] and Eq. (1). Such a method of calculation of the dimensionless growth rate affects the shape of the curve of the growth rate versus β.
5. Results and discussion 5.1. The comparison of growth rates calculated with the use of different methods In previous section, we concluded that the average value of the growth rate can be calculated only approximately if the growth time is not defined explicitly. In order to compare results obtained using different methods of calculation of growth rate, we found average growth rates R using Eq. (22a) as well as r using Eq. (22b). The rate R was calculated on the basis of the arithmetic mean of NiF , NiS sampled every preset small number of events. In extreme cases, the values of the rate R were negative, while the solid block number was observed to increase, even when the numbers of instantaneous values NiF , NiS used were over a million. Fig. 1 shows that the greatest differences between the rates R and r correspond to the simulations of the growth of crystal face without defects, at the lowest supersaturations, at which mononuclear growth mechanism operates. Under these conditions, the variations in values of NiF and NiS show some correlation with
Fig. 1. The plots of quotient of the average growth rates calculated in different ways. The growth rate R was calculated using Eq. (18), S with N F i , N i obtained by averaging the instantaneous values determined every preset small number of events, while the growth rate r was calculated according to Eq. (21). Both plots are obtained for a perfect face at α = 6 and Xs = 1.5a. (a) surface array size 120 × 120; (b) surface array size 30 × 30.
the time necessary for the preset number of events to occur, and therefore it seems that during the averaging of instantaneous values of NiF and NiS , if they are determined every preset number of events, an excessive pile-up of the samples with the lowest or greatest values may occur. The discrepancy between the results given by both methods of the average growth rate determination decreases with increasing size of the sur-
256
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
face array (cf. Figs. 1(a) and 1(b)) and with increasing β. When the growth rate is determined with the use of explicit growth time given by Eq. (20), we observe smaller scatter in the points obtained for the lowest supersaturations used in the plots r (β). Our calculations show that, in the original model of Gilmer and Bennema, it is possible to achieve better accordance between growth rates R and r by the way of sampling the values NiF and NiS every preset number of repetitions of the algorithm main loop. However, this method depends on specific features of this particular algorithm. 5.2. Growth of perfect crystal face 5.2.1. Mononuclear and multiple nucleation growth mechanisms It is well known that, at high α (α > 3.5), perfect crystal faces without steps generated by spiral dislocations grow by two-dimensional (2D) nucleation mechanism, if the supersaturation exceeds a critical value (e.g., Refs. [20–23]). Therefore, for 2D nucleation growth mechanism, supersaturations slightly higher than the critical one will be called “low supersaturations” in the remaining part of Section 5.2. Close to this critical value, small changes in β cause great relative changes in growth rate (see Table 1). Our investigation shows that, at low supersaturations and high values of α (α 4), the size of the surface array strongly influences the simulation results (cf. Table 1). The MC simulations reveal that, at α 4, crystal growth may occur either by formation of one 2D nucleus (Fig. 2(a)) spreading across the surface (mononuclear model [22,23]) or by formation of many 2D nuclei (multiple 2D nucleation) spreading on the surface according to birth and spread (B + S) model [23], depending on the size of the surface array and supersaturation value (see Table 2). At higher β, nuclei are sometimes formed on the spreading nuclei and then growth occurs by nucleus above nucleus (NaN) growth mechanism (Fig. 2(b)). The influence of the surface diffusion on the observed growth mechanism is negligible. According to the mononuclear growth model [22, 23], nucleus born on the surface spreads out and covers the entire surface before other critical size nuclei are born on this surface. Therefore, it is possible for the mononuclear growth mechanism to operate [20,22]
Table 1 The dependence of the average crystal growth rate r on β for α = 6, Xs = 1a and two different array sizes β
r Array 30 × 30
Array 120 × 120
0.35
0.00002
0.00038
0.38
0.00011
0.00133
0.40
0.00015
0.00176
0.42
0.00020
0.00309
0.45
0.00029
0.00463
0.50
0.00251
0.00711
0.55
0.00325
0.00997
0.60
0.00746
0.01298
0.65
0.00757
0.01732
0.70
0.02026
0.02463
if supersaturation and the area of the surface are sufficiently small. In crystal growth practice, this mechanism is hardly ever observed, because areas of real crystal surfaces are much too large. In our MC simulations performed for relatively low supersaturations and small array sizes, we did observe (Fig. 2(a)) this growth mechanism. With decreasing value of β, the number of coexisting nuclei decreases sharply. For the lowest values of β at which the face can still grow (at given α > 4), we observed single nucleus spreading across whole surface. The values of β, at which growth proceeds by mononuclear mechanism, decrease with increasing size of the surface array (cf. Table 2). It was shown [20,22] that mononuclear growth mechanism operates if the area S of crystal face fulfill the condition: S 1/2 < Lcr = (νv /J )1/3,
(24)
where Lcr is the critical length of face side, below which growth by mononuclear growth mechanism can occur, νv is the spread velocity of nucleated steps [22, 23] νv = D1 σ,
(25)
and J denotes the rate of 2D nucleus formation [22, 23] π γ 2 1/2 J = D2 β exp − . (26) β kT
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
257
Fig. 2. Growth by 2D nucleation mechanism at perfect crystal face (surface array size 120 × 120). (a) Mononuclear growth mechanism, α = 5, β = 0.2, Xs = 0; (b) NaN mechanism, α = 6, β = 0.7, Xs = 0. Table 2 The characteristic values β ∗ for which the transition from mononuclear to B + S growth mechanism is observed. The results are obtained for diffusion distance Xs = 1a α
Size of surface array 20 × 20
30 × 30
60 × 60
120 × 120
4
0.07
0.06 . . . 0.07
0.05
0.04 . . . 0.05
5
0.41
0.35
0.25
0.20
6
0.77
7
240 × 240
0.16
0.60
0.49
0.42
0.38
0.90
0.78
0.70
0.60
Here D1 and D2 are parameters independent of σ , β and α. The edge free energy γ is related to the parameter α [24]: γ α ≈ . (27) kT 4 After substituting Eqs. (25), (26) and (27) into Eq. (24) we obtain
S 1/2 < Lcr =
D1 D2
1/3
πα 2 . (28) σ 1/3 β −1/6 exp + 48β
Eq. (28) reveals that Lcr is very sensitive to changes in values of α and β, especially at low values of β. For example, it follows from Eq. (28) that at β = 0.1, the change in α from 4 to 5 causes an increase in Lcr by a factor of 2.7 × 103 . Furthermore, the change
258
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
in β from 0.04 to 0.05 causes reduction in Lcr by a factor 181, at α = 4. In our MC simulations, we did not observe such a strong influence of α and β on Lcr (cf. Table 2), at low β. At higher β (β > 0.5), as Lcr given by Eq. (28) is not so sensitive to changes in β and α, MC simulation results are in better accordance with Eq. (28). Mononuclear model predicts a linear dependence of the growth rate on the total area of crystal surface [23]. We found that, at given value of β, the growth rate increases with increasing area of the surface array, but it is rather difficult to describe this dependence by a linear function. However, we should remember that the theory assumes an infinite spread velocity of nucleus. Such an assumption is very close to the case of very small both the surface array size ( 60 × 60) and supersaturation values. Fig. 3 illustrates the dependence of the square root of the average growth rate r on the length of the array side. This figure shows that the dependence of the growth rate on the surface area is nearly linear at low supersaturations and for small surface arrays, where mononuclear growth mechanism operates and the time necessary to fill the whole layer by spreading nucleus is negligible. In Fig. 3, the curve obtained for β = 0.38 corresponds to mononuclear growth mechanism for the surface array sizes below 240 × 240 (cf. Table 2). It is seen in this figure (curve for β = 0.38) that, for the surface arrays of sizes over 120 × 120, the growth rate increases very slowly with increasing array size.
Fig. 3. The dependence of dimensionless growth rate r on the size of the surface array, for a perfect (001) face growing by 2D nucleation, α = 6, Xs = 1a.
In this range of array sizes, the spread time of nucleus cannot be neglected as it is assumed in mononuclear model. For the smallest arrays, the points corresponding to β = 0.38 (Fig. 3) show that, √ with increasing length of the surface array side, r increases a little slower than proportionally. This small discrepancy between mononuclear model and MC simulation results also appears for various values of mean diffusion distance Xs and is a result of too small surface array size, which is comparable to the critical size of twodimensional nuclei. Moreover, in the MC simulations, we used cyclic connection of the surface array sides, which is not used in theoretical model. For values of β high enough for the growth mechanism described by B + S model to operate, the growth rate increases slightly with increasing array size, tending to a constant value. However, it should be pointed out that the characteristic value of β (denoted β ∗ ) corresponding to the transition from one to the other growth mechanism increases sharply with increasing value of the parameter α (cf. Table 2). Although mononuclear and B + S growth mechanisms can be distinguished directly by observations of crystal surfaces, the plots of the number of grown layers versus growth time t (given by Eq. (20)) may also be useful. For the lowest values of β still enabling crystal growth, the curves have a form of irregular sharp steps (Fig. 4(a)), where the horizontal lines correspond to waiting for the formation of 2D nucleus and the nearly vertical lines correspond to rapid spreading of the nucleus. For slightly higher values of β (but still lower than the characteristic value β ∗ ) and for sufficiently large size of the surface array, the spread time may become visible in a plot (e.g., Fig. 4(b)). If the growth proceeds by B + S mechanism, the plotted curve is a diagonal line with random fluctuations (Fig. 4(c)). 5.2.2. Influence of surface diffusion on growth rate and surface morphology In the case of vapor growth, surface diffusion plays an essential role and it considerably increases the nucleation rate J (even by a factor of 102 ÷ 103 [25]). Consequently, crystal growth rate also increases considerably (e.g., Ref. [1]). It is known that in the case of solution growth, the influence of surface diffusion on the nucleation rate is not so pronounced [22]. In our simulations performed for
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
259
Fig. 5. The dependence of the dimensionless growth rate r on mean surface diffusion distance Xs plotted for a perfect face, α = 4, surface array size 120 × 120 and several values of β: (a) 0.05 β 0.10; (b) 0.5 β 0.8.
Fig. 4. The number of grown layers versus dimensionless growth time, for α = 6, Xs = 1a, surface array size 120 × 120.
α = 4 and several values of β. It is seen in this figure that the growth rate r increases moderately with increasing Xs /a. This effect can be explained by an increase in the spread velocity ν of 2D nuclei. According to B + S model [22,23] the normal growth rate R can be written: R = dν 2/3 J 1/3 ,
solution growth, we did not notice an increase in the nucleation rate with increasing mean diffusion distance Xs . Fig. 5 illustrates the dependence of dimensionless growth rate r on the mean diffusion distance Xs for
(29)
where d is the height of a nucleus and J is the nucleation rate given by Eq. (26). Since our investigation reveals that addition of growth units to the critical nucleus occurs by both direct impingement from the solution and surface diffusion, we assume that the total
260
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
spread velocity ν of a nucleus may be expressed as the sum of spread velocities νv and νs corresponding to direct incorporation and surface diffusion, respectively: ν = νv + νs .
(30)
The spread velocity νs , proportional to both Xs [23] and νv , can be expressed by Xs νv , (31) a where νv is given by Eq. (25) and coefficient A is assumed to be dependent on α and β. Upon substituting Eqs. (30) and (31) into Eq. (29), we find Xs 2/3 R = R0 1 + A , (32) a
νs = A
where R0 is the growth rate at Xs = 0. Therefore, we fitted our results (Fig. 5) according to the following equation Xs 2/3 r = r0 1 + A , (33) a and using the least-square method, we found coefficients A and r0 . The calculations performed for α = 4 reveal that for the lowest value of β (β = 0.07), coefficient A is equal to 0.39 and it decreases considerably with increasing β (even to 0.03 for β = 0.8). This means that at higher β, the spread velocity is determined by direct incorporation of growth units from the solution rather than by surface diffusion. At higher β, edges of nuclei are rough which should cause an increase in the spread velocity ν. On the other hand, however, surface diffusion reduces edge roughening.
These conflicting tendencies may cause reduction in the coefficient A at higher values of β. Although smoothing the edges of nuclei by surface diffusion was clearly visible, smoothing the crystal surface was not noticed. Similar calculations performed for α = 5 and α = 6 show that the coefficient A slightly increases with increasing α. This result indicates that the importance of surface diffusion in solution growth of perfect crystal face is greater for smooth surfaces, which seems to be natural. 5.3. Surfaces with a step We performed simulations for the (001) face containing a single step of one growth unit height. The step distance was assumed to be independent of supersaturation and equal to the length of the surface array side. For every array size and for various values of the parameter α 4, there is a range of supersaturations for which combined growth mechanism is observed: the 2D nucleation and the movement of the growth step (Fig. 6). This combined growth mechanism was also observed in MC simulations performed for vapor growth [1,19] and it is observed experimentally (e.g., Refs. [26–28]) for solution growth. The values of β corresponding to the combined growth mechanism decrease with increasing array size (and hence with increasing interstep distance). Therefore, the observed growth mechanism depends not only on α and β but also on the assumed interstep distance, especially at small sizes of the surface array.
Fig. 6. Combined growth mechanism: 2D nucleation and the growth by step movement, α = 5, β = 0.25, Xs = 1a, surface array size 120 × 120.
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
At high supersaturations, where 2D nucleation growth mechanism dominates, the dependence of R on β is clearly stronger than linear. The critical value of the parameter β, corresponding to the change from a dominant spiral growth mechanism in a dominant 2D nucleation mechanism increases sharply with increasing parameter α (Fig. 7). For α = 5 and the surface array side 60 × 60 (Fig. 7(a)), the critical value β is approximately equal to 0.25, while for α = 8, it is higher than 0.8 and hence Fig. 7(b) shows the dependence r (β) for dominating spiral growth mechanism. According to the surface diffusion theory of Burton, Cabrera and Frank (BCF) [15,23], for low supersaturations, at which only spiral growth mechanism operates, the dependence of the growth rate on β is linear if constant (i.e. independent of β ) interstep distance is assumed. The results of our MC simulations are in
261
good agreement with the theoretical predictions even for the lowest values of β, although some small deviation of the r = f (β) dependence from linearity appears in most simulations (Fig. 7(a) for β < 0.20 and Fig. 7(b)). This deviation is dependent on the mean diffusion distance Xs and the value of parameter α and it is a result of superposition of some effects influencing the growth rate: • Increasing β stimulates kinetic roughening of the step edge. This effect causes that, with increasing β, the average growth rate r increases a little faster than linearly (Fig. 7(b), curve Xs = 0). • Smoothing the step edge by surface diffusion. This effect is most pronounced for not very large values of β. With increasing mean diffusion distance Xs , the diffusion smoothing reduces the importance of kinetic roughening of the step edge, and in consequence, dependence r on β becomes more linear (Fig. 7(b)). The smoothing of the surface of step terraces by surface diffusion we observed, was not so distinct as the one reported in Refs. [1,5] concerning vapor growth. • With increasing β the number of blocks being added directly to the step edges increases faster than the number of blocks being added at any other location at the surface and subsequently transported to the edges. Therefore, the relative changes in growth rate dependence on diffusion distance Xs are the greatest for low values of β. Our simulations performed for solution growth show that in general, surface diffusion increases the growth rate, but this influence is clearly weaker than in the case of simulations performed for vapor growth [1, 19]. Moreover, for vapor growth, the critical value of β (corresponding to the change of the growth mechanism) strongly depends [19] on surface diffusion. According to our results, this dependence is insignificant in the case of solution growth (Fig. 7(a)).
6. Conclusions
Fig. 7. Plots of dependence of the growth rate r of the crystal face with single growth step on the parameter β for Xs = 3a (series ), Xs = 1a ( ) and Xs = 0 ( ). Surface array size is 60 × 60. (a) α = 5; (b) α = 8.
The following conclusions can be drawn from the present study: (1) At high supersaturations, the average crystal growth rate can be calculated with the use of time measure approximated as a number of events [18], but for a perfect face growing at low supersatura-
262
(2)
(3)
(4)
(5)
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
tions, this method gives large scatter in results and in extreme cases, the sign of the average growth rate may even be inconsistent with the observed growth of the crystal face. The growth time defined by Eq. (20) and the average growth rate (21), based on this definition of time, can be used for all values of α and β. The characteristic value β ∗ (and hence supersaturation σ ∗ ), at which a transition from mononuclear to B + S growth mechanism occurs, decreases with increasing size of the array representing crystal surface. According to the theory [22], small changes in β and α values should cause enormous changes in the critical length Lcr of the face, especially for not very high values of α (e.g., α = 4). In our MC simulations, we did not observe such a strong influence of α and β on Lcr . The mononuclear model predicts linear dependence of growth rate on total area of crystal surface but, according to our results, this dependence is not linear. We have found that, for large sizes of the surface array and values of β near β ∗ , the spreading time of 2D nucleus cannot be neglected as it is assumed in the theory. For higher values of β, where B + S mechanism operates, the average crystal growth rate slightly increases with increasing array size, tending to a constant value. Our MC simulations results indicate that for solution growth, surface diffusion affects the spread velocity of nuclei rather than the nucleation rate. The influence of surface diffusion on the growth rate is not so strong as in the case of vapor growth and it decreases with increasing parameter β. MC simulations performed for a face with single step show that the type of the dominating growth mechanism (2D nucleation or the growth step movement) depends considerably on parameters α and β and also on the size of the surface array. In contrast to the results reported for vapor growth, surface diffusion does not cause change of the growth mechanism.
Appendix A. Simulation algorithm In Section 3, we described the main assumptions for our implementation of fast algorithm of MC simulations. Here we discuss in detail the realization of this algorithm in the case of block growth of the (001) face. In this algorithm, three kinds of elementary events occurring on the crystal surface are taken into account: addition of a block to the surface (creation), removal of a block (annihilation) and surface diffusion. The frequencies of the events are defined as follows: (1) The frequency of any event on the phase boundary in a location with coordinates (x, y): − + kj+(x,y) + Kxy = ki(x,y)
4
ki(x,y),j (x ,y ) , (A.1)
q=1
where ki− is the annihilation frequency in a location with i solid lateral neighbors, kj+ is the creation frequency in a location with j solid lateral neighbors, kij is the surface diffusion frequency occurring from a location with i lateral neighbors onto a location with j lateral neighbors, q is the diffusion direction, and (x , y ) denote coordinates of a diffusion target location calculated as x = x + sin(πq/2) and y = y + cos(πq/2). (2) The frequency of any event at interface in column x of growth units: Kxy . (A.2) Kx = y
(3) The frequency of any event occurring in any location is defined as a sum of frequencies Kxy over all elements of the crystal surface array K= Kxy . (A.3) x
y
If we divide the linear search into 2 levels, first finding the column x of the crystal surface array, designated by hxy , and next finding the row y, the simulation will be realized according to the following algorithm: 1. Preparation of the crystal surface array hxy corresponding to the boundary of the solid and fluid phases, preparation of the arrays of frequencies ki− , ki+ , kij , Kxy , Kx , calculation of the sum K and setting the time counter t to zero.
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
2. Generation of a random number RND in the interval 0 . . . K with uniform distribution. 3. Finding the column number: x := 0 until (x < xmax and Kx RND) RND := RND − Kx x := x + 1
References [1] [2] [3] [4] [5] [6] [7]
4. Finding the row number: y := 0 until (y < ymax and Kxy RND) RND := RND − Kxy
[8] [9] [10]
y := y + 1 5. Selection of the event:
[11]
+ if (RND < ki(x,y) ) ⇒ creation
[12]
+ − + ki(x,y) ) otherwise, if (RND < ki(x,y)
[13] [14] [15]
⇒ annihilation otherwise ⇒ diffusion. The diffusion direction q is being selected as a number of one of the intervals with lengths + ki(x,y),j (x ,y ) , containing the number RND − ki(x,y) − − ki(x,y) . 6. Incrementing and/or decrementing elements of the crystal surface array hxy (and optionally hx y ) according to the event selected in point 5. 7. Incrementing the time counter: t := t + 1/K. 8. Updating the arrays Kxy , Kx and the sum K: for the coordinates (a, b) equal to (x, y), and (x , y ) in the case of diffusion, for the closest neighboring locations and for the locations neighboring to the closest neighbors the following operations should be performed: K := K − Kab , Ka := Ka − Kab , calculate new frequency Kab , K := K + Kab , Ka := Ka + Kab . 9. jump to point 2.
263
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
G.H. Gilmer, P. Bennema, J. App. Phys. 43 (1972) 1347. G.H. Gilmer, P. Bennema, J. Crystal Growth 13/14 (1972) 148. M. Kotrla, Comput. Phys. Commun. 97 (1996) 82. A.C. Levi, M. Kotrla, J. Phys.: Condens. Matter. 9 (1997) 299. R.-F. Xiao, J.I.D. Alexander, F. Rosenberger, J. Crystal Growth 109 (1991) 43. L. Guang-Zhao, J.P. van der Eerden, P. Bennema, J. Crystal Growth 58 (1982) 152. J. Karniewicz, A. Brozi, M. Rak, P. Pietrucha, SPIE Proc. 2373 (1995) 150. W.J.P. van Enckevort, A.C.J.F. van der Berg, K.B.G. Kreuwel, A.J. Derksen, M.S. Couto, J. Crystal Growth 166 (1996) 156. W.J.P. van Enckevort, A.C.J.F. van der Berg, J. Crystal Growth 183 (1998) 441. J.I.D. Alexander, in: J.P. van der Eerden, O.S.L. Bruinsma (Eds.), Science and Technology of Crystal Growth, Kluwer Academic Publishers, Dordrecht, 1995, p. 81. P. Bennema, J.P. van der Eerden, J. Crystal Growth 42 (1977) 201. P. Bennema, in: D.T.J. Hurle (Ed.), Handbook of Crystal Growth, Vol. 1a, North-Holland, Amsterdam, 1993, p. 477. M. Izdebski, M. Rak, A. Brozi, SPIE Proc. 3178 (1996) 135. F.L. Binsbergen, Kolloid Z. 237 (1970) 289. P. Bennema, G.H. Gilmer, in: P. Hartman (Ed.), Crystal Growth: An Introduction, North-Holland Publ., Amsterdam, 1973, p. 263. I. Prigogine, R. Defay, Chemical Thermodynamics, Longmans Green and Co., 1954. I. Prigogine, A. Bellemans, V. Mathot, The Molecular Theory of Solutions, North-Holland Publ., Amsterdam, 1957. F.F. Abraham, M. White, J. App. Phys. 41 (1970) 1841. P. Bennema, J. Crystal Growth 69 (1984) 182. A.I. Malkin, A.A. Chernov, I.V. Alexeev, J. Crystal Growth 97 (1989) 765. A.A. Chernov, Prog. Crystal Growth and Characterization 26 (1993) 121. A.A. Chernov, Crystal Growth, Crystallization Processes, Modern Crystallography, Vol. III, Springer, Berlin, 1984. M. Ohara, R.C. Reid, Modeling Crystal Growth Rates from Solution, Prentice-Hall Inc., Englewood Cliffs, NJ, 1973. P. Bennema, J. Boon, C. van Leeuwen, G.H. Gilmer, Krist. Tech. 27 (1973) 659. G.H. Pound, M.T. Simnad, L. Yang, J. Chem. Phys. 22 (1954) 1215. J.J. De Yoreo, T.A. Land, B. Dair, Phys. Rev. Lett. 73 (1994) 838. C.M. Pina, D. Bosbach, M. Prieto, A. Putnis, J. Crystal Growth 187 (1998) 119. C.M. Pina, U. Becker, P. Risthous, D. Bosbach, A. Putnis, Nature 395 (1998) 483.
Kinetic Monte Carlo study of crystal growth from solution Mirosława Rak ∗ , Marek Izdebski, Andrzej Brozi Institute of Physics, Technical University of Łód´z, ul. Wólcza´nska 219, PL-93-005 Łód´z, Poland Received 7 June 1999; received in revised form 26 March 2001
Abstract A study of both surface micromorphology and growth kinetics of the (001) face of a crystal growing from solution is performed on the basis of Monte Carlo simulations taking into account surface diffusion. Special attention is paid to growth simulations performed for layer-by-layer growth mechanism. An influence of the size of two-dimensional array describing crystal surface on simulation results is also discussed. The presented results show that the size of the surface array affects both growth mechanism and growth rate. The simulation results are confronted with predictions of crystal growth theories. The influence of surface diffusion on growth process is compared with some results reported for vapor growth, and it is found that for solution growth, the influence is much weaker. An implementation of fast simulation algorithm with explicit growth time is presented. It is revealed that calculations of the average crystal growth rate using time proportional to the number of events fail in the case of the perfect face growing at low supersaturations. 2001 Elsevier Science B.V. All rights reserved. PACS: 81.10.Aj; 81.10.Dn; 02.70.Lq; 07.05.Tp Keywords: Monte Carlo simulations; Crystal growth models; Solution growth; Layer-by-layer growth mechanism; Surface diffusion; Two-dimensional nucleation
1. Introduction Monte Carlo (MC) simulations are a convenient method for investigation of both surface micromorphology [1,2] and growth kinetics (e.g., Refs. [1–4]) of crystals. The model of Gilmer and Bennema [1,2] takes special position among numerous works in this field. It is based on the kinetic Ising’s model with the (001) face of Kossel crystal and on the solid on solid assumption. The model proposed in Refs. [1,2] is still used as a basis for construction of new ones (e.g., Refs. [5–9]) and its results are treated as a standard other works are compared to [10–12]. * Corresponding author.
E-mail address: [email protected] (M. Rak).
Although the statistical mechanical background of the simulations of crystal growth from solution were presented if Ref. [11], most of the published MC simulation results concern vapor growth, e.g., Refs. [1– 6,10]. The works concerning solution growth (e.g., Refs. [7–9,11,13]) are not numerous and they rather do not contain results obtained for growth conditions close to those used in practice of single crystal growth, i.e. for low supersaturations. The model proposed by Gilmer and Bennema was also designed for simulations of vapor growth, however, after modification of the formulae describing frequencies of elementary events (creation, annihilation and surface diffusion), it may be used for simulations of solution growth [7]. The aim of this paper is to study, using MC simulations, the growth kinetics and micromorphology of
0010-4655/01/$ – see front matter 2001 Elsevier Science B.V. All rights reserved. PII: S 0 0 1 0 - 4 6 5 5 ( 0 1 ) 0 0 2 3 8 - 7
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
the (001) face of a crystal growing from solution at supersaturations varying in wide range. The simulations are performed for both perfect crystal face and crystal face with steps. Special attention is paid to the growth by two-dimensional nucleation at supersaturations as low as possible. In our work, we also analyse the influence of the size of the two-dimensional array describing crystal surface (further on called the surface array) on the simulation results. The influence of surface diffusion on growth kinetics is also studied. Performing growth simulations for perfect crystal face growing at low supersaturations requires generation of enormous number of events. For this reason, we have applied a simulation algorithm which is much more efficient than the one used by Gilmer and Bennema [1,2]. In this algorithm, growth time is defined explicitly, and therefore the average growth rate can be calculated as a quotient of the normal growth of the face and the time of the growth. The measure of elapsed time used in Refs. [1,2] is not applicable in our model.
2. The model 2.1. Statement of the problem In this study, we simulate the growth of the (001) face of a simple cubic crystal. The crystal is a tetragonal Kossel crystal, thus we take into account only first neighbor interactions. The solid on solid (SOS) condition, which excludes solid overhangs, fluid inclusions in crystal and vacancies, is imposed. A two-dimensional (2D) square array represents a crystal surface. Growth units in the crystal-solution interface are represented by solid tetragonal blocks. Three kinds of elementary events are possible to occur on the crystal surface, namely: (1) creation (addition of a solid block to the surface) which represents adsorption of a growth unit at crystal surface; (2) annihilation (subtraction of a solid block from the surface) which represents desorption of a growth unit from crystal surface; (3) surface diffusion (subtraction of a surface solid block from a site of the surface array and immediate addition of it on one of the four nearest neighbor surface positions).
251
In the investigation presented here, we do not take into account the solvation and desolvation phenomena because their influence on growth kinetics and surface morphology was analyzed in our earlier paper [13]. It should be noted that the above presented assumptions, used in this model, are coincident with those used in Refs. [1,11], where they were discussed in detail. 2.2. Creation and annihilation frequencies The exchange of solid blocks between the solid and fluid phases is governed by the frequencies of creation and annihilation. The frequencies are equal to the transition probabilities per unit time for creation or annihilation and are related to the energy change of the system due to creation or annihilation. In Refs. [1, 2], the frequencies of creation ki+ and annihilation ki− were assumed to be independent of the position on the surface and supersaturation, respectively. In the present study concerning solution growth, we assume that the frequencies of creation ki+ and annihilation ki− are given by the formulae proposed by Binsbergen [1, 14]: β α ki+ = ft · exp − (2 − i) + , (1) 4 2 β α , (2) ki− = ft · exp (2 − i) − 4 2 where i is the number of solid lateral neighbors (i = 0, 1, 2, 3, 4). We do not consider vertical bonds because the SOS condition is imposed. The frequency ft denotes the number of trials to occur per unit time. It depends on the free energy Gη of activation of viscous flow [14]: Gη kT exp − ft = , (3) h kT where h is the Planck constant, k is the Boltzmann constant and T is temperature. The parameter α, occurring in Eqs. (1) and (2), describes thermal roughening of crystal surface and is defined [1,2,12,15] as 2(Φss + Φff ) − 4Φsf , (4) kT where Φss is the bond energy of a solid–solid nearest neighbor pair, Φsf is the energy of solid–fluid pair, Φff is the interaction energy between the average α=
252
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
contents of the neighboring blocks of fluid. Expression for the parameter α for a real solid–solution system is derived in Ref. [11]. Dimensionless parameter β, occurring in Eqs. (1) and (2), is sometimes called the kinetic roughening coefficient [12] and represents driving force for crystallization. It is proportional to the difference µ of chemical potentials [16,17] between the solute µfs and the solid phase µs [11]: β=
µ µfs − µs = , kT kT
(5)
where it is assumed that µs = µfs eq (corresponding to equilibrium, i.e. saturated solution) [11]. For perfect [16,17] and ideal [11,16,17] solutions it was derived (e.g., Ref. [11]) that µ = ln(1 + σ ), (6) kT where σ denotes the relative supersaturation. However, Eq. (6) is not valid for non ideal solutions since in such a case, chemical potentials depend not only on concentration of solute molecules but also on activity coefficients [11,16,17] which describe mutual interaction of solute particles. It should be pointed out that in our simulations, µ is assumed to be constant during growth process. To keep this quantity constant requires either an additional reservoir of solid particles in the crystallizer or evaporating solvent from the solution to the gas reservoir. It can be easily seen that the frequencies given by Eqs. (1) and (2) satisfy the condition ki+ α (7) = exp (i − 2) + β . 2 ki− On the other hand, Eq. (7) relating creation and annihilation frequencies was obtained [11] as a result of statistical mechanical arguments and the principle of microscopic reversibility which was supposed to be also valid for non-equilibrium situations. The fact, whether µ is produced by fluctuations or by growth conditions in a crystallizer does not make any difference [11]. It should be mentioned that Eq. (7) was also derived by Gilmer and Bennema [1]. However, the authors did not use statistical mechanical arguments. They considered [1] the principle of microscopic reversibility for equilibrium ( µ = 0) and took into
account the fact that in equilibrium, frequencies of creation and annihilation at kink sites (i = 2) must equal each other. In this way it was found that for equilibrium: ki+ α (i − 2) . (8) = exp 2 ki− For the case of non-equilibrium, according to rate theory it was shown [1] that the ratio ki+ /ki− increases by a factor exp(+β), which leads to Eq. (7). Eq. (7) enables a simple kinetic interpretation of driving force for crystallization. For a kink site (i = 2), we obtain k2+ /k2− = exp(β), which shows how unbalanced is the creation/annihilation of kink sites. Thus parameter β characterizes the growth kinetic of a 2D layer. It should be pointed out that our simulations concern layer-by-layer growth mechanism occurring at very low supersaturations, which means that in this paper, we study near-equilibrium growth. The kinetic roughening coefficient β is varied from the lowest value at which growth is observed to 0.8, and the parameter α is varied from 4 to 8. 2.3. Surface diffusion frequency In our simulations a surface diffusion jump is realized in a way similar to that presented in Ref. [1]. Therefore, the frequency kij of surface diffusion jump from a site with i lateral solid neighbors to a site with j lateral solid neighbors is calculated using the method presented in Ref. [1]. Since surface diffusion is represented by the combination of subtracting a solid block from a site and adding it on one of the four neighbor sites, the frequency kij can be assumed [1] to be equal to: kij = bki− kj+ ,
(9)
where b is a constant with the dimension of time. For i = j = 0 it follows from Eq. (9) that k00 (10) b = − +. k0 k0 Since an isolated (i = 0) adsorbed growth unit jumps to a new isolated (j = 0) site with frequency k00 , the average time τsdiff between surface diffusion jumps is given by: 1 . (11) τsdiff = k00
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
Additionally, the frequency k0− of annihilation from isolated position can be assumed [1] to be equal to the reciprocal of the mean lifetime τdeads of such a growth unit on the surface 1 (12) τdeads = − . k0 On the other hand, the mean diffusion distance Xs is expressed (e.g., Refs. [1,15]) by Xs = Ds τdeads , (13) where the surface diffusion coefficient is given by the Einstein equation (e.g., Refs. [1,15]): Ds =
a2
, (14) τsdiff where a denotes the lattice spacing. Combining Eqs. (9)–(14) we obtain 2 − + Xs ki kj . (15) kij = a k0+ After substituting Eqs. (1) and (2) into Eq. (15), we find [7] 2 β α Xs kij = ft exp (2 + j − i) − , (16) a 4 2 which differs from that used in Ref. [1] for vapor growth. In our simulations we consider values of Xs /a from 0 to 3. It should be pointed out that Eq. (16) satisfies the principle of microscopic reversibility as kij α (17) = exp (j − i) . kj i 2 In our simulations, the number j of solid lateral neighbors of a diffusing block at the target location is calculated after the removal of the block from its starting position. It seems natural, although the reverse possibility is also worth considering. The calculation of the neighbors number at the target position before the removal of the diffusing block enables preference of diffusion jumps along the shortest path, i.e. without change of the layer.
3. Method of simulation In the Monte Carlo method, selection of both the event to perform and its surface coordinates is
253
realized by comparison of random numbers with some number intervals. In the algorithm presented in Refs. [1,2] the lengths of the intervals are not varying during the simulation. The adaptation of the probability distribution of the currently generated events to the current state of the crystal surface is obtained by introducing another random number selection deciding whether to realize or abandon the event just chosen. The efficiency of the Gilmer and Bennema algorithm [1,2] depends strongly on the value of the parameter α. The ratio of the number of generated events to the number of the algorithm loop repetitions depends on the selection of formulae for the frequencies of elementary events, though in general, it decreases quickly with increasing value of parameter α. Apart from this, if we intend to take into account the influence of many factors (e.g., mechanical stresses, solvation, more than one kind of growth units) on the frequencies of elementary events, then the algorithm includes additional random number selections, which may cause abandoning the event. Every subsequent selection lowers the algorithm efficiency in geometrical progression. The algorithm used in our simulations provides both high effectiveness (independent of the parameter α and the solution supersaturation) and ability to take into account new parameters and phenomena, with only linear increase in simulation cost. In order to achieve this we introduced an array containing the sums of frequencies of all possible individual events for every possible surface coordinates (x, y). After every event some elements of the array situated near the site of the last event must be updated. Selection of the coordinates, where some event is to occur, is realized by selecting one of the number intervals (the widths of which are stored in an array) with the use of random number generator with uniform distribution. Then the configuration of the selected location is checked and one of the possible events is chosen with probability proportional to the frequency of this event in the location. Inspection of subsequent frequency array elements as a means for selection of the event coordinates is a very ineffective solution. On the other hand, the use of a binary tree structure is not a good solution either. The binary tree enables fast finding of coordinates, but its updating after every event is too time-consuming. The best results are obtained by dividing the linear search
254
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
into 2 to 4 levels, in which the coordinate ranges are being contracted. Every level of the search has its corresponding array storing sums of frequencies of events in the surface areas which may be chosen in this level. The implementation of the algorithm described above is presented in Appendix A. Detailed discussion of the efficiency of some MC algorithms together with ideas for its improvement may be found in Refs. [3,4]. The algorithm can be extended to take into account other kinds of elementary events with linearly increasing simulation cost. For example, the results of simulations taking into account solvation and desolvation phenomena were presented in Ref. [13]. In the performed simulations, the main loop of the algorithm was repeated until 11 layers grew on the crystal face. As the simulations started from perfectly smooth surface, the growth of the first layer was not taken into account in the calculations of the growth rate. In this way, the influence of the transition state occurring at the beginning of the simulation is limited.
algorithm of Gilmer and Bennema. The number of creations is a good measure of elapsed time only in the case of creation frequency independent of the number of lateral neighbors. In Section 5.1, we will discuss influence of some approximated measures of time on calculated average growth rates. In our simulations, time tn passing between the nth, and the next event is defined as an reciprocal of the sum K of the frequencies of all possible events in all sites, calculated for the surface configuration determined after the n events. The sum K is defined by Eq. (A.3) in Appendix A. Hence we find: tn =
1 . K(n)
(19)
Consequently, the time t passing in the system of the crystal since the beginning of the simulation is a sum of all elementary times tn calculated after occurrence of every event: t= tn . (20) n
4. Normal crystal growth rate In the model of Gilmer and Bennema, the average growth rate of a crystal face is calculated [1,2,10] on the basis of the difference between the averaged fluxes of the added and removed blocks (I + and I − ): d + (I − I − ) N 4 d + F = ki N i − ki− N Si , N
R =
(18)
i=0
where N Fi and N Si are the average numbers of fluid and solid blocks at the interface, respectively, having i solid lateral neighbors, d is the thickness of a one block layer and N is the number of all blocks in one layer. The average values N Fi , N Si should be found as an arithmetic mean of the instantaneous values NiF , NiS calculated every preset constant time period. However, the time measure used in the model presented in Ref. [2] causes the mean values to be approximate. It is possible to use the number of events [18], the number of creations [2,9] or the number of repetitions of the main algorithm loop as a measure of elapsed time. The latter method is suitable only for the original
The algorithm described in Appendix A contains the sum of frequencies K of all the possible events, which is updated after every occurring event. Therefore, the additional cost of calculation of the growth time defined by Eqs. (19) and (20) is negligible. Using the defined time, values of N Fi and N Si can be calculated according to the definition of the mean value. However, it is equivalent but simpler to calculate the average crystal growth rate as a quotient of normal growth of crystal face and the growth time: d N , (21) tN where N is an increment in the solid block number during time t. The increment N is taken from the results of simulation. The growth rate of the face depends directly on the thickness d of one layer of growth units and indirectly on the free energy Gη of viscous flow activation appearing in the definition of the frequency ft (Eq. (3)). As the dependence of the growth rate on the coefficient β at various values of parameters α and Xs is of particular interest to us, we introduce dimensionless average growth rates defined as: r=
R =
R dft
(22a)
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
255
and r =
r , dft
(22b)
where R and r are given by Eqs. (18) and (21), respectively. The dimensionless growth rate r is proportional to the growth rate expressed in physical units, e.g., in cm/s, if the activation energy Gη of viscous flow and temperature (and hence α) are constant, while the parameter β (corresponding to the supersaturation of the solution — Eq. (6)) may vary as it does not affect the proportionality coefficient between the two growth rates. In earlier papers [1,2,12,19], the dimensionless growth rate was calculated using the expression: R =
R . dk0+
(23)
It should be noted, however, that the frequency k0+ depends on the parameter β, and this dependence appears in both the vapor growth model [1,2] and Eq. (1). Such a method of calculation of the dimensionless growth rate affects the shape of the curve of the growth rate versus β.
5. Results and discussion 5.1. The comparison of growth rates calculated with the use of different methods In previous section, we concluded that the average value of the growth rate can be calculated only approximately if the growth time is not defined explicitly. In order to compare results obtained using different methods of calculation of growth rate, we found average growth rates R using Eq. (22a) as well as r using Eq. (22b). The rate R was calculated on the basis of the arithmetic mean of NiF , NiS sampled every preset small number of events. In extreme cases, the values of the rate R were negative, while the solid block number was observed to increase, even when the numbers of instantaneous values NiF , NiS used were over a million. Fig. 1 shows that the greatest differences between the rates R and r correspond to the simulations of the growth of crystal face without defects, at the lowest supersaturations, at which mononuclear growth mechanism operates. Under these conditions, the variations in values of NiF and NiS show some correlation with
Fig. 1. The plots of quotient of the average growth rates calculated in different ways. The growth rate R was calculated using Eq. (18), S with N F i , N i obtained by averaging the instantaneous values determined every preset small number of events, while the growth rate r was calculated according to Eq. (21). Both plots are obtained for a perfect face at α = 6 and Xs = 1.5a. (a) surface array size 120 × 120; (b) surface array size 30 × 30.
the time necessary for the preset number of events to occur, and therefore it seems that during the averaging of instantaneous values of NiF and NiS , if they are determined every preset number of events, an excessive pile-up of the samples with the lowest or greatest values may occur. The discrepancy between the results given by both methods of the average growth rate determination decreases with increasing size of the sur-
256
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
face array (cf. Figs. 1(a) and 1(b)) and with increasing β. When the growth rate is determined with the use of explicit growth time given by Eq. (20), we observe smaller scatter in the points obtained for the lowest supersaturations used in the plots r (β). Our calculations show that, in the original model of Gilmer and Bennema, it is possible to achieve better accordance between growth rates R and r by the way of sampling the values NiF and NiS every preset number of repetitions of the algorithm main loop. However, this method depends on specific features of this particular algorithm. 5.2. Growth of perfect crystal face 5.2.1. Mononuclear and multiple nucleation growth mechanisms It is well known that, at high α (α > 3.5), perfect crystal faces without steps generated by spiral dislocations grow by two-dimensional (2D) nucleation mechanism, if the supersaturation exceeds a critical value (e.g., Refs. [20–23]). Therefore, for 2D nucleation growth mechanism, supersaturations slightly higher than the critical one will be called “low supersaturations” in the remaining part of Section 5.2. Close to this critical value, small changes in β cause great relative changes in growth rate (see Table 1). Our investigation shows that, at low supersaturations and high values of α (α 4), the size of the surface array strongly influences the simulation results (cf. Table 1). The MC simulations reveal that, at α 4, crystal growth may occur either by formation of one 2D nucleus (Fig. 2(a)) spreading across the surface (mononuclear model [22,23]) or by formation of many 2D nuclei (multiple 2D nucleation) spreading on the surface according to birth and spread (B + S) model [23], depending on the size of the surface array and supersaturation value (see Table 2). At higher β, nuclei are sometimes formed on the spreading nuclei and then growth occurs by nucleus above nucleus (NaN) growth mechanism (Fig. 2(b)). The influence of the surface diffusion on the observed growth mechanism is negligible. According to the mononuclear growth model [22, 23], nucleus born on the surface spreads out and covers the entire surface before other critical size nuclei are born on this surface. Therefore, it is possible for the mononuclear growth mechanism to operate [20,22]
Table 1 The dependence of the average crystal growth rate r on β for α = 6, Xs = 1a and two different array sizes β
r Array 30 × 30
Array 120 × 120
0.35
0.00002
0.00038
0.38
0.00011
0.00133
0.40
0.00015
0.00176
0.42
0.00020
0.00309
0.45
0.00029
0.00463
0.50
0.00251
0.00711
0.55
0.00325
0.00997
0.60
0.00746
0.01298
0.65
0.00757
0.01732
0.70
0.02026
0.02463
if supersaturation and the area of the surface are sufficiently small. In crystal growth practice, this mechanism is hardly ever observed, because areas of real crystal surfaces are much too large. In our MC simulations performed for relatively low supersaturations and small array sizes, we did observe (Fig. 2(a)) this growth mechanism. With decreasing value of β, the number of coexisting nuclei decreases sharply. For the lowest values of β at which the face can still grow (at given α > 4), we observed single nucleus spreading across whole surface. The values of β, at which growth proceeds by mononuclear mechanism, decrease with increasing size of the surface array (cf. Table 2). It was shown [20,22] that mononuclear growth mechanism operates if the area S of crystal face fulfill the condition: S 1/2 < Lcr = (νv /J )1/3,
(24)
where Lcr is the critical length of face side, below which growth by mononuclear growth mechanism can occur, νv is the spread velocity of nucleated steps [22, 23] νv = D1 σ,
(25)
and J denotes the rate of 2D nucleus formation [22, 23] π γ 2 1/2 J = D2 β exp − . (26) β kT
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
257
Fig. 2. Growth by 2D nucleation mechanism at perfect crystal face (surface array size 120 × 120). (a) Mononuclear growth mechanism, α = 5, β = 0.2, Xs = 0; (b) NaN mechanism, α = 6, β = 0.7, Xs = 0. Table 2 The characteristic values β ∗ for which the transition from mononuclear to B + S growth mechanism is observed. The results are obtained for diffusion distance Xs = 1a α
Size of surface array 20 × 20
30 × 30
60 × 60
120 × 120
4
0.07
0.06 . . . 0.07
0.05
0.04 . . . 0.05
5
0.41
0.35
0.25
0.20
6
0.77
7
240 × 240
0.16
0.60
0.49
0.42
0.38
0.90
0.78
0.70
0.60
Here D1 and D2 are parameters independent of σ , β and α. The edge free energy γ is related to the parameter α [24]: γ α ≈ . (27) kT 4 After substituting Eqs. (25), (26) and (27) into Eq. (24) we obtain
S 1/2 < Lcr =
D1 D2
1/3
πα 2 . (28) σ 1/3 β −1/6 exp + 48β
Eq. (28) reveals that Lcr is very sensitive to changes in values of α and β, especially at low values of β. For example, it follows from Eq. (28) that at β = 0.1, the change in α from 4 to 5 causes an increase in Lcr by a factor of 2.7 × 103 . Furthermore, the change
258
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
in β from 0.04 to 0.05 causes reduction in Lcr by a factor 181, at α = 4. In our MC simulations, we did not observe such a strong influence of α and β on Lcr (cf. Table 2), at low β. At higher β (β > 0.5), as Lcr given by Eq. (28) is not so sensitive to changes in β and α, MC simulation results are in better accordance with Eq. (28). Mononuclear model predicts a linear dependence of the growth rate on the total area of crystal surface [23]. We found that, at given value of β, the growth rate increases with increasing area of the surface array, but it is rather difficult to describe this dependence by a linear function. However, we should remember that the theory assumes an infinite spread velocity of nucleus. Such an assumption is very close to the case of very small both the surface array size ( 60 × 60) and supersaturation values. Fig. 3 illustrates the dependence of the square root of the average growth rate r on the length of the array side. This figure shows that the dependence of the growth rate on the surface area is nearly linear at low supersaturations and for small surface arrays, where mononuclear growth mechanism operates and the time necessary to fill the whole layer by spreading nucleus is negligible. In Fig. 3, the curve obtained for β = 0.38 corresponds to mononuclear growth mechanism for the surface array sizes below 240 × 240 (cf. Table 2). It is seen in this figure (curve for β = 0.38) that, for the surface arrays of sizes over 120 × 120, the growth rate increases very slowly with increasing array size.
Fig. 3. The dependence of dimensionless growth rate r on the size of the surface array, for a perfect (001) face growing by 2D nucleation, α = 6, Xs = 1a.
In this range of array sizes, the spread time of nucleus cannot be neglected as it is assumed in mononuclear model. For the smallest arrays, the points corresponding to β = 0.38 (Fig. 3) show that, √ with increasing length of the surface array side, r increases a little slower than proportionally. This small discrepancy between mononuclear model and MC simulation results also appears for various values of mean diffusion distance Xs and is a result of too small surface array size, which is comparable to the critical size of twodimensional nuclei. Moreover, in the MC simulations, we used cyclic connection of the surface array sides, which is not used in theoretical model. For values of β high enough for the growth mechanism described by B + S model to operate, the growth rate increases slightly with increasing array size, tending to a constant value. However, it should be pointed out that the characteristic value of β (denoted β ∗ ) corresponding to the transition from one to the other growth mechanism increases sharply with increasing value of the parameter α (cf. Table 2). Although mononuclear and B + S growth mechanisms can be distinguished directly by observations of crystal surfaces, the plots of the number of grown layers versus growth time t (given by Eq. (20)) may also be useful. For the lowest values of β still enabling crystal growth, the curves have a form of irregular sharp steps (Fig. 4(a)), where the horizontal lines correspond to waiting for the formation of 2D nucleus and the nearly vertical lines correspond to rapid spreading of the nucleus. For slightly higher values of β (but still lower than the characteristic value β ∗ ) and for sufficiently large size of the surface array, the spread time may become visible in a plot (e.g., Fig. 4(b)). If the growth proceeds by B + S mechanism, the plotted curve is a diagonal line with random fluctuations (Fig. 4(c)). 5.2.2. Influence of surface diffusion on growth rate and surface morphology In the case of vapor growth, surface diffusion plays an essential role and it considerably increases the nucleation rate J (even by a factor of 102 ÷ 103 [25]). Consequently, crystal growth rate also increases considerably (e.g., Ref. [1]). It is known that in the case of solution growth, the influence of surface diffusion on the nucleation rate is not so pronounced [22]. In our simulations performed for
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
259
Fig. 5. The dependence of the dimensionless growth rate r on mean surface diffusion distance Xs plotted for a perfect face, α = 4, surface array size 120 × 120 and several values of β: (a) 0.05 β 0.10; (b) 0.5 β 0.8.
Fig. 4. The number of grown layers versus dimensionless growth time, for α = 6, Xs = 1a, surface array size 120 × 120.
α = 4 and several values of β. It is seen in this figure that the growth rate r increases moderately with increasing Xs /a. This effect can be explained by an increase in the spread velocity ν of 2D nuclei. According to B + S model [22,23] the normal growth rate R can be written: R = dν 2/3 J 1/3 ,
solution growth, we did not notice an increase in the nucleation rate with increasing mean diffusion distance Xs . Fig. 5 illustrates the dependence of dimensionless growth rate r on the mean diffusion distance Xs for
(29)
where d is the height of a nucleus and J is the nucleation rate given by Eq. (26). Since our investigation reveals that addition of growth units to the critical nucleus occurs by both direct impingement from the solution and surface diffusion, we assume that the total
260
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
spread velocity ν of a nucleus may be expressed as the sum of spread velocities νv and νs corresponding to direct incorporation and surface diffusion, respectively: ν = νv + νs .
(30)
The spread velocity νs , proportional to both Xs [23] and νv , can be expressed by Xs νv , (31) a where νv is given by Eq. (25) and coefficient A is assumed to be dependent on α and β. Upon substituting Eqs. (30) and (31) into Eq. (29), we find Xs 2/3 R = R0 1 + A , (32) a
νs = A
where R0 is the growth rate at Xs = 0. Therefore, we fitted our results (Fig. 5) according to the following equation Xs 2/3 r = r0 1 + A , (33) a and using the least-square method, we found coefficients A and r0 . The calculations performed for α = 4 reveal that for the lowest value of β (β = 0.07), coefficient A is equal to 0.39 and it decreases considerably with increasing β (even to 0.03 for β = 0.8). This means that at higher β, the spread velocity is determined by direct incorporation of growth units from the solution rather than by surface diffusion. At higher β, edges of nuclei are rough which should cause an increase in the spread velocity ν. On the other hand, however, surface diffusion reduces edge roughening.
These conflicting tendencies may cause reduction in the coefficient A at higher values of β. Although smoothing the edges of nuclei by surface diffusion was clearly visible, smoothing the crystal surface was not noticed. Similar calculations performed for α = 5 and α = 6 show that the coefficient A slightly increases with increasing α. This result indicates that the importance of surface diffusion in solution growth of perfect crystal face is greater for smooth surfaces, which seems to be natural. 5.3. Surfaces with a step We performed simulations for the (001) face containing a single step of one growth unit height. The step distance was assumed to be independent of supersaturation and equal to the length of the surface array side. For every array size and for various values of the parameter α 4, there is a range of supersaturations for which combined growth mechanism is observed: the 2D nucleation and the movement of the growth step (Fig. 6). This combined growth mechanism was also observed in MC simulations performed for vapor growth [1,19] and it is observed experimentally (e.g., Refs. [26–28]) for solution growth. The values of β corresponding to the combined growth mechanism decrease with increasing array size (and hence with increasing interstep distance). Therefore, the observed growth mechanism depends not only on α and β but also on the assumed interstep distance, especially at small sizes of the surface array.
Fig. 6. Combined growth mechanism: 2D nucleation and the growth by step movement, α = 5, β = 0.25, Xs = 1a, surface array size 120 × 120.
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
At high supersaturations, where 2D nucleation growth mechanism dominates, the dependence of R on β is clearly stronger than linear. The critical value of the parameter β, corresponding to the change from a dominant spiral growth mechanism in a dominant 2D nucleation mechanism increases sharply with increasing parameter α (Fig. 7). For α = 5 and the surface array side 60 × 60 (Fig. 7(a)), the critical value β is approximately equal to 0.25, while for α = 8, it is higher than 0.8 and hence Fig. 7(b) shows the dependence r (β) for dominating spiral growth mechanism. According to the surface diffusion theory of Burton, Cabrera and Frank (BCF) [15,23], for low supersaturations, at which only spiral growth mechanism operates, the dependence of the growth rate on β is linear if constant (i.e. independent of β ) interstep distance is assumed. The results of our MC simulations are in
261
good agreement with the theoretical predictions even for the lowest values of β, although some small deviation of the r = f (β) dependence from linearity appears in most simulations (Fig. 7(a) for β < 0.20 and Fig. 7(b)). This deviation is dependent on the mean diffusion distance Xs and the value of parameter α and it is a result of superposition of some effects influencing the growth rate: • Increasing β stimulates kinetic roughening of the step edge. This effect causes that, with increasing β, the average growth rate r increases a little faster than linearly (Fig. 7(b), curve Xs = 0). • Smoothing the step edge by surface diffusion. This effect is most pronounced for not very large values of β. With increasing mean diffusion distance Xs , the diffusion smoothing reduces the importance of kinetic roughening of the step edge, and in consequence, dependence r on β becomes more linear (Fig. 7(b)). The smoothing of the surface of step terraces by surface diffusion we observed, was not so distinct as the one reported in Refs. [1,5] concerning vapor growth. • With increasing β the number of blocks being added directly to the step edges increases faster than the number of blocks being added at any other location at the surface and subsequently transported to the edges. Therefore, the relative changes in growth rate dependence on diffusion distance Xs are the greatest for low values of β. Our simulations performed for solution growth show that in general, surface diffusion increases the growth rate, but this influence is clearly weaker than in the case of simulations performed for vapor growth [1, 19]. Moreover, for vapor growth, the critical value of β (corresponding to the change of the growth mechanism) strongly depends [19] on surface diffusion. According to our results, this dependence is insignificant in the case of solution growth (Fig. 7(a)).
6. Conclusions
Fig. 7. Plots of dependence of the growth rate r of the crystal face with single growth step on the parameter β for Xs = 3a (series ), Xs = 1a ( ) and Xs = 0 ( ). Surface array size is 60 × 60. (a) α = 5; (b) α = 8.
The following conclusions can be drawn from the present study: (1) At high supersaturations, the average crystal growth rate can be calculated with the use of time measure approximated as a number of events [18], but for a perfect face growing at low supersatura-
262
(2)
(3)
(4)
(5)
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
tions, this method gives large scatter in results and in extreme cases, the sign of the average growth rate may even be inconsistent with the observed growth of the crystal face. The growth time defined by Eq. (20) and the average growth rate (21), based on this definition of time, can be used for all values of α and β. The characteristic value β ∗ (and hence supersaturation σ ∗ ), at which a transition from mononuclear to B + S growth mechanism occurs, decreases with increasing size of the array representing crystal surface. According to the theory [22], small changes in β and α values should cause enormous changes in the critical length Lcr of the face, especially for not very high values of α (e.g., α = 4). In our MC simulations, we did not observe such a strong influence of α and β on Lcr . The mononuclear model predicts linear dependence of growth rate on total area of crystal surface but, according to our results, this dependence is not linear. We have found that, for large sizes of the surface array and values of β near β ∗ , the spreading time of 2D nucleus cannot be neglected as it is assumed in the theory. For higher values of β, where B + S mechanism operates, the average crystal growth rate slightly increases with increasing array size, tending to a constant value. Our MC simulations results indicate that for solution growth, surface diffusion affects the spread velocity of nuclei rather than the nucleation rate. The influence of surface diffusion on the growth rate is not so strong as in the case of vapor growth and it decreases with increasing parameter β. MC simulations performed for a face with single step show that the type of the dominating growth mechanism (2D nucleation or the growth step movement) depends considerably on parameters α and β and also on the size of the surface array. In contrast to the results reported for vapor growth, surface diffusion does not cause change of the growth mechanism.
Appendix A. Simulation algorithm In Section 3, we described the main assumptions for our implementation of fast algorithm of MC simulations. Here we discuss in detail the realization of this algorithm in the case of block growth of the (001) face. In this algorithm, three kinds of elementary events occurring on the crystal surface are taken into account: addition of a block to the surface (creation), removal of a block (annihilation) and surface diffusion. The frequencies of the events are defined as follows: (1) The frequency of any event on the phase boundary in a location with coordinates (x, y): − + kj+(x,y) + Kxy = ki(x,y)
4
ki(x,y),j (x ,y ) , (A.1)
q=1
where ki− is the annihilation frequency in a location with i solid lateral neighbors, kj+ is the creation frequency in a location with j solid lateral neighbors, kij is the surface diffusion frequency occurring from a location with i lateral neighbors onto a location with j lateral neighbors, q is the diffusion direction, and (x , y ) denote coordinates of a diffusion target location calculated as x = x + sin(πq/2) and y = y + cos(πq/2). (2) The frequency of any event at interface in column x of growth units: Kxy . (A.2) Kx = y
(3) The frequency of any event occurring in any location is defined as a sum of frequencies Kxy over all elements of the crystal surface array K= Kxy . (A.3) x
y
If we divide the linear search into 2 levels, first finding the column x of the crystal surface array, designated by hxy , and next finding the row y, the simulation will be realized according to the following algorithm: 1. Preparation of the crystal surface array hxy corresponding to the boundary of the solid and fluid phases, preparation of the arrays of frequencies ki− , ki+ , kij , Kxy , Kx , calculation of the sum K and setting the time counter t to zero.
M. Rak et al. / Computer Physics Communications 138 (2001) 250–263
2. Generation of a random number RND in the interval 0 . . . K with uniform distribution. 3. Finding the column number: x := 0 until (x < xmax and Kx RND) RND := RND − Kx x := x + 1
References [1] [2] [3] [4] [5] [6] [7]
4. Finding the row number: y := 0 until (y < ymax and Kxy RND) RND := RND − Kxy
[8] [9] [10]
y := y + 1 5. Selection of the event:
[11]
+ if (RND < ki(x,y) ) ⇒ creation
[12]
+ − + ki(x,y) ) otherwise, if (RND < ki(x,y)
[13] [14] [15]
⇒ annihilation otherwise ⇒ diffusion. The diffusion direction q is being selected as a number of one of the intervals with lengths + ki(x,y),j (x ,y ) , containing the number RND − ki(x,y) − − ki(x,y) . 6. Incrementing and/or decrementing elements of the crystal surface array hxy (and optionally hx y ) according to the event selected in point 5. 7. Incrementing the time counter: t := t + 1/K. 8. Updating the arrays Kxy , Kx and the sum K: for the coordinates (a, b) equal to (x, y), and (x , y ) in the case of diffusion, for the closest neighboring locations and for the locations neighboring to the closest neighbors the following operations should be performed: K := K − Kab , Ka := Ka − Kab , calculate new frequency Kab , K := K + Kab , Ka := Ka + Kab . 9. jump to point 2.
263
[16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28]
G.H. Gilmer, P. Bennema, J. App. Phys. 43 (1972) 1347. G.H. Gilmer, P. Bennema, J. Crystal Growth 13/14 (1972) 148. M. Kotrla, Comput. Phys. Commun. 97 (1996) 82. A.C. Levi, M. Kotrla, J. Phys.: Condens. Matter. 9 (1997) 299. R.-F. Xiao, J.I.D. Alexander, F. Rosenberger, J. Crystal Growth 109 (1991) 43. L. Guang-Zhao, J.P. van der Eerden, P. Bennema, J. Crystal Growth 58 (1982) 152. J. Karniewicz, A. Brozi, M. Rak, P. Pietrucha, SPIE Proc. 2373 (1995) 150. W.J.P. van Enckevort, A.C.J.F. van der Berg, K.B.G. Kreuwel, A.J. Derksen, M.S. Couto, J. Crystal Growth 166 (1996) 156. W.J.P. van Enckevort, A.C.J.F. van der Berg, J. Crystal Growth 183 (1998) 441. J.I.D. Alexander, in: J.P. van der Eerden, O.S.L. Bruinsma (Eds.), Science and Technology of Crystal Growth, Kluwer Academic Publishers, Dordrecht, 1995, p. 81. P. Bennema, J.P. van der Eerden, J. Crystal Growth 42 (1977) 201. P. Bennema, in: D.T.J. Hurle (Ed.), Handbook of Crystal Growth, Vol. 1a, North-Holland, Amsterdam, 1993, p. 477. M. Izdebski, M. Rak, A. Brozi, SPIE Proc. 3178 (1996) 135. F.L. Binsbergen, Kolloid Z. 237 (1970) 289. P. Bennema, G.H. Gilmer, in: P. Hartman (Ed.), Crystal Growth: An Introduction, North-Holland Publ., Amsterdam, 1973, p. 263. I. Prigogine, R. Defay, Chemical Thermodynamics, Longmans Green and Co., 1954. I. Prigogine, A. Bellemans, V. Mathot, The Molecular Theory of Solutions, North-Holland Publ., Amsterdam, 1957. F.F. Abraham, M. White, J. App. Phys. 41 (1970) 1841. P. Bennema, J. Crystal Growth 69 (1984) 182. A.I. Malkin, A.A. Chernov, I.V. Alexeev, J. Crystal Growth 97 (1989) 765. A.A. Chernov, Prog. Crystal Growth and Characterization 26 (1993) 121. A.A. Chernov, Crystal Growth, Crystallization Processes, Modern Crystallography, Vol. III, Springer, Berlin, 1984. M. Ohara, R.C. Reid, Modeling Crystal Growth Rates from Solution, Prentice-Hall Inc., Englewood Cliffs, NJ, 1973. P. Bennema, J. Boon, C. van Leeuwen, G.H. Gilmer, Krist. Tech. 27 (1973) 659. G.H. Pound, M.T. Simnad, L. Yang, J. Chem. Phys. 22 (1954) 1215. J.J. De Yoreo, T.A. Land, B. Dair, Phys. Rev. Lett. 73 (1994) 838. C.M. Pina, D. Bosbach, M. Prieto, A. Putnis, J. Crystal Growth 187 (1998) 119. C.M. Pina, U. Becker, P. Risthous, D. Bosbach, A. Putnis, Nature 395 (1998) 483.