Nucleation near the eutectic point in a Potts-lattice gas model , , Vishal Agarwal and Baron Peters

Citation: The Journal of Chemical Physics 140, 084111 (2014); doi: 10.1063/1.4865338 View online: http://dx.doi.org/10.1063/1.4865338 View Table of Contents: http://aip.scitation.org/toc/jcp/140/8 Published by the American Institute of Physics

Articles you may be interested in Nucleation in a Potts lattice gas model of crystallization from solution The Journal of Chemical Physics 131, 184101 (2009); 10.1063/1.3250934

THE JOURNAL OF CHEMICAL PHYSICS 140, 084111 (2014)

Nucleation near the eutectic point in a Potts-lattice gas model Vishal Agarwal1,a) and Baron Peters1,2,b) 1 2

Department of Chemical Engineering, University of California, Santa Barbara, California 93106, USA Department of Chemistry, University of California, Santa Barbara, California 93106, USA

(Received 21 October 2013; accepted 24 January 2014; published online 27 February 2014) We use the Potts-lattice gas model to study nucleation at and near the eutectic composition. We use rare-event methods to compute the free energy landscape for the competing nucleation products, and short trajectories at the barrier top to obtain prefactors. We introduce a procedure to tune the frequency of semigrand Monte Carlo moves so that the dynamics of a small closed system roughly resemble those of an infinite system. The non-dimensionalized nucleation rates follow trends as predicted by the classical nucleation theory. Finally, we develop corrections that convert free energy surfaces from closed (canonical) simulations into free energy surfaces from open (semigrand) simulations. The new corrections extend earlier corrections to now address situations like nucleation at the eutectic point where two products nucleate competitively. © 2014 AIP Publishing LLC. [http://dx.doi.org/10.1063/1.4865338] I. INTRODUCTION

Nucleation is the process by which a more stable phase forms from a metastable parent phase.1–7 Although nucleation has been studied for over 100 years, understanding nucleation mechanism(s) and kinetics still poses a significant challenge. Single component nucleation processes have been the focus of many computational and theoretical studies.8–13 More recently some interest has shifted towards solute precipitate nucleation.14–21 For most compositions, liquid to solid nucleation yields a single solid precipitate when the temperature is lowered through the liquidus. Nucleation at the eutectic composition, where two liquidus curves cross, occurs differently. Lowering the temperature through the eutectic point causes either of the two stable solids to nucleate.6, 7 Interestingly, nucleation and subsequent growth of the first eutectic nucleation product induces the nucleation of the second. Once initiated, the two phases grow together into exquisite microphase-separated microstructures.6, 22 Eutectic microstructures have a number of useful properties. First, the liquid-to-solid transformation of a eutectic mixture occurs with no macroscopic phase separation. As a consequence, mass transfer limitations to re-melting are minimal. Second, eutectic mixtures have low melting points relative to mixtures of other compositions. These properties are exploited for a number of applications. For example, a waterethylene glycol eutectic mixture is used as an antifreeze in cars,23 eutectic alloys are used for soldering,24 sodium and potassium eutectic alloys are used as a coolant in nuclear reactors,25 some printer inks have a eutectic composition,26 eutectic salt mixtures are used for storing thermal energy,27, 28

a) [email protected] b) [email protected]

0021-9606/2014/140(8)/084111/11/$30.00

and a chloride/urea eutectic mixture is used for synthesizing porous materials.29 Recently, Ganagalla and Punnathanam30 applied umbrella sampling to compute the nucleation barriers for a eutectic binary mixture of hard spheres. They found an increase in the nucleation barriers as the fluid composition approaches the eutectic. Interestingly, they found the critical nucleus having different composition than the pure solid phases indicating a fundamentally different mechanism at the eutectic point. Classical nucleation theory (CNT)31–37 provides a qualitative understanding of nucleation, but CNT assumes macroscopic properties in the interior and surface of the nucleus. Additionally, CNT assumes spherical nuclei4, 5 and ignores potentially important concentration fluctuations in the neighborhood of the nucleus.38 These problematic assumptions cause rates from CNT to deviate from experimental rates by many orders of magnitude.2, 4, 5 A more general route to the nucleation rate is to compute mean first passage times (MFPTs). The MFPT is the average time required for the system to leave the metastable state for the first time.39–41 Nucleation rates can then be simply obtained from the inverse of MFPTs. A number of computational studies suggest that rates obtained from MFPTs confirm qualitative predictions of the CNT while not relying on the more problematic assumptions.8, 42–44 There are two ways to compute MFPTs: one is the direct route that uses brute force trajectories and the other uses the statistical mechanics of rare events. Wedekind et al.43 computed MFPTs by tracking the largest nucleus as a function of time along thousands of MD trajectories. They showed how a critical cluster size and the nucleation rate could be identified directly from the obtained MFPT data. The method of Wedekind et al. still relies on brute force MD simulations, and therefore is limited to large undercoolings. For more moderate undercoolings, the nucleation barrier becomes large and a brute force trajectory in a small simulation volume will spend an extremely long time in the metastable state. Therefore, at

140, 084111-1

© 2014 AIP Publishing LLC

084111-2

V. Agarwal and B. Peters

moderate undercoolings, rare-event approaches become the most feasible way to obtain MFPTs and rates. Prior to the onset of growth, supersaturation of the liquid in experiments is virtually constant. A small nucleus in a large system at constant supersaturation exchanges solutes freely (and perhaps slowly) with its surroundings. Its thermodynamics are best modeled in open (grand or semigrand) ensembles. In a closed system, all nuclei including pre-critical nuclei must deplete the surrounding solute concentration to grow.45–48 The depletion effect is potentially important because nucleation rates and even mechanisms are particularly sensitive to supersaturation. To overcome this problem, one could use larger simulations,48 but large calculations are more computationally expensive. Alternatively, simulations with controlled chemical potentials require a large number of successful particle insertion moves to reach and sample compositional equilibrium. For condensed phases there are typically no openings for insertion of new particles. Each insertion becomes a non-trivial calculation, and simulations requiring many thousands of insertions are often not feasible. Therefore, a correction to estimate the open system results from those of a closed system simulation is required. Reguera et al.45 derived expressions for vapor-liquid nucleation in a small finite system. Grossier and Veesler47 obtained simple expression for correcting supersaturation when a single phase is nucleating from a metastable liquid in a small finite volume. In these previously considered finite systems, growth of the post-critical nucleus depletes the driving force and equilibrium is reached at a finite nucleus size. Previous corrections for finite size effects in simulations of nucleation in small closed systems concern depleted concentrations of the solute. However, in a eutectic mixture, nucleation of one product concentrates the other and induces its nucleation. We must therefore consider that precipitation and depletion of component-1 enriches the local concentration of component-2 which can also solidify. This work develops computational methods for computing solute precipitate nucleation barriers and rates in small systems with corrections designed to approximately recover the properties of infinitely large systems. We consider both thermodynamic and dynamical differences between small finite volume systems and infinitely large systems. Both open and closed ensembles can be simulated in the Potts-lattice gas (PLG) model, so we use this simple model to illustrate and test the new strategies. To approximate the dynamics of an open system, we show how semigrand Monte Carlo (MC) move frequencies can be tuned to approximate the natural dynamics of concentration fluctuations by diffusion. To estimate the free energy barriers in an open system using data from closed simulations in a finite volume, we generalize previous theories for nucleation in finite volumes. The remainder of this article is organized as follows: in Sec. II we present model details, simulation details, phase diagram calculations, and free energy surfaces. In Sec. III, we provide a novel methodology to model open system dynamics using finite semigrand MC simulations. In Sec. IV, we provide details of the methodology to obtain nucleation rates based on MFPTs. Finally, in Sec. V we derive simple corrections to the free energy surface

J. Chem. Phys. 140, 084111 (2014)

obtained from closed system simulations and compare it to the free energy surface from the open system. II. METHODOLOGY A. Model

Nucleation was studied using the hybrid PLG model.49, 50 The interactions in the PLG model depend on the species at neighboring sites, and also on the local lattice vector orientation at the neighboring sites. For a k = 2 component system, the particles in the PLG model interact with the following hamiltonian:49 ⎤ ⎡ 2   ⎣ [δm(i),k δm(j ),k (Gk + δs(i),s(j ) Ak )]⎦ , (1) H (r ) = − k=1

i,j 

where i, j represents the sum over the nearest neighbors, m(i) is the species at lattice site i (i.e., for a two component system m(i) = 1 or 2), and s(i) has been interpreted as the local lattice vector orientation at lattice site i.49, 51, 52 When the system is a liquid, there are no aligned lattice vectors. As the system becomes solid, large domains of sites will appear with coherent lattice orientations. If the neighboring sites in the PLG hamiltonian have the same identity k, the system is stabilized by Gk , and is additionally stabilized by Ak if the orientations are also same. The parameters Gk control the solubility of the components in each other, while parameters Ak control the melting/freezing transitions of the pure components. All simulations were carried out using two components (denoted as “1” and “2”), and each site having six possible orientations. Qualitatively different phase diagrams can be obtained by varying values of Gk and Ak as shown by Duff and Peters.49 For the present study, we used the PLG with the following parameter values: G1 = G2 = 0 and A1 = A2 = 1.275. For these choices of the PLG parameters, the system exhibits a single eutectic at x1 = 0.5 and kB T = 1.365 as explained below. B. Simulation details

We performed MC simulations on a simple cubic lattice periodically repeated in all the three dimensions. Random numbers were generated using the fast Mersenne Twister pseudo-random number generator based on 128-bit operations, which produces statistically good equi-distribution over long periods.53 A total of four types of moves were attempted in this work: (a) nearest-neighbor swap or Kawasaki move54 (i.e., selecting a set of two neighboring sites at random and swapping them), (b) particle orientation change (i.e., changing orientation of a randomly selected site with the new orientation chosen from a uniform distribution of all the orientations), (c) non-local particle swap (i.e., randomly choosing two different lattice sites and swapping them), and (d) particle identity change (i.e., changing identity of a randomly selected site). The simulations are performed with a fixed total number of particles, so particle identity changes are used instead of particle insertions. In all cases orientations at the lattice sites before and after the MC move were preserved unless otherwise noted. Different combination of move types were used

084111-3

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

C. Phase diagram

FIG. 1. Plot of the fugacity ratio vs mole fraction at kB T = 1.32 and kB T = 1.29. A linear fit is also shown.

for different calculations as explained below. The frequency of each move is also detailed below. All the above moves, except for the particle identity swap, were accepted or rejected based on the Boltzmann factor of the configurational change.55 Particle identity exchange from component 1 to component 2 is accepted with a probability Pacc (1 → 2) given by56

  f1 H(1→2) , (2) Pacc (1 → 2) = min 1, exp − f2 kB T where f1 /f2 is the fugacity ratio of component 1 to component 2, and H(1→2) is the energy change for the MC move. The fugacity ratio determines the supersaturation of the components and is a function of both the temperature and mole fraction. For a particular temperature, we determined the relationship between the fugacity ratio and the composition of the metastable liquid by the following procedure: 1. The system was equilibrated at the target temperature for 50 000 sweeps at a constant fugacity ratio. These calculations were performed with an initial configuration taken as the equilibrated structure at kB T = 1.45. 2. We then computed average mole fraction using another 20 000 sweeps for each value of the fugacity ratio. This procedure was repeated for different values of the fugacity ratio to obtain a plot of the fugacity ratio vs mole fraction at a constant temperature. Figure 1 shows a plot of the fugacity ratio vs average mole fraction at a constant value of kB T = 1.32 and kB T = 1.29. 3. We then fitted the fugacity ratio vs mole fraction data to a linear model. The above procedure was repeated at each temperature for which nucleation is studied in this work. We then used these linear parameters to choose appropriate fugacity ratios for the nucleation conditions shown on the temperaturecomposition phase diagram. Nucleation was studied in the temperature (kB T) range of 1.29-1.32 from the metastable liquid having mole fractions (x1 ) in the range of 0.5-0.54. We used simulation boxes of sizes 16 × 16 × 16 sites and 18 × 18 × 18 sites (for under-cooling of kB T > 0.07 and kB T ≤ 0.07, respectively) with periodic boundary conditions.

The details for obtaining the phase diagram are similar to that used in the previous work.49 For completeness, here we briefly summarize the procedure for obtaining the phase diagram. Phase diagram calculations were performed using an elongated periodic box of dimensions 300 × 25 × 25 to minimize the effects of the phase boundary. To obtain an equilibrium configuration, MC simulations were performed in the canonical ensemble with particle orientation changes and particle swaps attempted with equal probability. The moves were accepted or rejected based on the Metropolis criteria.55 The system was first equilibrated from an initially pure solid. Then the two phases were identified by computing the local and medium-range order parameters. The local order parameter is defined as49 M 1  1 θ (i) = + δm(i),m(j ) (δm(i),1 − δm(i),2 ), 2 2M j =1

(3)

where M is the sum of nearest, next-nearest, and next-nextnearest neighbors (for 3-dimensions M = 6 + 12 + 8 = 26 for a cubic lattice), θ (i) = 1 is for pure component 1, and θ (i) = 0 is for pure component 2. We note here that this local order parameter is different from the one used for identifying nucleus size in the free energy calculations as discussed below. The distinction between the two phases was further refined by computing the medium-range order parameter defined as57 ⎞ ⎛ M  1 ⎝ (i) = θ (j )⎠ . (4) θ (i) + 1+M j =1 The dividing value distinguishing the two phases was then identified from the distribution of the medium-range bond order parameters as shown in Figure 2. This dividing value was further used to determine the composition of the two phases as outlined in Duff and Peters.49 Because the local composition fluctuations are large near the eutectic, the gap in  is less pronounced than in previously studied systems.49 For the system studied in this work, the eutectic point lies at x1 = 0.5 and kB T = 1.365 as seen in Figure 3.

FIG. 2. Distributions of the medium-range bond order parameter () for the two phases at temperatures above (kB T = 1.37) and below (kB T = 1.36) the eutectic point. The composition boundary (shown as grey dots) is determined and used for assigning particles in the two phases.

084111-4

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

FIG. 3. Phase diagram for the PLG parameters G1 = G2 = 0 and A1 = A2 = 1.275 in 3-dimensions. The eutectic point, i.e., the intersection of the liquidus lines is approximately at x1 = 0.5 and kB T = 1.365. Filled diamonds represent the points at which homogeneous crystal nucleation was studied.

D. Cluster size: Solid nucleus identification

Computing nucleation rates requires identification of the cluster size or the nucleation coordinate along which the metastable liquid converts to the solid phase. A simple choice is to cluster contiguously neighboring sites having the same identity and same orientation as done in the previous work.58 However, at high concentrations near eutectic (mole fraction ∼0.5), clusters constructed in this way percolate throughout the simulation box.59, 60 Instead, we defined a local bond order parameter as follows: M  (i) = [δm(i),m(j ) δs(i),s(j ) ],

(5)

j =1

where the summation is over nearest, next-nearest, and nextnext-nearest neighbors (i.e., M = 26 on a cubic lattice). The term in the bracket is unity only when the particle i and particle j have the same identity and same orientation. Thus, the maximum value of the bond order parameter is 26. In order to distinguish the “solid-like” cluster from the metastable liquid, we obtained the distribution of the bond order parameter for a metastable liquid phase, for a solid phase, and for two solid phases in contact with each other. Following ten Wolde et al.,61 we considered contiguous particles i and j of the same identity and orientation to be connected if the order parameter for each of them exceeded a certain threshold value. We then compare distributions of the order parameter of different phases to determine the threshold value. Equilibration for these calculations was achieved in 50 000 MC sweeps using particle orientation changes and particle identity changes, as explained above, attempted with equal probability. The distributions of (i) were then obtained from another 20 000 MC sweeps with data collected at every sweep. Histograms of (i) for the three systems at kB T = 1.29 are shown in Figure 4. Comparing distributions for the liquid and single solid phase suggests that a threshold of the order parameter in between 15 and 20 would capture 100% of both the liquid and solid phases. However, we found that choosing a threshold parameter greater than 15 causes problems for

FIG. 4. Distributions of (i) for the metastable liquid (x1 = 0.5), one solid phase and two solid phases (x1 = 0.5) for the PLG model in 3-dimensions at kB T = 1.29. The distributions were computed for 20 000 MC sweeps with data collected at every sweep.

umbrella sampling at the post critical sizes. The system forms a network of domains interconnected by narrow bridges that cut across solid-solid interfaces. The narrow bridges are not recognized by the (i) > 15 criterion. Effectively, the whole system grows to form one large pseudo-grain despite the umbrella sampling restraint on cluster size. To overcome this problem we plotted the distribution of the order parameter for two solid phases present simultaneously in the same box, thus creating an interface and causing the order parameter distribution to increase for moderate values of (i). (i) ≤ 10 still includes essentially all of the liquid phase. We thus considered neighboring particles i and j of the same identity and orientation to be connected and belonging to a solid-like cluster if (i) and (j) both exceeded a value of 10. Jungblut et al.62 have outlined a method based on likelihood maximization63, 64 to more systematically optimize clustering thresholds.

E. Free energy surface

The free energy surface was obtained by monitoring the size of the largest cluster for the 1-rich and 2-rich phases (delg lg noted as “n1 ” and “n2 ,” respectively). The free energy as a function of cluster size in CNT reflects the population of clusters of size n (denoted as “FCNT (n)”). There are subtle but important differences between the CNT free energy and the Landau free energy surface as a function of the largest cluster size.65 The Landau free energy is F (n) = −kB T ln[p(nlg )], where n = nlg is the largest cluster in the simulation box. The CNT free energy, with the self consistency term,66 is F (n) − F (1) = −kB T ln[pop(n)/pop(1)] where pop(n) is the population of n-sized nuclei in the simulation box. The two definitions of the free energy are different: typically there are many clusters that are smaller than the largest cluster. The two free energy profiles become similar for large cluster sizes which are extremely rare and, by Poisson statistics, are even less likely to appear in twins, triplets, etc.17, 65 An important point is that the Landau free energy should be used with the MFPT formula, while the CNT free energy should be used with the CNT rate expression.17 In this work we use the MFPT and the Landau free energy to avoid making the assumptions associated with the CNT.

084111-5

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

The Landau free energy has a local minimum at some non-zero nucleus size corresponding to the typical maximum cluster size in the metastable state. The Landau free energy surfaces, both two- and one-dimensional, in this work are shifted so that the metastable minimum has a free energy of zero. We show below, for kB T = 1.29, x1 = 0.5, and a simulation box of size 16 × 16 × 16, the reference free energy lg minimum computed using the largest cluster size is at n1 lg = 12 and n2 = 12. In what follows, we drop superscript “lg” and all reference to the nucleus size is corresponding to the largest one. For nucleation near a eutectic point, there is a largest cluster for each of the two possible solid phases, i.e., n1 and n2 . Free energy can be computed by sampling the probability, P(n1 , n2 ), of observing the system with given largest nucleus sizes, n1 and n2 , F (n1 , n2 ) = −kB T ln[P (n1 , n2 )] + constant.

(6)

To reduce the computational cost we computed P(n1 , n2 ) by the following procedure. We first applied umbrella sampling67, 68 to compute P(n1 + n2 ). For umbrella sampling, we used hard wall restraints along with partially overlapping windows of size six for the variable n1 + n2 . We then computed the joint probability distribution P(n1 , n2 ) using P (n1 , n2 ) = P [n1 , (n1 +n2 )] = P [n1 |(n1 +n2 )] · P (n1 + n2 ), (7) where P[n1 |(n1 + n2 )] is the probability of finding a nucleus of the largest size n1 for a fixed value of n1 + n2 . Equation (7) shows how the joint probability can be extracted from the single variable umbrella sampling data. The first equality in Eq. (7) comes from the fact that for a given n1 , the value of n1 + n2 determines a unique value for n2 . To reduce the computational expense of computing the order parameter at every step, we performed hybrid MC-MC.69, 70 In this methodology, nucleus size is computed after k steps and these k steps are accepted only if the order parameter falls in the umbrella sampling window. Probability distributions were obtained for 100 000 MC sweeps after an equilibration run of 50 000 MC sweeps. The free energy for nucleation at kB T = 1.29 and x1 = 0.5 is shown in Figure 5. A snapshot of simulation cell for the nucleus sizes of n1 = 104 and n2 = 45 is shown in Figure 6. We note here that the free energy data have been symmetrized (i.e., P (n1 , n2 ) = [P˜ (n1 , n2 ) + P˜ (n2 , n1 )]/2, where P˜ is the raw data). Because we are using the same PLG parameters for both species, the probability of forming the 1-rich and 2-rich phase for x1 = 0.5 is essentially equal. Therefore, sampling points can be increased, and hence the computational cost can be reduced, by symmetrizing the data. We note here that free energy data were symmetrized only to obtain the twodimensional free energy surface for the symmetric x1 = 0.5 case. The saddle point for nucleation of both the 1-rich and 2-rich phases lies at n = 300 due to symmetry. Note that the channels for nucleation are very narrow. They are separated by a high “forbidden” region. Therefore, nucleation of the 1rich and 2-rich phases can be considered as separate competing pathways. Before computing rates, we must consider the dynamics of the nucleus size variables.

FIG. 5. Free energy surface as a function of the largest nucleus size of components 1 and 2 at kB T = 1.29 and x1 = 0.5. The nucleus size is defined as the number of contiguous particles of the same identity and orientation each having (i) > 10. The order parameter is defined by Eq. (5).

III. DYNAMICS IN AN OPEN SYSTEM

For semigrand canonical MC simulations, the dynamics of the overall simulation box composition depends on the frequency of the identity swaps. In this section, we tune the particle identity swap frequency to approximate the dynamics of real composition fluctuations in an open system. Note that the particle identity swap frequency has no effect on the free energy barriers. The compositional dynamics can, in principle, influence the nucleation prefactor. We expect this effect to be small and therefore we did not analyze the relationship between the prefactor and identity swap frequency. A more important reason to simulate as nearly as possible the true compositional dynamics pertains to some important mechanistic artifacts that might otherwise emerge. If identity swaps are

FIG. 6. Snapshot of 16 × 16 × 16 simulation cell during an umbrella sampling run. Nuclei in this snapshot are of sizes n1 = 104 and n2 = 45. For clarity, molecules in the metastable liquid are shown with smaller spheres. Different orientations are represented with different colors.

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

R 2 (τ ) = 6Ds τ,

(8)

where τ is the number of MC sweeps. For these calculations, we performed MC simulations with the nearest-neighbor swap and orientation change moves attempted with equal probability. The system was first quenched from kB T = 1.45 to the target temperature and equilibrated for 1000 MC sweeps. Mean square displacement was then computed for another 20 000 MC sweeps. A plot of mean square displacement vs MC sweep computed at kB T = 1.29 is shown in Figure 7. By applying a linear fit, we extracted a self diffusion constant of Ds = 0.06 sites2 /sweep at kB T = 1.29 and x1 = 0.5. Now we consider the dynamics of composition fluctuations. Let N1 be the number of particles of component one in the simulation box at time t. The compositional relaxation time by semigrand moves can be measured by computing the N1 autocorrelation function and fitting to an exponential decay as

 δN1 (τ )δN1 (0) τ , (9) exp − δN1 (0)2  τrelax where δN1 = N1 − N1 , N1  is the average concentration of N1 in the metastable liquid, and τ relax is the relaxation time. The fluctuations in N1 will decay faster if identity swap moves are attempted more frequently as shown in Figure 8. One can, therefore, find an identity swap frequency at which the relax-

4 2

not sufficiently frequent, nucleus growth will be effectively arrested by solute depletion effects in the early post-critical stages. If identity swaps are too frequent, then a less intuitive, but still serious mechanistic artifact may occur. Berezhkovskii and Zitserman71, 72 describe a saddle point avoidance effect that could be important, particularly for nucleation near the eutectic in binary systems. If the interfacial free energy between the two solid phases at the eutectic point is small, then the forbidden region in the free energy landscape of Figure 5 will also be stabilized. If composition fluctuations are sufficiently slow, it may then become faster to nucleate via the diagonal n1 ≈ n2 , along which the overall composition of the fluid need not change. The critical nucleus would then resemble a Janus particle, already containing a single layer of the typical eutectic microstructure. These two possible mechanistic changes explain why choosing a realistic compositional relaxation time can be important. The interfacial free energy between the two solids in our study is prohibitively large, giving two narrow channels for which saddle avoidance is unimportant and for which rates can be computed separately. However, post-critical growth can still be artificially limited if the identity swaps are too infrequent. We now present a methodology to tune identity swap frequencies so that the compositional dynamics of a small system roughly correspond to those in the same volume if it were embedded within an infinite system. In a real system with a given rate of solute self diffusion, a solute concentration excess spread over a L × L × L region will decay on a timescale τ ∼ L2 /Ds . Here we are ignoring any differences between the self and transport diffusivity. For a particular temperature, we first computed the self diffusion coefficient (denoted as “Ds ”) from mean square displacement (denoted as “R2 (τ )”) using the Einstein relation

Mean Square Displacement (

084111-6

3

6Ds 2

= 0.36147 τ

2

1

0

0

1

2

3

4

5

6

7

8

9

10

MC Sweep No. (τ) FIG. 7. Mean square displacement as a function of MC sweeps at kB T = 1.29 and x1 = 0.5. The slope of the linear fit is used to calculate the self diffusion coefficient.

ation time becomes equal to L2 /Ds . This is the swap frequency at which the simulation cell will roughly behave as though part of an infinitely large system. The N1 autocorrelation function at kB T = 1.29 and x1 = 0.5 for different swap frequencies is shown in Figure 8. As shown in Figure 8, at 3 identity swaps/sweep the relaxation time becomes equal to the time required to diffuse a cell length L. Identity swap frequency was recalibrated in this way for every nucleation condition. IV. NUCLEATION RATES

Non-equilibrium probabilistic evolution of the nucleus size variables can be studied using the Zeldovich-Frenkel equation33, 37, 73, 74 or the essentially equivalent Smoluchowski equation.17, 75 Induction times can be obtained from the latter by computing the MFPT,39, 40, 43   exp[βF (n)] dn, (10) τMF P T = exp[−βF (n)]dn D(n‡ ) ∪ ∩ where ∪ and ∩ represents integration over the metastable basin and barrier top, respectively. D(n‡ ) is the attachment/

FIG. 8. N1 autocorrelation function plotted as a function of MC sweeps and identity swap frequency at x1 = 0.5 and kB T = 1.29. L2 /Ds is the approximate relaxation time expected when relaxation occurs by true diffusion in an infinitely large system.

084111-7

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

detachment frequency at the barrier top. Nucleation rates (“Jnxn ”) can further be obtained by inverse of the MFPTs and the simulation box volume V as Jnxn = (τMF P T V )−1 .

(11)

In principle, one should use Langer’s multidimensional form of Kramer’s theory to compute rates.76–78 Here, as mentioned in Sec. II E, the channels for nucleation are very narrow, a one-dimensional treatment is sufficient, and the rate through each channel can be computed separately. To obtain the free energy profile for nucleation through the channel leading to phase 1, we integrate out n2 from the channel leading to the 1-rich phase. Formally, we project the channel leading to phase-1 onto variable n1 : F (n1 )  n∗ = −kT ln 0 2 exp[−βF (n1 , n2 )]dn2 . The integration is over values of n2 between 0 and a cutoff n∗2 = 50. This calculation gives the free energy as a function of n1 within an assumption that nucleation of the 2-rich phase has not yet occurred. Similarly, a separate barrier F(n2 ) for nucleation of phase 2 can be computed as long as nucleation of the 1-rich phase has not yet occurred. This is done by integrating over n1 at each value of n2 in the channel leading to 2-rich phase. Again, a cutoff n∗2 = 50 is needed to exclude states where the 1-rich phase has already nucleated. To compute nucleation rates from the F(n1 ) and F(n2 ) profiles, we have fitted the uppermost ∼3kB T of the barrier to a quadratic expression. The dynamical information, i.e., the attachment/detachment frequency (“D(n‡ )”) was estimated by the procedure outlined by Im and Roux,79 δn(τ )2  , D(n‡ ) = 2τ

(12)

where δn(τ ) = n(τ ) − n(τ ). We computed D(n‡ ) by launching trajectories from the ensemble of configurations at the barrier top. The length of the trajectories was short enough to neglect the curvature effects of the free energy surface at the barrier top.58 Averages were computed using 100 trajectories, each 1 sweep long, from each of the 100 sample configurations drawn from the ensemble of 500 configurations that were generated by umbrella sampling. A plot of D(n‡ ) studied at different temperatures and compositions is shown in Figure 9. As expected, for the same mole fraction in the

metastable liquid, D(n‡ ) increases with the increase in tem2/3 perature. At one point in Figure 9, D(n‡ )/(n‡ Ds ) decreases with temperature. The supplementary material shows that the 2/3 decrease actually emerges from the normalization n‡ Ds .80 Attachment/detachment frequency for nucleation is generally by activated hopping in solid solutions, giving the nucleation prefactor an Arrhenius dependence. The nucleation barrier itself decreases with temperature due to increased supercooling. The combined temperature dependence of the attachment/detachment frequency and the nucleation barrier usually causes nucleation rates to have a maximum as a function of temperature. This maximum is often called the “nose” in the rate of nucleation.6, 7 The present model does not take activated hopping into account. The nose can be artificially built into simulations, by associating each sweep with a time interval that depends on temperature by an Arrhenius law. We opted not to build an Arrhenius law into our diffusivities because it would require an arbitrary choice of the activation energy for diffusion. This would in turn arbitrarily influence the temperature corresponding to the “nose.” Our simulations thus do not show the nose at all. Although absolute rates from the CNT are rarely accurate, CNT often does accurately predict trends in the dependence on supercooling and supersaturation. For the present model of eutectic nucleation we find that the MFPT rates closely follow the expected dependence on supercooling from the CNT. The rate of nucleation obtained from the CNT is given by JCNT = ρ1 D(n‡ )Ze−FCNT (n‡ )/kB T ,

(13)

where ρ 1 is the monomer concentration and FCNT (n‡ ) is the (CNT) free energy barrier for nucleation. Z is the so-called Zeldovich factor which takes into account the curvature of the barrier top.3 According to CNT the free energy barrier is dependent on the temperature as81 Tm2 FCNT (n‡ ) ∝ , kB T T (T )2

(14)

where Tm is the freezing temperature and T is the degree of supercooling. Alternatively, this means

 Tm2 JCNT ∝− . (15) ln ρ1 D(n‡ ) T (T )2 A plot of J /L3 D(n‡ ) vs Tm2 /T (T )2 is shown in Figure 10. The results approximately collapse to a single line confirming the trend predicted by the CNT. We note here that the MFPT results make no assumptions on the nucleus shape, capillarity, or interfacial energy as in the CNT. V. FINITE SIZE EFFECTS ON NUCLEATION FREE ENERGY

FIG. 9. Non-dimentionalized attachment/detachment frequency at the barrier top. Nucleus sizes at the barrier top are shown in parentheses.

In Sec. II E, we obtained nucleation free energy barriers by semigrand simulation of a small open system. In Sec. IV we used the free energy barriers to obtain nucleation rates. As mentioned in Introduction, it is often not feasible to directly simulate open ensembles in atomistic simulations. Although, finite size effects on nucleation of single nucleus have been

084111-8

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

where μ is the chemical potential, n is the number of molecules in the nucleus, σ is the surface tension, A is the surface area of the nucleus, and superscripts and subscripts (m) (α) have been defined above. Also, μ(α) 1 = (μ1 − μ1 ) and (β) (β) (m) μ2 = (μ2 − μ2 ) are the driving forces for the nucleation of the α-phase and β-phase, respectively. The thermodynamic driving force for non-ideal liquid solutions can be written as82

 γi xi , (17) μi = kB T ln γi∗ xi∗

FIG. 10. Variation of the nucleation rate as a function of degree of supercooling.

well documented,45–48 finite size effects on competitive nucleation of two phases require some discussion. For example, near the eutectic point nucleation of one phase increases supersaturation of the other thus inducing its nucleation. We begin by deriving thermodynamic equations governing the formation of the two solid phases in a homogeneous liquid in a closed system. Consider a hypothetical system consisting of a metastable liquid phase (denoted with superscript “m”), a nucleus of the 1-rich phase (denoted with superscript “α”), and a nucleus of the 2-rich phase (denoted with superscript “β”) as shown in Figure 11. In what follows, component 1 is denoted with subscript “1,” and component 2 with subscript “2.” The total volume and temperature of the system is constant during the nucleation process. It should be noted that, CNT is often derived using Euler’s equations which are valid because the solute chemical potential in the metastable fluid is assumed to be independent of the nucleus size. In a finite system, when chemical potential in the metastable fluid is a function of the nucleus size, one needs to start with the differential form. Nucleus formation depletes the supersaturation of the environment, so differential growth of the nucleus causes a differential change in the solute chemical potential. The free energy change due to small changes in size for nuclei of two types is given by (β)

(β)

(β)

(β)

(α) (α) (α) dF = −μ(α) 1 dn1 − μ2 dn2 + σ1 dA1 + σ2 dA2 , (16)

where xi and xi∗ are the mole fractions of the ith component in the metastable liquid and reference metastable liquid phase, respectively. γ i and γi∗ are the activity coefficients of the ith component in the metastable phase and reference metastable liquid phase, respectively. The activity coefficient accounts for non-idealities in the chemical potentials for components of a fluid mixture. For an ideal solution, the activity coefficients approach 1. Note that we have defined the driving force according to the metastable fluid properties and properties of the bulk stable solid. We explain below that this could lead to some errors while correcting the free energy surface. Let us now define a difference function, δF, as   (β)  (β)  − Fconf n(α) (18) δF = F∞ n(α) 1 , n2 1 , n2 , (β)

where F∞ (n(α) 1 , n2 ) is the free energy for formation of the (β) nucleus sizes n(α) 1 and n2 when the total system volume is infinity (i.e., supersaturation is constant throughout the nucleation process, which is one of the assumptions of the CNT). (β) Fconf (n(α) 1 , n2 ) denotes the free energy of formation in a closed system when changes in the supersaturation cannot be neglected. Further, substituting Eqs. (16) and (17) in Eq. (18), we obtain



  γ1 x1 γ2 x2 dn1 + kB T ln dn2 , d(δF ) = kB T ln γ1,0 x1,0 γ2,0 x2,0 (19) where subscript “0” refers to the initial values (note that the reference values in Eq. (17) cancel out). We note here that in deriving this equation we have assumed contribution from surface energy will nearly cancel for the same size nucleus in the difference function. Activity coefficients can be approximated using the twosuffix Margules equation83 ln γ1 =

A (1 − x1 )2 , kB T

(20)

ln γ2 =

A (1 − x2 )2 , kB T

(21)

where A is independent of mole fractions but weakly dependent on temperature and is determined empirically in practice. We estimated A from the fugacity ratio vs mole fraction data as shown in Figure 1 from the semigrand simulations (for details see Sec. II B). The fugacity ratio for a liquid state can be written as82 FIG. 11. Schematic figure of the thermodynamic system composed of a metastable liquid, nucleus of 1-rich phase, and nucleus of 2-rich phase.

f1 γ1 x1 f1,0 = , f2 γ2 x2 f2,0

(22)

084111-9

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014)

FIG. 12. Least square fit to the fugacity ratio simulation data using Eq. (23) at kB T = 1.29.

where f1,0 and f2,0 are the fugacity of the pure components. Because we are using a symmetric system, f1,0 = f2,0 . Taking logarithms on both sides and using Eq. (20) yields



 f1 x1 = A(1 − x1 )2 − Ax12 + ln . (23) ln f2 1 − x1 We then performed nonlinear least square fit to the fugacity ratio vs mole fraction data and obtained a value of A = 1.536. The least square fit is shown in Figure 12. The mole fractions can be written as N1,0 − n(α) N (m) 1 , (24) x1 =  (m) 1 (m)  =  (β)  N1 + N2 N1,0 + N2,0 − n(α) 1 − n2 where Ni represents the number of atoms of the ith component, and superscripts and subscripts are defined above. We assume here that phase α pure 1-component and β pure (β) 2-component. Substituting ξ1 = n(α) 1 /N1,0 and ξ2 = n2 /N2,0 into Eq. (24) yields 

1 − ξ1 . (25) x1 = x1,0 1 − ξ1 x1,0 − ξ2 x2,0 Similarly, the mole fraction of the component 2 in the metastable liquid can be written as 

1 − ξ2 x2 = x2,0 . (26) 1 − ξ1 x1,0 − ξ2 x2,0

FIG. 13. Correction to the free energy at x1, 0 = 0.5 and kB T = 1.29.

ξ 2 = 0. On substituting these values into Eq. (27) and expanding the logarithms, the free energy correction (δF) vanishes. The free energy correction is plotted as a function of the sizes of two competing nuclei in Figure 13. As expected, the correction is always negative indicating that the closed-system free energy barrier is always greater than the open-system free energy barrier.45, 47 Also, the correction is zero when the nuclei of both types grow simultaneously. Next, we applied this correction to the closed system free energy surface with A/kB T = 1.536 as shown in Figure 14. We find here that one requires A/kB T = 1.2 to perfectly match the estimated open system free energy from the closed system free energy. The discrepancy may stem from overestimation of the effective driving force for nucleation. We have derived the correction δF using bulk properties for the stable solid phase. Figure 15 shows how the nucleus might instead initially adopt the solid composition that is in equilibrium with the surrounding fluid.84 The intermediate solid phase is a metastable extension of the solidus, and therefore it is potentially metastable towards enrichment to the stable equilibrium solid. If the gap between the metastable solidus and the stable solvus curve is significant, it can lead to a smaller effective driving force.

Substituting Eqs. (20), (25), and (26) in Eq. (18) and integrating, we get δF kB T (N1,0 + N2,0 ) = [−x1,0 (1 − ξ1 )ln(1 − ξ1 ) − x2,0 (1 − ξ2 )ln(1 − ξ2 ) +(1 − ξ1 x1,0 − ξ2 x2,0 )ln(1 − ξ1 x1,0 − ξ2 x2,0 )]  2 2 x2,0 ξ2 x1,0 (1 − ξ1 )2 x1,0 ξ1 x2,0 A + + kB T 1 − ξ1 x1,0 (1 − ξ1 x1,0 − ξ2 x2,0 )(1 − ξ1 x1,0 )  2 2 − x2,0 ξ2 x1,0 −x1,0 ξ1 x2,0 .

(27)

δF has been non-dimensionalized by kB T times the initial finite amount of the metastable liquid. For an infinite system ξ 1 ,

FIG. 14. Open and closed system free energies obtained by integrating the free energy below n2 = 50 along n2 -axis. Also shown is the closed system free energy with corrections.

084111-10

V. Agarwal and B. Peters

J. Chem. Phys. 140, 084111 (2014) 4 K.

FIG. 15. Schematic showing the eutectic point and the metastable extension of the solidus. The latter identifies the composition of nuclei that are at equilibrium with the metastable fluid.

VI. CONCLUSIONS

We have investigated nucleation at several temperatures and compositions near the eutectic point in a Potts-lattice Gas model. Nucleation rates were computed from MFPT. The MFPT at each condition and for each of the two possible nucleation products was obtained from the free energy as a function of the largest nucleus size and the frequency of attachment/detachment from the critical nucleus. Free energy barriers were computed directly by umbrella sampling in the semigrand canonical ensemble, and also by applying a correction for finite size effects to free energies from umbrella sampling in a closed (canonical) ensemble. Additionally, we tuned the frequency of species identity swaps so that composition fluctuations in the finite periodic cell decay at the rate which is expected for diffusion in and out of the simulation subvolume in an infinitely large system. These new procedures allow us to approximate the thermodynamic and dynamical properties of an open system in an infinite environment using results from simulation in small, closed systems. The corrections should facilitate atomistic simulations of solute precipitate nucleation where open systems and large scale simulations are not feasible. For nucleation near the eutectic point in the Potts-lattice gas model, the temperature and composition dependence of computed nucleation rates followed qualitative trends predicted by the classical nucleation theory. ACKNOWLEDGMENTS

V. Agarwal was supported by the MRSEC of the National Science Foundation under Award No. DMR1121053. B. Peters was supported by the National Science Foundation CAREER Award No. CHE0955502. We thank Dr. Nathan Duff, Dr. Brandon C. Knott, the IRG3 team, and the reviewers for helpful suggestions. We acknowledge support from the Center for Scientific Computing from the CNSI, MRL: an NSF MRSEC (DMR-1121053), NSF CNS-0960316 and Hewlett Packard. 1 R. J. Davey, S. L. M. Schroeder, and J. H. ter Horst, Angew. Chem. Int. Ed.

(English) 52, 2166 (2013). 2 D. Erdemir, A. Y. Lee, and A. S. Myerson, Acc. Chem. Res. 42, 621 (2009). 3 F.

F. Abraham, Homogeneous Nucleation Theory (Academic Press, New York, 1974).

F. Kelton and A. L. Greer, Nucleation in Condensed Matter: Applications in Materials and Biology (Elsevier Ltd., Amsterdam, 2010). 5 D. Kashchiev, Nucleation: Basic Theory with Applications (ButterworthHeinemann, Oxford, MA, 2000). 6 W. Kurz and D. J. Fisher, Fundamentals of Solidification (Trans Tech Publications, 1989). 7 D. A. Porter and K. E. Easterling, Phase Transformations in Metals and Alloys (Stanley Thornes Ltd., 2001). 8 A. Pérez and A. Rubio, J. Chem. Phys. 135, 244505 (2011). 9 D. Moroni, P. ten Wolde, and P. Bolhuis, Phys. Rev. Lett. 94, 235703 (2005). 10 S. Auer and D. Frenkel, Nature 409, 1020 (2001). 11 S. Auer and D. Frenkel, J. Chem. Phys. 120, 3015 (2004). 12 W. Lechner, C. Dellago, and P. G. Bolhuis, Phys. Rev. Lett. 106, 085701 (2011). 13 D. Kashchiev and A. Firoozabadi, J. Cryst. Growth 243, 476 (2002). 14 J. Anwar and D. Zahn, Angew. Chem. Int. Ed. Engl. 50, 1996 (2011). 15 P. R. T. Wolde and D. Frenkel, Science 277, 1975 (1997). 16 E. Clouet, in Fundamentals of Modeling for Metals Processing, edited by D. U. Furrer and S. L. Semiatin (ASM International, 2010), Vol. 22, pp. 203–219. 17 V. Agarwal and B. Peters, “Solute precipitate nucleation: A review of theory and simulation advances,” Adv. Chem. Phys. (in press). 18 R. Demichelis, P. Raiteri, J. D. Gale, D. Quigley, and D. Gebauer, Nat. Commun. 2, 590 (2011). 19 S. Punnathanam and P. A. Monson, J. Chem. Phys. 125, 024508 (2006). 20 R. Cabriolu, D. Kashchiev, and S. Auer, J. Chem. Phys. 133, 225101 (2010). 21 D. Kashchiev and G. M. van Rosmalen, Cryst. Res. Technol. 38, 555 (2003). 22 J. W. Christian, The Theory of Transformation in Metals and Alloys (Parts I and II) (Elsevier Ltd., Oxford, 2002). 23 J. B. Ott, J. R. Goates, and J. D. Lamb, J. Chem. Thermodyn. 4, 123 (1972). 24 S. Kang, R. Rai, and S. Purushothaman, J. Electron. Mater. 25, 1113 (1996). 25 V. Subbotin, M. Arnol’dov, F. Kozlov, and A. Shimkevich, At. Energ. 92, 29 (2002). 26 N. Davies, B. Nicholas, A. Davies, and M. Nicholas, “Patent: Hot melt ink for jet printing – contg. image former and fusible eutectic comps, giving sharp images esp. on plastics,” U.S. patent WO9404618-A1 (1994). 27 G. J. Janz, Handbook of Molten Salts (Academic Press Inc., New York, 1967). 28 M. K. Rathod and J. Banerjee, Renewable Sustainable Energy Rev. 18, 246 (2013). 29 E. Cooper, C. Andrews, and P. Wheatley, Nature 430, 1012 (2004). 30 S. R. Ganagalla and S. N. Punnathanam, J. Chem. Phys. 138, 174503 (2013). 31 R. Becker and W. Döring, Ann. Phys. (Berlin, Ger.) 416, 719 (1935). 32 L. Farkas, Z. Phys. Chem. 175, 236 (1927). 33 J. Frenkel, J. Chem. Phys. 7, 200 (1939). 34 J. Gibbs, Trans. Conn. Acad. 3, 108 (1878). 35 D. Turnbull and J. Fisher, J. Chem. Phys. 17, 71 (1949). 36 M. Volmer and A. Weber, Z. Phys. Chem. 119, 277 (1926). 37 J. Zeldovich, Acta Physicochim. URSS 18, 1 (1943). 38 B. Peters, J. Chem. Phys. 135, 044107 (2011). 39 R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford University Press, New York, 2001). 40 A. Nitzan, Chemical Dynamics in Condensed Phases: Relaxation, Transfer, and Reactions in Condensed Molecular Systems (Oxford University Press, 2006). 41 P. Hänggi, P. Talkner, and M. Borkovec, Rev. Mod. Phys. 62, 251 (1990). 42 J. Wedekind, G. Chkonia, J. Wölk, R. Strey, and D. Reguera, J. Chem. Phys. 131, 114506 (2009). 43 J. Wedekind, R. Strey, and D. Reguera, J. Chem. Phys. 126, 134103 (2007). 44 G. Chkonia, J. Wölk, R. Strey, J. Wedekind, and D. Reguera, J. Chem. Phys. 130, 064505 (2009). 45 D. Reguera, R. K. Bowles, Y. Djikaev, and H. Reiss, J. Chem. Phys. 118, 340 (2003). 46 J. W. Schmelzer and A. S. Abyzov, J. Non-Cryst. Solids 384, 2–7 (2014). 47 R. Grossier and S. Veesler, Cryst. Growth Des. 9, 1917 (2009). 48 J. Wedekind, D. Reguera, and R. Strey, J. Chem. Phys. 125, 214505 (2006). 49 N. Duff and B. Peters, J. Chem. Phys. 131, 184101 (2009). 50 J. Ashkin and E. Teller, Phys. Rev. 64, 178 (1943). 51 D. Bellucci, V. Cannillo, and A. Sola, Ceram. Int. 36, 1983 (2010).

084111-11 52 R.

V. Agarwal and B. Peters

LeSar, Introduction to Computational Materials Science (Cambridge University Press, Cambridge, 2013). 53 M. Saito, and M. Matsumoto, in Monte Carlo and Quasi-Monte Carlo Methods, edited by A. Keller, S. Heinrich, and H. Niederreiter (SpringerVerlag Berlin Heidelberg, 2008), pp. 607–622. 54 K. Kawasaki, Phys. Rev. 145, 224 (1966). 55 N. Metropolis, A. Rosenbluth, M. Rosenbluth, and A. H. Teller, J. Chem. Phys. 21, 1087 (1953). 56 D. A. Kofke and E. D. Glandt, Mol. Phys. 64, 1105 (1988). 57 W. Lechner and C. Dellago, J. Chem. Phys. 129, 114707 (2008). 58 B. C. Knott, N. Duff, M. F. Doherty, and B. Peters, J. Chem. Phys. 131, 224112 (2009). 59 S. Hayward, D. W. Heermann, and K. Binder, J. Stat. Phys. 49, 1053 (1987). 60 J. W. Esam, in Phase Transitions and Critical Phenomena, edited by C. Domb and M. S. Green (Academic Press Inc., New York, 1972), pp. 197–269. 61 P. ten Wolde, M. RuizMontero, and D. Frenkel, J. Chem. Phys. 104, 9932 (1996). 62 S. Jungblut, A. Singraber, and C. Dellago, Mol. Phys. 111, 3527–3533 (2013). 63 B. Peters, Molec. Sim. 36, 1265 (2010). 64 G. Beckham and B. Peters, J. Phys. Chem. Lett. 2, 1133 (2011). 65 L. Maibaum, Phys. Rev. Lett. 101, 019601 (2008). 66 S. L. Girshick and C.-P. Chiu, J. Chem. Phys. 93, 1273 (1990). 67 G. Torrie and J. Valleau, Chem. Phys. Lett. 28, 578 (1974). 68 D. Chandler, Introduction to Modern Statistical Mechanics (Oxford University Press, Inc., 1987), pp. 168–175.

J. Chem. Phys. 140, 084111 (2014) 69 S.

Duane, A. Kennedy, B. Pendleton, and D. Roweth, Phys. Lett. B 195, 216 (1987). 70 D. Frenkel and B. Smit, Understanding Molecular Simulation: From Algorithms to Applications (Academic Press Inc., San Diego, CA, 2002). 71 L. Berezhkovskii and V. Zitserman, J. Chem. Phys. 102, 3331 (1995). 72 A. Berezhkovskii, L. Berezhkovskii, and V. Zitzerman, Chem. Phys. 130, 55 (1989). 73 D. Reguera, J. Rub, and A. Pérez-Madrid, Physica A 259, 10 (1998). 74 D. Reguera, J. M. Rub, and A. Perez-Madrid, J. Chem. Phys. 109, 5987 (1998). 75 M. von Smoluchowski, Ann. Phys. 326, 756 (1906). 76 H. A. Kramers, Physica 7, 284 (1940). 77 J. Langer, Ann. Phys. 54, 258 (1969). 78 B. Peters, J. Chem. Phys. 131, 244103 (2009). 79 W. Im and B. Roux, J. Mol. Biol. 319, 1177 (2002). 80 See supplementary material at http://dx.doi.org/10.1063/1.4865338 for additional details on temperature dependence of the attachment/detachment frequency. 81 P. Yi and G. C. Rutledge, Annu. Rev. Chem. Biomol. Eng. 3, 157 (2012). 82 J. Prausnitz, R. Lichtenthaler, and E. de Azevedo, Molecular Thermodynamics of Fluid-Phase Equilibria, 3rd ed. (Prentice-Hall Inc., New Jersey, 1999). 83 M. Margules, Sitzungsber. Akad. Wiss. Wien, Math.-Naturwiss. Kl., Abt. 2B 104, 1243 (1895). 84 R. Ni, F. Smallenburg, L. Filion, and M. Dijkstra, Mol. Phys. 109, 1213 (2011).

Nucleation near the eutectic point in a Potts-lattice gas model.

We use the Potts-lattice gas model to study nucleation at and near the eutectic composition. We use rare-event methods to compute the free energy land...
2MB Sizes 2 Downloads 3 Views