Numerical test of finite-size scaling predictions for the droplet condensation-evaporation transition
Andreas Nußbaumer, Johannes Zierenberg, Elmar Bittner, Wolfhard Janke
NNumerical test of finite-size scaling predictions forthe droplet condensation-evaporation transition
Andreas Nußbaumer , Johannes Zierenberg , Elmar Bittner andWolfhard Janke Institut f¨ur Physik, Johannes Gutenberg Universit¨at Mainz, Staudinger Weg 7, D-55128Mainz, Germany Institut f¨ur Theoretische Physik, Universit¨at Leipzig, Postfach 100 920, D-04009 Leipzig,Germany Institut f¨ur Theoretische Physik, Universit¨at Heidelberg, Philosophenweg 16, D-69120Heidelberg, GermanyE-mail: [email protected]
Abstract.
We numerically study the finite-size droplet condensation-evaporation transitionin two dimensions. We consider and compare two orthogonal approaches, namely at fixedtemperature and at fixed density, making use of parallel multicanonical simulations. Theequivalence between Ising model and lattice gas allows us to compare to analytical predictions.We recover the known background density (at fixed temperature) and transition temperature(at fixed density) in the thermodynamic limit and compare our finite-size deviations to thepredicted leading-order finite-size corrections.
1. Introduction
Droplet formation is an essential process in nature with a variety of analogues in biologicalsystems and material science. Of course, this generally involves non-equilibrium processes, e.g.,the formation of nucleation prerequisites from local fluctuations while the surrounding gas actsas a density bath. Here, we consider instead a canonical setup in a finite system of size V with fixed temperature T and particle number N . If the gas is supersaturated, a variationof N or T results in the formation of equilibrium droplets [1, 2, 3, 4, 5]. In simple terms,the particle excess is subdivided to form a single macroscopic droplet in equilibrium with thesurrounding vapor plus remaining excess. The resulting theory has been supported by numerouscomputational studies at fixed temperature, including the two-dimensional lattice gas [6, 7, 8]and three-dimensional Lennard-Jones gas [9, 10, 11]. The orthogonal approach at fixed densityhas received less attention [12], but recently enabled us to come closer to the asymptotic scalingregime for two- and three-dimensional lattice gas and three-dimensional Lennard-Jones gas [13].In the following, we will compare the leading-order scaling corrections for a two-dimensionallattice gas at fixed temperature and fixed density from analytical predictions with numericalresults.
2. Model
We consider a lattice gas in d = 2 dimensions. Excluded volume is modeled as lattice sites beingeither occupied by exactly one particle or empty, n i = { , } . Short-range interaction is included a r X i v : . [ c ond - m a t . s t a t - m ec h ] M a y y nearest-neighbor interaction ( (cid:104) i, j (cid:105) ). The Hamiltonian is then H = − (cid:88) (cid:104) i.j (cid:105) n i n j , (1)which is equivalent to an Ising model at fixed magnetization with T Is = 4 T and coupling constant J = 1 [14]. This originates in the identification of the spin state s i = 2 n i −
1, which allows oneto rewrite the Ising Hamiltonian H Is = − (cid:88) (cid:104) i,j (cid:105) s i s j = 4 H − d ( V − N ) , (2)considering that, on a simple hypercubic lattice, the sum over nearest neighbors yields d contributions and N = (cid:80) n i . The equivalence is then established by equating the Boltzmannfactors exp( − β H ), where β = ( k B T ) − . Evaluating β Is H Is = β H , where the constant shiftis neglected due to physical unimportance, results in a simple rescaling of the temperature.When making use of the Ising equivalence all energy-related observables have to be rescaledcorrespondingly. The canonical Ising model is then equivalent to a grand-canonical lattice gas.Exploiting the equivalence to the 2D Ising model, the critical point is located at β c = 4 β Is c =2 ln (cid:0) √ (cid:1) ≈ .
763 or k B T c ≈ . m is observed and described by the Onsager-Yang solution [15, 16]: m ( β Is ) = (cid:2) − sinh − (2 β Is ) (cid:3) / . (3)This is directly related to the (grand-canonical) equilibrium background density ρ = (1 − m ) / χ = βκ = ˆ κ [8]and may be evaluated from sufficiently long series expansions (see, e.g., Ref. [17] and referencestherein), where χ ( β Is ) = β Is n (cid:88) i =0 c i u i with u = 12 sinh(2 β Is ) , (4)and c = { , , , , , , , , , , , , , ... } , hereconsidered up to the 300 th term. The equilibrium shape of a 2D Ising droplet is describedby the Wulff plot (or shape), given by [18] W = 4( β Is ) (cid:90) β Is σ d x cosh − (cid:20) cosh (2 β Is )sinh(2 β Is ) − cosh( x ) (cid:21) , (5)where σ = 2 + ln[tanh( β Is )] /β Is and cosh − is referring to the inverse hyperbolic cosine. Thiswill be relevant for the surface free energy of a (Wulff shaped) droplet of unit volume τ IsW = 2 √ W .Being energy-related, the interface tension gets converted as τ W = τ IsW /
3. Method
In order to obtain numerical data at fixed temperature and at fixed density, we employed twodifferent kinds of simulations. In both cases the underlying algorithm is the multicanonicalmethod [19, 20, 21, 22]. The condensation transition is a first-order phase transition for whichthis method is well suited because it potentially allows one to overcome barriers in the freeenergy. The principle idea is to replace the Boltzmann weight exp( − βE ) by an a priori unknownweight function W ( E ), which is iteratively adapted in order to yield a flat histogram over a The coefficients were obtained from [17]. igure 1.
Snapshots of a two-dimensional lattice gas at fixed density above (left) and below(right) the condensation-evaporation temperature, showing a homogeneous gas phase and adroplet in equilibrium with surrounding vapor, respectively. Both plots show N = 2500 particleson a lattice of linear size L = 500, i.e., ρ = 0 . c it holds h ( E min ) /h ( E max ) > c . In part, we further make use of aparallel implementation [23, 24], which exploits the fact that in each iteration the histogram isan estimate of the probability distribution belonging to the current weight function. This allowsone to distribute the sampling to independent Markov chains and to obtain a joint estimate as asimple sum of individual histograms. The procedure scales very well for the problem at hand [8].The involved Monte Carlo updates include particle displacements to free nearest neighborsand random jumps to free sites. This corresponds to local and global Kawasaki updates. Errorsare obtained using the Jackknife method [25] and standard error propagation. For further detailsof the employed methods we refer to Refs. [6, 7, 8, 13].
4. Theory of leading-order correction
A natural approach to particle condensation/evaporation is the consideration of a grand-canonical ensemble at fixed temperature, where the system relaxes to an equilibrium backgroundcontribution N = ρ V . Below the corresponding critical temperature the system is dominatedeither by particles (fluid branch) or by void space (gas branch) and ρ (cid:54) = 1 /
2. Let us nowconsider the dilute gas branch well below the critical temperature and fix
N > N . This resultsin the canonical ensemble of a supersaturated gas with particle excess δN = N − N . Initially,this excess goes into the gas phase while for sufficiently large excess droplet formation occurs.In equilibrium droplet formation, the probability for intermediate-sized droplets was shownto vanish [3] and the scenario reduces to a homogeneous gas phase and an inhomogeneousphase of a droplet in equilibrium with surrounding vapor. This may be considered as theinterplay of entropy maximization by fluctuations in the gas phase and energy minimizationby forming a droplet [3, 4], see also Fig. 1. At fixed temperature it is possible to consider(fixed) thermal fluctuations and relate them to infinite-size temperature-dependent quantities,like the isothermal compressibility ˆ κ and the normalized surface free energy τ W , see Sec. 2.The free energy may then be approximated by a contribution F fluc from the fluctuation ofarticle excess δN and a contribution F drop from the single macroscopic droplet of size V D : F fluc = ( δN ) κV and F drop = τ W ( V D ) d − d . (6)These contributions are idealized with possible sources of corrections in both the Gaussianapproximation and the droplet shape for finite systems.In the two-phase scenario, the particle excess may be decomposed into the excess inside thedroplet δN D and the excess in the fluctuating phase δN F , i.e., δN = N − N = δN D + δN F .Linking the droplet size to the particle excess inside the droplet, one expects δN D = ( ρ L − ρ ) V D ,where ρ L and ρ are the background liquid and gas density, respectively. Then, one mayintroduce a scalar fraction λ = δN D /δN, (7)such that δN D = λδN and δN F = (1 − λ ) δN . The total free energy F = F drop + F fluc becomes F = τ W (cid:18) λδNρ L − ρ (cid:19) d − d + (1 − λ ) ( δN ) κV = τ W (cid:18) δNρ L − ρ (cid:19) d − d (cid:16) λ d − d + ∆(1 − λ ) (cid:17) , (8)with a dimensionless “density” parameter∆ = ( ρ L − ρ ) d − d κτ W ( δN ) d +1 d V = ( ρ L − ρ ) d − d κτ W ( ρ − ρ ) d +1 d V d . (9)At fixed temperature ρ L , ρ , ˆ κ, τ W are constants and ∆ may be interpreted as an unusual density.In principle, all constants may be estimated and for the present case, equivalent to the 2D Isingmodel, the parameters are even known exactly or with very high precision.This leading-order formulation allows one to obtain the fraction of particles inside the largestdroplet λ as a function of ∆ in the limit of large systems, by minimizing Eq. (8) with respectto λ . It turns out (for details see Refs. [3, 4]) that there exists a constant ∆ c below whichno droplet forms ( λ = 0) and above which a single macroscopic droplet exists with non-trivial λ > λ c : ∆ c = 1 d (cid:18) d + 12 (cid:19) d +1 d = 0 . ... and λ c = 2 d + 1 = 2 / . (10)The result λ (∆) describes the expectation value of the equilibrium droplet size in the limitof large systems without any free parameter. As mentioned before, this already includes theleading-order finite-size corrections for idealized assumptions. In fact, Eq. (9) may be rewrittenat each finite-size transition density ρ c , where ∆( ρ c ) = ∆ c . For a lattice gas model with particle-hole symmetry, ρ L = 1 − ρ , this yields to leading order ρ c = ρ + (cid:32) κτ W ∆ c (1 − ρ ) d − d (cid:33) dd +1 V − d +1 , or m c = m − m (cid:18) χτ IsW ∆ c m (cid:19) dd +1 V − d +1 , (11)where the second formulation is in terms of the Ising model at fixed magnetization. This is thenotation of Biskup et al. [3, 4], but is in quantitative agreement with the (independent) resultof Neuhaus and Hager [2]. For the 2D Ising model, they find the same leading scaling behavior∆ m ( L ) = A cond L − / , with the amplitude A cond = 0 . ... for β Is = 0 . A = 0 . ... for the corresponding constants. For β Is = 0 . m (cid:39) . χ (cid:39) . τ IsW (cid:39) . . .
25 0 . .
35 0 . .
45 0 . . − a ρ / T a T / ρ ρ ( T ) T = T Is / a T a ρ Figure 2.
Leading-order normalized finite-size correction amplitude at fixed temperature T andfixed density ρ . The densities ρ ( T ) on the lower x-axis are calculated from the Onsager-Yangsolution Eq. (3). The data points are obtained by numerically evaluating Eq. (14) and Eq. (15).Dashed lines indicate the parameters considered below.In the following, we will numerically review the leading-order finite-size scaling predictions.We consider directly the scaling of the factual transition density (or magnetization) at fixedtemperature [6, 7] and compare to the scaling of the transition temperature at fixed density [13].Albeit the possibility of logarithmic corrections, we consider empirical higher-order correctionsas powers of the leading term if necessary. Our fit ansatz for the leading-order correction is ρ c = ρ + a T V − / + O (cid:16) V − / (cid:17) at fixed T , (12) T c = T + a ρ V − / + O (cid:16) V − / (cid:17) at fixed ρ. (13)By comparing to Eq. (11), a T can be related to ∆ c . Similarly, a relation for a ρ followsfrom a Taylor expansion around T of a reformulated Eq. (9), namely ∆ / V − / = f ( ρ, T ).From the equivalence to the Ising model, we identify f ( ρ, T ) with f ( m, T Is ) = ( m ( T Is ) − m )[ m ( T Is ) / / [ χ ( T Is ) τ IsW ( T Is )] − / (see Ref. [13] for details), which leads to∆ c = a / T m / /χτ IsW at fixed
T , (14)∆ c = (cid:2) a ρ f (cid:48) ( m, T Is0 ) (cid:3) / at fixed ρ. (15)This may be evaluated numerically considering Eqs. (3)–(5) and (10) to yield estimates of theleading-order correction. Figure 2 shows results for a range of temperatures and accordingdensities, related by the Onsager-Yang solution Eq. (3). The relative leading-order correctionsat fixed temperature are almost an order of magnitude larger than at fixed density. However, inthe practical finite-size scaling analysis higher-order corrections are of more significance, whichcannot be estimated by the current theory.
5. Results
We will start the discussion with our findings at fixed temperature k B T = 0 . < ∆ c0 . . .
81 0 0 . . . λ ∆ analyti L = L = L = L = L = . . . . . . .
018 0 0 .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . ρ ρ c V − / datahigher-order (cid:28)tanalyti al predi tion Figure 3.
Droplet formation at fixed temperature k B T = 0 .
375 [6, 7]. (left) Fraction of excessversus rescaled dimensionless density ∆. (right) Finite-size scaling of the transition density ρ c .The black line shows the predicted leading-order scaling and the blue line is a higher-order fit tothe data. ρ is the analytically known thermodynamic limit. The dotted vertical line indicatesthe end of the fitting range ( L = 80).the system is in the gaseous phase and for ∆ > ∆ c a single macroscopic droplet forms inaccordance with theory. In fact, the results are obtained in the equivalent formulation of anIsing model at fixed magnetization, but here discussed in the generic formulation of a latticegas. In order to ensure ergodic sampling, we augmented the Kawasaki Monte Carlo updatewith a multicanonical scheme to locally sample a flat histogram including transition states andensuring a good sampling of the energy probability distribution up to suppressions of 30 ordersof magnitude. The fraction of excess in the largest droplet, Eq. (7), is obtained by a two-stepprocess. Firstly, all clusters in the system are determined, where we define a cluster as alike sitesconnected via the nearest-neighbor property. The largest cluster is the background (or gaseousphase), while the second largest cluster forms the droplet we are interested in. Then, the volumeof the droplet is identified as all sites confined by the boundary, including the enclosed holes.The fraction λ is then the ratio of the droplet volume and the expected equilibrium volume offull excess V δ = ( N − ρ L ) / (1 − ρ ). The result of this procedure is shown in Fig. 3 (left).We also show the limiting case that can be found as the solution of Eq. (8) as solid black line.Every data point is the average of 10 Monte Carlo sweeps ( L updates). At the infinite-sizetransition density, the analytical prediction is λ c = 2 /
3, see Eq. (10). Hence, we estimate thefinite-size transition point as the density ρ c for which λ ( ρ c , L ) = 2 /
3. Technically, we do a linearinterpolation of the data points in Fig. 3 (left) and search for the intersection.The resulting scaling of the transition density is shown in Fig. 3 (right) in dependence on V − / as expected from Eq. (12). The thermodynamic limit is given by the Onsager solution, ρ (cid:39) . k B T = 0 .
375 that ∆ c (cid:39) . ≈ (6 . a T ) / and thus a T ≈ . shown in the figure as a black line. The data does not allow for a qualitativelysatisfying leading-order fit, but the largest system size is already close to the predicted finite-sizedeviation. Considering the next higher empirical correction, i.e., ρ c = ρ + a T V − / + b T V − / ,allows for a decent fit which yields for L >
80 the result ρ = 0 . a T = 0 . Q ≈ .
02. The estimated limit is in good agreement with the For k B T = k B T Is / .
375 we obtain m (cid:39) . ρ (cid:39) . χ = ˆ κ (cid:39) . τ IsW (cid:39) . τ W (cid:39) . − − . .
32 0 .
34 0 .
36 0 .
38 0 . .
42 0 .
44 0 . C V T L = L = L = L = L = L = . . .
01 0 .
02 0 .
03 0 .
04 0 .
05 0 .
06 0 . T T c V − / dataleading-order (cid:28)tanalyti al predi tion Figure 4.
Droplet formation at fixed density ρ = 10 − [13]. (left) Specific heat with exemplarydata points that indicate the size of the error. (right) Finite-size scaling of the transitiontemperature. The black line shows the numerical evaluation of Eq. (9) at fixed ρ and the blueline is a leading-order fit. T is the analytically known thermodynamic limit. The dotted verticalline indicates the end of the fitting range ( L = 400).analytical prediction ρ marked by the arrow. Using Eq. (15) with error propagation yields theestimate ∆ c = 0 . , /
3] strongly influences the size of the higher-ordercorrections.We next turn to the orthogonal setup of a fixed density ρ = 10 − with varyingtemperature [13]. Here, the multicanonical method is more straightforward and for each systemsize we may perform a single (yet parallel) simulation with up to 128 cores and 1 . × measurements in the final production run. This yields a full temperature range of expectationvalues and the transition temperature may be determined precisely as the peak of the specificheat C V = k B β (cid:0) (cid:104) E (cid:105) − (cid:104) E (cid:105) (cid:1) /N , see Fig. 4 (left). For details see Ref. [13], from which werecapture the data in order to compare to the leading-order results at fixed temperature. For ρ = 10 − , the numerical evaluation of the Onsager solution Eq. (3) yields T (cid:39) . c (cid:39) . ≈ ( − . a ρ ) / and thus a ρ ≈ − . L (cid:39) L ≥
400 yields T = 0 . a ρ = − . Q ≈ .
25. The limit is in good agreement with the analytical value.Using Eq. (15) with error propagation yields the estimate ∆ c = 0 . . Conclusions We numerically investigated the leading-order finite-size scaling corrections of the two-dimensional droplet condensation-evaporation transition at fixed temperature and density,respectively. In both cases we could recover the analytically known thermodynamic limits,in part by considering empirical higher-order corrections. The leading-order corrections werefound to be slightly smaller than but consistent with the analytical predictions [3, 4]. This showsthat the established theory is able to qualitatively describe the finite-size deviations; in case ofthe two-dimensional lattice gas already for system sizes L (cid:39) Acknowledgments
The project was funded by the European Union and the Free State of Saxony. This work hasbeen partially supported by the DFG (Grant No. JA 483/31-1), the Leipzig Graduate School“BuildMoNa”, and the Deutsch-Franz¨osische Hochschule DFH-UFA (Grant No. CDFA-02-07).The authors gratefully acknowledge the computing time provided by the John von NeumannInstitute for Computing (NIC) on the supercomputer JUROPA at J¨ulich Supercomputing Centre(JSC).
References [1] Binder K and Kalos M H 1980
J. Stat. Phys. J. Stat. Phys.
Europhys. Lett. J. Stat. Phys.
Physica A
Europhys. Lett. Phys. Rev. E J. Phys: Conf. Ser.
J. Chem. Phys.
J. Chem. Phys.
Phys. Rev. E Physica A
Phys. Rev. E Phys. Rev. Nuovo Cimento (Suppl.) Phys. Rev. J. Phys. A: Math.Theor. J. Phys. A: Math. Gen. Phys. Lett. B
Phys. Rev. Lett. Int. J. Mod. Phys. C Physica A
Comput. Phys. Commun.