- Email: [email protected]
Nucleus size and Zeldovich factor in two-dimensional nucleation at the Kossel crystal (0 0 1) surface
Surface Science 574 (2005) 77–88 www.elsevier.com/locate/susc
Nucleus size and Zeldovich factor in two-dimensional nucleation at the Kossel crystal (0 0 1) surface J.H. ter Horst *, P.J. Jansens Laboratory for Process Equipment, Delft University of Technology, Leeghwaterstraat 44, 2628 CA Delft, The Netherlands Received 19 May 2004; accepted for publication 14 October 2004 Available online 2 November 2004
Abstract A Monte Carlo simulation method to determine the growth probability of two-dimensional clusters on the Kossel (0 0 1) surface is described. From this growth probability the nucleus size and the Zeldovich factor for two-dimensional nucleation were obtained. For low supersaturations the Gibbs–Thomson equation describes the supersaturation dependency of the nucleus size very well. At higher supersaturations a deviation from the predicted nucleus sizes is observed. This is mainly due to the size dependence of the apparent specific step free energy and/or the shape factor of the nucleus. A small systematic deviation between the theoretical and actual growth probability is observed and explained. 2004 Elsevier B.V. All rights reserved. Keywords: Nucleation; Crystallization; Computer simulations; Monte Carlo simulations
1. Introduction It is generally not possible to predict sufficiently accurate nucleus sizes and nucleation rates for complex crystalline structures using purely theoretical considerations. For these complex crystal structures accurate step free energy (or in the case
* Corresponding author. Tel.: +31 15 278 6678; fax: +31 15 278 6975. E-mail address: [email protected] (J.H. ter Horst).
of 3D nucleation the interfacial energy) predictions and shape factor predictions are needed in order to calculate the nucleus size and nucleation rate correctly. Molecular simulations could aid to improve the accuracy in predictions by simulating the nucleation process. Molecular simulations can furthermore be used to improve the fundamental understanding of the nucleation process. In this paper a newly proposed general method to determine the nucleus size and the Zeldovich factor will be validated by using it on a simple model system. The 2D cluster formation on the Kossel crystal (0 0 1) surface was chosen as the
0039-6028/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.susc.2004.10.020
78
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
model system. For the Kossel model system relatively accurate values of the shape factor and the specific step free energy can be determined. As a first objective, this paper is to demonstrate that this method determines the nucleus size and Zeldovich factor. This paper is furthermore meant to find the origins of the difference between molecular simulation results and Classical Nucleation Theory (CNT) predictions thus studying the limitations of the CNT. The formation of a two-dimensional cluster consisting of n growth units (GU) results in (1) a free energy gain due to the supersaturation in the ambient phase and (2) an excess of free energy due to the creation of a (closed) step. The free energy gain is linearly related to the cluster size n and the supersaturation while the excess free energy can be seen as a total step free energy that is a function of the length L of the cluster periphery and an apparent specific step free energy j. The work W(n) to form a n-sized cluster then is written as [1]: W ðnÞ ¼ n Dl þ Lj ¼ n Dl þ cajn1=2
ð1Þ
where c is a shape factor, a is the distance between nearest neighbours, and Dl is the supersaturation, which is the difference in chemical potential between the ambient phase and the crystalline phase. The j is apparent because the actual specific step free energy may not be isotropic. At large cluster sizes the energy gain is dominant while at small sizes the total step free energy is larger. Therefore a cluster size exists for which the work to form it is maximal. This cluster size is the nucleus size n = n* and the corresponding energy barrier is the nucleation work W*. The nucleus is in unstable equilibrium with the ambient phase. Clusters smaller than the nucleus size tend to decay while clusters larger than the nucleus size tend to grow out. In CNT it is assumed that the apparent specific step free energy and the shape factor are independent of the size n and the supersaturation Dl and that they have the values j = j1 and c = c1 of an equilibrium shaped and macroscopically large 2D nucleus. Setting the derivative of the work W(n) to the cluster size n to zero results in the 2D variant of the classical Gibbs–Thomson equation [1] with aid of which the nucleus size nCNT can be predicted:
nCNT ¼
c21 a2 j21 4 Dl2
ð2Þ
The assumptions on the shape factor and the apparent step free energy make Eq. (2) approximate. It is not clear down to what nucleus sizes and up to what supersaturations these approximations are acceptable. One of the purposes of the simulation study described in this paper is to check the validity of Eq. (2) by comparing the nucleus size from simulations and the nucleus size determined from Eq. (2). First a recently proposed equation for cluster growth probabilities [2] is introduced (Section 2). These cluster growth probabilities can be determined with aid of computer simulation techniques. This approach is followed for 2D nucleation on the Kossel crystal (0 0 1) surface by performing kinetic Monte Carlo simulations (Section 3). Theoretical estimates for the shape factor c1 and the apparent specific step free energy j1 to be used in the CNT predictions are given and discussed in Section 4. The cluster growth probabilities from simulations result in the nucleus size and the Zeldovich factor at a certain supersaturation. In Section 5 the simulation results are presented and discussed after which in Section 6 two important aspects of the growth probability method and the kinetic theory of nucleation are discussed.
2. Growth probability of n-sized clusters Nucleation is known to occur by formation of molecular (or atomic) clusters of various sizes. The cluster size changes randomly as a result of successive attachments and detachments of single growth units (GU: molecules, atoms, ions) to and from the cluster. The cluster of size n = n* is the nucleus. It is that particular cluster to and from which molecules are attached and detached with equal frequency. Clusters of size n < n* tend to decay, because per unit time less molecules are attached to than detached from them. On the contrary, the clusters of size n > n* are characterized by attachment frequencies greater than the detachment ones and for that reason these clusters tend to grow up to macroscopic sizes. As mole-
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
cular attachments and detachments are random events, a given n-sized cluster can grow and reach a macroscopic size only with a certain growth probability P(n). Using the kinetic theory of nucleation it turns out that the growth probability P(n) at the expense of certain approximations but without loss of generality can be expressed in terms of a numerical factor b and the nucleus size n* [2]: P ðnÞ ¼ 12½1 þ erfðbðn n ÞÞ
ð3Þ
where erf is the error function. The growth probability expression is valid for any kind of one-component nucleation (two- and three-dimensional, homogeneous and heterogeneous nucleation of vapours, liquids or solids) [2]. The numerical factor b is related to the width D* of the nucleus region and to the Zeldovich factor z [1]: b¼
p1=2 ¼ p1=2 z D
ð4Þ
The Zeldovich factor is an important parameter in the pre-exponential kinetic factor of the nucleation rate equation. Outside the nucleus region for n < n* D*/2 the growth probability P(n) approaches 0 while for n > n* + D*/2 the value of P(n) approaches 1. This paper will show that the growth probability expression is valuable to determine nucleation characteristics with aid of computer simulations. The growth probability as a function of initial cluster size at a certain supersaturation can be determined with aid of computer simulations. From the simulated growth probability P(n) as a function of the cluster size at a certain supersaturation the nucleus size n* and the numerical factor b can be determined through Eq. (3). This has been done for a liquid droplet in a supersaturated vapour Lennard–Jones compound [2]. In this paper it is focussed on the 2D cluster formation on the Kossel (0 0 1) surface as a first step towards 3D nucleation simulations of e.g., polymorphic crystalline compounds. Such 3D nucleation simulations lead to the experimentally advantageous prediction of the supersaturations under which concomitant nucleation can be avoided.
79
3. Simulation details It was chosen to study the growth probability of 2D clusters on the Kossel (0 0 1) surface. A Kossel crystal has an idealized simple cubic crystal structure with a lattice constant a. Both crystal and ambient phase, either vapour, melt or solution, are at temperature T. The ambient phase has the same lattice structure as the crystalline compound. The bond energy ess refers to the interaction between growth units (GU) in the crystalline compound while the energy ell refers to the interaction between ambient phase units. The bond energy esl is the interaction energy between a GU and an ambient phase unit. Now the overall lateral bond strength / can be defined [3]: / ess þ ell 2esl ¼ kT 2kT
ð5Þ
This is the energy loss per growth unit upon the breaking of a solid–solid bond which coincides with the breaking of an ambient phase-ambient phase bond and the creation of 2 solid-ambient phase bonds. The growth probability of a mono-atomically thick 2D cluster consisting of n GU on such a Kossel (0 0 1) surface can be studied with the n-fold way Monte Carlo (MC) algorithm [4]. Values of / = 1.33, 2 and 5kT were chosen in the simulations. An increase in bond strength / is equivalent to a decrease in temperature. At T = 0 K or an infinitely large bond strength the (0 0 1) surface is atomically smooth. Lowering the bond strength or increasing the temperature will lead to an increase of adatoms, vacancies and kinks in the steps. At a relatively low bond strength of / = 1.33kT (corresponding to a relatively high temperature) a macroscopically large nucleus on the (0 0 1) surface is circular while at a relatively high bond strength of / = 5kT (relatively low temperature) the macroscopically large nucleus has a shape between a square and a circle (see Section 4). Each MC step a GU is either created or annihilated while surface diffusion is neglected. There are eight types of transitions: creation and annihilation of a solid GU with either 1, 2, 3 or 4 lateral solid neighbours (the creation of GU without
80
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
lateral neighbours was not considered in the simulations: only the cluster grows or decays in the simulation). Then k þ i can be defined as the creation rate of GU with i neighbours and k i as the annihilation rate of GU with i neighbours. To ensure microscopic reversibility [5], also known as the detailed balance condition in molecular simulation theory [6], the quotient of the rates of these transitions must be equal to the quotient of the corresponding Boltzmann factors: kþ U ðiþ Þ U ði Þ i ¼ exp kT k i / Dl ¼ exp 2ði 2Þ þ ð6Þ kT kT Here U(i+) is the energy of the system after a GU with i neighbours is created, U(i) is the energy of the system after a GU with i neighbours is removed, U(i+) U(i) is the energy difference between a filled and empty solid position with i neighbours and Dl is the supersaturation. Microscopic reversibility leaves either k þ i or k i free to choose in a simulation. The creation rate þ kþ can be chosen to be independent of the i ¼ k number of neighbours i: Dl k þ ¼ k eq exp ð7Þ kT The rate keq is the single GU addition rate per site at equilibrium between the crystal solid phase and the ambient phase. Physically this corresponds to a random impingement of GU on the cluster. From Eqs. (6) and (7) it then follows that the annihilation rate k i for a GU is independent of the supersaturation and given by: / k i ¼ k eq exp 2ði 2Þ ð8Þ kT The defined rates ensure that at a supersaturation of Dl = 0 the Kossel (0 0 1) surface is in equilibrium and thus does not grow or dissolve. The transition rates k+ and k i are proportional to the transition probabilities p+ and p i to be used in the MC simulation. This MC scheme is also known as the random rain model [7]. Only the ambient phase changes as compared to equilibrium conditions when a supersaturation is applied,
typically the case for solution and vapour growth [7]. Due to the definition of the system of a developing two-dimensional cluster there are a limited number of sites N where a GU can be either created or annihilated. Each of these N positions corresponds to a translation probability (p+ for a crea tion position and p or p for an 1 , p2 , p3 4 annihilation position). Each MC step a GU is either created or annihilated which means that the summation over all N translation probabilities is unity. Through a random number in the range [0, 1) the translation site is picked from these N sites. The growth probability P(n) of an n-sized cluster is determined by performing a number M of MC simulations (usually M = 300–10 000) starting with the same initial cluster size n. A minimum nmin P 1 and a maximum nmax cluster size were chosen for which the values P(nmin) = 0 and P(nmax) = 1 were sufficiently accurate approximations. The simulation was stopped only when the cluster size reached either n = nmin or n = nmax. The number M+ of simulations resulting in a cluster of size n = nmax gives the growth probability: P ðnÞ ¼
Mþ M
ð9Þ
Rectangular initial cluster shapes with minimum step energy were used. The initial clusters used in the MC simulations all had a rectangular shape and a size of 2 · 1, 2 · 2, 3 · 2, 3 · 3, 4 · 3, 4 · 4, 5 · 4, . . ., 25 · 24, 25 · 25. In addition, at / = 2kT and Dl = 0.3kT apart from the 10 · 10 sized and thus square shaped initial cluster, elongated initial shapes of 4 · 25 and 5 · 20 sized initial clusters were used. The elongated initial shapes did not have a significant effect on the calculated growth probability of that cluster size n = 100. This indicates that the shape equilibration process is generally much faster than the growth and decay processes, validating the choice of as square as possible rectangular initial clusters. Furthermore in case of a small nucleus size the shape factor is mainly determined by the lowest energy cluster configuration (square) because of the relatively large energy difference between the low and higher cluster configurations for small clusters.
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
4. The nucleus size using the classical nucleation theory The nucleus size is that cluster size for which the attachment frequency and the detachment frequency of growth units are equal. The nucleus is therefore in (unstable) equilibrium with the supersaturated ambient phase and will thus have the equilibrium shape. In the case of a macroscopically large nucleus on the Kossel crystal (0 0 1) surface the exact form of the specific step free energy j(//kT, h) is known [8]: ajð/=kT ,hÞ ¼ j cos hjsinh1 ðaj cos hjÞ kT ð10Þ þ j sin hjsinh1 ðaj sin hjÞ where " #1=2 2 1 b2 a¼ ð11Þ b 1 þ ðsin2 2h þ b2 cos2 2hÞ1=2 and
/ 2 sinh kT b¼ ð12Þ / 2 cosh kT Eq. (10) reflects that for anisotropic lattices the actual specific step free energy j = j(//kT, h) changes with the bond strength //kT and the orientation angle h from the [1 0 0] direction. For a specific bond strength a polar diagram of the actual specific step free energy j(//kT, h) can be made. The equilibrium shape is that shape for which the perpendicular distance of all steps to the nucleus centre is proportional to the actual specific step free energy j(//kT, h) of that step (Gibbs–Curie–Wulff rule) [9]. Since the actual specific step free energy j(//kT, h) is known from Eq. (10) the equilibrium shape can be constructed which was done numerically. From the equilibrium shape construction then the surface area (which is a relative nucleus size n1 ) and the periphery length (which is a relative nucleus periphery length L1 ) are known. Then a good approximation of the shape factor c1 of the nucleus to be used in the Gibbs–Thomson equation is given by: L ð13Þ c1 ¼ p1 ffiffiffiffiffiffi a n1
81
Here for L1 the periphery length of the calculated equilibrium shape and for n1 the surface area of the calculated equilibrium shape is taken. The shape factor is calculated to be c1 = 3.545, 3.554 and 3.681 for respectively / = 1.33, 2 and 5kT. At / = 1.33kT the nucleus is near perfectly spherical (c1 2p1/2). The shape factor at a bond strength of / = 5kT indicates a rather spherical nucleus (c1 = 3.681 compared to c = 4 for a square shape at an infinitely large bond strength or a temperatures of 0 K). In classical nucleation theory the total step free energy of a nucleus is approximated by the value 1=2 c1 aj1 n1 . The total step free energy of the equilibrium shape can be calculated because at each infinitely small length dL of the periphery of the equilibrium shape the specific step free energy j(//kT, L) is known from the polar diagram and the equilibrium shape construction. The total step free energy of a nucleus in classical nucleation theory and of a macroscopically large nucleus at equilibrium are equal if for n1 the surface area of the calculated equilibrium shape is taken: Z L 1=2 c1 aj1 n1 ¼ jð/=kT ,LÞ dL ð14Þ 0
In Eq. (14) only the apparent specific step free energy j1 of the nucleus is unknown and thus can be determined, which was done numerically. In the case of bond strength / = 1.33, 2 and 5kT this leads to an apparent specific step free energy of respectively j1 = 0.793, 1.768 and 5.244kT/a. These values can be compared to the specific step free energy of an infinitely long step [10] j = 0.788, 1.728 and 4.993kT/a for a bond strength of respectively / = 1.33, 2 and 5kT. The difference between the apparent specific step free energy of the nucleus and the specific step free energy of an infinitely long step reflects the effect of the corners of the nucleus. A perfectly spherical or perfectly square shaped nucleus has an apparent specific step free energy equal to the specific step free energy of an infinitely long step. A shape between these two extremes results in an apparent specific step free energy which is slightly larger than the specific step free energy of an infinitely large step, e.g., j1 = 5.244kT/a compared to j = 4.993kT/a
82
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
Table 1 The shape factor c1, and the specific step free energy j1 as a function of the bond strength //kT calculated with the method described in Section 4 //kT
c1
aj1/kT
1.33 2 5
3.545 3.554 3.681
0.793 1.768 5.244
for / = 5kT due to the higher specific step free energy j(//kT, L) at the corners. In Table 1 the obtained values for the shape factor c1 and step free energy j1 are summarized. These values can be used to determine the nucleus size nCNT as a function of the supersaturation with aid of the Gibbs–Thomson equation. It should be noted however that there are some differences between the conditions of the simulations and the theory used to calculate c1 and j1. First of all, the values in Table 1 relate to macroscopically large nuclei while in the simulations the clusters have a finite size. It is expected that the finite size of the nucleus increases the shape factor c compared to c1 because less nucleus configurations are possible and the relative difference between the lowest and the higher energy configurations of the nucleus is larger. The nucleus shape therefore becomes more square-like at smaller sizes. The same argument leads to a lower entropy of the system from which it follows that the apparent specific step free energy j1 increases. A cluster would have a larger total step free energy in the simulations than is predicted and therefore the nucleus size n* will be larger than predicted with the classical nucleation theory. It is expected that this effect becomes larger at smaller nucleus sizes. The second difference between simulations and theory is the prevailing supersaturation in the simulations. The supersaturation decreases the specific step free energy [11] and therefore possibly also the apparent specific step free energy and the shape factor of the nucleus. The third difference between simulation and theory is that within the simulations there is the possibility that vacancies in the interior of the nucleus appear. This is an entropy effect causing a roughening of the surface. This would mean that the total step free energy is lower. During the
simulations such vacancies were rarely observed. It can be said that this may cause a difference for the value from simulations with the theoretical value in the extrapolation towards macroscopically large nuclei (zero supersaturation). Even at zero supersaturation on a flat surface vacancies in that surface will be created and annihilated continuously. This effect would be higher for lower bond strengths. It is therefore expected that for large nuclei and small supersaturations (where the effect of the first two differences is not present) still a small difference between the predicted nucleus size and the nucleus size from simulations may be present. This difference then is larger for smaller bond strengths.
5. The nucleus size from the cluster growth probabilities A cluster much larger than the nucleus size always survives and grows out to become a new growth layer while a cluster much smaller than the nucleus size has an infinitely small growth probability. A cluster of a size within the nucleus region has a certain probability 0 < P(n) < 1: While some clusters of a certain size grow out to macroscopic sizes, others will decay. This is observed in subsequent MC simulations of clusters with an initial size n. Fig. 1 shows the development of the cluster size during three simulations for an initial cluster size of n = 144 at a supersaturation of Dl = 0.275kT and a bond energy / = 2kT. One simulation results in a decaying cluster to a size of nmin = 5, while in two simulations the cluster grows to a size nmax = 300. In 700 subsequent simulations at a supersaturation of Dl = 0.275kT there were 368 simulations resulting in cluster growth to a value of nmax = 300 while 332 simulations resulted in a decay to nmin = 5. Because P(nmin) = 0 and P(nmax) = 1 in this case the growth probability P(n = 144) = 368/700 = 0.526. 5.1. Cluster growth probabilities Using a certain supersaturation and bond energy the growth probabilities for a large number of initial cluster sizes were determined. Fig. 2
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
83
300 250 200 n 150 100 50 0 0
10000
20000
30000
40000
simulation steps Fig. 1. The development of the cluster size of an initial cluster n = 144 at a supersaturation of Dl = 0.275kT in three MC simulations (diamonds, squares and triangles). The developing cluster size n is shown with an interval of 100 simulation steps.
tion of the nucleus size n* and the numerical factor b. At Dl = 0.275kT and a bond strength / = 2kT a nucleus size of n* = 142.3 ± 0.4 and a numerical factor b = (22.9 ± 0.5) · 103 were determined. From Eq. (4) it then follows that the Zeldovich factor z = (12.9 ± 0.3) · 103 and the nucleus region D* = 77.5.
1 0.8 0.8 0.275
0.6
0.175
P (n ) 0.4 0.2
5.2. Nucleus size
0 1
100
200
300
400
500
600
n
Fig. 2. The growth probability P(n) of a cluster at a supersaturation of Dl = 0.175kT (diamonds), 0.275kT (squares) and 0.8kT (triangles) for bond strength / = 2kT. The lines through the points are the best fits of Eq. (3) to the simulation data.
shows the growth probabilities obtained at a supersaturation of Dl = 0.175, 0.275 and 0.8kT for a bond strength of / = 2kT. For a supersaturation Dl = 0.275kT small clusters n < 70 decay and large clusters n > 240 grow. In between there is the nucleus region where a fraction of the clusters decay and another fraction grows out. The growth probability P(n) as a function of the initial cluster size n is described by Eq. (3). As shown by the solid lines in Fig. 2 Eq. (3) describes the behaviour of P(n) as a function of the initial cluster size very well. Therefore the fit of Eq. (3) to the growth probabilities P(n), determined from the simulations, results in an accurate determina-
The described technique was used to determine the nucleus size n* and Zeldovich factor z at a large number of supersaturations for the bond strengths / = 1.33, 2 and 5kT. All bond strength values are chosen to be higher than / = 0.88kT below which the surface is reported to be rough. The simulations are therefore performed below the thermal roughening temperature. Fig. 3 shows the nucleus size n* as a function of the supersaturation. The standard deviation in the nucleus size n* obtained from the growth probabilities is less than 1% at lower supersaturations and increases going to smaller nucleus size. In all events the standard deviation was smaller than the size of the symbols in Fig. 3. In Fig. 3 the solid lines are the predictions with CNT (Eq. (2)) using j1 and c1 from Table 1 for the three broken bond values. For low supersaturations the Gibbs–Thomson equation predicts the nucleus size very well. In all cases however
84
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88 1000
1 5 0.9
(n*CNT/n*)
1/2
100 5
n* 10
2
0.8
2
0.7 1.33
1.33
0.6
1 0
1
2
3
4
5
6
∆ µ /kT
Fig. 3. The nucleus size n* as a function of the supersaturation Dl/kT for the bond strengths / = 1.33 (diamonds), 2 (squares) and 5kT (triangles). The solid lines are drawn with aid of the Gibbs–Thomson equation (Eq. (2)) and the values in Table 1.
the Gibbs–Thomson equation predicts a lower value than obtained from the simulations. To investigate the difference in predicted nucleus size nCNT and simulated nucleus size n* the value of (nCNT =n Þ1=2 is defined. From Eq. (2) it follows that in the case that the shape factor c and the apparent specific step free energy j of the nucleus are independent of the nucleus size n*: 1=2 nCNT c1 j1 ð15Þ ¼ n cj The value of ðnCNT =n Þ1=2 would be unity over the whole range of supersaturations when c and j were independent of the nucleus size n* and the nucleus size prediction is correct. It would give a constant value other than unity when the assumption of the n*-independent c and j were valid while c1 and/or j1 were inaccurately predicted. The value of 1=2 ðnCNT =n Þ would not be constant with the supersaturation if the assumption of the n*-independent c and j does not hold. Since the simulations are performed on discrete clusters the assumption of the n*-independent c and j can be expected not to hold at higher supersaturations where the nucleus size is small. 1=2 Fig. 4 shows the value of ðnCNT =n Þ as a function of the supersaturation. It can immediately be concluded that the apparent specific step free energy and/or the shape factor change with supersat1=2 uration because of the change of ðnCNT =n Þ with supersaturation. The change is largest for / =
0.5 0
1
2
3
4
5
6
∆ µ /κΤ
Fig. 4. The value of ðnCNT =n Þ1=2 as a function of the supersaturation Dl/kT. The slope indicates that the apparent specific step free energy j of the nucleus and/or the shape factor c are nucleus size dependent. The effect is least pronounced at a bond energy / = 5kT. Extrapolation to zero supersaturation shows a lower value than unity with the lowest value for / = 1.33kT. This indicates that the theoretical apparent specific step free energy and/or shape factor of the nucleus may be slightly higher than the actual apparent specific step free energy and/or shape factor.
1.33kT while at / = 5kT the change is just visible in Fig. 4. Because of the gradual change in the 1=2 it seems that the apparent value of ðnCNT =n Þ step free energy and/or the shape factor also gradually change with supersaturation, even at nucleus sizes as large as n* = 500. The value of ðnCNT =n Þ1=2 is lower than 1 because the Gibbs–Thomson equation in all cases predicts a lower value for the nucleus size than that obtained from the simulations: the nucleus size is underestimated by the Gibbs–Thomson equa1=2 tion. The value of ðnCNT =n Þ increases towards a value near 1 when the supersaturation is decreased. This means that the predicted value nCNT approaches the nucleus size n*: the prediction is better at lower supersaturations because there the nucleus size is larger and the values of the shape factor and apparent specific step free energy approach the theoretical values. It seems that an extrapolation towards zero supersaturation in Fig. 4 does not result in a value of unity for ðnCNT =n Þ1=2 , especially in the case of a bond strength of / = 1.33kT. This indicates that the apparent specific step free energy j1 and/or the shape factor c1 are slightly underestimated for all bond strengths. This may be due to the fact
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
that in the simulations vacancies are allowed to appear in the clusters (see Section 4). 5.3. Zeldovich factor The Zeldovich factor obtained from the growth probability simulations can be compared to the Zeldovich factor predicted with aid of the classical nucleation theory. In case of 2D heterogeneous nucleation the Zeldovich factor zCNT from classical nucleation theory is given by [1]: 1=2 Dl=kT zCNT ¼ ð16Þ 4pnCNT Fig. 5 shows the Zeldovich factor z determined from the numerical factor b obtained from the simulations. The solid lines are the predicted Zeldovich factors zCNT given by Eq. (16) for all three broken bond energies. Apart from the higher supersaturations the Zeldovich factor zCNT lies within the significance of z. At higher supersaturations at the lower bond strengths first there seems to be a sudden rise after which the Zeldovich factor z levels off. Since the numerical factor b can only have a value between 0 and 1 this levelling off is expected to occur below a value of z 0.56. At these high supersaturations the nucleus region D* 2 is very narrow and therefore the Zeldovich factors above about z = 0.25 have to be taken with care. In these cases the nucleus region could not be sampled accurately to obtain a
0.6 0.5 2 0.4 z 0.3 5 0.2
1.33
0.1 0 0
1
2
3
4
5
6
∆ µ /kT
Fig. 5. The Zeldovich factor z determined from the numerical factor b obtained from the simulations for / = 1.33 (diamonds), 2 (squares) and 5kT (triangles). The lines give the predicted Zeldovich factors from Eq. (17).
85
sufficiently accurate Zeldovich factor since only 1 or 2 initial cluster sizes fall in this range. While the standard deviation of the nucleus size n* determined from the simulation data is less than 1%, the standard deviation of the Zeldovich factor is larger and increasing with supersaturation. This will be discussed in the next section.
6. Discussion 6.1. A cluster of size n = 1 The growth probability P(n) is the fraction of clusters of size n that grow out to become a crystal monolayer before they decay. Similarly a decay probability Q(n) can be defined which is the fraction of clusters of size n that decay to a size of n = 1 before they are able to grow out. The decay probability Q(n) is given to be [2]: QðnÞ ¼ 1 P ðnÞ ¼
X ðnÞ CðnÞ
ð17Þ
where X(n) is the actual concentration of n-sized clusters. In a supersaturated ambient phase C(n) is an imaginary cluster concentration which would set up if truly stable equilibrium were possible. The C(n) is generally called the equilibrium concentration of n-sized clusters. A further discussion on the equilibrium concentration C and actual concentration X of clusters falls outside the scope of this paper and can be found elsewhere [1]. Eq. (17) can be used for n = 1 and thus to calculate P(n = 1) if X(n = 1) and C(n = 1) are known. An exact expression for the right hand side of Eq. (17) valid for n P 2 is known [1]. A basic assumption in kinetic theory of nucleation is responsible for the value of X(n = 1). This basic assumption is that the actual concentration X(n = 1) of clusters of size 1 is equal to the equilibrium concentration C(n = 1) of clusters of size 1: C(n = 1) = X(n = 1). This leads with aid of Eq. (17) to the assumption that the growth probability P(n = 1) = 0: a cluster of 1 GU has no chance to survive and grow out to become a macroscopically large cluster. In a supersaturated ambient phase the probability to form a macroscopically large cluster is
86
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
nonzero. As this cluster comprises of GU this means that the growth probability P(n = 1) is nonzero: a cluster of 1 GU has a certain (small) chance to survive and grow out to become a macroscopically large cluster. This indicates that the basic assumption in kinetic nucleation theory must therefore be considered only as an approximation. In the kinetic Monte Carlo simulations (and in experiments as well) at sufficiently large supersaturations the growth probability P(n = 1) of a cluster of 1 GU becomes larger than 0. A cluster of 1 GU has only one possible shape configuration. At the periphery of a cluster of 1 GU there are four positions where a GU with one neighbour can be added. Since the cluster contains 1 GU only this GU which has no neighbours can be removed. Assuming that the growth probability P(n = 2) = 1 for a cluster of size n = 2, the growth probability can be determined directly by P ðn ¼ 1Þ pðn ¼ 1Þ ¼ 4k þ =ð4k þ þ k i Þ where p(n) indicates the probability that a cluster of size n grows to a size of n + 1 rather than decay to a size of n 1. For the bond strength / = 2kT and a supersaturation of Dl = 3kT the probability then has the considerable value of P(n = 1) = 0.026 in stead of the assumed growth probability of 0. It is thus important to notice that the basic assumption in kinetic theory of nucleation is only valid at sufficiently low supersaturations and the kinetic theory of nucleation can only be applied at these supersaturations. One implication of this is that the kinetic theory of nucleation should be applied with great care near thermal and kinetic roughening. The basic assumption in kinetic theory furthermore implies that the growth probability method can be used to obtain accurate nucleus sizes and Zeldovich factors when the growth probability P(n = 1) is sufficiently low. With aid of the determined nucleus size n* and numerical factor b the growth probability P(n = 1) can be determined with aid of Eq. (3). For the bond strength / = 1.33kT at a supersaturation of Dl = 0.5kT the nucleus size n* = 12.9 and the numerical factor b = 0.118 lead to a growth probability for n = 1 of P(n = 1) = 0.023. For the bond strength / = 5kT the calculated growth probability for n = 1 is P(n = 1) = 0 down to a nucleus size of n* = 3.26
(Dl = 5.5kT). The assumption thus seems not to hold for / = 1.33kT at supersaturations Dl P 0.5kT while for / = 5kT it seems sufficiently accurate even down to a nucleus size of n* = 3.26. 6.2. The shape of the growth probability curve By fitting the growth probabilities P(n) from the simulations to Eq. (3) the nucleus size n* and Zeldovich factor z were obtained. It was noticed however that the accuracy in the Zeldovich factor did not increase when the number of simulations per growth probability was increased from 200 to 5000 or more. For instance for the bond strength / = 1.33kT the standard deviation is about 13.5% at higher supersaturation (Dl = 1kT, n* = 4.44, M = 1000). Increasing the number of simulations for Dl = 1kT to M = 5000 does not decrease the standard deviation. This indicates a small systematic deviation between the simulated growth probabilities and the fitted growth probabilities. This deviation can be captured in the value DP(n) which is the difference between the actual growth probability P(n) from simulations and the fitted growth probability Px(n) calculated with the values for n* and b from the fit of Eq. (3) to P(n). A positive value for DP(n) means that the actual growth probability P(n) is higher than Px(n). Fig. 6 shows the deviation DP(n) as a function of the growth probability P(n) for the bond strength / = 1.33kT at the supersaturations Dl = 0.25– 1kT. It should be noted that the small systematic deviation does not explain the difference between predicted and actual nucleus size discussed in Section 5. For comparison: at a bond strength of / = 1.33kT and a supersaturation of Dl = 0.25kT the nucleus size was determined to be n* = 43.1 ± 0.4 while the predicted nucleus size nCNT ¼ 31:2. As can be seen in Fig. 6 the deviation is positive around a growth probability of about 0.25 < P(n) < 0.75 and negative at growth probabilities going to 0 or 1. For the bond strength / = 2kT a similar figure was obtained with a maximum value of DP(n) which is lower than the maximum value of DP(n) at / = 1.33kT. The deviation was not larger than the uncertainty in case of the simulations
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88 0.1
0.05
∆ P(n)
0
-0.05
-0.1 0
0.25
0.5
0.75
1
87
abilities to compensate the effect of smaller actual growth probabilities at both the larger and the smaller cluster sizes. The fitting procedure therefore results in a slightly overestimated nucleus size. At a bond strength of / = 1.33kT and a supersaturation of Dl = 0.25kT the nucleus size was determined to be n* = 43.1 ± 0.4 from the fit while the simulated growth probability P(n = 42) = 0.495 seems to indicate that n* is only slightly higher than n = 42.
P(n)
Fig. 6. The difference DP(n) between actual and fitted growth probability for the bond energy / = 1.33kT and the supersaturations in the range of Dl = 0.25–1kT as a function of the actual growth probability P(n).
with a bond strength / = 5kT. It was furthermore noticed that the deviation maximum at about P(n) = 0.5 increases with supersaturation. The deviation is due to the fact that the shape of the growth probability curve from Eq. (3) is slightly different than the shape of the actual growth probability. By evaluation of the different approximations used in the kinetic theory of nucleation it was recognized that the deviation is caused by the approximation in the equilibrium cluster concentration C(n). The equilibrium concentration C(n) of clusters is approximated by [1]: 2 CðnÞ ¼ C exp b2 ðn n Þ ð18Þ The actual equilibrium concentration is larger than the approximated equilibrium concentration at cluster sizes below and smaller than the approximated equilibrium concentration at cluster sizes above the nucleus size n* [1]. This comes down to a smaller actual nucleus region below the nucleus size and a larger one above the nucleus size while the nucleus size is not affected. This means that the growth probability from the simulations is correct around the nucleus size but is overestimated at both smaller and larger cluster sizes. Since the nucleus size and Zeldovich factor are obtained by fitting the obtained growth probabilities over the whole range of cluster sizes to Eq. (3), the fitted growth probabilities around the nucleus size are a bit smaller than the simulated growth prob-
7. Conclusions The growth probability of 2D clusters can readily be determined from computer simulations. From the growth probabilities the nucleus size and the Zeldovich factor of 2D nuclei on the Kossel (1 0 0) surface were determined. At lower supersaturations the determined nucleus sizes approach the nucleus sizes determined with aid of the classical nucleation theory. There are indications that the theoretical values for the apparent specific step free energy and/or shape factor of the nucleus for these lower supersaturations are slightly higher than the ones extrapolated from the simulation data. This is most probably due to the neglect of the creation and annihilation of vacancies in the interior of the 2D clusters. At high supersaturations the nucleus becomes too small to justify the assumptions made to determine the apparent specific step free energy and the shape factor of the nucleus to be used in the classical nucleation theory. In the simulations the step free energy and/or the shape factor of the nucleus are a function of the nucleus size and supersaturation while this is not accounted for in Classical Nucleation Theory. The kinetic theory of nucleation and thus also the growth probability method assumes that the growth probability of a cluster of size 1 is 0. For small bond strengths and extreme supersaturations this is not the case in the simulations. A small systematic deviation between actual and theoretical growth probability is found which originates from the approximation used for the equilibrium cluster concentration. The proposed method will be useful in predicting 3D nucleus
88
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
size, Zeldovich factor and nucleation rate of more complex compounds like polymorphic compounds. Acknowledgment The authors would like to thank Dimo Kashchiev and Hugo Meekes for the stimulating discussions on this subject.
References [1] D. Kashchiev, Nucleation: Basic Theory with Applications, Butterworth-Heinemann, Oxford, 2000. [2] J.H. ter Horst, D. Kashchiev, J. Chem. Phys. 119 (4) (2003) 2241.
[3] D. Kashchiev, J.P. van der Eerden, C. van Leeuwen, J. Crystal Growth 40 (1977) 47. [4] A.B. Bortz, M.H. Kalos, J.L. Lebowitz, J. Comput. Phys. 17 (1975) 10. [5] P. Bennema, J.P. van der Eerden, J. Crystal Growth 42 (1977) 201. [6] D. Frenkel, B. Smit, Understanding Molecular Simulations: from Algorithms to Applications, second ed., Academic Press, San Diego, 2002. [7] H.M. Cuppen, H. Meekes, W.J.P. van Enckevort, P. Bennema, E. Vlieg, Surf. Sci. 525 (2003) 1. [8] C. Rottman, M. Wortis, Phys. Rev. B 24 (11) (1981) 6274. [9] A.A. Chernov, in: H.J. Queisser (Ed.), Springer Series in Solid State Sciences 36: Modern Crystallography III, crystal growth, part 1, Springer-Verlag, Berlin, 1984. [10] H.J. Leamy, G.H. Gilmer, K.A. Jackson, in: J.M. Blakely (Ed.), Surface Physics of Crystalline Materials, Academic Press, New York, 1975. [11] H.M. Cuppen, H. Meekes, W.J.P. van Enckevort, E. Vlieg, H.J.F. Knops, Phys. Rev. B 69 (2004) 245404.
Nucleus size and Zeldovich factor in two-dimensional nucleation at the Kossel crystal (0 0 1) surface J.H. ter Horst *, P.J. Jansens Laboratory for Process Equipment, Delft University of Technology, Leeghwaterstraat 44, 2628 CA Delft, The Netherlands Received 19 May 2004; accepted for publication 14 October 2004 Available online 2 November 2004
Abstract A Monte Carlo simulation method to determine the growth probability of two-dimensional clusters on the Kossel (0 0 1) surface is described. From this growth probability the nucleus size and the Zeldovich factor for two-dimensional nucleation were obtained. For low supersaturations the Gibbs–Thomson equation describes the supersaturation dependency of the nucleus size very well. At higher supersaturations a deviation from the predicted nucleus sizes is observed. This is mainly due to the size dependence of the apparent specific step free energy and/or the shape factor of the nucleus. A small systematic deviation between the theoretical and actual growth probability is observed and explained. 2004 Elsevier B.V. All rights reserved. Keywords: Nucleation; Crystallization; Computer simulations; Monte Carlo simulations
1. Introduction It is generally not possible to predict sufficiently accurate nucleus sizes and nucleation rates for complex crystalline structures using purely theoretical considerations. For these complex crystal structures accurate step free energy (or in the case
* Corresponding author. Tel.: +31 15 278 6678; fax: +31 15 278 6975. E-mail address: [email protected] (J.H. ter Horst).
of 3D nucleation the interfacial energy) predictions and shape factor predictions are needed in order to calculate the nucleus size and nucleation rate correctly. Molecular simulations could aid to improve the accuracy in predictions by simulating the nucleation process. Molecular simulations can furthermore be used to improve the fundamental understanding of the nucleation process. In this paper a newly proposed general method to determine the nucleus size and the Zeldovich factor will be validated by using it on a simple model system. The 2D cluster formation on the Kossel crystal (0 0 1) surface was chosen as the
0039-6028/$ - see front matter 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.susc.2004.10.020
78
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
model system. For the Kossel model system relatively accurate values of the shape factor and the specific step free energy can be determined. As a first objective, this paper is to demonstrate that this method determines the nucleus size and Zeldovich factor. This paper is furthermore meant to find the origins of the difference between molecular simulation results and Classical Nucleation Theory (CNT) predictions thus studying the limitations of the CNT. The formation of a two-dimensional cluster consisting of n growth units (GU) results in (1) a free energy gain due to the supersaturation in the ambient phase and (2) an excess of free energy due to the creation of a (closed) step. The free energy gain is linearly related to the cluster size n and the supersaturation while the excess free energy can be seen as a total step free energy that is a function of the length L of the cluster periphery and an apparent specific step free energy j. The work W(n) to form a n-sized cluster then is written as [1]: W ðnÞ ¼ n Dl þ Lj ¼ n Dl þ cajn1=2
ð1Þ
where c is a shape factor, a is the distance between nearest neighbours, and Dl is the supersaturation, which is the difference in chemical potential between the ambient phase and the crystalline phase. The j is apparent because the actual specific step free energy may not be isotropic. At large cluster sizes the energy gain is dominant while at small sizes the total step free energy is larger. Therefore a cluster size exists for which the work to form it is maximal. This cluster size is the nucleus size n = n* and the corresponding energy barrier is the nucleation work W*. The nucleus is in unstable equilibrium with the ambient phase. Clusters smaller than the nucleus size tend to decay while clusters larger than the nucleus size tend to grow out. In CNT it is assumed that the apparent specific step free energy and the shape factor are independent of the size n and the supersaturation Dl and that they have the values j = j1 and c = c1 of an equilibrium shaped and macroscopically large 2D nucleus. Setting the derivative of the work W(n) to the cluster size n to zero results in the 2D variant of the classical Gibbs–Thomson equation [1] with aid of which the nucleus size nCNT can be predicted:
nCNT ¼
c21 a2 j21 4 Dl2
ð2Þ
The assumptions on the shape factor and the apparent step free energy make Eq. (2) approximate. It is not clear down to what nucleus sizes and up to what supersaturations these approximations are acceptable. One of the purposes of the simulation study described in this paper is to check the validity of Eq. (2) by comparing the nucleus size from simulations and the nucleus size determined from Eq. (2). First a recently proposed equation for cluster growth probabilities [2] is introduced (Section 2). These cluster growth probabilities can be determined with aid of computer simulation techniques. This approach is followed for 2D nucleation on the Kossel crystal (0 0 1) surface by performing kinetic Monte Carlo simulations (Section 3). Theoretical estimates for the shape factor c1 and the apparent specific step free energy j1 to be used in the CNT predictions are given and discussed in Section 4. The cluster growth probabilities from simulations result in the nucleus size and the Zeldovich factor at a certain supersaturation. In Section 5 the simulation results are presented and discussed after which in Section 6 two important aspects of the growth probability method and the kinetic theory of nucleation are discussed.
2. Growth probability of n-sized clusters Nucleation is known to occur by formation of molecular (or atomic) clusters of various sizes. The cluster size changes randomly as a result of successive attachments and detachments of single growth units (GU: molecules, atoms, ions) to and from the cluster. The cluster of size n = n* is the nucleus. It is that particular cluster to and from which molecules are attached and detached with equal frequency. Clusters of size n < n* tend to decay, because per unit time less molecules are attached to than detached from them. On the contrary, the clusters of size n > n* are characterized by attachment frequencies greater than the detachment ones and for that reason these clusters tend to grow up to macroscopic sizes. As mole-
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
cular attachments and detachments are random events, a given n-sized cluster can grow and reach a macroscopic size only with a certain growth probability P(n). Using the kinetic theory of nucleation it turns out that the growth probability P(n) at the expense of certain approximations but without loss of generality can be expressed in terms of a numerical factor b and the nucleus size n* [2]: P ðnÞ ¼ 12½1 þ erfðbðn n ÞÞ
ð3Þ
where erf is the error function. The growth probability expression is valid for any kind of one-component nucleation (two- and three-dimensional, homogeneous and heterogeneous nucleation of vapours, liquids or solids) [2]. The numerical factor b is related to the width D* of the nucleus region and to the Zeldovich factor z [1]: b¼
p1=2 ¼ p1=2 z D
ð4Þ
The Zeldovich factor is an important parameter in the pre-exponential kinetic factor of the nucleation rate equation. Outside the nucleus region for n < n* D*/2 the growth probability P(n) approaches 0 while for n > n* + D*/2 the value of P(n) approaches 1. This paper will show that the growth probability expression is valuable to determine nucleation characteristics with aid of computer simulations. The growth probability as a function of initial cluster size at a certain supersaturation can be determined with aid of computer simulations. From the simulated growth probability P(n) as a function of the cluster size at a certain supersaturation the nucleus size n* and the numerical factor b can be determined through Eq. (3). This has been done for a liquid droplet in a supersaturated vapour Lennard–Jones compound [2]. In this paper it is focussed on the 2D cluster formation on the Kossel (0 0 1) surface as a first step towards 3D nucleation simulations of e.g., polymorphic crystalline compounds. Such 3D nucleation simulations lead to the experimentally advantageous prediction of the supersaturations under which concomitant nucleation can be avoided.
79
3. Simulation details It was chosen to study the growth probability of 2D clusters on the Kossel (0 0 1) surface. A Kossel crystal has an idealized simple cubic crystal structure with a lattice constant a. Both crystal and ambient phase, either vapour, melt or solution, are at temperature T. The ambient phase has the same lattice structure as the crystalline compound. The bond energy ess refers to the interaction between growth units (GU) in the crystalline compound while the energy ell refers to the interaction between ambient phase units. The bond energy esl is the interaction energy between a GU and an ambient phase unit. Now the overall lateral bond strength / can be defined [3]: / ess þ ell 2esl ¼ kT 2kT
ð5Þ
This is the energy loss per growth unit upon the breaking of a solid–solid bond which coincides with the breaking of an ambient phase-ambient phase bond and the creation of 2 solid-ambient phase bonds. The growth probability of a mono-atomically thick 2D cluster consisting of n GU on such a Kossel (0 0 1) surface can be studied with the n-fold way Monte Carlo (MC) algorithm [4]. Values of / = 1.33, 2 and 5kT were chosen in the simulations. An increase in bond strength / is equivalent to a decrease in temperature. At T = 0 K or an infinitely large bond strength the (0 0 1) surface is atomically smooth. Lowering the bond strength or increasing the temperature will lead to an increase of adatoms, vacancies and kinks in the steps. At a relatively low bond strength of / = 1.33kT (corresponding to a relatively high temperature) a macroscopically large nucleus on the (0 0 1) surface is circular while at a relatively high bond strength of / = 5kT (relatively low temperature) the macroscopically large nucleus has a shape between a square and a circle (see Section 4). Each MC step a GU is either created or annihilated while surface diffusion is neglected. There are eight types of transitions: creation and annihilation of a solid GU with either 1, 2, 3 or 4 lateral solid neighbours (the creation of GU without
80
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
lateral neighbours was not considered in the simulations: only the cluster grows or decays in the simulation). Then k þ i can be defined as the creation rate of GU with i neighbours and k i as the annihilation rate of GU with i neighbours. To ensure microscopic reversibility [5], also known as the detailed balance condition in molecular simulation theory [6], the quotient of the rates of these transitions must be equal to the quotient of the corresponding Boltzmann factors: kþ U ðiþ Þ U ði Þ i ¼ exp kT k i / Dl ¼ exp 2ði 2Þ þ ð6Þ kT kT Here U(i+) is the energy of the system after a GU with i neighbours is created, U(i) is the energy of the system after a GU with i neighbours is removed, U(i+) U(i) is the energy difference between a filled and empty solid position with i neighbours and Dl is the supersaturation. Microscopic reversibility leaves either k þ i or k i free to choose in a simulation. The creation rate þ kþ can be chosen to be independent of the i ¼ k number of neighbours i: Dl k þ ¼ k eq exp ð7Þ kT The rate keq is the single GU addition rate per site at equilibrium between the crystal solid phase and the ambient phase. Physically this corresponds to a random impingement of GU on the cluster. From Eqs. (6) and (7) it then follows that the annihilation rate k i for a GU is independent of the supersaturation and given by: / k i ¼ k eq exp 2ði 2Þ ð8Þ kT The defined rates ensure that at a supersaturation of Dl = 0 the Kossel (0 0 1) surface is in equilibrium and thus does not grow or dissolve. The transition rates k+ and k i are proportional to the transition probabilities p+ and p i to be used in the MC simulation. This MC scheme is also known as the random rain model [7]. Only the ambient phase changes as compared to equilibrium conditions when a supersaturation is applied,
typically the case for solution and vapour growth [7]. Due to the definition of the system of a developing two-dimensional cluster there are a limited number of sites N where a GU can be either created or annihilated. Each of these N positions corresponds to a translation probability (p+ for a crea tion position and p or p for an 1 , p2 , p3 4 annihilation position). Each MC step a GU is either created or annihilated which means that the summation over all N translation probabilities is unity. Through a random number in the range [0, 1) the translation site is picked from these N sites. The growth probability P(n) of an n-sized cluster is determined by performing a number M of MC simulations (usually M = 300–10 000) starting with the same initial cluster size n. A minimum nmin P 1 and a maximum nmax cluster size were chosen for which the values P(nmin) = 0 and P(nmax) = 1 were sufficiently accurate approximations. The simulation was stopped only when the cluster size reached either n = nmin or n = nmax. The number M+ of simulations resulting in a cluster of size n = nmax gives the growth probability: P ðnÞ ¼
Mþ M
ð9Þ
Rectangular initial cluster shapes with minimum step energy were used. The initial clusters used in the MC simulations all had a rectangular shape and a size of 2 · 1, 2 · 2, 3 · 2, 3 · 3, 4 · 3, 4 · 4, 5 · 4, . . ., 25 · 24, 25 · 25. In addition, at / = 2kT and Dl = 0.3kT apart from the 10 · 10 sized and thus square shaped initial cluster, elongated initial shapes of 4 · 25 and 5 · 20 sized initial clusters were used. The elongated initial shapes did not have a significant effect on the calculated growth probability of that cluster size n = 100. This indicates that the shape equilibration process is generally much faster than the growth and decay processes, validating the choice of as square as possible rectangular initial clusters. Furthermore in case of a small nucleus size the shape factor is mainly determined by the lowest energy cluster configuration (square) because of the relatively large energy difference between the low and higher cluster configurations for small clusters.
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
4. The nucleus size using the classical nucleation theory The nucleus size is that cluster size for which the attachment frequency and the detachment frequency of growth units are equal. The nucleus is therefore in (unstable) equilibrium with the supersaturated ambient phase and will thus have the equilibrium shape. In the case of a macroscopically large nucleus on the Kossel crystal (0 0 1) surface the exact form of the specific step free energy j(//kT, h) is known [8]: ajð/=kT ,hÞ ¼ j cos hjsinh1 ðaj cos hjÞ kT ð10Þ þ j sin hjsinh1 ðaj sin hjÞ where " #1=2 2 1 b2 a¼ ð11Þ b 1 þ ðsin2 2h þ b2 cos2 2hÞ1=2 and
/ 2 sinh kT b¼ ð12Þ / 2 cosh kT Eq. (10) reflects that for anisotropic lattices the actual specific step free energy j = j(//kT, h) changes with the bond strength //kT and the orientation angle h from the [1 0 0] direction. For a specific bond strength a polar diagram of the actual specific step free energy j(//kT, h) can be made. The equilibrium shape is that shape for which the perpendicular distance of all steps to the nucleus centre is proportional to the actual specific step free energy j(//kT, h) of that step (Gibbs–Curie–Wulff rule) [9]. Since the actual specific step free energy j(//kT, h) is known from Eq. (10) the equilibrium shape can be constructed which was done numerically. From the equilibrium shape construction then the surface area (which is a relative nucleus size n1 ) and the periphery length (which is a relative nucleus periphery length L1 ) are known. Then a good approximation of the shape factor c1 of the nucleus to be used in the Gibbs–Thomson equation is given by: L ð13Þ c1 ¼ p1 ffiffiffiffiffiffi a n1
81
Here for L1 the periphery length of the calculated equilibrium shape and for n1 the surface area of the calculated equilibrium shape is taken. The shape factor is calculated to be c1 = 3.545, 3.554 and 3.681 for respectively / = 1.33, 2 and 5kT. At / = 1.33kT the nucleus is near perfectly spherical (c1 2p1/2). The shape factor at a bond strength of / = 5kT indicates a rather spherical nucleus (c1 = 3.681 compared to c = 4 for a square shape at an infinitely large bond strength or a temperatures of 0 K). In classical nucleation theory the total step free energy of a nucleus is approximated by the value 1=2 c1 aj1 n1 . The total step free energy of the equilibrium shape can be calculated because at each infinitely small length dL of the periphery of the equilibrium shape the specific step free energy j(//kT, L) is known from the polar diagram and the equilibrium shape construction. The total step free energy of a nucleus in classical nucleation theory and of a macroscopically large nucleus at equilibrium are equal if for n1 the surface area of the calculated equilibrium shape is taken: Z L 1=2 c1 aj1 n1 ¼ jð/=kT ,LÞ dL ð14Þ 0
In Eq. (14) only the apparent specific step free energy j1 of the nucleus is unknown and thus can be determined, which was done numerically. In the case of bond strength / = 1.33, 2 and 5kT this leads to an apparent specific step free energy of respectively j1 = 0.793, 1.768 and 5.244kT/a. These values can be compared to the specific step free energy of an infinitely long step [10] j = 0.788, 1.728 and 4.993kT/a for a bond strength of respectively / = 1.33, 2 and 5kT. The difference between the apparent specific step free energy of the nucleus and the specific step free energy of an infinitely long step reflects the effect of the corners of the nucleus. A perfectly spherical or perfectly square shaped nucleus has an apparent specific step free energy equal to the specific step free energy of an infinitely long step. A shape between these two extremes results in an apparent specific step free energy which is slightly larger than the specific step free energy of an infinitely large step, e.g., j1 = 5.244kT/a compared to j = 4.993kT/a
82
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
Table 1 The shape factor c1, and the specific step free energy j1 as a function of the bond strength //kT calculated with the method described in Section 4 //kT
c1
aj1/kT
1.33 2 5
3.545 3.554 3.681
0.793 1.768 5.244
for / = 5kT due to the higher specific step free energy j(//kT, L) at the corners. In Table 1 the obtained values for the shape factor c1 and step free energy j1 are summarized. These values can be used to determine the nucleus size nCNT as a function of the supersaturation with aid of the Gibbs–Thomson equation. It should be noted however that there are some differences between the conditions of the simulations and the theory used to calculate c1 and j1. First of all, the values in Table 1 relate to macroscopically large nuclei while in the simulations the clusters have a finite size. It is expected that the finite size of the nucleus increases the shape factor c compared to c1 because less nucleus configurations are possible and the relative difference between the lowest and the higher energy configurations of the nucleus is larger. The nucleus shape therefore becomes more square-like at smaller sizes. The same argument leads to a lower entropy of the system from which it follows that the apparent specific step free energy j1 increases. A cluster would have a larger total step free energy in the simulations than is predicted and therefore the nucleus size n* will be larger than predicted with the classical nucleation theory. It is expected that this effect becomes larger at smaller nucleus sizes. The second difference between simulations and theory is the prevailing supersaturation in the simulations. The supersaturation decreases the specific step free energy [11] and therefore possibly also the apparent specific step free energy and the shape factor of the nucleus. The third difference between simulation and theory is that within the simulations there is the possibility that vacancies in the interior of the nucleus appear. This is an entropy effect causing a roughening of the surface. This would mean that the total step free energy is lower. During the
simulations such vacancies were rarely observed. It can be said that this may cause a difference for the value from simulations with the theoretical value in the extrapolation towards macroscopically large nuclei (zero supersaturation). Even at zero supersaturation on a flat surface vacancies in that surface will be created and annihilated continuously. This effect would be higher for lower bond strengths. It is therefore expected that for large nuclei and small supersaturations (where the effect of the first two differences is not present) still a small difference between the predicted nucleus size and the nucleus size from simulations may be present. This difference then is larger for smaller bond strengths.
5. The nucleus size from the cluster growth probabilities A cluster much larger than the nucleus size always survives and grows out to become a new growth layer while a cluster much smaller than the nucleus size has an infinitely small growth probability. A cluster of a size within the nucleus region has a certain probability 0 < P(n) < 1: While some clusters of a certain size grow out to macroscopic sizes, others will decay. This is observed in subsequent MC simulations of clusters with an initial size n. Fig. 1 shows the development of the cluster size during three simulations for an initial cluster size of n = 144 at a supersaturation of Dl = 0.275kT and a bond energy / = 2kT. One simulation results in a decaying cluster to a size of nmin = 5, while in two simulations the cluster grows to a size nmax = 300. In 700 subsequent simulations at a supersaturation of Dl = 0.275kT there were 368 simulations resulting in cluster growth to a value of nmax = 300 while 332 simulations resulted in a decay to nmin = 5. Because P(nmin) = 0 and P(nmax) = 1 in this case the growth probability P(n = 144) = 368/700 = 0.526. 5.1. Cluster growth probabilities Using a certain supersaturation and bond energy the growth probabilities for a large number of initial cluster sizes were determined. Fig. 2
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
83
300 250 200 n 150 100 50 0 0
10000
20000
30000
40000
simulation steps Fig. 1. The development of the cluster size of an initial cluster n = 144 at a supersaturation of Dl = 0.275kT in three MC simulations (diamonds, squares and triangles). The developing cluster size n is shown with an interval of 100 simulation steps.
tion of the nucleus size n* and the numerical factor b. At Dl = 0.275kT and a bond strength / = 2kT a nucleus size of n* = 142.3 ± 0.4 and a numerical factor b = (22.9 ± 0.5) · 103 were determined. From Eq. (4) it then follows that the Zeldovich factor z = (12.9 ± 0.3) · 103 and the nucleus region D* = 77.5.
1 0.8 0.8 0.275
0.6
0.175
P (n ) 0.4 0.2
5.2. Nucleus size
0 1
100
200
300
400
500
600
n
Fig. 2. The growth probability P(n) of a cluster at a supersaturation of Dl = 0.175kT (diamonds), 0.275kT (squares) and 0.8kT (triangles) for bond strength / = 2kT. The lines through the points are the best fits of Eq. (3) to the simulation data.
shows the growth probabilities obtained at a supersaturation of Dl = 0.175, 0.275 and 0.8kT for a bond strength of / = 2kT. For a supersaturation Dl = 0.275kT small clusters n < 70 decay and large clusters n > 240 grow. In between there is the nucleus region where a fraction of the clusters decay and another fraction grows out. The growth probability P(n) as a function of the initial cluster size n is described by Eq. (3). As shown by the solid lines in Fig. 2 Eq. (3) describes the behaviour of P(n) as a function of the initial cluster size very well. Therefore the fit of Eq. (3) to the growth probabilities P(n), determined from the simulations, results in an accurate determina-
The described technique was used to determine the nucleus size n* and Zeldovich factor z at a large number of supersaturations for the bond strengths / = 1.33, 2 and 5kT. All bond strength values are chosen to be higher than / = 0.88kT below which the surface is reported to be rough. The simulations are therefore performed below the thermal roughening temperature. Fig. 3 shows the nucleus size n* as a function of the supersaturation. The standard deviation in the nucleus size n* obtained from the growth probabilities is less than 1% at lower supersaturations and increases going to smaller nucleus size. In all events the standard deviation was smaller than the size of the symbols in Fig. 3. In Fig. 3 the solid lines are the predictions with CNT (Eq. (2)) using j1 and c1 from Table 1 for the three broken bond values. For low supersaturations the Gibbs–Thomson equation predicts the nucleus size very well. In all cases however
84
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88 1000
1 5 0.9
(n*CNT/n*)
1/2
100 5
n* 10
2
0.8
2
0.7 1.33
1.33
0.6
1 0
1
2
3
4
5
6
∆ µ /kT
Fig. 3. The nucleus size n* as a function of the supersaturation Dl/kT for the bond strengths / = 1.33 (diamonds), 2 (squares) and 5kT (triangles). The solid lines are drawn with aid of the Gibbs–Thomson equation (Eq. (2)) and the values in Table 1.
the Gibbs–Thomson equation predicts a lower value than obtained from the simulations. To investigate the difference in predicted nucleus size nCNT and simulated nucleus size n* the value of (nCNT =n Þ1=2 is defined. From Eq. (2) it follows that in the case that the shape factor c and the apparent specific step free energy j of the nucleus are independent of the nucleus size n*: 1=2 nCNT c1 j1 ð15Þ ¼ n cj The value of ðnCNT =n Þ1=2 would be unity over the whole range of supersaturations when c and j were independent of the nucleus size n* and the nucleus size prediction is correct. It would give a constant value other than unity when the assumption of the n*-independent c and j were valid while c1 and/or j1 were inaccurately predicted. The value of 1=2 ðnCNT =n Þ would not be constant with the supersaturation if the assumption of the n*-independent c and j does not hold. Since the simulations are performed on discrete clusters the assumption of the n*-independent c and j can be expected not to hold at higher supersaturations where the nucleus size is small. 1=2 Fig. 4 shows the value of ðnCNT =n Þ as a function of the supersaturation. It can immediately be concluded that the apparent specific step free energy and/or the shape factor change with supersat1=2 uration because of the change of ðnCNT =n Þ with supersaturation. The change is largest for / =
0.5 0
1
2
3
4
5
6
∆ µ /κΤ
Fig. 4. The value of ðnCNT =n Þ1=2 as a function of the supersaturation Dl/kT. The slope indicates that the apparent specific step free energy j of the nucleus and/or the shape factor c are nucleus size dependent. The effect is least pronounced at a bond energy / = 5kT. Extrapolation to zero supersaturation shows a lower value than unity with the lowest value for / = 1.33kT. This indicates that the theoretical apparent specific step free energy and/or shape factor of the nucleus may be slightly higher than the actual apparent specific step free energy and/or shape factor.
1.33kT while at / = 5kT the change is just visible in Fig. 4. Because of the gradual change in the 1=2 it seems that the apparent value of ðnCNT =n Þ step free energy and/or the shape factor also gradually change with supersaturation, even at nucleus sizes as large as n* = 500. The value of ðnCNT =n Þ1=2 is lower than 1 because the Gibbs–Thomson equation in all cases predicts a lower value for the nucleus size than that obtained from the simulations: the nucleus size is underestimated by the Gibbs–Thomson equa1=2 tion. The value of ðnCNT =n Þ increases towards a value near 1 when the supersaturation is decreased. This means that the predicted value nCNT approaches the nucleus size n*: the prediction is better at lower supersaturations because there the nucleus size is larger and the values of the shape factor and apparent specific step free energy approach the theoretical values. It seems that an extrapolation towards zero supersaturation in Fig. 4 does not result in a value of unity for ðnCNT =n Þ1=2 , especially in the case of a bond strength of / = 1.33kT. This indicates that the apparent specific step free energy j1 and/or the shape factor c1 are slightly underestimated for all bond strengths. This may be due to the fact
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
that in the simulations vacancies are allowed to appear in the clusters (see Section 4). 5.3. Zeldovich factor The Zeldovich factor obtained from the growth probability simulations can be compared to the Zeldovich factor predicted with aid of the classical nucleation theory. In case of 2D heterogeneous nucleation the Zeldovich factor zCNT from classical nucleation theory is given by [1]: 1=2 Dl=kT zCNT ¼ ð16Þ 4pnCNT Fig. 5 shows the Zeldovich factor z determined from the numerical factor b obtained from the simulations. The solid lines are the predicted Zeldovich factors zCNT given by Eq. (16) for all three broken bond energies. Apart from the higher supersaturations the Zeldovich factor zCNT lies within the significance of z. At higher supersaturations at the lower bond strengths first there seems to be a sudden rise after which the Zeldovich factor z levels off. Since the numerical factor b can only have a value between 0 and 1 this levelling off is expected to occur below a value of z 0.56. At these high supersaturations the nucleus region D* 2 is very narrow and therefore the Zeldovich factors above about z = 0.25 have to be taken with care. In these cases the nucleus region could not be sampled accurately to obtain a
0.6 0.5 2 0.4 z 0.3 5 0.2
1.33
0.1 0 0
1
2
3
4
5
6
∆ µ /kT
Fig. 5. The Zeldovich factor z determined from the numerical factor b obtained from the simulations for / = 1.33 (diamonds), 2 (squares) and 5kT (triangles). The lines give the predicted Zeldovich factors from Eq. (17).
85
sufficiently accurate Zeldovich factor since only 1 or 2 initial cluster sizes fall in this range. While the standard deviation of the nucleus size n* determined from the simulation data is less than 1%, the standard deviation of the Zeldovich factor is larger and increasing with supersaturation. This will be discussed in the next section.
6. Discussion 6.1. A cluster of size n = 1 The growth probability P(n) is the fraction of clusters of size n that grow out to become a crystal monolayer before they decay. Similarly a decay probability Q(n) can be defined which is the fraction of clusters of size n that decay to a size of n = 1 before they are able to grow out. The decay probability Q(n) is given to be [2]: QðnÞ ¼ 1 P ðnÞ ¼
X ðnÞ CðnÞ
ð17Þ
where X(n) is the actual concentration of n-sized clusters. In a supersaturated ambient phase C(n) is an imaginary cluster concentration which would set up if truly stable equilibrium were possible. The C(n) is generally called the equilibrium concentration of n-sized clusters. A further discussion on the equilibrium concentration C and actual concentration X of clusters falls outside the scope of this paper and can be found elsewhere [1]. Eq. (17) can be used for n = 1 and thus to calculate P(n = 1) if X(n = 1) and C(n = 1) are known. An exact expression for the right hand side of Eq. (17) valid for n P 2 is known [1]. A basic assumption in kinetic theory of nucleation is responsible for the value of X(n = 1). This basic assumption is that the actual concentration X(n = 1) of clusters of size 1 is equal to the equilibrium concentration C(n = 1) of clusters of size 1: C(n = 1) = X(n = 1). This leads with aid of Eq. (17) to the assumption that the growth probability P(n = 1) = 0: a cluster of 1 GU has no chance to survive and grow out to become a macroscopically large cluster. In a supersaturated ambient phase the probability to form a macroscopically large cluster is
86
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
nonzero. As this cluster comprises of GU this means that the growth probability P(n = 1) is nonzero: a cluster of 1 GU has a certain (small) chance to survive and grow out to become a macroscopically large cluster. This indicates that the basic assumption in kinetic nucleation theory must therefore be considered only as an approximation. In the kinetic Monte Carlo simulations (and in experiments as well) at sufficiently large supersaturations the growth probability P(n = 1) of a cluster of 1 GU becomes larger than 0. A cluster of 1 GU has only one possible shape configuration. At the periphery of a cluster of 1 GU there are four positions where a GU with one neighbour can be added. Since the cluster contains 1 GU only this GU which has no neighbours can be removed. Assuming that the growth probability P(n = 2) = 1 for a cluster of size n = 2, the growth probability can be determined directly by P ðn ¼ 1Þ pðn ¼ 1Þ ¼ 4k þ =ð4k þ þ k i Þ where p(n) indicates the probability that a cluster of size n grows to a size of n + 1 rather than decay to a size of n 1. For the bond strength / = 2kT and a supersaturation of Dl = 3kT the probability then has the considerable value of P(n = 1) = 0.026 in stead of the assumed growth probability of 0. It is thus important to notice that the basic assumption in kinetic theory of nucleation is only valid at sufficiently low supersaturations and the kinetic theory of nucleation can only be applied at these supersaturations. One implication of this is that the kinetic theory of nucleation should be applied with great care near thermal and kinetic roughening. The basic assumption in kinetic theory furthermore implies that the growth probability method can be used to obtain accurate nucleus sizes and Zeldovich factors when the growth probability P(n = 1) is sufficiently low. With aid of the determined nucleus size n* and numerical factor b the growth probability P(n = 1) can be determined with aid of Eq. (3). For the bond strength / = 1.33kT at a supersaturation of Dl = 0.5kT the nucleus size n* = 12.9 and the numerical factor b = 0.118 lead to a growth probability for n = 1 of P(n = 1) = 0.023. For the bond strength / = 5kT the calculated growth probability for n = 1 is P(n = 1) = 0 down to a nucleus size of n* = 3.26
(Dl = 5.5kT). The assumption thus seems not to hold for / = 1.33kT at supersaturations Dl P 0.5kT while for / = 5kT it seems sufficiently accurate even down to a nucleus size of n* = 3.26. 6.2. The shape of the growth probability curve By fitting the growth probabilities P(n) from the simulations to Eq. (3) the nucleus size n* and Zeldovich factor z were obtained. It was noticed however that the accuracy in the Zeldovich factor did not increase when the number of simulations per growth probability was increased from 200 to 5000 or more. For instance for the bond strength / = 1.33kT the standard deviation is about 13.5% at higher supersaturation (Dl = 1kT, n* = 4.44, M = 1000). Increasing the number of simulations for Dl = 1kT to M = 5000 does not decrease the standard deviation. This indicates a small systematic deviation between the simulated growth probabilities and the fitted growth probabilities. This deviation can be captured in the value DP(n) which is the difference between the actual growth probability P(n) from simulations and the fitted growth probability Px(n) calculated with the values for n* and b from the fit of Eq. (3) to P(n). A positive value for DP(n) means that the actual growth probability P(n) is higher than Px(n). Fig. 6 shows the deviation DP(n) as a function of the growth probability P(n) for the bond strength / = 1.33kT at the supersaturations Dl = 0.25– 1kT. It should be noted that the small systematic deviation does not explain the difference between predicted and actual nucleus size discussed in Section 5. For comparison: at a bond strength of / = 1.33kT and a supersaturation of Dl = 0.25kT the nucleus size was determined to be n* = 43.1 ± 0.4 while the predicted nucleus size nCNT ¼ 31:2. As can be seen in Fig. 6 the deviation is positive around a growth probability of about 0.25 < P(n) < 0.75 and negative at growth probabilities going to 0 or 1. For the bond strength / = 2kT a similar figure was obtained with a maximum value of DP(n) which is lower than the maximum value of DP(n) at / = 1.33kT. The deviation was not larger than the uncertainty in case of the simulations
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88 0.1
0.05
∆ P(n)
0
-0.05
-0.1 0
0.25
0.5
0.75
1
87
abilities to compensate the effect of smaller actual growth probabilities at both the larger and the smaller cluster sizes. The fitting procedure therefore results in a slightly overestimated nucleus size. At a bond strength of / = 1.33kT and a supersaturation of Dl = 0.25kT the nucleus size was determined to be n* = 43.1 ± 0.4 from the fit while the simulated growth probability P(n = 42) = 0.495 seems to indicate that n* is only slightly higher than n = 42.
P(n)
Fig. 6. The difference DP(n) between actual and fitted growth probability for the bond energy / = 1.33kT and the supersaturations in the range of Dl = 0.25–1kT as a function of the actual growth probability P(n).
with a bond strength / = 5kT. It was furthermore noticed that the deviation maximum at about P(n) = 0.5 increases with supersaturation. The deviation is due to the fact that the shape of the growth probability curve from Eq. (3) is slightly different than the shape of the actual growth probability. By evaluation of the different approximations used in the kinetic theory of nucleation it was recognized that the deviation is caused by the approximation in the equilibrium cluster concentration C(n). The equilibrium concentration C(n) of clusters is approximated by [1]: 2 CðnÞ ¼ C exp b2 ðn n Þ ð18Þ The actual equilibrium concentration is larger than the approximated equilibrium concentration at cluster sizes below and smaller than the approximated equilibrium concentration at cluster sizes above the nucleus size n* [1]. This comes down to a smaller actual nucleus region below the nucleus size and a larger one above the nucleus size while the nucleus size is not affected. This means that the growth probability from the simulations is correct around the nucleus size but is overestimated at both smaller and larger cluster sizes. Since the nucleus size and Zeldovich factor are obtained by fitting the obtained growth probabilities over the whole range of cluster sizes to Eq. (3), the fitted growth probabilities around the nucleus size are a bit smaller than the simulated growth prob-
7. Conclusions The growth probability of 2D clusters can readily be determined from computer simulations. From the growth probabilities the nucleus size and the Zeldovich factor of 2D nuclei on the Kossel (1 0 0) surface were determined. At lower supersaturations the determined nucleus sizes approach the nucleus sizes determined with aid of the classical nucleation theory. There are indications that the theoretical values for the apparent specific step free energy and/or shape factor of the nucleus for these lower supersaturations are slightly higher than the ones extrapolated from the simulation data. This is most probably due to the neglect of the creation and annihilation of vacancies in the interior of the 2D clusters. At high supersaturations the nucleus becomes too small to justify the assumptions made to determine the apparent specific step free energy and the shape factor of the nucleus to be used in the classical nucleation theory. In the simulations the step free energy and/or the shape factor of the nucleus are a function of the nucleus size and supersaturation while this is not accounted for in Classical Nucleation Theory. The kinetic theory of nucleation and thus also the growth probability method assumes that the growth probability of a cluster of size 1 is 0. For small bond strengths and extreme supersaturations this is not the case in the simulations. A small systematic deviation between actual and theoretical growth probability is found which originates from the approximation used for the equilibrium cluster concentration. The proposed method will be useful in predicting 3D nucleus
88
J.H. ter Horst, P.J. Jansens / Surface Science 574 (2005) 77–88
size, Zeldovich factor and nucleation rate of more complex compounds like polymorphic compounds. Acknowledgment The authors would like to thank Dimo Kashchiev and Hugo Meekes for the stimulating discussions on this subject.
References [1] D. Kashchiev, Nucleation: Basic Theory with Applications, Butterworth-Heinemann, Oxford, 2000. [2] J.H. ter Horst, D. Kashchiev, J. Chem. Phys. 119 (4) (2003) 2241.
[3] D. Kashchiev, J.P. van der Eerden, C. van Leeuwen, J. Crystal Growth 40 (1977) 47. [4] A.B. Bortz, M.H. Kalos, J.L. Lebowitz, J. Comput. Phys. 17 (1975) 10. [5] P. Bennema, J.P. van der Eerden, J. Crystal Growth 42 (1977) 201. [6] D. Frenkel, B. Smit, Understanding Molecular Simulations: from Algorithms to Applications, second ed., Academic Press, San Diego, 2002. [7] H.M. Cuppen, H. Meekes, W.J.P. van Enckevort, P. Bennema, E. Vlieg, Surf. Sci. 525 (2003) 1. [8] C. Rottman, M. Wortis, Phys. Rev. B 24 (11) (1981) 6274. [9] A.A. Chernov, in: H.J. Queisser (Ed.), Springer Series in Solid State Sciences 36: Modern Crystallography III, crystal growth, part 1, Springer-Verlag, Berlin, 1984. [10] H.J. Leamy, G.H. Gilmer, K.A. Jackson, in: J.M. Blakely (Ed.), Surface Physics of Crystalline Materials, Academic Press, New York, 1975. [11] H.M. Cuppen, H. Meekes, W.J.P. van Enckevort, E. Vlieg, H.J.F. Knops, Phys. Rev. B 69 (2004) 245404.