On the Definition of a Monte Carlo Model for Binary Crystal Growth

782

J. Phys. Chem. B 2007, 111, 782-791

On the Definition of a Monte Carlo Model for Binary Crystal Growth J. H. Los,* W. J. P. van Enckevort, H. Meekes, and E. Vlieg IMM Solid State Chemistry, Radboud UniVersity of Nijmegen, ToernooiVeld, 6525 ED Nijmegen, The Netherlands ReceiVed: September 19, 2006; In Final Form: NoVember 20, 2006

We show that consistency of the transition probabilities in a lattice Monte Carlo (MC) model for binary crystal growth with the thermodynamic properties of a system does not guarantee the MC simulations near equilibrium to be in agreement with the thermodynamic equilibrium phase diagram for that system. The deviations remain small for systems with small bond energies, but they can increase significantly for systems with large melting entropy, typical for molecular systems. These deviations are attributed to the surface kinetics, which is responsible for a metastable zone below the liquidus line where no growth occurs, even in the absence of a 2D nucleation barrier. Here we propose an extension of the MC model that introduces a freedom of choice in the transition probabilities while staying within the thermodynamic constraints. This freedom can be used to eliminate the discrepancy between the MC simulations and the thermodynamic equilibrium phase diagram. Agreement is achieved for that choice of the transition probabilities yielding the fastest decrease of the free energy (i.e., largest growth rate) of the system at a temperature slightly below the equilibrium temperature. An analytical model is developed, which reproduces quite well the MC results, enabling a straightforward determination of the optimal set of transition probabilities. Application of both the MC and analytical model to conditions well away from equilibrium, giving rise to kinetic phase diagrams, shows that the effect of kinetics on segregation is even stronger than that predicted by previous models.

1. Introduction It is well-known that equilibrium is often not achieved in nature. Experimentally, at best one can reach equilibrium within a well-defined, finite, and closed system. Of more general applicability are the assumptions that the system will strive for equilibrium and that there is equilibrium locally on a sufficiently small length scale. These are basic concepts in nonequilibrium thermodynamics. A pronounced example of systems where large deviations from equilibrium can be expected are systems involving mixed solid phases, here meant as multicomponent solid solutions, due to the (usually) very low diffusion rate in the solid phase. This hasbeendemonstratedbothexperimentally1-5 andtheoretically.6-15 It certainly holds for molecular systems. In atomic systems, like metals, the diffusion rate may still be appreciable enabling relaxation to the equilibrium state within shorter time scales, especially at high temperature. In this work we think of molecular systems and we will neglect the diffusion in the solid phase completely. In that case the kinetics during crystallization is of decisive importance for the state of the crystallized product just after crystallization, i.e., its composition, possibly with composition gradients, and its stability. Usually equilibrium is not achieved, even if the crystallization is performed very slowly, as has recently been demonstrated for p-dichloro-/pdibromobenzene mixtures.5 The changing liquid composition during the crystallization process leads to a final solid phase built up of shells with a gradual change in composition, even if one assumes that the actual growth composition of the solid phase is in agreement with the thermodynamic equilibrium phase diagram (EPD). * Address correspondence to this author.

Usually it is assumed that at a crystallization temperature slightly below the thermodynamic equilibrium temperature Teq the growth composition is (approximately) in agreement with that prescribed by the EPD. In the present work we investigate this assumption by means of lattice Monte Carlo (MC) simulations, based on the Kossel model.16,17 If the assumption is right, then the growth rate in the MC model should vanish at Teq. Let TMC eq be the temperature where the growth rate vanishes, as actually “measured” in a MC simulation for a vicinal surface with a permanent availability of kinks to rule out the barrier for 2D nucleation. Then, for a thermodynamically consistent set of transition probabilities, to be defined for the MC model, one can be sure that TMC eq e Teq. A solid phase being formed from a liquid phase above its Teq would imply a macroscopic increase of the total free energy, which would be a severe violation of thermodynamics. On the other hand, it is possible to find that TMC eq < Teq, implying no growth within a certain temperature interval below Teq, due to some kind of kinetic barrier. This may be considered as a time scale problem. However, if the growth rate becomes very small at some temperature well below Teq this has important implications, comparable to those for 2D nucleation. Moreover, note that if TMC eq < Teq, then also the solid-phase composition just below TMC eq can be expected to be different from that according to the EPD. In this paper we will show that the before mentioned discrepancy between the results from MC simulations and the EPD can indeed occur, even with a thermodynamically consistent set of transition probabilities. This is due to the surface

10.1021/jp066133+ CCC: $37.00 © 2007 American Chemical Society Published on Web 01/09/2007

A Monte Carlo Model for Binary Crystal Growth kinetics. However, we also show that with a reasonable extension of the MC model it is possible to overcome this discrepancy. In previous work on MC simulation for binary crystal growth,18-22 the matching between the MC model used and the thermodynamic EPD of the issued system, often assumed to be ideal,20-22 was not investigated in detail. In other work12,13,15 TMC eq is fitted to Teq by a single parameter, the melting entropy ∆S. However, we found that this latter approach in general does not give the thermodynamic equilibrium solid-phase composition in the limit of vanishing undercooling, although the differences remain small for systems with not too large melting entropies, for example, semiconductor systems. An additional problem of this approach is that if the melting entropy of the mixture, including a possible excess contribution, is known from experiments, then this imposes a constraint on ∆S that takes away the possibility to adjust ∆S by fitting the thermodynamic Teq. By the melting entropy we mean here the change in total entropy on melting minus the change in the configuration entropy of mixing. Sometimes this is called the vibrational entropy change, although it actually contains also a contribution from the difference in positional entropy in the liquid and the solid phase.23 In our definition of the MC model we will assume that the melting energy and entropy of a mixture are known. They are uniquely related to the pure component melting energies and entropies, and to the excess energy and entropy. Each of these thermodynamic properties are experimentally accessible, e.g., from calorimetric measurements,24-26 or alternatively, they could be obtained from molecular simulation methods.27 Analytical models which are able to reproduce the MC results are quite useful because of the gain in calculation efficiency. And quite often they reveal important mechanisms, providing better understanding. In previous work we presented two analytical models for kinetic segregation, enabling the calculation of kinetic (or nonequilibrium) phase diagrams (KPDs) within certain approximations. One is rooted in merely bulk thermodynamic considerations, the so-called linear kinetic segregation (LKS) model.10 The other so-called mean field kink site kinetic segregation (MFKKS) model28 is more kinetic, based on the kinetics for kink sites, assuming a mean field occupation of the solid neighbors of a particle at a kink site. Both models reproduce the trends in the MC results, but the MFKKS model is generally more accurate. However, the latter model becomes less accurate for nonideal systems at conditions close to equilibrium. Therefore, we here present an extension of this model that goes beyond mean field, which it also suitable for nonideal systems. Away from equilibrium, transport limitations in the liquid phase will start to play a role. In that case the analytical descriptions should be extended with additional equations, relating the surface liquid properties to those in the bulk. These equations have been given elsewhere,29 where they were applied in combination with the LKS model. In section 2, a first MC lattice model for binary growth is defined. In section 3, the EPDs according to this MC model are compared with the thermodynamic EPDs for several model systems. A new analytical kinetic segregation model, supporting the results of section 3, is developed in section 4. In section 5 we propose an extended MC model. Both EPDs and KPDs according to the extended MC model are compared with various analytical models in sections 5 and 6. We conclude with a summary in section 7.

J. Phys. Chem. B, Vol. 111, No. 4, 2007 783 2. Monte Carlo Simulation Model The here described MC simulation model for the growth of a binary crystal described here is an extension of the Kossel model16,17 for pure crystal growth, which is a lattice model with simple cubic symmetry. In its binary version each configuration is represented by a 3D integer matrix H(x,y,z) for integer coordinates x, y, and z, with 0 < x e Nx, 0 < y e Ny, and 0 < z e Nz. The z-direction is chosen perpendicular to the growth surface. Each entry of H has the value 0, 1, and 2 indicating a site occupied by a liquid particle and sites occupied by a solid particle 1 and 2, respectively. Thus, only the configurations of the solid phase are stored. The liquid phase is represented by its given composition through the mole fractions (xl1, xl2) at the liquid-solid interface. During a simulation particles of type 1 and 2 can attach and detach to/from the surface according to the transition probabilities for these events, and respecting the solid-on-solid condition.16 Adopting the random rain model, i.e., assuming a constant attachment rate independent of the surface site where it attaches, the attachment probability for component i ()1, 2) is given by + l + l l R+ i ) Ki xi ) Ki,0 γi xi

(1)

in accordance with chemical reaction rate theory, implying an + l attachment rate for component i equal to K+ i ) Ki,0 γi, where + Ki,0 is a kinetic constant for particle i and where γli is the activity coefficient of component i in the liquid phase. The detachment rate for a particle of component i is fixed by imposing detailed balance, which is assumed to hold for each of the possible events as determined by the nearest neighbor configurations of the particle i. It states Ki,ev

K+ i

(

) exp -

)

∆F ˜ i,ev kBT

(2)

where ∆F ˜ i,ev is the change in the Helmholtz free energy for a selected event (ev), i.e., not including the configurational entropy contribution, which is already implicitly taken into account by the MC procedure. We have

∆F ˜ i,ev ) ∆Ui,ev - T∆S˜ i,ev

(3)

where ∆Ui,ev and ∆S˜ i,ev are the changes in energy and vibrational entropy respectively for the given event. For a given detachment of a solid particle i at the surface, ∆Ui,ev is the broken bond energy given by 2

∆Ui,ev )

∑ j)1

2

ss nss ij ∆φij +

nslij ∆φslij ∑ j)1

(4)

where nss ij is the number of solid-solid ij-bonds being transformed into liquid-solid ij-bonds and nslij is the number of solid-liquid ij-bonds being transformed into liquid-liquid ijss ls bonds during the detachment, and where ∆φss ij ) φij - φij and sl sl ll ∆φij ) φij - φij are the corresponding bond energy changes, the bond energy, defined as a positive energy, with φPP′ ij between a particle i in phase P (l or s) and a particle j in phase P′ (l or s). Following previous work,12,15 we will assume the liquid-solid interactions to be equal to the liquid-liquid interactions, implying ∆φslij ) 0, so that 2

∆Ui,ev )

nss ∑ ij ∆φij j)1

(5)

784 J. Phys. Chem. B, Vol. 111, No. 4, 2007

Los et al.

where we omitted the superscript ss in the bond energy changes for convenience. To simulate the binary crystal growth for a given mixture these bond energy changes ∆φij have to be linked to the thermodynamic properties of the given system. In many cases the energy in phase P can be reasonably well described by

UP ) NP(xP1 UP1 + xP2 UP2 + xP1 xP2 uE,P 0 )

(6)

where NP is the total number of particles in phase P and where the excess energy is described by one single parameter uE,P 0 , which represents the lowest order term in a Redlich-Kister expression30 for the excess energy. The most simple approximation to obtain the bond energies within phase P from the thermodynamic data is to use a mean field (MF) approximation for the total energy reading:

1 P PP P P PP P PP UPMF ) - ZNP(xP1 (xP1 φPP 11 + x2 φ12 ) + x2 (x1 φ21 + x2 φ22 )) 2 1 P PP P P E,P ) - ZNP(xP1 φPP 11 + x2 φ22 + 2x1 x2 φ ) 2

(7)

PP where we used xP1 + xP2 ) 1 and φPP 12 ) φ21 and where Z ) 6 E,P is the number of nearest neighbors and φ is an excess bond energy defined as

1 PP PP φE,P ) φPP 12 - (φ11 + φ22 ) 2

(8)

Equating eqs 6 and 7 directly leads to the bond energies in terms of the thermodynamic properties of the system:

φPP ii

2UPi )(i ) 1, 2) Z

φPP 12

and

UP1 + UP2 + uE,P 0 )Z (9)

The excess energy in the liquid phase is included in the activity coefficient γli appearing in the attachment probabilities (eq 1). On average, nonideality in the liquid phase gives a contribution of kBT ln (γli) to F ˜ i,ev, which results in a factor 1/γli at the righthand side of eq 2. But this factor cancels against the same factor + l on the left-hand side of eq 2, using K+ i ) Ki,0 γi, so that eq 2 can be rewritten as

(

+ ) Ki,0 exp Ki,ev

)

∆F ˜ ′i,ev kBT

(10)

where now F ˜ ′i,ev does not contain the excess energy in the liquid phase anymore. Within this approach the ∆φij in eq 5 have to be taken as

∆φii )

2∆Ui (i ) 1, 2) Z

and

∆φ12 )

∆U1 + ∆U2 - uE,s 0 Z (11)

where ∆Ui ) Uli - Usi (i ) 1, 2) are the pure component melting energies. Alternatively, the liquid-phase excess energy can be included E,I E,s in the ∆φij, replacing -uE,s 0 by u0 - u0 in eq 11, but then the + + l attachment rate should be taken equal to Ki,0 instead of Ki,0 γ i. However, this would be in (non-severe) disagreement with chemical reaction rate theory. If we assume that the excess vibrational entropy is zero, then the total gain in the vibrational part of the entropy after dissolving a crystal with N1 particles 1 and N2 particles 2 should

be equal to

∆S˜ ) N1∆S1 + N2∆S2

(12)

where ∆Si ) ∆Ui/Ti and Ti (i ) 1, 2) are the pure component melting entropies and temperatures of the two components. Assuming the vibrational entropy change ∆S˜ i,ev associated with the detachment of a particle i to be independent of the surface site, eq 12 implies

∆S˜ i,ev ) ∆Si

(13)

This is the most simple choice for ∆S˜ i,ev being consistent with the given thermodynamic property of zero excess vibrational entropy, and which at least gives matching between MC and thermodynamic equilibrium in the pure component limits. For nonideal systems the mean field approximation (eq 7) is not exact, and consequently also the value for ∆φ12 by eq 11 is not exact. Nevertheless, we will use eq 11 as an exact relation to specify ∆φij indirectly by giving ∆Ui (i ) 1, 2) and uE,s 0 . E,s Equation 6 is then not exact for the given uE,s , unless u ) 0. 0 0 Instead, for the liquid phase we do not need the link with a lattice and we will assume that eq 6 is “exact” as an empirical expression. Accordingly, for the examples below with a nonideal regular (i.e., zero excess entropy) liquid solution, γli in eq 1 is given by eq 46 in Appendix A. In the presentation of the results we will use dimensionless model parameters, including the energy parameters

i )

∆Ui kBT2

and

E,P )

uE,P 0 kBT2

(14)

the pure component temperatures Θi and the undercooling ∆Θ, defined as

Θi )

Ti T2

and

∆Θ ) Θeq - Θ

(15)

implying Θ2 ) 1, where Θ ) T/T2 and Θeq ) Teq/T2. Correspondingly, dimensionless pure component melting entropies are defined by σi ) ∆Si/kB ) i/Θi. All the relations involving the “normal” thermodynamic quantities hold also for these dimensionless quantities, but with kB ) 1. 3. Confrontation of Equilibrium Phase Diagrams In general, for the model described in section 2 with given parameters i, E,P, and Θ2 the thermodynamic EPD cannot be calculated exactly, unless uE,s 0 ) 0. This problem is equivalent to the 3D Ising problem, which has not been solved analytically so far. In Appendix A, we give the basic ingredients for calculating approximate EPDs within two approximations: the mean field approximation (MFA) and the quasichemical approximation (QCA). To make our point clear, we first consider a system with ideal E,s mixing properties, i.e., uE,l 0 ) u0 ) 0. In that case, it costs no energy to create two 12-bonds out of a 11- and a 22-bond and it holds that ∆φ12 ) (∆φ11 + ∆φ22)/2, as follows directly from eq 11. At a fixed composition each configuration within the lattice representation of the bulk phase has exactly the same energy so that the equilibrium distribution of the two particles over the lattice points will be completely random. This means that the configurational entropy is ideal and exactly known. Consequently, the thermodynamic EPD is exactly known in this case and can be calculated as outlined in Appendix A. Moreover,

A Monte Carlo Model for Binary Crystal Growth

J. Phys. Chem. B, Vol. 111, No. 4, 2007 785

Figure 1. Snapshot of the surface of a binary Kossel crystal with a diagonal step achieved by tilted periodic boundary conditions (eq 16).

the bond energies from eq 11 give an exact representation of the system. Therefore, we might expect that the EPD resulting from MC simulations will be in agreement with the thermodynamic EPD. To determine the EPD according to the MC model, denoted as the MCEPD from now on, we proceeded as follows. To construct the MC equilibrium liquidus we have to determine the temperature where the growth rate vanishes for a set of liquid compositions. Once this MC “equilibrium” temperature, TMC eq , has been determined, the corresponding point on the MC “equilibrium” solidus is determined by the composition of the solid phase that is formed at a temperature slightly below TMC eq . To exclude inhibition of growth by 2D nucleation barriers, the growth simulations were performed on a surface with a permanent availability of kink sites. In fact, to optimize the statistics we included a whole row of kink sites, realized by a diagonal step (see Figure 1) via the tilted periodic boundary conditions:

Figure 2. Equilibrium phase diagram from MC simulations (filled circles) with ∆S˜ i,ev ) ∆Si (eq 13), compared with the thermodynamic E,s equilibrium phase diagram (full lines) for ideal systems (uE,l 0 ) u0 ) 0) with different values of 1 and 2, as indicated in the graphs, with Θ1 ) 0.8 (upper graphs) and 0.9 (lower graphs). The dotted lines guide the eye for the MC results.

H(Nx+1,y,z) ) H(1,y,z-1) H(0,y,z) ) H(Nx,y,z+1) H(x,Ny+1,z) ) H(x,1,z-1) H(x,0,z) ) H(x,Ny,z+1)

(16)

For given liquid composition, the temperature was chosen somewhat below a first estimate of the equilibrium temperature and the simulation was started. If no growth, or even dissolution, occurred, the temperature was lowered until growth was observed. Then we increased and decreased the temperature repeatedly, reducing the temperature interval between no-growth and growth temperatures to find the closest possible estimate of the “equilibrium” temperature. Each time the simulation for a new temperature was started with the sample from the latest run showing a positive growth rate. All MC results presented in this work were obtained with a surface area equal to (Nx,Ny) ) (50,100). In several cases we performed simulations with larger surfaces, including (100,100)- and (100,200)-surfaces, but these did not show any change in the results. The kinetic constants K+ i (i ) 1, 2) were both taken equal to one, adopting some suitable time unit. Indeed, + for isomorphous components one expects K+ 1,0 ) K2,0. However, a different choice should not affect the “equilibrium” properties. We verified this for several cases taking K+ 1,0 ) + + or 2K ) K . 2K+ 2,0 1,0 2,0 A comparison of MCEPDs with the corresponding thermodynamic EPDs is shown in Figure 2 for ideal systems with several values of the interaction parameters i as indicated in the graphs. Clearly, the thermodynamic EPD is well reproduced by the MC model when 1 ) 2 (left-hand graphs), i.e.,when ∆φ11 ) ∆φ22, but this is not the case when 1 and 2 are different (left-hand graphs). Furthermore, the agreement between the EPDs remains reasonable for the case of Figure 2b where the bond energies are relatively small. However, large deviations occur for the case of Figure 2d where the bond energies are

Figure 3. The sticking fraction as a function of the relative undercooling ∆Θ for the four systems of Figure 2 with xl2 ) 0.25, as indicated by the letters a, b, c, and d. The inset is a zoom-in for small ∆Θ.

large. We note that the relative difference (2 - 1)/(2 + 1) has been taken to be the same in the cases of Figure 2b and Figure 2d. A measure for the growth rate is the sticking fraction, defined as

R ) R1 + R2 )

N+ 1 - N1

N+

+

N+ 2 - N2

N+

(17)

+ where Ri ) (N+ i - Ni )/N is the partial sticking fraction for + component i, with Ni the number of attachments and Ni the number of detachments for that component and N+ ) N+ 1 + N+ 2 the total number of attachments. Figure 3 shows the sticking fraction as a function of ∆θ for the same four systems as in Figure 2, and for a liquid composition equal to xl2 ) 0.25. The inset clearly shows that the sticking fraction for the system of Figure 1d vanishes rather rigorously at approximately ∆θ = 0.012, i.e., significantly above the ∆θ ) 0 expected thermodynamically. The results in Figures 2 and 3 illustrate a general trend in our MC simulations. Unless 1 ) 2, the MCEPD is not in agreement with the thermodynamic EPD. The discrepancy increases for increasing difference between 1 and 2 and also for increasing av ) (1 + 2)/2. Apparently, in these cases there is a barrier for growth, a kind of metastable zone, not attributed to 2D nucleation but due to the chemical potentials of the particles at the surface being different from those in the bulk crystal. In the MC model there is neither diffusion in the solid phase nor direct exchange of particles between the liquid and the solid bulk phase, which would enable the relaxation to the


Los et al.

true equilibrium state. This more or less correctly mimics the reality of vanishing or extremely slow diffusion in the solid phase. The system can (and will) only try to minimize free energy at the surface, where the liquid and solid phase are in contact. But this does not necessarily lead to the situation of thermodynamic equilibrium between the bulk phases, as is shown here.

approximation. The net formation rate of ij-bonds in the direction of the step edge is given by l sk sk Rij ) K+ ˜i xi xj - pijK ij xi

where K ˜ij is the average detachment rate for a kink site particle of type i with a neighbor of type j along the step edge

4. Kink Site Kinetic Segregation Model

K ˜ij )

As has been demonstrated by the MC simulations in the previous section, for large interaction parameters i the surface kinetics become more and more important, overruling thermodynamics even at growth conditions close to equilibrium. An analytical approach to this problem should be based on the surface kinetics. Such an approach is the so-called mean field kink site kinetic segregation (MFKKS) model,28 which is based on the kinetics at kink sites and which assumes a mean field environment of a particle at a kink site. In general, the MFKKS model gives good agreement with MC results, but it becomes inaccurate for nonideal systems at conditions close to equilibrium, mainly due to the mean field approximation. Therefore, here we give an extension of this model beyond mean field, which we will call the extended kink site kinetic segregation (EKKS) model. In the EKKS model, as in the MFKKS model, the net incorporation rate Ri of component i is approximated by

Ri )

Nk(K+ i

xli

-

K ˜i

xsk i )

K ˜i )

2

∑ pijpikpilK-ijkl

with pij the probability that the neighbor along the step edge of a kink particle of type i is of type j, and likewise for pik and pil for the other two neighbors k and l, the one within the upper terrace and the other below in the bulk respectively. Kijkl is the detachment rate for the particle i with neighbors j, k, and l, which according to eqs 3 and 5 equals

(

)

∆φij + ∆φik + ∆φil - T∆S˜ i,ev kBT

(20)

with ∆S˜ i,ev given by eq 13. At growth conditions the composition of the growing solid, xsi , follows from

xs1 xs2

)

l sk R1 K + ˜1 x1 - K 1 x1 ) + l R2 K x - K ˜ - xsk 2

2

2

pikpilK∑ ijkl k,l)1

(23)

l sk sk ˜K+ i xi xj ) pijK ij xi

(24)

Dividing eq 24 for ij ) i2 by that for ij ) i1 one obtains - sk ˜ i1 xsk pi2 K x2 2 ) - sk ) pi1 K ˜ x Q xsk i2

1

(25)

i 1

where we defined Qi as

Qi )

K ˜ i2 K ˜ i1

[

) exp -

]

(∆φi2 - ∆φi1) kBTeq

(26)

with use of eqs 20 and 23. Equation 25 together with the requirement pi1 + pi2 ) 1 gives

pi1 )

xsk 1 Qi

and

sk xsk 1 Qi + x2

pi2 )

xsk 2 sk xsk 1 Qi + x2

(i ) 1, 2)

(19)

j,k,l)1

+ Kijkl ) Ki,0 exp -

2

The first term in eq 22 describes the creation of an ij-bond by the attachment of a particle i at a empty kink site with a neighbor j along the edge, and the second term describes the elimination of an ij-bond by the detachment of a particle i for such a configuration. At equilibrium Rij ) 0 for any pair ij, implying

(18)

where Nk is the average number of kink sites on the surface and xsk i is the average occupation (in mole fraction) of kink sites by component i. K ˜i is the average detachment rate of a particle i at a kink site, given by

(22)

(21)

2

The difference with the MFKKS model is that in the latter model pij is assumed to be equal to the solid bulk mole fraction xsj , and likewise pik ) xsk and pil ) xsl . For the EKKS model we need additional relations for the nearest neighbor correlations pij. We assume these correlations to be equal in each direction for the isotropic Kossel model considered here. From MC simulations we found this assumption to be valid in very good

(27)

Assuming that the average occupation of a neighbor site of a kink particle is equal to the bulk solid composition provides one additional, independent relation: sk s xsk 1 p11 + x2 p21 ) x1

(28)

sk s The second relation xsk 1 p12 + x2 p22 ) x2 is not independent if s s sk sk we impose x1 + x2 ) 1 and x1 + x2 ) 1. For a given liquid composition, the “equilibrium” temperature according to the EKKS model, which we will denote as TKM eq (KM stands for kink model), and the corresponding “equilibrium” solid composition can be computed by solving eq 18 for R1 ) R2 ) 0 and eq 28, inserting K ˜i from eq 19 and pij from eq 27. The eqs sk 18 and 28 together with xs1 + xs2 ) 1 and xsk 1 + x2 ) 1 provide s sk five equations to solve the five unknowns: x1, xs2, xsk 1 , x2 , and KM Teq . At nonequilibrium (growth) conditions eqs 24 and 27 do not hold, but in that case we can derive an alternative relation for pij. Let xsk ij be the probability that the atom at the kink and its neighbor along the step edge are of type i and j, respectively. sk sk We expect the ratio xi2 /xi1 to be equal to Ri2/Ri1, which sk sk together with xij ) xi pij and eq 22 leads to

l sk - sk ˜ i2 xi pi2 K+ i xi x2 - pi2K ) + l sk pi1 K x x - p K ˜ - xsk i

i

1

i1 i1

i

(i ) 1, 2)

(29)


J. Phys. Chem. B, Vol. 111, No. 4, 2007 787 this should not change the MC results near equilibrium too much, since the relevant exchange between liquid and solid is dominated by the exchanges at the kink sites. Therefore, we here introduce a more fundamental variation in the transition probabilities, maintaining the physically plausible random rain model. Our variation of the MC model consists of permitting the vibrational entropy change for events to be environment dependent. As an introduction to this model for a regular solution, we take

∆S˜ ev,1 ) ∆S1 + xs2∆S′12

Figure 4. MCEPDs (symbols and dotted lines) for different choices of ∆S˜ i,ev, compared with the thermodynamic EPD (full lines) and the KMEPD (dashed line) for an ideal system with 1 ) 40 and 2 ) 46. The MCEPD and KMEPD in graph a are for ∆S˜ i,ev ) ∆Si, whereas those in graph b are for ∆S˜ i,ev according to eq 34 with ∆Sii,av ) 0 (circles) and ∆Sii,av ) ∆Sav/6 (crosses) and with ∆S22 11 from eq 42 giving the highest TMC eq in both cases. In graph b the KMEPD and EPD completely overlap.

Inserting the K ˜ij from eq 23 and substituting pi2 ) 1 - pi1, eq 29 leaves a third-order polynomial in pi1 to be solved. Then, for a given liquid composition and temperature T < TKM eq the growth composition is calculated by solving eqs 21 and 28 sk together with xs1 + xs2 ) 1 and xsk ˜1 + x2 ) 1, including K i from eq 19, and pij obtained from eq 29 with pi2 ) 1 - pi1. These equations uniquely fix the four unknowns xsi , xsk i (i ) 1, 2). In Figure 4a the “EPD” according to the EKKS model, which we will call KMEPD, is compared with the MCEPD and the thermodynamic EPD for an ideal system with 1 ) 40, 2 ) 46, and Θ1 ) 0.9. For such relatively large interaction parameters, typically occurring for molecular systems (e.g., fats systems), especially the difference between the MC “equilibrium” solidus and the thermodynamic equilibrium solidus becomes quite large. The solidus according to the EKKS model shows a much better agreement with the MC solidus, showing the relevance of a kinetic description. We also observe that MC TKM eq is somewhat below Teq . In fact, the EKKS model is an analytical description of a restricted MC model (only growth at kink sites), with less freedom for the system to find its lowest free energy state than the MC model. Therefore, we expect that MC TKM eq e Teq e Teq and this is what we see in Figure 4a. 5. Extended Monte Carlo Model The large deviations from the thermodynamic EPDs, as shown in the previous sections, prompt a reconsideration of the MC model, which after all, being a lattice model, is a rather crude approximation of reality. Although we believe that the shown kinetic effects can certainly play a role in practice, one can expect that nature’s ability to find thermodynamic equilibrium is larger than that predicted by the type of MC simulations used here. The deviations are a consequence of the kinetics, and the kinetics depend on the transition probabilities. A first question to be addressed is whether these transition probabilities are uniquely defined within the constraints imposed by the given thermodynamic properties. This is not the case. First of all, the detailed balance eq 2 only gives the ratio of the detachment and attachment probability for a given event, leaving the freedom to choose a different kinetic constant for each of the possible events,31,32 instead of the random rain model. However,

∆S˜ ev,2 ) ∆S2 - xs1∆S′12 (30)

and

instead of eq 13, where ∆S˜ ′12 is a mixing vibrational entropy contribution. Then, the total melting entropy for a crystal of Ns ) Ns1 + Ns2 particles is simply evaluated, using that Nsi ) xsi Ns, as

∆S˜ ) Ns1∆S˜ ev,1 + Ns2∆S˜ ev,2 ) Ns1∆S1 + Ns2∆S2

(31)

as it should be for a mixture with zero excess vibrational entropy. Thus, the expressions 30 for ∆S˜ i,ev are consistent with the thermodynamic property of a regular solution for any value of ∆S˜ ′12. For the implementation in the MC program the ∆S˜ i,ev values in eq 30 are inconvenient, because the probabilities now depend on the composition of the solid phase, which in principle is not a priori known. This means that the determination of the solidphase composition has to be performed self-consistently. Alternative expressions for ∆S˜ i,ev, not requiring selfconsistency, can be formulated by introducing a nearest neighbor dependence in the vibrational entropy. In fact, the environment dependence of the vibrational entropy of a growth unit will primarily be determined by its nearest neighbors. We start by writing the vibrational entropy SPi of a single particle i within a bulk phase P as

SPi

)

P S0,i

+

1

2

PP nPP ∑ ij Sij 2 j)1

(32)

where the factor 1/2 comes from the fact that SPP ij is a pair PP P contribution with SPP 12 ) S21 , and where S0,i is an on-site term. For the vibrational entropy of a solid particle at the surface, Ssurf i , we write s ) S0,i + Ssurf i

1

2

1

2

nslij Sslij ∑nssij Sssij + 2 ∑ 2 j)1 j)1

(33)

Assuming again Sslij ) Sllij, similar as for the bond energies, the change in the vibrational entropy associated with the detachment of a particle i becomes 2

∆S˜ i,ev ) ∆S0,i +

nss ∑ ij ∆Sij j)1

(34)

where ∆S0,i ) Sl0,i - Ss0,i and ∆Sij ) Sllij - Sss ij , so that ∆S12 ) ∆S21. Equation 34 is physically quite reasonable. For positive ∆Sij, the larger the number of broken bonds the larger the entropy change. In fact, it brings into account that a weakly (with few bonds) bonded particle at the surface is still somewhat liquidlike with a higher vibrational entropy than a bulk solid particle. The on-site term ∆S0,i can be interpreted as the change in


Los et al.

positional entropy between a “free”, liquid particle and a particle that is bonded to the surface. According to eq 34 the total gain in vibrational entropy after dissolution of a macroscopic amount of crystal bulk phase is given by 2

∆S˜ )

ss ss Nsi ∆S0,i + Nss ∑ 11∆S11 + N22∆S22 + N12∆S12 i)1

(35)

where Nss ij is the total number of ij-bonds in the crystal. For the pure crystals ∆S˜ should be equal to the pure component melting entropy, implying two relations to be satisfied, namely:

1 ∆Si ) ∆S0,i + Z∆Sii 2

(i ) 1, 2)

(36)

The particle numbers Nsi and the bond numbers Nss ij are not independent, but satisfy the exact relations

1 1 s ss ss s ss Nss 11 ) (ZN1 - N12); N22 ) (ZN2 - N12) 2 2

(37)

∆Sii,av in fractions of ∆Sav ) (∆S1 + ∆S2)/2. Reasonable values for ∆Sii,av for Z ) 6 satisfy 0 < ∆Sii,av < ∆Sav/3 roughly. For ∆Sii,av > ∆Sav/3, the pair free energies ∆φij - T∆Sij will tend to zero or become even negative: e.g., for ∆Sii,av ) ∆Sav/3 and ∆S11 = ∆S22 we have T∆Sii = T∆Sav/3 = ∆Uav/3 = φii at temperatures not too far from the melting point. The above statements are illustrated in Figure 4b, where the MCEPDs for various choices of the parameters ∆S22 11 and ∆Sii,av are compared with the thermodynamic EPDs for the ideal system of Figure 4a. For two choices of ∆Sii,av, and ∆S22 11 adjusted by maximizing TMC eq , the MCEPD are both in almost perfect agreement with the thermodynamic EPD. Thus, the MCEPD is independent of ∆Sii,av. Moreover, in both cases the MC ∆S22 11-values yielding maximal Teq are equal within the error margins. For the KMEPD according to the EKKS model it can be shown analytically that it is independent of ∆Sii,av. The EKKS model considers only kink site events for which it holds that ss ˜ i,ev in eq 20 now becomes nss ii ) 3 - nij and the ∆S

∆S˜ i,ev ) ∆S0,i +

∑ nssik∆Sik ) ∆S0,i + 3∆Sii + nssij (∆Sij -

k)i,j

Substitution of eq 37 into eq 35 and using eq 36 leads to

1 E ss ∆Sii) ) ∆Si + nss ij (∆Sjj - ∆Sii) + nij ∆S ) 2 1 E ss ss ∆Si + nss ij ∆Sii + nij ∆S (40) 2

2

∆S˜ )

E Nsi ∆Si + Nss ∑ 12∆S i)1

(38)

where we have defined the excess melting entropy per 12-bond, ∆SE, as

1 ∆SE ) ∆S12 - (∆S11 + ∆S22) 2

(39)

Thus, for a system with zero excess melting entropy thermodynamic consistency is achieved by taking ∆SE ) 0 and consequently ∆S12 ) (∆S11 + ∆S22)/2. Indeed, we took ∆SE ) -SE,s ) 0 for all model systems in this work. However, also for any nonzero ∆SE, which we assume to be given as a thermodynamic property, consistency is guaranteed by taking ∆S12 in agreement with eq 39, after specifying ∆S11 and ∆S22. These latter two parameters remain as free parameters specifying ∆S˜ i,ev together with the thermodynamic constraints, and introduce a certain indefiniteness in ∆Si,ev, and thus in the transition probabilities. HoweVer, from simulations it turns out that the MCEPD, within the error margin of the MC method, depends only on the difference parameter ∆S22 11 ) ∆S22 - ∆S11 and that a practically perfect agreement with the thermodynamic MC EPD is obtained for that choice of ∆S22 11 for which Teq is maximal. Thus, in this way the MCEPD can be fitted to the thermodynamic EPD by adjusting ∆S22 11. This still leaves one free parameter, which we choose to be ∆Sii,av ) (∆S11 + ∆S22)/2. This second parameter has an effect similar to a scaling of Jackson’s R-factor,33,34 which determines the surface roughness. Roughly speaking, for a large ∆Sii,av, the total pair free energies ∆φij - T∆Sij become small, and weak pair interactions imply a low roughening temperature. In principle, ∆Sii,av should somehow be fitted to experimental data, such as observed roughness or 2D or 3D nucleation rates, involving also the lattice symmetry,35 or it should be estimated via free energy calculations of surface configurations. For the reasons just mentioned it makes sense to use the two parameters ∆S22 11 and ∆Sii,av, instead of ∆S11 and ∆S22. Giving these two parameters uniquely fixes the quantities ∆S0,i and ∆Sij in eq 34 via eqs 36 and 39. We will express both ∆S22 11 and

for i ) 1, 2 and j * i, where we have used eq 39 and where ∆Sjjii ) ∆Sjj - ∆Sii. Equation 40 immediately shows that ∆S˜ i,ev, and thus the KMEPD, depends on ∆Sjjii (apart from the thermodynamic properties) and not on ∆Sii,av, just as in the MC 22 simulations. If we maximize TKM eq as a function of ∆S11, then the maximum occurs at the same value as that for which TMC eq is maximal in the MC simulation, both for ideal and nonideal systems. This means that we can use the EKKS model to obtain the parameter S22 11 required for the MC model to be in agreement with thermodynamic equilibrium. In addition, at this 22 sk optimal value for ∆S22 11, say ∆S11,m, the kink composition xi according to the EKKS model appears to be exactly equal to the solid bulk composition xsi . This is not so surprising, because in that case the configurational entropy in the occupation of kink sites by the two particles is equal to the ideal configurational entropy in the bulk, which makes the change in the chemical potential for a kink site event equal the bulk chemical potential difference, so that equilibrium for kink site events implies equilibrium between the bulk phases. The s 22 observation that xsk i ) xi facilitates the determination of ∆S11,m from the EKKS model is as follows. For an ideal system it holds that pij ) xsj , as predicted by the EKKS model for zero excess energy parameters. In that case, s eqs 25 and 26, with xsk i ) xi , give rise to K ˜ i2

[

) 1 w exp -

K ˜ i1

]

(∆φi2 - ∆φi1) - Teq(∆Si2 - ∆Si1) )1 kBTeq (41)

which together with eqs 11 and 39 for zero excess parameters yields jj ) ∆Sii,m

(∆φjj - ∆φii) ii ) -∆Sjj,m Teq

(j * i)

(42)


Figure 5. MCEPDs (symbols and dotted lines) for two nonideal systems with 1 ) 2 ) 40, E,l and E,s as indicated in the graphs, and different choices of ∆Si,ev. The thermodynamic EPD is calculated within the two approximations outlined in Appendix A, namely the MFA (dashed lines) and the QCA (full lines). Also the KMEPD for maximal TKM eq is drawn (dashed dotted lines), but it overlaps completely with liquid-solid equilibrium lines of the EPD within the QCA. The open circles are the MC results for ∆Si,ev ) ∆Si, whereas the filled circles and crosses are the MC results for ∆Si,ev according to eq 34 with ∆Sii,av ) 0 and ∆Sii,av ) ∆Sm,av/6, respectively, and with the ∆S22 11 yielding 22 the highest TMC eq . Note that the value of ∆S11 is not a constant but depends on the liquid composition. s For a nonideal system, substituting xsk i ) xi in the EKKS jj model equations and considering ∆Sii,m as an unknown leaves jj , which four unknowns, namely, xsi (i ) 1, 2), Teq, and ∆Sii,m can be solved from eq 18 for Ri ) 0 (i ) 1, 2), eq 28, and xs1 + xs2 ) 1, and using pij obtained by solving eq 29 with pi2 ) 1 - pi1. As an illustration for a nonideal system, Figure 5 shows a comparison of phase diagrams for a eutectic system of two components with equal pure component thermodynamic properties, as for racemic mixtures, with a solid-phase excess energy parameter yielding phase separation but with some mixing in the solid phase, for an ideal (upper graph) and a nonideal liquid phase (lower graph). The results in Figure 5b show that taking into account nonideality in the liquid phase via the activity coefficients in the attachment probabilities (eq 1) is appropriate. In both graphs there is very good agreement between the MCEPD and the KMEPD for maximal TKM eq and the thermodynamic EPD, calculated within the quasichemical approximation (QCA) (see Appendix A), which should be very close to the true EPD in this case. For the simple choice with ∆Si,ev ) ∆Si the MC simulations predict much more mixing in the solid phase. As for the ideal systems, the MC equilibrium properties are unaffected by variations in ∆Sii,av. is It is worth noting that the KMEDP for maximal TMC/KM eq practically equal to the EPD within the QCA. This suggests that these two approximations to equilibrium are equivalent. Although the mathematical expressions are completely different, both models include a nearest neighbor correlation beyond mean field, neglecting correlations beyond nearest neighbors. Thus, the EKKS model provides a good approximation of the EPD. At the same time it provides the value for ∆Sjjii yielding optimal agreement between this analytical EPD and the MCEPD.

6. Confrontation of Kinetic Phase Diagrams The logical extension of the previous approach of determining ∆Sjjii by maximizing TMC eq to nonequilibrium (growth) condi-


Figure 6. KPDs (symbols) for undercoolings ∆Θ ) 0.03 (a) and 0.06 (b) from MC simulations compared with those according to the EKKS model (full lines) for the ideal system of Figure 2d. The open circles represent the input set of liquid compositions in the MC simulations. In both models ∆S˜ i,ev was taken according to eq 34 with the ∆S22 11 yielding maximal growth rate and with, for the MC simulations, ∆Sii,av ) 0 (filled circles) and ∆Sii,av ) ∆Sm,av/6 (crosses). The dashed-dotted lines represent the results according to the LKS model. As a reference, also the thermodynamic EPD (dashed line) is drawn.

tions is to impose the growth rate (sticking fraction) to be maximal. Actually, the latter constraint covers both nonequilibrium and equilibrium, considering equilibrium as the limit of vanishing undercooling. Also in this case, to obtain the optimal value for ∆Sjjii, i.e., the one yielding maximal growth rate, one can use the EKKS model for nonequilibrium since the optimal ∆Sjjii according to this model is practically equal to that for the MC model for any value of ∆Sii,av. Here it should be noted that in the MC model the growth rate depends on ∆Sii,av. Large ∆Sii,av-values yield a high kink density, giving large growth rates. But for any given ∆Sii,av-value the maximal growth rate in the MC simulations occurs at practically the same ∆Sjjii-value, which is equal to the value that maximizes the growth rate per kink site in the EKKS model. Note that in the EKKS model we do not need to calculate the kink density in eq 18 since it cancels on the right-hand side of eq 21. For the maximal growth rate solutions according to the EKKS model the composition at the kink site is not equal to the bulk composition as previously for equilibrium, so that the above simplification for finding the maximal TKM eq cannot be applied in this case. Instead, we have to perform the maximization using an appropriate numerical method. Figure 6 shows kinetic phase diagrams (KPDs) for two different undercoolings resulting from MC simulations (MCKPDs) compared to those according to the EKKS model for the ideal system of Figure 2d, with ∆S22 11 yielding the maximal growth rate, being the same in both models. The figure also shows the prediction by the linear kinetic segregation (LKS) model10 (see Appendix B). In the MC simulations two values for ∆Sii,av were used, namely, ∆Sii,av ) 0 and ∆Sii,av ) ∆Sii,av/ 6. As for the MCEPDs, also the MCKPDs are hardly affected by a change in ∆Sii,av. For both undercoolings the KPDs according to the MC and the EKKS model match almost perfectly, whereas the LKS model, although giving the right trends, shows an appreciable discrepancy. Similar performances of the analytical models are observed for the nonideal systems of Figure 5, as shown in Figure 7. For the relatively small undercooling ∆Θ ) 0.0025, the KPD shows two kinetic solidus branches, corresponding to a 1-rich and a 2-rich growing phase. At the equimolar liquid composition xli


Figure 7. KPDs (symbols) for undercoolings ∆Θ ) 0.0025 (in both graphs) and 0.02 (in both graphs) from MC simulations compared with those according to the EKKS model (full lines) for the nonideal systems of Figure 5a,b. The open circles represent the input set of liquid compositions for the MC simulations. ∆S˜ i,ev was taken according to eq 34 with the ∆S22 11 yielding the maximal growth rate and with, in the MC simulations, ∆Sii,av ) 0 (filled circles) and ∆Sii,av ) ∆Sm,av/6 (crosses). The latter results (crosses) are for a surface with a step parallel to [100]. The dashed lines represent the EPD and the dashed-dotted lines represent the KPD according to the LKS model, both calculated within the quasichemical approximation (QCA).

) 0.5 the EKKS model gives two solutions with local minima for the growth rate, represented by the endpoints of the two solidus branches. Due to the symmetry of the system these two minimal growth rates are equal. For the same reason of symmetry, the MC simulation can be restricted to compositions from xl2 ) 0 to 0.5. Starting from the left pure component side, for each run xl2 was increased and the simulation was started using the binary Kossel crystal from the previous run as the initial configuration. For ∆Θ ) 0.0025, arriving at xl2 ) 0.5 we find the steady-state growth with the solid-phase composition that is closest to the composition from the previous run, corresponding almost perfectly to the end point of the left-hand solidus branch for the EKKS model. Apparently, the phase on the right-hand branch with equal growth rate does not nucleate in the given case. Also in this nonideal case the MCKPD is not affected by a change in ∆Sii,aV. As another variation, Figure 7 includes results for a surface with a step along the [100] direction. Now the formation of kinks is required in the growth process, but as we see this does not change the segregation, showing the transferability of the EKKS model to surfaces with no “a priori” availability of kink sites. 7. Summary After having shown the impact of the surface kinetics on the liquid-solid phase behavior (segregation) for binary solid solution growth, even at near-equilibrium conditions and especially for high melting entropy systems, we have addressed the problem of how to define the transition probabilities in a lattice Monte Carlo (MC) model for such systems. Basically, there exists a certain indefiniteness in the chemical potential change of an individual component during its transition from the liquid to the solid phase (or visa versa). Only the total sum of these chemical potential changes over many events and both components during growth or dissolution of a mixed crystal is fixed by thermodynamic constraints, as imposed by the given pure component and excess thermodynamic properties. Within this context, we formulated an extended MC model where this arbitrariness is built in by introducing a neighbor

Los et al. dependence in the vibrational entropy change ∆S˜ i,ev during the transition of a component i. In this model ∆S˜ i,ev is split up into an on-site term ∆S0,i, which reflects the positional entropy change, and pair contributions ∆Sij describing the vibrational entropy change. Within the thermodynamic constraints this introduces two extra parameters in the transition probabilities, which we have taken as ∆S22 11 ) ∆S22 - ∆S11 and ∆Sii,av ) (∆S11 + ∆S22)/2. It turns out that the MC “equilibrium” and kinetic phase diagrams (MCEPDs and MCKPDs) depend, within the error margins, only on ∆S22 11. The arbitrariness in the choice can be overcome by imposing the maximal growth of ∆S22 11 rate, i.e., the fastest decrease in the total free energy of the system. This additional constraint leads to a “perfect” agreement between the MCEPD and the thermodynamic EPD and is therefore a plausible way to fix ∆S22 11. The other parameter, ∆Sii,av, has influence on the roughening temperature, the growth rate, and the width of the metastable zone. If no further experimental or theoretical information on these properties is available, then a standard choice would be ∆Sii,av ) 0. However, we expect that ∆Sii,av > 0 in general, which remains to be investigated in further detail. In principle a similar extension (with just one parameter ∆S11) can also be applied for pure systems, in order to correct for the extremely high undercoolings required in MC simulations to enable 2D nucleation growth for large ∆S (molecular) systems.31 We also developed an analytical kinetic model based on integration at kink sites, which includes a nearest neighbor correlation beyond the mean field approximation. This so-called extended kink site kinetic segregation (EKKS) model, which includes the previous extension of the MC model, has its own “EPD”, called KMEPD, which analytically depends only on ∆S22 11, supporting the MC results. It allows for a straightforward determination of the optimal ∆S22 11-value, which is shown to be the same as in the MC model. The EKKS model is able to predict the liquid-solid equilibrium lines (not known analytically for nonideal lattice-based systems) with the same accuracy as the quasichemical approximation (QCA). KPDs for different systems, both ideal and nonideal systems, according to the EKKS model and a previously presented model,10 the so-called linear kinetic segregation (LKS) model (see Appendix B), have been calculated and compared to those obtained from the MC model. The KPDs from the EKKS model reproduce the MC results well for any choice of the system parameters and the ∆S22 11-value, although the best agreement occurs for the ∆S22 11-value giving maximal growth rate. Compared to these results for the optimal ∆S22 11-value, the performance of the LKS model within the QCA, which is completely independent of ∆S22 11, is quite reasonable but clearly less accurate, especially for systems with a large difference in pure component melting heats. Essentially, the reduction of the segregation for increasing undercooling or, more generally, the role of kinetics as predicted by the MC and the EKKS models is even stronger than that predicted by the LKS model. Acknowledgment. This work has been sponsored by Stichting Technische Wetenschappen (STW), The Netherlands. Appendix A: Equilibrium Two-phase equilibrium lines in an EPD follow from imposing equality of the component’s chemical potentials in the two phases for each component, leading to s1 s2 s2 γli xli ) γsi xsi Ci(T) and γs1 i xi ) γi xi

(i ) 1, 2) (43)



for liquid-solid and solid-solid (s1-s2) equilibrium, respectively, where Ci(T) is a pure component thermodynamic factor and the activity coefficients γPi are related to the total excess free energy FE,P, given through

Ci(T) ) exp

[

( )]

∆Ui 1 1 kB Ti T

and kBT ln (γPi ) )

( ) ∂FE,P ∂NPi

P,T,NjP

(44) respectively, with NPi the amount of component i in phase P. For the liquid phase we will assume a regular solution, implying36,37

γli ) exp

( ) (xPj )2uE,l 0 kBT

(j * i)

(45)

For the solid phase, within the lattice model described in section 2, FE,s is not known exactly. Here we give two analytical approximations. Mean Field Approximation (MFA). Within the MFA, FE,s ) ZNsxs1 xs2τE,s, yielding

γsi ) exp

(

)

(xsj )2ZτE,s kBT

(j * i)

(46)

where we used xsi ) Nsi /Ns and Ns ) Ns1 + Ns2 and where τE,s ) φE,s - TSE,s is the excess free energy parameter, with SE,s ) ss ss -∆SE ) (∆Sss 11 + ∆S22)/2 - ∆S12. Quasichemical Approximation (QCA). The QCA38,39 includes a nearest neighbor correlation beyond mean field. Within the QCA, FE,s is given by

(

[ ( ) ( ) ( )])

xs11 1 FE,s ) ZNs 2xs1 xs2τE,s + kBT (xs1)2 ln s 2 + 2 (x1) (xs2)2 ln

xs22

(xs2)2

+ 2xs1 xs2 ln

xs12

2xs1 xs2

(47)

where xsij are the equilibrium fractions of ij-bonds given by

4xs1 xs2

xs12 ) 1+

x ( ( ) )

and

τE,s 1 + 4 exp - 1 xs1 xs2 kBT 1 xsii ) xsi - xs12 2

(i ) 1, 2) (48)

Away from the critical temperature for phase separation both models give remarkably good results, when compared with numerical results from MC simulation for the Ising problem, but the accuracy of the QCA remains much better when approaching the critical temperature. Appendix B: LKS Model The linear kinetic segregation (LKS) model10 is based on a nonequilibrium detailed balance equation for bulk chemical potential differences. For given liquid composition and temperature, the solid composition follows from

xs1 xs2

)

l l s s K+ 1,0 (γ1 x1 - γ1 x1C1(T)) l l s s K+ 2,0 (γ2 x2 - γ2 x2C2(T))

(49)

and xs1 + xs2 ) 1, where Ci(T) and γPi are defined by eq 44, with FE,s either in the MFA or the QCA. References and Notes (1) Reynolds, P. A. Mol. Phys. 1975, 29, 519. (2) Campisano, S. U.; Floti, G.; Baeri, P.; Grimaldi, M. G.; Rimini, E. Appl. Phys. Lett. 1980, 37, 719. (3) Baeri, P.; Poate, J. M.; Campisano, S. U.; Floti, G.; Rimini, E.; Cullis, A. G. Appl. Phys. Lett. 1980, 37, 912. (4) White, C. W.; Wilson, S. R.; Appleton, B. R.; Young, F. W., Jr. J. Appl. Phys. 1980, 51, 738. (5) Matovic, M.; van Miltenburg, J. C.; Los, J. H. CALPHAD 2006, 30, 209. (6) Aziz, M. J. J. Appl. Phys. 1982, 53, 1158. (7) Aziz, M. J. J. Appl. Phys. Lett. 1983, 43, 552. (8) Chvoj, Z. Cryst. Res. Technol. 1986, 21, 1003. (9) Jackson, K. A.; Gilmer, G. H.; Temkin, D. E. Phys. ReV. Lett. 1995, 75, 2530. (10) Los, J. H.; van Enckevort, W. J. P.; Vlieg E.; Flo¨ter, E. J. Phys. Chem. B 2002, 106, 7321. (11) Los, J. H.; van Enckevort, W. J. P.; Vlieg, E.; Flo¨ter, E.; Gandolfo, F. G. J. Phys. Chem. B 2002, 106, 7331. (12) Jackson, K. A.; Gilmer, G. H.; Temkin, D. E.; Beatty, K. M. J. Cryst. Growth 1996, 163, 461. (13) Jackson, K. A.; Beatty, K. M.; Gudgel, K. A. J. Cryst. Growth 2004, 271, 481. (14) Beatty, K. M.; Jackson, K. A. J. Cryst. Growth 1997, 174, 28. (15) Beatty, K. M.; Jackson, K. A. J. Cryst. Growth 2004, 271, 495. (16) Gilmer, G. H.; Bennema, P. J. Cryst. Growth 1972, 13, 148. (17) Gilmer, G. H.; Bennema, P. J. Appl. Phys. 1972, 43, 1347. (18) Cherepanova, T. A.; van der Eerden J. P.; Bennema, P. J. Cryst. Growth 1978, 44, 537. (19) Cherepanova, T. A.; Dzelme, J. B. Cryst. Res. Technol. 1981, 16, 399. (20) Matsumoto, N.; Kitamura, M. J. Cryst. Growth 2001, 222, 667. (21) Matsumoto, N.; Kitamura, M. J. Cryst. Growth 2002, 237 (1), 51. (22) Kitamura, M.; Matsumoto, N. J. Cryst. Growth 2004, 260, 243. (23) Ubbelohde, A. R. Melting and Crystal Structure; Clarendon Press: Oxford, UK, 1965. (24) Bouwstra, J. A.; Brouwer, N.; van Genderen, A. C. G.; Oonk, H. A. J. Thermochim. Acta 1980, 38, 97. (25) Bouwstra, J. A.; Oonk, H. A. J. CALPHAD 1982, 6, 11. (26) van der Linde, P. R.; Bolech, M.; den Besten, R.; Verdonk, M. L.; van Miltenburg, J. C.; Oonk, H. A. J. J. Chem. Thermodyn. 2002, 34, 613. (27) Frenkel, D.; Smit, B. Understanding Molecular Simulations; Academic Press: San Diego, CA, 1996. (28) Los, J. H.; van den Heuvel, M.; van Enckevort, W. J. P.; Vlieg, E.; Oonk, H. A. J.; Matovic, M.; van Miltenburg, J. C. CALPHAD 2006, 30, 216. (29) Los, J. H.; Matovic, M. J. Phys. Chem. B 2005, 30, 14632. (30) Redlich, O.; Kister, A. T. Ind. Eng. Chem. 1948, 40, 345. (31) Boerrigter, S. X. M.; Josten, G. P. H.; van de Streek, J.; Hollander, F. F. A.; Los, J. H.; Cuppen, H. M.; Bennema, P.; Meekes, H. J. Phys. Chem. A 2004, 108, 5894. (32) Cuppen, H. M.; Meekes, H.; van Enckevort, W. J. P.; Bennema, P.; Vlieg, E. Surf. Sci. 2003, 525, 1. (33) Jackson, K. A. Mater. Sci. Eng. 1984, 65, 7. (34) Beatty, K. M.; Jackson, K. A. J. Cryst. Growth 2000, 211, 13. (35) van Enckevort, W. J. P.; van der Eerden, J. P. J. Cryst. Growth 1979, 47, 501. (36) Hildebrand, J. H. J. Am. Chem. Soc. 1929, 51, 66. (37) Guggenheim, E. A. Proc. R. Soc. London, Ser. A 1935, A148, 304. (38) Rushbrooke, G. Proc. R. Soc. London, Ser. A 1938, 166, 296. (39) Guggenheim, E. A. Proc. R. Soc. London, Ser. A 1938, 169, 134.

On the Definition of a Monte Carlo Model for Binary Crystal Growth

Recommend Documents