Monte Carlo Tests of Nucleation Concepts in the Lattice Gas Model
aa r X i v : . [ c ond - m a t . s o f t ] M a y Monte Carlo Tests of Nucleation Concepts in the Lattice Gas Model
Fabian Schmitz, Peter Virnau and Kurt Binder
Institute of Physics, Johannes Gutenberg-Universit¨at Mainz, Germany
The conventional theory of homogeneous and heterogeneous nucleation in a supersaturated vaporis tested by Monte Carlo simulations of the lattice gas (Ising) model with nearest-neighbor attractiveinteractions on the simple cubic lattice. The theory considers the nucleation process as a slow (quasi-static) cluster (droplet) growth over a free energy barrier ∆ F ∗ , constructed in terms of a balance ofsurface and bulk term of a “critical droplet” of radius R ∗ , implying that the rates of droplet growthand shrinking essentially balance each other for droplet radius R = R ∗ . For heterogeneous nucleationat surfaces, the barrier is reduced by a factor depending on the contact angle. Using the definition of“physical” clusters based on the Fortuin-Kasteleyn mapping, the time-dependence of the cluster sizedistribution is studied for “quenching experiments” in the kinetic Ising model, and the cluster size ℓ ∗ where the cluster growth rate changes sign is estimated. These studies of nucleation kinetics arecompared to studies where the relation between cluster size and supersaturation is estimated fromequilibrium simulations of phase coexistence between droplet and vapor in the canonical ensemble.The chemical potential is estimated from a lattice version of the Widom particle insertion method.For large droplets it is shown that the “physical clusters” have a volume consistent with the estimatesfrom the lever rule. “Geometrical clusters” (defined such that each site belonging to the cluster isoccupied and has at least one occupied neighbor site) yield valid results only for temperatures lessthan 60% of the critical temperature, where the cluster shape is non-spherical. We show how thechemical potential can be used to numerically estimate ∆ F ∗ also for non-spherical cluster shapes. I. INTRODUCTION
Since the theory of nucleation phenomena was introduced a long time ago [1–3], the question under which conditionsthe “conventional theory” of nucleation is accurate has been debated (see e.g. [4–24]) and this debate continues untiltoday. For the simplest case of homogeneous nucleation (by statistical fluctuations in the bulk) of a one-componentliquid droplet from the vapor, the basic statement of the theory is that under typical conditions nucleation processesare rare events, where a free energy barrier ∆ F ∗ very much larger than the thermal energy k B T is overcome, andhence the nucleation rate is given by an Arrhenius law, j = ω exp( − ∆ F ∗ /k B T ) . (1)Here j is the number of nuclei, i.e. droplets that have much larger radii R than the critical radius R ∗ associated withthe free energy barrier ∆ F ∗ of the saddle point in configuration space, that are formed per unit volume and unit time; ω is a kinetic prefactor. Now ∆ F ∗ is estimated from the standard assumption that the formation free energy of adroplet of radius R can be written as a sum of a volume term ( ∝ πR / ∝ πR ), i.e.∆ F ( R ) = − πR µ ( ρ ℓ − ρ v ) + 4 πR γ vℓ . (2)Since the liquid droplet can freely exchange particles with the surrounding vapor, it is natural to describe its ther-modynamic potential choosing the chemical potential µ and temperature T as variables, and expand the differencein thermodynamic potentials of liquid and vapor at the coexistence curve, ∆ µ = µ − µ coex , ρ v and ρ ℓ denoting thedensities of the coexisting vapor ( v ) and liquid ( ℓ ) phases. According to the capillarity approximation, the curvaturedependence of the interfacial tension γ vℓ is neglected, γ vℓ is taken for a macroscopic and flat vapor-liquid interface.Then the critical radius R ∗ follows from ∂ ∆ F ( R ) ∂R (cid:12)(cid:12)(cid:12)(cid:12) R ∗ = 0 , R ∗ = 2 γ vℓ ∆ µ ( ρ ℓ − ρ v ) , (3)and the associated free energy barrier is ∆ F ∗ hom = 4 π R ∗ ) γ vℓ . (4)However, since typically ∆ F ∗ hom is less than 100 k B T , the critical droplet is a nanoscale object, and thus the treatmentEqs. (1)-(4) is questionable. Experiments (e.g. [25–27]) were not able to yield clear-cut results on the validity ofEqs. (1)-(4), and how to improve this simple approach: critical droplets are rare phenomena, typically one observesonly the combined effect of nucleation and growth; also the results are often “contaminated” by heterogeneousnucleation events due to ions, dust, etc. [28–32], and since j varies rapidly with the supersaturation, only a smallwindow of parameters is suitable for investigation. Therefore this problem has been very attractive, in principle, forthe study via computer simulation. However, despite numerous attempts (e.g. [8–10, 16, 21–24, 33]), this approach isalso hampered by two principal difficulties:(i) Computer simulations can often only study a small number of decades in time, [9, 33], which in typical casescorrespond to small barriers (∆ F ∗ ≤ k B T ) rather than the larger ones which are of more interest in thecontext of experiments.(ii) On the atomistic scale, it is a difficult and not generally solved problem to decide which particles belong toa droplet and which particles belong to its environment; the vapor-liquid interface is diffuse and fluctuating[34, 35].For these reasons, many of the available simulation studies have addressed nucleation in the simplistic Ising (latticegas) model, [9, 10, 18, 22, 33, 35–56], first of all since it can be very efficiently simulated, and secondly because one candefine more precisely what is meant by a “cluster”. Associating Ising spins σ i = +1 at a lattice site i with a particle, σ i = − pp ( T ) = 1 − exp( − J/k B T ) , (5) J being the Ising model exchange constant.Due to Eq. (5), the “physical clusters” defined in this way are typically smaller than the geometrical clusters, andtheir percolation point can be shown to coincide with the critical point [59–61]. A geometrical cluster hence cancontain several physical clusters. Note that to apply Eq. (5), random numbers are used, and hence physical clustersare not deterministically defined from the spin configuration, but rather have some stochastic character. This presentsa slight difficulty in using physical clusters in the study of cluster dynamics.While Eq. (5) has been used in the context of simulations of critical phenomena in the Ising model, applying veryefficient Swendsen-Wang [60] and Wolff [68] simulation algorithms, this result has almost always been ignored in thecontext of simulations of nucleation phenomena [18, 22, 47–52]. While it is allright to ignore the difference betweengeometrical and physical clusters in the limit T → p ( T ) → H acts [55, 56]. First of all, in this way also a systematic investigation of heterogeneous nucleation atplanar walls becomes feasible; secondly, due to the reduction of the barrier ∆ F ∗ het in comparison to ∆ F ∗ hom ; nucleationfor reasonably large values of R ∗ becomes accessible to study.In Sec. II, we consider the equilibrium of the lattice gas model for ρ v < ρ < ρ ℓ in systems in a L × L × L geometrywith periodic boundary conditions, to show that physical clusters do occupy precisely the volume predicted by thelever rule analysis, [21, 23, 24, 55] as they should when the thermodynamic limit is approached. We present evidencethat physical clusters are correctly identified by both the lever rule method and the approach based on the “atomistic”identification of clusters based on Eq. (5) at all temperatures, from zero temperature up to the critical temperature T c . In contrast, Eq.(3), which implies a spherical droplet shape, is found to work only at temperatures distinctlyabove the interface roughening transition temperature T R [71, 72], even for very large radii R . We attribute thesediscrepancies to the fact that due to the anisotropy of the interface tension for our lattice model pronounced deviationsof the average droplet shape from a sphere occur [73–76], presenting data on the shape of large droplets. In Sec. III,we describe our results on the dynamics of the droplet size distribution and on the attempt to find R ∗ from the sizewhere growth and shrinking processes of clusters are balanced. This study is also carried out for systems with a freesurface, for which the contact angles for various values of the surface field have been estimated previously [55, 56],since in this case much lower barriers (for large clusters) result, which is crucial for making this study feasible withmanageable effort. However, the radii R ∗ predicted from this analysis of kinetics show slight deviations from the radii R ∗ predicted from ∆ F ∗ het . Possible reasons for this discrepancy will be discussed. Finally, Sec. IV summarizes ourconclusions. II. MICROSCOPICALLY DEFINED CLUSTERS VERSUS MACROSCOPIC DOMAINS IN THERMALEQUILIBRIUM
As is obvious from Eq. (2) and the reasoning behind it, this approach is adequate when one deals with the descriptionof macroscopically large domains in equilibrium with a surrounding bulk phase. However, one needs to find anextension of the concept that can be applied also to nanoscopically small droplets, “clusters” in the lattice gas modelthat contain perhaps only of the order of 100 fluid particles. In this section, we want to confirm the idea that one mustuse the concept of “physical clusters” based on Eq. (5) for this purpose, rather than the “geometrical clusters” that areso widely used when the lattice gas model is used to test nucleation theory concepts. While the geometrical clustersare appropriate if one works at extremely low temperatures where the clusters basically have the shape of small cubes[52], this region clearly is inappropriate when one has the application for vapor-to-liquid nucleation in mind, wheredroplets are spherical, and their interfaces are rough and fluctuating rather than smooth planar facets. In fact, manystudies of nucleation in the lattice gas model have been made in d = 3 dimensions at temperatures near T /T c = 0 . T R /T c ≈ . k B T c /J = 4 . k B T R /J ≈ .
44, it is clear thattemperatures much closer to T c must be studied to render the assumption of a spherical droplet shape accurate. Infact, this assumption of a spherical droplet shape is accurate when the difference between the interfacial stiffness [80]and the interfacial free energy becomes negligibly small. Numerical studies of Hasenbusch and Pinn [81] indicate thatthis is only the case for k B T /J ≥ .
9. As a consequence, it is clear that most of the existing studies of nucleationphenomena in the Ising model, that were based on geometrically defined clusters, and had to be done at much lowertemperatures, are inconclusive: the deviation of the average droplet shape from a sphere enhances the surface termin Eq. (2); but the fluctuation corrections discovered for small droplets by the “lever rule method” [21, 23, 24, 55]show that γ vℓ ( R ) < γ vℓ ( ∞ ) for small R and hence the surface term in Eq. (2) is decreased. Thus, it hardly can be asurprise that some of the studies concluded that the nucleation barriers predicted by classical nucleation theory andthe capillarity approximation (that ignores the R -dependence of γ vℓ ( R )) are too high, and others concluded they aretoo low, or even reported good agreement. We take the latter finding as indication that the two opposing effects haveaccidentally more or less canceled each other.This problem is the motivation for the present section, which attempts to show that “physical clusters” based onEq. (5) are appropriate to identify clusters in the lattice gas model, irrespective of temperature and cluster size, andare equivalent to the droplets of the “lever rule method”, for large enough droplets.Fig. 1 recalls this approach: one samples for a system of volume V = L × L × L the effective thermodynamicpotential f L ( T, ρ ) per lattice site as a function of density ρ , f L ( T, ρ ) = [ F ( N, V, T ) − F ( N = V ρ ℓ , V, T )] /V , (6)where N is the number of occupied lattice sites ( ρ i = (1 + σ i ) / σ i = ± i ). Since phase coexistence between bulk liquid (at density ρ ℓ ) and vapor (at density ρ v ) occurs at a chemicalpotential µ coex that corresponds to the “field” (in magnetic notation) H = 0, ρ ℓ and ρ v are simply related to thespontaneous magnetization of the Ising ferromagnet as ρ ℓ = (1 + m sp ) / ρ v = (1 − m sp ) /
2, and H translates into µ via H = ( µ − µ coex ) / µ/
2. The accurate sampling of f L ( T, ρ ) for large L is a nontrivial task, it requires the useof advanced methods such as “multicanonical Monte Carlo” [93] or “Wang Landau sampling” [94, 95] or “successiveumbrella sampling” [90, 91], see [83, 84] for background on such techniques. From Eq. (6), one defines a chemicalpotential function as a derivative, ˜ µ ( N, V, T ) = ∂F ( N, V, T ) ∂N (cid:12)(cid:12)(cid:12)(cid:12) V,T , (7)∆ µ L ( T, ρ ) = ˜ µ ( N, V, T ) − µ coex . (8)One recognizes that the isotherms ∆ µ L ( T, ρ ) vs. ρ exhibit a loop: the homogeneous vapor remains stable also forsome region where µ > µ coex , until a peak occurs, which indicates the “droplet evaporation/condensation transition”[21, 85, 92]: in the first regime where ∆ µ L ( T, ρ ) decreases with ρ , a (more or less spherical or cubical) droplet coexistswith surrounding vapor. Here, we are not interested in the further transitions that one can recognize from this curve,where the droplet changes shape from spherical to cylindrical, or to a slab configuration, etc. [21, 23, 24]. Instead, we -0.4-0.3-0.2-0.100.10.20.30.4 c h e m i ca l po t e n ti a l ∆ µ L ( T , ρ ) / k B T ρ fr ee e n e r gy p e r p a r ti c l e f L ( T , ρ ) / k B T ρρ ’ ρ ’’ FIG. 1: Illustration of the numerical procedure for determining the density triplets ρ ′ , ρ and ρ ′′ and the associated values f L ( T, ρ ′ ), f L ( T, ρ ) and f L ( T, ρ ′′ ) for a given choice of T and L , according to the “lever rule method”. The density ρ must bechosen such that it is not too close to the peak of the curve ∆ µ L /k B T versus ρ , which is due to the “droplet evaporation-condensation transition”, because all microstates that are sampled over should contain a droplet in the system. It must alsobe chosen not too close to the first kink in the curve ∆ µ L /k B T versus ρ , which is due to the transition of the droplet to acylindrical shape (stabilized by the periodic boundary condition in the direction of the cylinders). Apart from these constraints,the choice of ρ is arbitrary, and studying different choices of ρ provides useful consistency checks on the results. The densities ρ ′ and ρ ′′ then can be read off from the curve ∆ µ L /k B T versus ρ , since it is required that ∆ µ L ( T, ρ ′ ) = ∆ µ L ( T, ρ ) = ∆ µ L ( T, ρ ′′ ).The corresponding values of the free energy densities then can be read off from the lower part of the figure, and using Eqs. (9),(10) with N exc = 0, the surface free energy F S is extracted. The actual data shown here refer to the case k B T /J = 3 . L = 20. emphasize the key idea of the “lever rule method”: one can identify a range of choices for the chemical potential µ where three states of the finite system can exist in equilibrium with the same chemical potential, namely a homogeneousvapor at density ρ ′ > ρ v , a homogeneous liquid at density ρ ′′ > ρ ℓ , and a state where two-phase coexistence betweenthe droplet and surrounding vapor occurs. Since the vapor in this case exists at the same chemical potential as thepure vapor, it must be of the same physical nature as the state with density ρ ′ , and similarly, the liquid in the dropletcan be identified with the liquid at ρ ′′ . Making now use of the fact that for large enough systems a system can besuitably decomposed into independent subsystems, we write for the free energy, with V ′′ the volume taken by thedroplet, V ′ = V − V ′′ , V f L ( T, ρ ) = V ′ f L ( T, ρ ′ ) + V ′′ f L ( T, ρ ′′ ) + F S , (9)where the free energy densities f L ( T, ρ ′ ) and f L ( T, ρ ′′ ) are explicitly known, and also f L ( T, ρ ) is known: thus, whenthe droplet volume V ′′ is known, the surface free energy F S of the droplet, which is defined via Eq. (9), is determined.A similar decomposition can readily be written down for the particle number, N = V ρ = V ′ ρ ′ + V ′′ ρ ′′ + N exc , (10)where we have allowed for an excess number N exc of particles, to be associated with the interface. If we consider thedefinition of an “equimolar dividing surface” [34], N exc = 0, and then reading off ρ ′ and ρ ′′ from the construction inFig. 1 we see that Eq. (10) readily yields V ′ and V ′′ for the considered density ρ , and via Eq. (9) we can immediatelyextract F S from the data. Note that these arguments do not invoke the assumption that the dividing surface needsto be a sphere. If one makes the assumption, V ′′ = 4 πR /
3, and then one can write also F S = 4 πR γ vℓ ( R ). Thus,it is assumed that all interactions of particles inside the droplet (volume region V ′′ ) with particles inside the vapor(volume region V ′ ) are restricted to the interfacial region, and hence can be accounted for by their contribution tothe surface free energy F S . (a) den s i t y ρ ′ chemical potential ∆µ /k B Tk B T/J=3.0, gcank B T/J=4.0, gcank B T/J=3.0, cank B T/J=4.0, can (b) den s i t y ρ ′′ chemical potential ∆µ /k B Tk B T/J=3.0, gcank B T/J=4.0, gcank B T/J=3.0, cank B T/J=4.0, can
FIG. 2: (Color online) (a) Plot of ρ ′ ( T, ∆ µ ) versus ∆ µ/k B T , for several temperatures k B T /J as indicated. The symbolsrepresent data recorded in the grand-canonical ensemble of the lattice gas, while the curves are obtained in the canonicalensemble, recording ∆ µ = ∆ µ ( T, ρ ) via the lattice version of the Widom particle insertion method [55, 87, 88]. The chosenlattice size was L = 60. (b) Same as (a), but for ρ ′′ ( T, ∆ µ ). Both methods agree perfectly. However, the method defined via Eqs. (9), (10), and illustrated in Fig. 1 becomes difficult to apply for very largedroplets, because the sampling of f L ( T, ρ ) then becomes unreliable or would require an unaffordable effort. Themethod is also difficult to apply for small droplets, because then one must use relatively small simulation boxes toensure the stability of the inhomogeneous state where a droplet coexists with surrounding vapor [37]. Thus, it isvery desirable to complement the approach by a more “microscopic” identification of droplets, and this is possiblevia Eq. (5). In particular, it has been shown that apart from finite size effects (see e.g. [61]) that for L → ∞ the spontaneous magnetization of the Ising ferromagnet m sp coincides with the percolation probability P , whichis defined [86] as the fraction of sites belonging to the largest “physical cluster” in the system. When we henceanalyze a configuration at a density ρ , where (cf. Fig. 1) a large cluster is present in the system, using Eq. (5) todefine clusters the largest cluster will include ℓ sites, which we hence can associate with its (total) magnetization M = m ′′ V ′′ where m ′′ is then the magnetization per site (note that the Ising magnet/lattice gas isomorphism impliesthat ρ ′′ = (1 + m ′′ ) / V ′′ from a “measurement” of the averagesize h ℓ i of the largest cluster in the system via V ′′ = h ℓ i m ′′ . (11)Note that this volume in general differs from the volume of a geometrical cluster: If h ℓ geom i counts all occupied sitesbelonging to the geometrical cluster, and noting that the density in a large geometrical cluster is just the bulk density,namely ρ ′′ = (1 + m ′′ ) /
2, the volume taken by the geometrical cluster is given by V ′′ geom = h ℓ geom i ρ ′′ = 2 h ℓ geom i m ′′ . (12)In the limit L → ∞ , where also h ℓ i and V ′′ get macroscopically large, m ′′ tends to m sp , while for finite L it is clearthat m ′′ slightly exceeds m sp (and ρ ′′ exceeds ρ ℓ , see Fig. 1). However, recording the relation ρ = ρ ( T, ∆ µ ) [or theequivalent relation m = m ( T, H ) of the Ising ferromagnet] very precisely is an easy task, see Fig. 2, since both statesat densities ρ ′ , ρ ′′ (Fig. 1) are homogeneous, not affected by heterophase fluctuations, and since the temperaturesstudied are still well below T c , statistical fluctuations are small, and finite size effects are negligible. Fig. 2 presentsrepresentative results for both ρ ′ ( T, ∆ µ ) and ρ ′′ ( T, ∆ µ ) versus ∆ µ . Note we also have used the lattice version of theWidom particle insertion method [55, 87, 88] to record the inverse function ∆ µ = ∆ µ ( T, ρ ) from simulations in thecanonical ensemble, where ρ was chosen as the independent control variable. The perfect agreement between bothapproaches not only serves as a test of the accuracy and correctness of our numerical procedures, but also showsthat for the chosen temperatures and linear dimensions finite size effects on states in “pure” phases are completelynegligible, since finite size effects are known to differ in the two ensembles [83, 84], but are not detected here at all. (a) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length L ρ =0.08 ρ =0.07 ρ =0.06 ρ =0.05 ρ =0.04 ρ =0.03 (b) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length L ρ =0.12 ρ =0.11 ρ =0.10 ρ =0.09 ρ =0.08 ρ =0.07 ρ =0.06 ρ =0.05 (c) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length L ρ =0.12 ρ =0.11 ρ =0.10 ρ =0.09 ρ =0.08 ρ =0.07 (d) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length L ρ =0.20 ρ =0.19 ρ =0.18 ρ =0.17 FIG. 3: (Color online) Plot of the ratio V ′′ /V ′′ lr of the droplet volume as obtained from the “Coniglio-Klein-Swendsen-Wang”[59, 60] cluster definition, using Eq. (5) to define physical clusters, and their volume V ′′ (Eq. (11)), and from the lever rule(Eq. (13)), as a function of 1 /L , for four temperatures: k B T /J = 2 . . . . To record V ′′ as defined in Eq. (11) from the simulations, we performed simulations in the canonical ensemble,choosing various values of ρ , and equilibrate a large droplet coexisting with surrounding vapor. The initial stateis then chosen putting a droplet with the size predicted by the lever rule (and density ρ = (1 + m sp ) /
2) intothe simulation box, which then is carefully equilibrated. The standard method to simulate the Ising model withconserved magnetization (which corresponds to the lattice gas model in the canonical
N V T ensemble) is the “spinexchange algorithm” [83, 84]. However, the standard nearest neighbor exchange method implies that any local excessof magnetization (or density, respectively) can only relax diffusively, and the resulting “hydrodynamic slowing down”[83, 84] hampers the fast approach towards thermal equilibrium. We thus used instead a single spin-flip algorithm inwhich the total magnetization is restricted to two neighboring values { M − , M } : So if the system has magnetization M , flipping a down-spin (which would mean a transition M → M + 2) is automatically rejected, and if the system isin the state M −
2, flipping of up-spins is forbidden. For large systems ( L → ∞ ) any corrections to a strictly canonicalsimulation at a magnetization per spin m = M/L are of order 1 /L and hence negligible.Choosing very large systems (up to L = 160, rather than L = 20 as used in Fig. 1) the ratio of V ′′ as found fromEq. (11) and V ′′ lr from the lever rule (with the assumption N exc ≡ V ′′ lr = V ρ − ρ ′ ρ ′′ − ρ ′ , (13)is plotted vs. 1 /L in Fig. 3 for several temperatures and various densities ρ . These data show that V ′′ /V ′′ lr → L → ∞ , irrespective of temperature and density (in the density region where a droplet not affected by the periodicboundary conditions is present, as explained in Fig. 1). The fact that V ′′ /V ′′ lr extrapolates to unity for L → ∞ notprecisely, but only within some error, is due to the fact that statistical errors affect both the estimation of h ℓ i and of V ′′ lr (via errors in the estimation of ∆ µ and hence ρ ′ ). As expected, using the proper definition of physical clusters,one can work at arbitrary temperatures, both at T < T R (a), T slightly above T R (b) or T rather close to T c (d).While in cases (a) the simple geometrical cluster definition would also work, since essentially all bonds are “active”( p ( T ) is almost unity { Eq. (5) } ), and the non-spherical shape of the clusters does not matter in this context. But thegeometrical cluster definition would clearly break down in case (d) due to the proximity of the percolation transitionthat occurs for geometrical clusters at a density not much larger than those included in Fig. 3(d) [57, 64, 65]. Onthe other hand, we note from the fact that there always occurs asymptotically in the relation V ′′ /V ′′ lr a correction oforder 1 /L , that for physical clusters the assumption N exc = 0, that is often (but not always [24]) made in the leverrule method, does not hold: i.e., when we assume that N exc is proportional to the surface area of the droplet, we canwrite N exc = C ( V ′′ ) / ρ exc , (14)where C is a geometrical factor ( C = (36 π ) / for a spherical droplet, C = 6 for a cube), and ρ exc is an excess densityof the particles due to the interface of the droplet. Using V ′′ lr from Eq. (13) as a first-order estimate in Eq. (14), onereadily finds that the term N exc /V in Eq. (10) yields a 1 /L correction, N exc V = CL (cid:18) ρ − ρ ′ ρ ′′ − ρ ′ (cid:19) / ρ exc , (15)and hence Eq. (10) would yield instead of Eq. (13) V ′′ lr (corrected) V ′′ lr = 1 − CL (cid:18) ρ − ρ ′ ρ ′′ − ρ ′ (cid:19) / ρ exc ρ − ρ ′ , L → ∞ = 1 − Cρ exc L ( ρ ′′ − ρ ′ ) − / ( ρ − ρ ′ ) − / (16)which is qualitatively in accord with Fig. 3. In order to test Eq. (16) and obtain estimates for the temperaturedependence of ρ exc , Fig. 4 plots our numerical results for the ratios V ′′ /V ′′ lr versus ( ρ ′′ − ρ ′ ) − / ( ρ − ρ ′ ) − / /L . We seea very good data collapse at straight lines going through unity at the ordinate within numerical error at all studiedtemperatures. (a) d r op l e t v o l u m e r a t i o V ’’ / V l r ’’ ( ρ ’’ - ρ ’) -2/3 ( ρ - ρ ’) -1/3 /Lk B T/J=2.0k B T/J=2.6k B T/J=3.0 (b) d r op l e t v o l u m e r a t i o V ’’ / V l r ’’ ( ρ ’’ - ρ ’) -2/3 ( ρ - ρ ’) -1/3 /Lk B T/J=3.0k B T/J=4.0k B T/J=4.3
FIG. 4: (Color online) Plot of V ′′ /V ′′ lr versus ( ρ ′′ − ρ ′ ) − / ( ρ − ρ ′ ) − / /L at low temperatures (a) and at temperatures closer to T c (b), including all densities ρ that were analyzed. Straight line fits to the data are included. Note that the different symbolscorrespond to different choices of ρ in Fig. 3. Since the data in Fig. 3 suggest that for “physical droplets” in the lattice gas model an appreciable “interfacialadsorption” (as expressed in ρ exc in Eq. (14)) occurs, it is of interest to not only study the average volume of the (a) den s i t y ρ Radius R L=60L=80L=100L=130L=160R lr R’’R
Geom (b) den s i t y ρ Radius R L=60L=80L=100L=130L=160R lr R’’R
Geom (c) den s i t y ρ Radius R L=60L=80L=100L=130L=160R lr R’’R
Geom (d) a v e r age s qua r ed i n t e r f a c i a l w i d t h < w > logarithm of droplet radius log(R)T=3T=4T=4.3 FIG. 5: (Color online) Radial density profiles ρ ( R ) plotted vs. R at k B T /J = 3 . T = 4 . k B T /J = 4 . L × L × L with L = 60 , , , ,
160 (from left to right). The vapor density and the liquid density are independentof L , as expected. Vertical straight lines indicate the radii that follow (assuming R = (3 V ′′ / π ) / , i.e. a spherical shape) fromEq. (11) with V ′′ = h ℓ i /m ′′ or from the simple geometrical cluster definition, Eq. (12), where V ′′ geom = 2 h ℓ geom i / (1 + m ′′ ) orfrom the lever rule Eq. (13). A horizontal line at ρ = 0 . R → ∞ the density profile should becomesymmetric with respect to this line, due to the spin reversal symmetry of the Ising model, and then the correct cluster definitionmust yield a cluster radius compatible with this inflection point of the profile. All data are obtained from averages over severalhundred independent droplet configurations. The dotted vertical lines are the radii predicted from the lever rule. Part (d) givesa plot of the average squared interfacial width h w i vs. ln R (in this plot, R is taken from the condition ρ ( R ) = 0 . droplets h ℓ i but also their radial density profile (Fig. 5). It is seen that for radii R which are in the range from 10 to20 lattice constants (corresponding to droplet volumes in the range from about 4500 to 36000, so these are alreadyclusters of a mesoscopic size, with a huge nucleation barrier far beyond observation in a simulation of nucleationevents or in experiments) the profiles are very broad, and their width increases slightly with increasing droplet size.For comparison, the prediction for the cluster radius resulting from the widely used standard geometrical definition ofclusters is also included: it happens that this geometrical radius is still fairly close to the correct radius at k B T /J = 3 . T c the geometrical cluster definition yields completely unreasonable results, due to the onset ofpercolation phenomena, and for k B T /J = 4 . .
3, the geometrical radii are indeed unreasonably large.It is interesting to note (Fig. 5d) that the width of the interfacial profile increases with increasing R . This phe-nomenon is well-known for planar interfaces and attributed to capillary waves. Since for a large droplet the surfaceis locally planar, most of the capillary wave spectrum is not affected by the interface curvature. Thus we may, as afirst approximation, take over the result for the broadening of a planar interface of linear dimension L , replacing L by the droplet radius R [98–100] h w i = w + 14 e γ ln (cid:18) Rλ min (cid:19) (17)where w is the “intrinsic width” of the interfacial profile, e γ the “interfacial stiffness” (note that a factor 1 /k B T is absorbed in its definition) and λ min a short wavelength cutoff, whose precise value is not known. Near T c theinterfacial stiffness coincides with the interfacial free energy, while e γ → ∞ at T → T R , and then the capillary wavebroadening disappears. If one accepts the above formula, and uses the data of Hasenbusch and Pinn [81], one predictsfor the slope [4 e γ ] − ≈ . k B T /J = 4 .
0) or 7 .
98 ( k B T /J = 4 . R term are actually somewhat smaller, namely about 2.13 at k B T /J = 4 . .
83 at k B T /J = 4 .
3, but of thesame order of magnitude.We deliberately do not discuss the radial droplet density profiles for k B T /J = 2 . k B T /J = 2 .
6, however, sinceat these temperatures the droplet shape shows distinct deviations from the spherical shape. This can be checkeddirectly by recording contours of constant density in slices of width = 1 (taking into account 3 lattice planes throughthe droplet’s center of mass, parallel to the xy plane, the xz plane and the yz plane, respectively). Averaging thesedensity contours over several hundred statistically independent observations the plots shown in Fig. 6 are obtained.They are all taken at the same fixed density ρ = 0 .
33. Note that for this density, the stable state would be a cylindricaldroplet (stabilized by the periodic boundary conditions) but for such large systems ( L = 100) the compact dropletshapes shown here (chosen by an appropriate initial condition) are always perfectly metastable. At k B T /J = 2 . x = 0 or y = 0 is not evident, due to finite-size effects at the points where the facets join the roundsections, which replace the sharp edges of the cubes at nonzero temperature). At k B T /J = 2 .
6, where we exceed theroughening temperature slightly, there clearly occur no longer any facets, but the density contours in Fig. 6 are stilldistinctly non-circular: the diameter in diagonal direction clearly is about 7% larger than in the lattice directions.Even at k B T /J = 3 .
0, we find an enhancement of the diameter in diagonal direction of about 3%. At k B T /J = 4 . k B T /J = 4 . (a)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (b)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (c)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (d)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (e)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 (f)
10 20 30 40 50 60 70 80 90 100 10 20 30 40 50 60 70 80 90 100 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
FIG. 6: (Color online) Contours of constant density ρ ( x, y ) = ρ i in the planes x = 0, y = 0 and z = 0 cutting throughdroplets situated at the origin in a box of dimension L = 100 with periodic boundary conditions in every direction, for thecases k B T /J = 1 . k B T /J = 2 . k B T /J = 2 . k B T /J = 3 . k B T /J = 4 . k B T /J = 4 . ρ = 0 .
33. In addition to the color-coded average density per site, the contours for three densities ρ i = ( ρ v ( T ) + 0 . / , . , ( ρ ℓ ( T ) + 0 . / (a) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length LV ′′ V Geom V CNT (b) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length LV ′′ V Geom V CNT (c) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length LV ′′ V Geom V CNT (d) ∞
200 130 100 80 60 d r op l e t v o l u m e r a t i o V / V l r ′′ length LV ′′ V Geom V CNT
FIG. 7: (Color online) Droplet volume V ′′ relative to the lever rule estimate V ′′ lr plotted versus 1 /L for four temperatures k B T /J = 2 . . . . V ′′ geom = 2 h ℓ geom i / (1 + m ′′ ), the “physical cluster” definition V ′′ = h ℓ i /m ′′ { Eq. (11) } andthe result of classical nucleation theory V ′′ = 4 π ( R ∗ ) / R ∗ is computed from (3), using h ∆ µ i as “measured” in thesimulation by the lattice version of the Widom particle insertion method. At each temperature, several densities were used: ρ = 0 . , . , . . . , .
12 (a), 0 . , . , . . . , .
12 (b), 0 . , . , . . . , .
20 (c), 0 . , . , .
28 (d).
As a final part of our analysis of static properties of physical droplets in the Ising model, we exploit the fact that thechemical potential µ can be “measured” by the lattice version of the Widom particle insertion method [55, 87, 88] alsoin the state when the system is inhomogeneous, e.g. for the case of interest when a droplet coexists with surroundingvapor. Actually, already in [55] it was shown that µ actually stays spatially constant in such a situation. However,while the estimation of µ from Eq. (7) requires a very accurate estimation of the free energy density f L ( ρ, T ) { Eq. (6),Fig. 1 } , and in practice this works only for not so large linear dimension L (such as L = 20 in Fig. 1), the estimationof µ from the particle insertion method still works for volumes that are orders of magnitude larger. As a consequence,we can use Eq. (3) to test whether or not the actual droplet volumes V ′′ are compatible with conventional nucleationtheory for large droplets (assuming that the droplets are spherical) so R ∗ = (3 V ′′ / π ) / holds for a critical droplet.Fig. 7 presents a plot of the cluster volume V ′′ versus inverse linear dimension for several densities, at the temper-atures k B T /J = 2 .
6, 3 .
0, 4 . .
3, using the classical nucleation theory prediction V ′′ CNT = 4 π ( R ∗ ) / R ∗ given by Eq. (3) for comparison. For this purposes γ vℓ is taken from the work of Hasenbusch and Pinn [81, 82], andso there occur no unknown parameters whatsoever. The geometrical cluster volume V ′′ geom { Eq. (12) } is close to theestimate based on the Coniglio-Klein-Swendsen-Wang “physical cluster”-definition, Eq. (11), at k B T /J = 2 . k B T /J = 3 .
0, while at k B T /J = 4 . V ′′ geom then is systematically toohigh (in comparison with all other estimates, including V ′′ lr { Eq. (13) } , which is used as a convenient normalization).It is interesting to observe that the classical nucleation theory estimates based on Eq. (3) are systematically too lowfor T = 2 . T = 3 .
0, while for T ≥ . V ′′ we imply that2the cluster volume is spherical and hence we underestimate the surface area: the geometrical factor C ( T ) introducedin Eq. (14) increases from about 4.836 for the spherical shape near T c up to 6.0 for the cube, as the temperature islowered. Since the non-spherical droplet shapes at k B T /J = 2 . µ then yields atoo low radius R ∗ . The observation that the anisotropy of surface tension in the lattice gas model becomes noticeablefor k B T /J < . V ′′ /V ′′ lr = 1 for L → ∞ , since in this limit the lever rule { Eq. (13) } is trivially true with ρ ′ = ρ v and ρ ′′ = ρ ℓ . Themethod based on the definition of physical clusters, Eq. (11), indeed is nicely compatible with this expectation atall temperatures; although it is somewhat unsatisfactory (and unexpected) that at finite L there occurs a surfacecorrection due to the surface excess N exc noted in Eq. (14) (and discussed above). However, it is clear that the twoother methods do not give results that are correct for L → ∞ in general: while the method based on the geometriccluster definition still gives essentially correct results at k B T /J = 2 . T > T R , the volume of geometrical clusters issystematically too large. At k B T /J = 4 .
3, the error is as large as 60 to 80% even asymptotically, and for the clustersizes that were actually studied the overestimation actually is by a factor two to five (Fig. 7d)! In view of the fact thatthe geometric cluster definition must break down due to the percolation transition [57], this failure is not unexpected,but we are not aware that it ever has been quantified previously. In the regime from k B T /J = 2 . .
0, the error ofthe geometric cluster volume raises from a few percent to 15 to 40%.The results obtained from the classical nucleation theory via the “measurement” of the supersaturation ∆ µ , on theother hand, yield essentially the correct result at high temperatures ( k B T /J ≥ . L , and hence R ∗ , in the considered range. But it is remarkable thatfor k B T /J = 3 . V ′′ CNT /V ′′ lr ≈ .
94 (Fig. 7b) and for k B T /J = 2 . V ′′ CNT /V ′′ lr ≈ .
89 (Fig. 7a). Thisdiscrepancy becomes worse at lower temperatures (e.g. V ′′ CNT /V ′′ lr ≈ .
68 at k B T /J = 2 . V ′′ CNT /V ′′ lr is independent of R ∗ seems to be at variance with the finding of a curvature-dependent surface tension γ vℓ ( R ) due to Winter et al. [21, 23, 24, 55, 56, 96]. In fact, evidence was provided that γ vℓ ( R ) = γ vℓ ( ∞ )1 + 2( l/R ) (18)where l is a length proportional to the correlation length in the bulk [96]. Note that due to the spin reversal symmetryof the Ising model one can show [97] that a Tolman correction ( ∝ /R ) must be absent for R → ∞ . However, fromEqs. (2) and (18), it is straightforward to show that R ∗ ∆ µ ( ρ ℓ − ρ v ) = 2 γ ( ∞ ) 1 + 4( l/R ∗ ) (1 + 2( l/R ∗ ) ) = 2 γ ( ∞ ) (cid:2) − l/R ∗ ) + O (( l/R ∗ ) ) (cid:3) . (19)Hence for R ∗ ≫ l it is clear that the result for R ∗ is still given by Eq. (3), and for k B T /J ≤ . V ′′ CNT /V ′′ lr that are seen in Fig. 7 can be attributed fully to the deviation of the droplet shape from a perfect sphere,caused by the anisotropy of the interfacial free energy. In Fig. 8, we now present the ratio of the intercepts V ′′ lr /V ′′ CNT for R ∗ → ∞ as a function of temperature, since we know that { Eq. (3) } V ∗ CNT = (4 π/ γ vℓ ( ∞ ) / [∆ µ ( ρ ℓ − ρ v )] ,while Eq. (4) yielded ∆ F ∗ hom /k B T = (36 π ) / γ vℓ ( ∞ ) / V ∗ CNT ) / = m coex HV ∗ CNT for a spherical droplet. At T = 0,however, the droplet is a perfect cube, and for intermediate temperatures, its shape (for V ∗ → ∞ ) is given by theWulff construction [101] (and hence not explicitly known). However, for large droplet volume V we can write ingeneral ∆ F ( V ) = − m coex HV + γ e A e V − / V / , V → ∞ (20)when we have assumed that the droplets of different linear dimension R (for R → ∞ ) have the same shape at fixedtemperature, so we can write V = e V R for the droplet volume ( e V is then formally the volume for R = 1) and thesurface area is A = e AR . For instance, for a sphere we have e V sphere = 4 π/ e A sphere = 4 π , and for the cube e V cube = 1 and e A = 6. Minimizing ∆ F ( V ) with respect to V yields( V ∗ ) / = 13 γm coex H e A e V − / (21)3for a general shape, which is in between sphere and cube. The barrier then can be written as∆ F ∗ = 127 γ ( m coex H ) e A e V − = m coex HV ∗ . (22)Using now the fact that by choosing a particular large droplet volume in our simulation, H is automatically fixed forany volume V ∗ lr due to the thermal equilibrium situation constructed in our simulation. So it makes sense to estimatethe ratio of barriers as ∆ F ∗ ∆ F ∗ CNT = V ∗ lr V ∗ CNT (23)For T = 0, we know that γ is again the interface tension of the planar surface, also used in the classical nucleation theoryfor the spherical surface. Hence when we write the ratio of ∆ F ∗ / ∆ F ∗ CNT , using Eq. (22), the term γ / (27( m coex H ) )cancels, ∆ F ∗ T =0 ∆ F ∗ CNT = e A e V − e A e V − = 6 π . (24)This asymptotic value of V ∗ lr V ∗ CNT should be reached for T →
0, while V ∗ lr V ∗ CNT → T → T c . Our numerical results(Fig. 8) are compatible with this expectation. π v o l u m e r a t i o V l r ′′ / V CN T ′′ temperature k B T/J L=60L=80L=100T R T C FIG. 8: (Color online) Temperature variation of V ′′ lr /V ′′ CNT . The locations of the roughening temperature T R and the criticaltemperature T c are indicated by dotted vertical lines. As stated in the main text, the values of γ vl are taken from [82],using Monte Carlo results for temperatures above k B T /J = 2 . k B T /J = 2 .
0. Note that for temperatures smaller than k B T /J = 1 .
4, the magnetization at coexistence m coex is almost indistinguishable from its saturation value m coex = 1, and then our implementation of the Widom particleinsertion method cannot be applied to the chosen system sizes. III. TIME EVOLUTION OF THE DROPLET SIZE DISTRIBUTION AND DROPLET GROWTH RATES
We now consider time-dependence of the system where we start the system at time t = 0 in an equilibrated statein the vapor at µ = µ coex , but switch on at time t = 0 a chemical potential µ > µ coex (or equivalently, a positivemagnetic field H in the notation of Ising ferromagnet) at which the liquid is the stable phase. As mentioned in theintroduction, we consider heterogeneous in addition to homogeneous nucleation, choosing a L × L × D geometry withtwo walls at z = 1 and z = D , choosing surface fields H , H D such that H D favors the vapor but H (acting at thewall at z = 1) favors the liquid. The reason for this choice is, that by proper choice of H one can adjust the contactangle θ at which sessile wall-attached macroscopic droplets can occur. The barrier against heterogeneous nucleation∆ F ∗ het is predicted to be very much reduced, in comparison to the barrier ∆ F ∗ hom against homogeneous nucleation, ifthe contact angle is small, since [28, 29] ∆ F ∗ het = ∆ F ∗ hom f ( θ ) ,f ( θ ) = (1 − cos θ ) (2 + cos θ ) / . (25)As shown with the “lever rule” method [55, 56] indeed rather large wall attached droplets can be simulated inequilibrium with supersaturated vapor which have barriers of order ∆ F ∗ het ≈ k B T or so only, for suitable choices of H and H , and so a comparison with kinetic studies then seems reachable, and varying H over some range providesan additional variable to test the theory. (a)
100 1000 10000 100000 1e+06 1e+07 0.05 0.1 0.15 0.2 0.25 0.3 0.35 li f e t i m e τ M S chemical potential ∆µ /k B TH =0.0H =0.1H =0.2H =0.3H =0.4 (b)
100 1000 10000 100000 1e+06 1e+07 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 li f e t i m e τ M S chemical potential ∆µ /k B TH =0.0H =0.1H =0.2H =0.3H =0.4 FIG. 9: (Color online) Mean lifetime of metastable states in L × L × D Ising systems at k B T /J = 3 . k B T /J = 4 . µ/k B T , for various choices of H as indicated. This lifetime was measured as the firsttime, when the time-dependent magnetization m ( t ) becomes positive, as described in more detail in the main text. L = 60 and D = 30, with H D /J = − . x and y directions only. The average lifetimedecreases with higher surface fields H and highter chemical potentials ∆ µ/k B T . As a first step, preliminary runs were performed with the single spin flip Metropolis algorithm [83, 84] monitoringthe average lifetime of the metastable vapor. This was done using a large sample (10 ) of equilibrated initial statesat H = 0, where the considered field H (or chemical potential µ − µ coex , respectively) was then switched on and thetime recorded when the (initially negative) magnetization reaches the value m = 0 for the first time. Fig. 9 showsestimates for the resulting mean first passage times for a range of choices of H as a function of the field. When this“lifetime” of the state with m < , the system is rather unstable, nucleation occursfast and is followed by fast domain growth as well; such fast decays of unstable systems are not suitable for tests ofnucleation theory. It is seen, that in a rather narrow interval of fields H (for each value of H ) the lifetime increasesfrom 10 to 10 . Such parameter combinations ( H, H ) will be studied in the following only; if we would study caseswhere the lifetime is significantly larger than 10 , no critical droplet would be formed during affordable simulationtimes.Fig. 10 then shows typical time evolutions of the size distribution of n ℓ ( t ) of physical clusters of size ℓ , normalizingthem by the equilibrium cluster concentration n eq ℓ for H = 0. Here n eq ℓ is defined as the average number of physicalclusters per lattice site, and we have the sum rule ρ = N X ℓ =1 ℓn eq ℓ (26)5 (a) no r m a li z ed c l u s t e r c on c en t r a t i on n l ( t ) / n l , eq time in Monte Carlo steps per spin t (b) no r m a li z ed c l u s t e r c on c en t r a t i on n l ( t ) / n l , eq time in Monte Carlo steps per spin t FIG. 10: (Color online) Time dependence of the cluster concentration ratio n ℓ ( t ) /n eqℓ in equilibrium at H = 0, for the case k B T /J = 3, H /J = 0 .
4, and
H/J = 0 .
15 (a) and 0 .
17 (b). Curves show the choices ℓ = 40 , , . . . ,
120 (from bottom to top). since every site occupied by a particle must be part of some cluster. Of course, here we are only interested in nottoo small clusters, and hence Fig. 10 focuses on clusters with ℓ ≥
40. In the time evolution of n ℓ ( t ) /n eq ℓ we recognizethree regimes: for times of order t ≤ , n ℓ ( t ) /n eq ℓ is rapidly rising: this period of time corresponds to the relaxationfrom the initial state (where H = 0) towards the metastable state. In the latter, n ℓ ( t ) is almost constant for at leastone, or even several decades of time. Then a decay of these plateau values sets in, which is due to the fact that toomany much larger clusters have grown, the volume fraction of the system that is still in the metastable phase shrinks,and so less clusters of intermediate size (as studied in Fig. 10) are observed. This behavior is qualitatively similar toprevious studies (e.g. [9, 33]) which were based on the geometrical cluster definition, however. We also remark thatin the case of Fig. 10b the lifetime of the plateau extends up to t ≈ τ MS ≈ t distinctly less than τ MS should be analyzed.While some of the previous work on the studies of the kinetics of cluster growth in metastable Ising models (e.g.[9, 33]) tried to use directly n ℓ ( t ) to extract information on the validity of nucleation theory concepts, we here try toimplement a different concept. Namely, we follow the trajectories of individual (large) clusters with respect to theirsize in time, { ℓ i ( t ) } → { ℓ ′ i ( t + ∆ t ) } → . . . , where i is an index to label individual clusters. To ensure that the i ’thcluster at time t + ∆ t is actually a descendent of the i ’th cluster at the time t , we have to choose ∆ t small enough,and also record the location (center of gravity ~X i ( t )) and the components of the gyration radius R i,α ( t ) = ℓ i ( t ) − ℓ i ( t ) X k =1 x k,α − ℓ i ( t ) ( ℓ i ( t ) X k =1 x k,α ) / , (27)where x k,α is the α ’th Cartesian coordinate for the k ’th lattice site belonging to the cluster with label i , consistingof ℓ i ( t ) lattice sites. In each time step in which an analysis of the clusters is performed, the set of coordinates { ~X i ( t ) , ~R i ( t ) } is recorded.Note that there occurs the difficulty that the number of large clusters is not constant during the simulation: clustersform and decay or split into parts, and since we know that ℓ geom exceeds ℓ for each cluster, and the assignment of the“active bonds” according to Eq. (5) to identify from the geometrical cluster the associate physical clusters is a randomprocess, some random shift of ~X i ( t ) would occur even if we carry out two successive cluster identifications from thesame spin configuration (∆ t = 0). Of course, such shifts should be small in comparison with ~R i ( t ), but as ∆ t is chosennonzero it is clear that useful results are only obtained if ∆ t is small enough, and ℓ is large in comparison to clustersthat correspond to typical thermal fluctuations [58, 61]. Hence only clusters for which ℓ > ℓ min are considered (forthe temperature k B T /J = 3 . ℓ min = 10). So if by such criteria (for details see [89]) it is ensuredthat the i ’th cluster with size ℓ ′ i ( t + ∆ t ) is a descendent of the i ’th cluster with size ℓ i ( t ) at time t , we can define areaction rate Γ( ℓ ) as Γ( ℓ ) = (cid:28) ℓ ′ i ( t + ∆ t ) − ℓ i ( t )]∆ t (cid:29) ℓ i ( t ) (28)6Here the index ℓ i ( t ) stands for an average over a sampling of all cluster trajectories recorded in the simulation, anda smoothing procedure of the (otherwise too noisy) data with a triangular smoothing function [89] was applied. Ofcourse, in order to collect statistically significant data on Γ( ℓ ), it is necessary to perform many runs for each parametercombination ( T, H, H ) that is studied. We observe that the lifetime of the metastable stale in such runs is fluctuatingdramatically, and so it is necessary to choose the run time of each run individually, rather than the same for all runs.It was decided to stop each run automatically when the largest cluster size ℓ max i ( t ) = L /
20 was reached. Of course,then only clusters with ℓ ≪ L /
20 could be studied, to avoid artifacts caused by this cutoff. (a) -10-5 0 5 10 15 20 25 0 50 100 150 200 250 300 350 400 c l u s t e r r ea c t i on r a t e Γ ( l ) cluster size l/m ′′ H/J=0.43H/J=0.45H/J=0.49 (b) -15-10-5 0 5 10 0 50 100 150 200 250 300 350 400 c l u s t e r r ea c t i on r a t e Γ ( l ) cluster size l/m ′′ H/J=0.26H/J=0.28H/J=0.32
FIG. 11: (Color online) Cluster reaction rate Γ( ℓ ) versus cluster size ℓ for a bulk system at k B T /J = 3 . k B T /J = 3 . H /J = 0 . H are shown,increasing from bottom to top, for clusters grow more likely with larger field strength. Full curves are based on Eq. (28), whilebroken curves are based on the approximation based on the use of the largest cluster only. In the heterogeneous case (b), bothmethods agree very well at lower fields because there is only one larger cluster in the system at a time. While Eq. (28) is based on using all clusters ℓ i ( t ) > ℓ min at each time t , one can simplify matters by restrictingthe analysis only to the trajectory of the biggest cluster in the system [89]. When one does this, one ignores possibleproblems from the fact that from time to time the identity of the largest cluster changes. Fig. 11 shows now typicalresults for Γ( ℓ ), using both this latter approximation and the method based on Eq. (28). Both methods yield similartrends, although they differ somewhat in detail (particularly in the case of homogeneous nucleation in the bulk).What one expects theoretically for Γ( ℓ ) is a monotonous increase of Γ( ℓ ) with ℓ , where Γ( ℓ ) is negative for clusterssmaller than the critical cluster size ℓ ∗ while Γ( ℓ ) is positive for ℓ > ℓ ∗ . The method yields a second positive part for ℓ slightly larger than ℓ min . This is an artifact of ignoring clusters smaller than ℓ min , which vanishes for ℓ min → ℓ ) at the right hand side, and convert to the clustervolume V ∗ according to Eq. (11), we obtain the data shown in Fig. 12. It is seen that for a given value of the field H (or µ − µ coex , respectively) classical nucleation theory underestimates the volume of the critical cluster: in otherwords, a given volume V ∗ leads to a larger value of H . Since (according to classical nucleation theory and Eq. (4))the barriers scale as H − , this means when one studies nucleation barriers as function of cluster volume V ∗ or clusterradius R ∗ , one finds lower barriers in the simulation rather than predicted. Qualitatively, the data from the presentanalysis of cluster kinetics confirm the findings from the lever rule method of Winter et al. [55], as far as homogeneousnucleation is concerned (Fig. 12a), although some questions on systematic errors in both methods have not been fullysettled. Nevertheless, the qualitative agreement between these quite different approaches is satisfactory. For the caseof heterogeneous nucleation, however, for a given supersaturation the critical cluster volume predicted from clusterkinetics (Fig. 11) is distinctly larger than the corresponding results from the static methods. We have no explanationfor this discrepancy.7 (a)
80 90 100 120 140 160 180 200 220 240 260 0.24 0.25 0.26 0.27 0.28 0.29 0.3 0.31 0.32 0.33 0.34 c r i t i c a l c l u s t e r v o l u m e V * chemical potential ∆µ /k B T allbiggestlever ruleCNT (b)
100 200 300 400 500 600 700 800 0.12 0.15 0.18 0.21 0.24 0.27 c r i t i c a l c l u s t e r v o l u m e V * chemical potential ∆µ /k B TH /J=0.0H /J=0.1H /J=0.2H /J=0.3H /J=0.4CNTlever rule FIG. 12: (Color online) Log-log plot of the critical cluster volume V ∗ against the normalized chemical potential difference ∆ µ for k B T /J = 3 for homogeneous nucleation in the bulk (a) and heterogeneous nucleation at a wall with several choices of thesurface field H /J , as indicated. The broken straight line is the prediction of the classical nucleation theory, Eq. (3), amendedby the enhancement factor (1 .
064 at k B T /J = 3) taken from Fig. 8. The dotted lines are the corresponding results based onthe leverrule. Results for the method based on Eq. (28) are shown in both plots, the method based on the biggest cluster isonly drawn in (a), as both methods yield the same results in the heterogeneous case. In case (a) the lever rule data were takenfrom a system at linear dimension L = 15, and thus possibly affected by some finite size effects. Note also that in (b) theprediction of the classical nucleation theory and the leverrule results lie on top of each other. IV. CONCLUDING DISCUSSION
In the present work, we have studied aspects of nucleation theory by simulation of clusters and their dynamics, usingthe Ising (lattice gas) model on the simple cubic lattice. Both homogeneous nucleation and heterogeneous nucleationat planar walls (where a “surface field” acts) have been considered.Although many aspects of this problem have been studied before in works of various groups extending over severaldecades, most of the previous work is inconclusive since it relied on the use of the “geometric” cluster definition.We have given evidence that this geometric cluster definition does not yield correct results for large clusters at thetemperatures far above the roughening transition temperature where the clusters have spherical shape; at temper-atures below the roughening transition temperature the geometric clusters and the “physical clusters” are basicallyindistinguishable, but due to the pronounced anisotropy effects a simple analysis of nucleation phenomena is notpossible.However, in the limit of large droplet volumes V → ∞ , where one can neglect any corrections to the decompositionof the droplet formation free energy into the bulk term plus a surface correction, one can compute the nucleation freeenergy barrier ∆ F ∗ from measuring the excess chemical potential ∆ µ that is in equilibrium with a given V . Fig. 8shows the enhancement of ∆ F ∗ with respect to the standard result for spherical droplets (using Eqs. (3), (4)). Wethus show that in the Ising (lattice gas) model this enhancement gradually rises from unity as the temperature islowered from the critical temperature, reaches almost 10% at T /T c = 0 .
6, and rises steeply below the rougheningtemperature towards the low temperature limit 6 /π ≈ .
91 (Fig. 8). This enhancement reflects the consequences ofthe anisotropy of the interfacial free energy, such as the gradual crossover of the droplet shape from a sphere to acube (Fig. 6).On the other hand, we demonstrate that physical clusters do give consistent results, at least in the bulk whenone is concerned with homogeneous nucleation. We show that in the limit where the droplets get macroscopicallylarge, they converge against the simple lever rule predictions. However, we do find an (unexpected) surface excessin the particle number of such clusters also in this case. We also demonstrate the validity of the relation betweenchemical potential (of the supersaturated vapor) and the droplet radius that classical nucleation theory predicts forlarge droplets near the critical temperature. We also give evidence that the droplet-vapor interface is broadened dueto capillary waves; we remind the reader that mean-field type theories and density functional theories [11, 14] cannotinclude such capillary wave effects (which also should give rise to a correction term on the droplet formation freeenergy, not yet included in Eq. (2)).We would also like to stress that many of our considerations can be carried over to a study of clusters in d = 2dimensions, where a construction as in Fig. 1 also holds. However, we expect two distinctions: (i) The rougheningtransition temperature T R is zero, so the crossover of droplet shape from the circle to the square occurs without anysingularity even for arbitrarily large droplets. (ii) Percolation coincides with the critical point, but geometrical clustersstill are too large, and to describe nucleation, physical clusters defined via Eq. (5) should also be used. Of course,it would be very desirable to carry these considerations over to nucleation in off-lattice models of fluids. However, aprecise analogue of Eq. (5) is still not known, and hence other concepts to define physical clusters [58] need to beused, if one wishes to study nucleation near the critical point.In the second part we present a first study of the time evolution of the cluster population based on the “physicalcluster” definition. However, due to the large computer resources needed for this study, only data at a single tem-perature ( k B T /J = 3 .
0) are presented. In order to allow a comparison of this part of the study with our results onstatic properties of critical droplets, as studied in the first part of the paper, we use a criterion to estimate the criticaldroplet size from the balance between droplet growth and shrinking processes. In the case of homogeneous nucleation,the results obtained in this way are roughly compatible with the results obtained from the static lever rule method.Studying droplet volumes in the range from 100 to 200, clear deviations from the classical nucleation theory are seen,which can be attributed to a decrease of the nucleation barrier due to fluctuation effects. However, in the case ofheterogeneous nucleation, a rather large discrepancy between the results of the statics and dynamics of droplets isfound. This discrepancy is not understood yet, and must be left as a challenging problem for the future.
Acknowledgements : We thank D. Winter for providing us with the data from Ref. [55] that were included inFig. 12 for comparison. One of us (F. S.) thanks the Deutsche Forschungsgemeinschaft for partial support under grantNo VI 237/4-3 (SPP 1296).9 [1] M. Volmer and A. Weber, Z. phys. Chem. , 277 (1926)[2] R. Becker and W. D¨oring, Ann. Phys. , 719 (1935)[3] Yu. B. Zeldovich, Acta Physicochim. URSS , 1 (1943)[4] J. Feder, K. C. Russell, J. Lothe, and G. M. Pound, Adv. Phys. , 111 (1966)[5] H. Reiss, J. L. Katz, and E. R. Cohen, J. Stat. Phys. , 83 (1968)[6] J. S. Langer, Ann. Phys. (N.Y.) , 108 (1967); Ann. Phys. , 258 (1969)[7] Nucleation , edited A. C. Zettlemoyer (M. Dekker, New York, 1969)[8] F. F. Abraham,
Homogeneous Nucleation Theory (Academic, New York, 1974)[9] K. Binder and D. Stauffer, Adv. Phys. , 343 (1976)[10] K. Binder, Rep. Progr. Phys. , 783 (1987)[11] D. W. Oxtoby and R. Evans, J. Chem. Phys. , 7521 (1988)[12] A. Dillmann and G. E. A. Meier, Chem. Phys. Lett. , 71 (1989)[13] H. Reis, A. Tabazadeh, and J. Talbot, J. Chem. Phys. , 1266 (1990)[14] D. W. Oxtoby, J. Phys.: Condens. Matter , 7627 (1992)[15] A. Laaksonen, V. Talanquer, and D. W. Oxtoby, Annu. Rev. Phys. Chem. , 489 (1995)[16] P. R. ten Wolde and D. Frenkel, J. Chem. Phys. , 9901 (1998)[17] D. Kashchiev, Nucleation, Basic Theory with Applications (Butterworth-Heinemann, Oxford, 2000)[18] A. C. Pan and D. Chandler, J. Phys. Chem. B , 19681 (2004)[19] Nucleation,
Comptes Rendus Physique , Vol. 7 (2006), special issue, edited by S. Balibar and J. Villain.[20] K. Binder, in
Kinetics of Phase Transitions , edited by S. Puri and V. Wadhavan (CRC Press, Boca Raton, 2009), Chap.2[21] M. Schrader, P. Virnau, and K. Binder, Phys. Rev. E , 061104 (2009)[22] S. Ryu and W. Cai, Phys. Rev. E , 030601(R) (2010)[23] B. J. Block, S. K. Das, M. Oettel, P. Virnau and K. Binder, J. Chem. Phys. , 154702 (2010)[24] A. Troester, M. Oettel, B. Block, P. Virnau and K. Binder, J. Chem. Phys. , 064709 (2012)[25] Y. Vilsanen, R. Strey, and H. Reiss, J. Chem. Phys. , 4680 (1993)[26] A. Fladerer and R. Strey, J. Chem. Phys. , 164710 (2006)[27] K. Iland, J. W¨olk, R. Strey and D. Kashchiev, J. Chem. Phys. , 154506 (2007)[28] D. Turnbull, J. Appl. Phys. , 1022 (1950)[29] D. Turnbull, J. Chem. Phys. , 198 (1950)[30] A. W. Castleman, Adv. Colloid Interface Sci. , 73 (1979[31] H. Biloni, in Physical Metallurgy (R. W. Cahn and P. Haaasen, eds) Amsterdam, North-Holland, 1983) p. 477[32] J. Curtius, Comptes Rendus Phys. , 1027 (2006)[33] K. Binder and H. M¨uller-Krumbhaar, Phys. Rev. B , 2328 (1974)[34] J. S. Rowlinson and B. Widom, Molecular Theory of Capillarity (Oxford Univ. Press, Oxford, 1982)[35] K. Binder and D. Stauffer, J. Stat. Phys. , 49 (1972)[36] K. Binder and E. Stoll, Phys. Rev. Lett. , 47 (1973)[37] K. Binder and M. H. Kalos, J. Stat. Phys. , 363 (1980)[38] D. Stauffer, A. Coniglio, and D. W. Heermann, Phys. Rev. Lett. , 1299 (1982)[39] D. Stauffer, Int. J. Mod. Phys. C , 1071 (1992)[40] H. Furukawa and K. Binder, Phys. Rev. A , 556 (1982)[41] J. Marro and R. Toral, Physica , 563 (1983)[42] O. Penrose, J. Lebowitz, J. Marro, M. Kalos, and J. Tobochnik, J. Stat. Phys. , 399 (1984)[43] D. W. Heermann, A. Coniglio, W. Klein, and D. Stauffer, J. Stat. Phys. , 447 (1984)[44] H. Tomita and S. Miyashita, Phys. Rev. B , 8886 (1992)[45] P. A. Rikvold, H. Tomita, S. Miyashita, and S. W. Sides, Phys. Rev. E , 5080 (1994)[46] P. A. Rikvold, and B. M. Gorman, in Annual Reviews of Computational Physics, Vol. 1 (D. Stauffer, ed.) p. 149 (WorldScientific, Singapore, 1994)[47] M. Acharyya and D. Stauffer, Eur. Phys. J. B , 571 (1998)[48] H. Vehkamaki and I. J. Ford, Phys. Rev. E , 6483 (1999)[49] V. A. Shneidman, K. A. Jackson, and K. M. Beatty, Phys. Rev. B , 3579 (1999)[50] V. A. Shneidman, K. A. Jackson, and K. M. Beatty, J. Chem. Phys. , 6932 (1999)[51] R. A. Ramos, P. A. Rikvold, and M. A. Novotny, Phys. Rev. B , 9053 (1999)[52] S. Wonczak, R. Strey, and D. Stauffer, J. Chem. Phys. , 1976 (2000)[53] A. T. Bustillos, D. W. Heermann and C. E. Cordeiro, J. Chem. Phys. , 4864 (2004)[54] K. Brendel, G. T. Barkema, and H. van Beijeren, Phys. Rev. E , 031601 (2005)[55] D. Winter, P. Virnau and K. Binder, J. Phys.: Condens. Matter , 464118 (2009)[56] D. Winter, P. Virnau and K. Binder, Phys. Rev. Lett. , 225703 (2009)[57] H. M¨uller-Krumbhaar, Phys. Lett. , 27 (1974)[58] K. Binder, Ann. Phys. , 390 (1976)[59] A. Coniglio and W. Klein, J. Phys. A , 2775 (1980) [60] R. H. Swendsen and J. S. Wang, Phys. Rev. Lett. , 86 (1987)[61] M. D’Onorio De Meo, D. W. Heermann, and K. Binder, J. Stat. Phys. , 585 (1990)[62] A. Coniglio, J. Phys. A8 , 1773 (1975)[63] A. Coniglio, F. Peruggi, C. Nappi and L. Russo, J. Phys. A , 205 (1977)[64] K. Binder, Solid State Commun. , 191 (1980)[65] S. Hayward, D. W. Heermann, and K. Binder, J. Stat. Phys. , 1053 (1987)[66] P. W. Kasteleyn and C. M. Fortuin, J. Phys. Soc. Japan (Suppl) 11 (1969)[67] C. M. Fortuin and P. W. Kasteleyn, Physica , 536 (1972)[68] U. Wolff, Phys. Rev. Lett. , 361 (1989)[69] K. Kawasaki, in Phase Transitions and Critical Phenomena, Vol. 2 (C. Domb and M. S. Green, eds.) (Academic Press,London, 1972)[70] H. M¨uller-Krumbhaar and K. Binder, J. Stat. Phys. , 1 (1973)[71] J. D. Weeks, in Ordering in Strongly Fluctuating Condensed Matter Systems (T. Riste, ed.) p. 293 (Plenum Press, NewYork, 1980)[72] H. Van Beijeren and I. Nolden, in
Structure and Dynamics of Surfaces II (W. Schommers and P. Blanckenhagen, eds.) p.259 (Springer, Berlin 1987)[73] G. Wulff, Z. Kristallogr. Mineral. , 449 (1901)[74] C. Herring, Phys. Rev. , 87 (1951)[75] C. Rottman and M. Wortis, Phys. Rev. B , 6274 (1981)[76] M. Holzer, Phys. Rev. B , 10570 (1990)[77] K.K. Mon, S.Wansleben, D.P. Landau and K.Binder, Phys. Rev. B , 7089 (1989)[78] M. Hasenbusch and K. Pinn, J. Phys. A , 63 (1997)[79] A. M. Ferrenberg and D. P. Landau, Phys. Rev. B , 5081 (1991)[80] V. Privman, Phys. Rev. Lett. , 183 (1988)[81] M. Hasenbusch and K. Pinn, Physica A , 342 (1993)[82] M. Hasenbusch and K. Pinn, Physica A , 189 (1994)[83] D. P. Landau and K. Binder, A Guide to Monte Carlo Simulation in Statistical Physics, 3rd ed (Cambridge Univ. Press,Cambridge 2009)[84] K. Binder and D. W. Heermann,
Monte Carlo Simulation in Statistical Physics: An Introduction . 5th ed. (Springer, Berlin2010)[85] K. Binder, Physica A , 99 (2003)[86] D. Stauffer and A. Aharony,
Introduction to percolation theory (Taylor and Francis, London, 1994)[87] B. Widom, J. Chem. Phys. , 2808 (1963)[88] B. Widom, J. Phys. Chem. , 869 (1982)[89] F. Schmitz, diploma thesis (Johannes Gutenberg-Universit¨at Mainz, 2011, unpublished.)[90] P. Virnau, M. M¨uller, J. Chem. Phys. , 10925 (2004)[91] P. Virnau, M. M¨uller, L. G. MacDowell, K. Binder, J. Chem. Phys. , 2169 (2004)[92] L.G. MacDowell, P.Virnau, M.M¨uller, and K.Binder, J. Chem. Phys. 120, 5293 (2004)[93] B. A. Berg, T. Neuhaus, Phys. Rev. Lett. , 9 (1992)[94] F. Wang and D. P. Landau, Phys. Rev. Lett. , 2050 (2001)[95] F. Wang and D. P. Landau, Phys. Rev. E , 056101 (2001)[96] S. K. Das and K. Binder, Phys. Rev. Lett. , 235702 (2011)[97] M. P. A. Fisher and M. Wortis, Phys. Rev. B , 6252 (1984)[98] D. Jasnov, Rep. Prog. Phys. , 1059 (1984)[99] K. Binder, M. M¨uller, F. Schmid, and A. Werner, Adv. Colloid and Interface Science , 237 (2001)[100] K. Binder and M. M¨uller, Int. J. Mod. Phys. C , 1093 (2000)[101] G.V. Wulff, Z. Kristallogr. Mineral.34