Microcanonical Determination of Interface Tension

Microcanonical Determination of the Interface Tension of Flat and Curved Interfaces from Monte Carlo Simulations

Abstract

The investigation of phase coexistence in systems with multi-component order parameters in finite systems is discussed, and as a generic example, Monte Carlo simulations of the two-dimensional q-state Potts model (q=30) on square lattices () are presented. It is shown that the microcanonical ensemble is well-suited both to find the precise location of the first order phase transition and to obtain an accurate estimate for the interfacial free energy between coexisting ordered and disordered phases. For this purpose, a microcanonical version of the heatbath algorithm is implemented. The finite size behaviour of the loop in the curve describing the inverse temperature versus energy density is discussed, emphasizing that the extrema do not have the meaning of van der Waals-like ”spinodal points” separating metastable from unstable states, but rather describe the onset of heterophase states: droplet/bubble evaporation/condensation transitions. Thus all parts of these loops, including the parts that correspond to a negative specific heat, describe phase coexistence in full thermal equilibrium. However, the estimates for the curvature-dependent interface tension of the droplets and bubbles suffer from unexpected and unexplained large finite size effects which need further study.

pacs:
05.10.Ln,64.60.De,64.60.qe,64.75.Gh,68.03.Cd

1 Introduction

In the theory of first order phase transitions, a quantity of major interest is the interface tension between phases in metastable coexistence. In particular, in the study of nucleation phenomena one faces the problem to determine the curvature dependence of interface free energy/entropy between a nucleating droplet or bubble of the emerging stable and its surrounding metastable phase [1, 2, 3, 4, 5, 6, 7]. While an exact understanding of metastability in infinite systems is currently still lacking from the point of view of rigorous statistical mechanics [8, 9], one can nevertheless try to obtain valuable information from the study of phase separation in finite systems [10, 11, 12]. As analytical calculations are quite difficult for any nontrivial Hamiltonian, one resorts to simulations. For systems like Ising type spin models (lattice gases) [13], binary mixtures [14] and simple fluids [15, 14], that have been studied by simulation, a convenient scalar order parameter like the total magnetization or particle number, which is extensive and even additive under partitions of the total system volume into subvolumes, is available, and whose density serves to distinguish between the different phases. Phase coexistence in an equilibrated finite system is characterized by the identity of the corresponding conjugate intensive quantity, whose physical meaning is that of an applied magnetic field or chemical potential. Thus one can extract the order parameter bulk densities of the coexisting phases from analyzing free energies obtained from Monte Carlo [10, 11, 15, 14] or molecular dynamics [16, 17, 18] simulations. In particular, armed with an additive order parameter, one can determine the equimolar volumes of the coexisting phases by evaluating the condition of vanishing adsorption of this order parameter for the corresponding choice of dividing surface. From such data it is straightforward to compute the interface tension [15, 14].

However, problems arise in cases for which a scalar additive order parameter is not available. On the one hand, the different phases may be characterized by different values of a multi-component rather than a simple scalar quantity (e.g. the formal vector of particle numbers of different species in multi-component fluids). On the other hand, there are cases where an explicit microscopic expression for an order parameter is frankly not known like e.g. in the study of protein folding.

An extensive parameter that is always on hand and whose density – except for somewhat degenerate cases like hard spheres – may be utilized to distinguish different coexisting metastable and stable phases is the energy . But unlike in the case of magnetization or particle number, the notion of an “equi-energetic” surface seems to be completely counter-intuitive, as attributing zero energy to the interface is in conflict with the fundamental physical principle that interfaces are regions with energy densities that are usually considerably larger than bulk ones. Thus, it is not clear how a reasonable measure of the corresponding subvolumes of the phases can be obtained from a microcanonical analysis based solely on monitoring the energy density. The purpose of the present paper is to show how this can be done, and that energy may still be a “good” parameter to determine interface tensions.

2 Review of the Grand Canonical Route to the Interface Tension

We start by reviewing the traditional Gibbs dividing surface approach to the description of phase coexistence of two phases labelled and in a fluid [19, 20, 3, 21, 22]. Suppose first that the fluid has only a single chemical component. In spite of the fact that the problem of how to rigorously define coexistence of a droplet of phase in (unstable) equilibrium with the surrounding phase has not been solved up to date for the case of an infinite system, it makes sense to consider such a situation in a large but finite system, in which such an equilibrium may actually be stabilized by imposing appropriate thermodynamic boundary conditions. Phase separation into two connected regions is then usually detected by observing regions of different density by monitoring the average inhomogeneous density profile . Following Gibbs one may then choose an arbitrary dividing surface, defined as a levelling surface of zero thickness normal to the density gradient field, thus splitting the total volume into subvolumes

 V=Vα+Vβ. (1)

Once this separation is agreed upon, any other extensive observable can be split into homogeneous and so-called excess contributions as follows. If assumes homogeneous equilibrium densities in the phases , we set

 M=Mα+Mβ+Mx, (2)

where

 Mα≡Vαmα,Mβ≡Vβmβ. (3)

In a similar way, we construct the excess energy , entropy , Helmholtz free energy , grand potential , particle number , and so on from their homogeneous densities and the division (1). In a box of volume , planar interfaces will usually form parallel to one pair of limiting walls. The excess quantities defined above will then generally depend on the chosen position of the dividing surface with area . A notable exception is , as can be understood from the fact that the corresponding homogeneous densities are the negative pressures of both phases, which in case of a planar interface must agree by simple stability arguments. Thus, for a planar phase separation geometry, the interface tension

 σ=Ωx/A (4)

is well defined, regardless of any parallel shift of the dividing surface. In contrast, the adsorption

 Γ=Nx/A (5)

changes with such a parallel shift of . In turn, the condition of vanishing adsorption then uniquely fixes the position of and leads to an intuitively appealing definition of the “actual” position of the interface, known as the equimolar surface.

While it is straightforward to generalize the definition of the position of such an equimolar dividing surface from the planar to that of a spherical or otherwise curved case, the definition of the corresponding interface tension now requires considerably more care. For a spherical interface in let us agree to use to label the phase inside the spherical volume and to label the surrounding one outside the sphere. To stabilize a curved interface, the pressures on both sides of the interface must necessarily be different. In classical macroscopic physics this fact is encoded in the Laplace-Young (LY) equation

 Δp=pα−pβ=2σR, (6)

which was derived from a mechanical analysis of the surface tension at the beginning of the 19th century. When promoting the statistical mechanics definition (4) from the planar case to that of curved interfaces, since for , the grand potentials densities of both phases will also disagree and so

 Ωx=Ωx(R)=Vω−Vα(R)ωα−Vβ(R)ωβ (7)

is bound to pick up a dependence on the radius. Thus, for a curved interface the interface tension appearing in (6) will itself be -dependent. At this point, however, it is crucial to realize that within the Gibbs dividing surface approach the choice of the radius , being a purely theoretical construct, is in principle arbitrary, such that the classical relation (6) can only hold for a distinguished value of . Nevertheless, physically observable quantities should not depend on the position of the artificially introduced dividing surface. Indeed, differentiating the definition of the excess grand potential in d dimensions with respect to while keeping all other variables fixed, which is called a notional derivative and indicated by a bracket notation , one arrives at the generalized LY equation [19, 20]

 Δp=(d−1)σ(R)R+[dσ(R)dR]. (8)

The classical LY equation (6) may only be recovered from this equation for the special choice of radius , for which

 [dσ(R)dR]R=Rs≡0. (9)

To account for this fact, the corresponding dividing surface is known as the surface of tension. Thus, is a stationary point of under a notional variation of . In fact, it is not hard to see [19, 20] that is indeed a minimum of , which can explicitly be deduced from the general form of in or dimensions

 σ(R)σ(Rs)=1+⎧⎪ ⎪⎨⎪ ⎪⎩12(R−RsR)2RRs,d=213(R−RsR)2(Rs+2R)Rs,d=3, (10)

according to which is universally determined from knowledge of and . The difference between radius of the equimolar surface and , which has become famous under the name Tolman length, is known to be of molecular sizes. Of course, in principle one could also compute physical observables from e.g. the equimolar interface tension or any other choice of radius , since all physical information is encoded in any such choice. However, the choice is particularly convenient, as the condition (9) allows to condense many formulas to a considerable more compact and manageable form. To some extent this is also true for the choice , but the definition explicitly makes use of the particle number as an extensive parameter. For a multi-component fluid of different chemical components, which we may label by an index , this complicates manners in a considerable way. In fact, the particle number , adsorption , associated chemical potential are then both promoted from scalar quantities to formal vectors , , , where and so on. One now must deal with -component averaged density profiles . In general, such systems have complex phase diagrams, and if we concentrate on a particular transition, the levelling surfaces of different density components may yield different equimolar surfaces for each chosen component, relative to which the remaining components form inhomogeneous adsorption layers. In principle it may still be possible to eliminate clumsy adsorption terms from formulas to some extent by an ingenious choice of dividing surface, the so-called Koenig dividing surface [23]. However, to fulfil its defining condition, all components of the adsorption have to be balanced simultaneously with the multiple components of the associated chemical potential, which leads to profound numerical difficulties in the evaluation of simulation results.

To illustrate these calamities, let us review the practical steps to calculate the interface tension from Monte Carlo free energy simulation data for a one component fluid and compare them to the expected effort for a multi-component one. We choose a cubic simulation box of size with periodic boundary conditions. Carrying out grand canonical Monte Carlo simulations at some fixed temperature chosen somewhat lower than the critical temperature , one first has to determine the coexistence chemical potential . Practically, this is done by implementing the equal–weight rule [24] numerically by performing a histogram re-weighting [25] to the simulated grand canonical probability distribution of particle numbers and extrapolating the result to the thermodynamic limit . For a multiple component fluid, these steps are already quite involved, as the approximately Gaussian-shaped peaks of are defined on a multidimensional space of variables and in presence of possible low temperature phases one has to deal with of them in order to pin down the individual components of the vectorial chemical potential to their coexistence values .

Nevertheless, suppose that this task has been carried out successfully. In the single component case, one next needs to resolve the detailed structure of by e.g. Wang-Landau sampling, eliminating residual inaccuracies by a weighted Monte Carlo production run. In this way, one obtains a dimensionless excess free energy density

 −β^f(L)(T0,ρ)≡(1/V)lnPT0Vμ0(N), (11)

which depends on the scalar density (we have dropped a purely -dependent normalization part which is irrelevant for what follows; the usual thermodynamic notation should not lead to any confusion with the label for one of the two phases). At a first order phase transition, the finite size potential has a double-welled shape with a flat central plateau, and an analysis of its fine details reveals several distinct ranges of the total density between its two minima, for which one may observe phase separation in planar, cylindrical (in ) and spherical shapes induced by the imposed periodic boundary conditions. These regions are also detectable in the derived finite size canonical “excess” chemical potential

 ^μ(L)(T0,ρ)=⎛⎜⎝∂^f(L)(T0,ρ)∂ρ⎞⎟⎠T0. (12)

From the distorted “Van der Waals loop” shape of , one can infer that for small enough the equation generally allows for at least three roots corresponding to two bulk densities of coexisting phases at total density . This allows to introduce the total dimensionless grand potential density

 ω(L)(T0,μ)=ˆf(L)(T0,ρ(T0,μ))−μρ(L)(T0,μ) (13)

and the grand potential densities

 ω(L)α(T0,μ) = ˆf(L)(T0,ρα(T0,μ))−μρ(L)α(T0,μ), (14) ω(L)β(T0,μ) = ˆf(L)(T0,ρβ(T0,μ))−μρ(L)β(T0,μ), (15)

from which it is straightforward to compute for a given volume partition , and thus, using (4), the notionally -dependent interface tension . In detail, and are found by numerical minimization of with respect to , while can be calculated from the spherical equimolar volume determined by the lever rule

 Vα(Re)V=ρ−ραρβ−ρα,Vβ(Re)V=ρβ−ρρβ−ρα, (16)

for in and in , respectively. In passing we note that since , is exclusively determined by (cf. [15, 14]).

All this is very fine – however, in a multi-component setting the numerical effort to carry out these steps successfully is again forbidding. What is needed is an intrinsically scalar approach.

3 Microcanonical Approach

As mentioned above, a basic extensive scalar quantity that is always available is the total energy of the system. It is by now well known that while phase transitions are only well defined in an infinite system in the strict sense, they can nevertheless be studied conveniently by analyzing the so-called “convex intruder” in the microcanonical entropy of a finite version of the system [26, 27].

By its very definition, such a convex intruder in gives rise to a corresponding anomaly of the microcanonical inverse temperature and a possibly negative branch of the accompanying microcanonical specific heat with quasi-singularities at the intruder boundaries (for a schematic picture, see Figure 1 below). In a finite (or at least in some sense “small” [26]) system such an observation does not contradict thermodynamic stability requirements [28]. One now observes that up to a sign, the overall appearance of is quite similar to that of the excess chemical potential discussed above. Indeed, parallels between the former constructions and the ones carried below are not accidental.

Recording the energy probability distribution at fixed for different values of , we can locate the coexistence chemical potential by applying the equal weight rule [29, 24, 30] to the two separate approximately Gaussian-shaped peaks, as will be discussed below for the case of the Potts model. In any case, this may be feasible even in case of a multi-component order parameter, since the domain of is the scalar quantity . For the rest of the simulation, the parameters and remain fixed. It is then convenient to formally replace the energy by a grand-canonical version

 E−μ0N → E (17)

in the same way as one may formally replace the original Hamiltonian by a grand-canonical one in computing the grand-canonical partition function. For our purposes, the “background” homogeneous chemical potential can thus be eliminated from the list of thermodynamic variables. With fixed, we now record the number of states

 g(E,V)≡eS(E,V) (18)

with grand canonical energy in a flat histogram simulation followed by a weighted Monte Carlo production run. Once is known, we arrive at the grand canonical partition function in the form

 Z(β,V)=∑Ee−βE+S(E,V)=∑ee−V(βe+s(L)(e)), (19)

where we introduced the energy and entropy densities and , valid for arbitrary . Suppose that at a particular inverse temperatures the sum is dominated by its largest summand in the limit of large system size. The corresponding energy density is determined by the equation

 β≡(∂s(L)(e)∂e)V∣∣∣e=e(L)(β)=β(L)(e)|e=e(L)(β), (20)

which expresses the equality of inverse canonical temperature and microcanonical temperature for this special value of . One may then approximate by

 βω(L)(β,μ0)≈βe(L)(β)−s(L)(e(L)(β)), (21)

i.e. as a Legendre transform. This equation is quite similar to (13), playing the role of , that of and that of . Similar to the strategy employed above for the chemical potential, for a certain range of inverse temperatures the equation has (at least) three different solutions corresponding to the total energy density and those of the two coexisting phases . If we now let approach this interval around , the presence of the convex intruder in makes the above discrete saddle point approximation break down since not one but at least three “saddle points” are found, which in the non-degenerate case of three roots correspond to the dimensionless grand potential densities

 βω(L)(β) = βe−s(L)(e), (22) βω(L)α(β) = βeα−s(L)(eα), (23) βω(L)β(β) = βeβ−s(L)(eβ). (24)

At this stage, we have reached our goal announced above, as these quantities, once they are determined, allow to compute the excess grand potential and thus the interface tension from formula (7). By minimization of w.r.t  we obtain the radius of the surface of tension and the corresponding surface tension . A calculation of is, of course, not feasible in this approach.

4 The 2d Potts Model

We illustrate our strategy taking the example of the nearest neighbour Potts model in . The -state Potts model [31], whose Hamiltonian on a simple cubic 2d lattice of sites in zero external field is given by

 H[{s(x)}]=∑⟨xy⟩[1−δs(x),s(y)], (25)

where is ideally suited for this purpose for a number of reasons. (i) First of all, for the model undergoes a temperature-driven first order phase transition. (ii) Regarding this transition, a wealth of rigorous results [31, 32, 33, 34, 35] is available in the literature which can serve to benchmark our simulation results. For the first order phase transition temperature of a bulk system, one has the exact analytic expression

 1/T0=β0=ln(1+√q). (26)

In [36, 34] it was reported that, as was expected from general arguments, the inverse temperature at which the ratio of the two weights of the thermal energy probability distribution is just , agrees with the exact bulk value up to exponentially small corrections. Thus, serves as a convenient definition of a finite size transition temperature. Other rigorous results include the latent heat per volume [33], the limiting internal energy densities at and the difference of specific heats [32]. Furthermore, the reduced interface tension (i.e. the interface free energy density at multiplied by ) between the disordered and one of the ordered phases along the square lattice direction was rigorously determined [35] to be

 2σo/d=4∞∑n=0ln1+wn1−wn, (27)

where

 wn:=(√2cosh(n+1/2)π22v)−1 (28)

with

 v:=ln[12(√√q+2+√√q−2)]. (29)

Equally important for us is the fact that there is even a rigorous calculation of the full anisotropic interface tension available [37]. However, since these calculations are too involved to be reproduced here, we content ourselves with noting that at the resulting anisotropy for calculated from the formulae in [37] is vanishingly small. This happenstance is a very important prerequisite for any attempt to apply our evaluation strategy for the interface tension, which rests on a presupposed spherical symmetry of bubbles and droplets.

(iii) The state Potts model’s order parameter is not scalar, but has dimension . Since this is a central issue in the present context, let us briefly review its nuts and bolts. Guided by physical intuition, a scalar “order parameter” could be defined by the following reasoning [31, 29]. Let denote the number of spins of a given microstate with value , where . Let . Then

 M:=q⟨Nmax⟩−Nq−1. (30)

Obviously is confined to values , and for complete disorder, while for any perfectly ordered domain. In the Ising case , indeed corresponds to the modulus of the magnetization density of the system. Thus, if we break up the system volume into subvolumes and add up their different ’s, generally , i.e.  is not additive between subsystems.

On the other hand, in [38] Zia and Wallace construct a full -component order parameter. They introduce unit vectors in , such that the following relations are satisfied:

 e(a)e(b)=qδab−1q−1. (31)

Any set of such vectors defines a generalized tetrahedron in , the vectors pointing from the center to each corner. Trivially, the vectors cannot be linearly independent. Instead, they satisfy the geometrically evident sum rule

 ∑ae(a)=0. (32)

With the one-to-one correspondence understood, one associates a component order parameter

 M:=∑x∈Γe(x) ∈Rq−1 (33)

with each given microstate, which we will call the magnetization. In terms of the occupation numbers of the Potts spin states

 M=q∑a=1N(a)e(a). (34)

If one multiplies (34) with any of the unit vectors , it is easy to see that

 M(b):=e(b)M=qN(b)−Nq−1. (35)

Clearly, agrees with as defined in (30) provided . This clarifies the role of as well as the low temperature domain structure of the model. Namely, suppose that . Then, we can rewrite (34) as

 M=∑a≠b(N(b)−N(a))≥0(−e(a)). (36)

The domains are thus geometrically represented by convex cones enclosed by the set of vectors , and (35) gives the projection of onto the average direction , the symmetric group of permutations of numbers acts as the underlying symmetry.

At small values of , several very small fluctuation “clusters” of the same or competing spin values may coexist, causing the system to jump randomly from one domain cone to another. On the other hand, the parameter is always additive by construction, but does depend on the chosen direction in -space. Only for values of larger than a certain threshold do we find agreement of the order parameters computed from (30) and (35), since then the direction along which the projection from to occurs is uniquely determined by the value of the majority occupation number. A cluster decomposition performed during the course of a simulation can provide information in analyzing these fluctuations in order parameter topology.

Thus, the parameter (30) is only additive for magnetizations and both belong to a single common domain , and thus cannot be used meaningfully as a parameter in a Gibbs dividing surface construction. On the other hand, sampling the free energy as a function of the full order parameter is out of the question. In other words, we are exactly in the situation outlined earlier, and thus embark on a microcanonical strategy instead. Well, not quite. Actually, the situation for the Potts model is not as complicated as the one outlined in the introductory section. This is largely due to the fact that there is no need to determine a coexistence “chemical potential”, which would correspond to an external vectorial magnetic field. As the different -spins of the Potts model are all coupled in the same way, this external field is fixed to be exactly zero.

5 The Microcanonical Heat Bath Algorithm

Any successful design for a Monte Carlo algorithm devoted to the study of phase separation at phase transitions of pronounced first order character must address the phenomenon of exponentially diverging autocorrelation times

 τ(L,β)∼exp(2βσ∞Ld−1) (37)

accompanying first order phase transitions which known as hyper-critical slowing down [39] (in (37) a -dimensional cubic box of volume with periodic boundary conditions is assumed). A Wang-Landau type of algorithm is capable of overcoming large entropy or free energy barriers separating different stable or metastable phases [40] and is in principle straightforward to implement. In its original microcanonical version, the algorithm directly yields the density of states (or, for a discrete system, rather number of states)

 g(E)=∑νδ(Eν−E) (38)

of a system having microstates with energies (we are somewhat casual about the use of the Dirac delta function for simplicity, and we have put for the same reason) with high precision, which is related to the microcanonical entropy by . To control possible residual errors, one may thus determine an approximate microcanonical density of states by Wang-Landau simulations and perform subsequent biased Monte Carlo production runs with statistical weights based on this approximate density of states. In fact, knowledge of conveniently allows to determine a multitude of other -dependent observables at various temperatures simultaneously with high precision, as is explained in more detail below.

It remains to construct a suitable move set for a microcanonical Wang-Landau simulation scheme. Single -spin updates are simple to implement and may be a reasonable choice for Ising systems far from criticality, but are quite inefficient in exploring the regions of phase space of the large Potts model which are of interest for studying phase separation, namely those configurations, in which a single ordered domain of Potts spins of, say, value coexists with a disordered background. Indeed, suppose that during the course of the simulation, our random walk arrives at a particular configuration, in which almost all Potts spins agree with this condition, while just a few, say, are yet disordered. In such a microstate, chances are only that of the “missing” sites is indeed drawn and its spin be assigned the “right” value in creating the next trial configuration, thus matching the surrounding domain.

Within the canonical ensemble, it is well known that the heat bath algorithm [41, 42] is superior to the standard Metropolis scheme for high Potts models. In detail, let denote the total energy of the system in the microstate . Choose a random site and let denote the value of the Potts spin variable at this particular site. Furthermore, let denote the nearest neighbour sites of . Defining the local energy at site by

 Eq≡E(local)(q,x):=−z∑i=1δq,s(yi), (39)

one can split

 E(μ)=E(local)(qμ,x)+E(μ)rest (40)

and define a set of heat bath probabilities

 pq:=e−βEq∑q−1n=0e−βEn, (41)

which are manifestly independent of the value of the initial central Potts spin . The canonical heat bath algorithm amounts to choosing a new value for this spin with probability in every step. It is easy to see that the resulting algorithm satisfies detailed balance as well as ergodicity. In terms of generation and acceptance probabilities, we have and , i.e. the stochastic character only enters in the generation of configurations, which are, once generated, always accepted.

It is straightforward to translate these ideas from the canonical to the microcanonical setting. To illustrate the correspondence, let us denote the canonical Boltzmann weights by

 πν:=e−βE(ν)Z(β). (42)

Then the canonical heat bath probabilities can quite trivially be rewritten as

 pqν:=e−β(Eq+Erest)∑q−1n=0e−β(Eq+Erest)=πqν∑q−1n=0πn. (43)

Now, to translate the algorithm to the microcanonical ensemble, we simply replace . In terms of the microcanonical entropy , we can rewrite the above probabilities as

 pqν=1/g(Eν)∑q−1n=01/g(En)=e−S(Eν)∑q−1n=0e−S(En). (44)

Once we have determined the microcanonical entropy , it is in principle straightforward to obtain the free energy as a function of and some other (preferable scalar) observable , where denotes e.g. the magnetization , the projection of the order parameter along an arbitrary fixed internal direction , the size of the largest geometric or Swendsen-Wang cluster, all simultaneously computed for different temperatures from the underlying microstates visited during the course of a single microcanonical biased Monte Carlo simulation.

To sketch this procedure, we consider the constrained microcanonical density of states , which is formally written as

 g(E,o)=∑{s(x)}δ(o−O[{s(x)}])δ(E−H[{s(x)}]). (45)

The corresponding (conditional) probability to find the value of at total energy is

 P(o|E)=g(E,o)∑o′g(E,o′)=g(E,o)g(E)∝h(o,E), (46)

where denotes a two-dimensional histogram recorded during the course of the simulation. But, according to the rules of conditional probabilities, this precisely implies that the canonical probability to find the value of at inverse temperature is

 P(o|β)=∑EP(o|E)P(E|β)∝∑Eh(o,E)P(E|β). (47)

To obtain this probability at any given temperature from a microcanonical Monte Carlo simulation biased by the predetermined density of states , we thus only need to reweight the recorded histograms by the known function . The desired constrained free energy density

 f(β,o)=−1Nβln∑{s(x)}δ(o−O[{s(x)}])e−βH[{s(x)}] (48)

can be recast in a similar way, since

 Missing \left or extra \right (49) = ∑Ee−βE∑{s(x)}δ(o−O[{s(x)}])δ(E−H[{s(x)}]) = Z(β)∑EP(E|β)g(E,o)=Z(β)P(o|β),

and is thus (up to an unimportant constant) given by

 f(β,o)=−1NβlnP(o|β). (50)

6 Microcanonical Results

We have conducted a series of Landau-Wang simulations followed by weighted Monte Carlo production runs of a square lattice Potts model of size to determine the density of states and thus the entropy density and the microcanonical temperature . Since we are interested in the development of phase separation in small to finite system sizes, we carried out simulations for linear sizes . For the Potts model, the resulting signs of phase separation were not observed to be very pronounced. However, increasing to clearly revealed the expected convex intruder. But even then, from merely looking at the entropy density it is virtually impossible to detect the delicate features appearing at finite system sizes that we are interested in (cf. Figure 4).

However, numerically taking the derivative of with respect to , we obtain the microcanonical inverse temperature which provides a detailed view of the delicate substructures hidden in (cf. Figure 5). Inspection of the resulting curves gives a first hint on the quality of our simulation data, as one should take into account that numerically differentiating potentially noisy data should greatly magnify any statistical irregularities and errors in such data.

To resolve the convex intruder in the original entropy data , it turns out to be convenient to consider the auxiliary dimensionless function [43]

 Λ(L)(E,β)/N=λ(L)(e,β):=βe−s(L)(e). (51)

Of course, coincides up to a constant with the logarithm of the canonical energy probability function at inverse temperature . However, for our present purposes we may regard as a finite size “Landau potential”, i.e. an incomplete Legendre transform of the microcanonical entropy, and compare its features to those of canonical Landau potentials (cf. [44, 45]). Thus, let us tune the parameter to values near the bulk inverse transition temperature and analyze the resulting shape of as a function of . As expected, for such temperatures resembles a somewhat distorted double-well shape (cf. Figure 6) with two pronounced minima at energy densities and (the subscripts “c” and “v” correspond to “condensed” and “vapor”) separated by a large “Landau free energy barrier” with a “flat” central plateau of practically constant and vanishing or at least quite small slope. We identify and with the equilibrium energy densities of the bulk “condensed” (ordered) and “vapor” (disordered) phases, while the thermodynamics of their possible coexistence configurations in encoded in the features of the potential well between them. The flat central region of the potential, which signals phase separated configurations with a slab-like interface geometry [44], is, of course, also reflected in the central linear section of at the level of (cf. Figure 5).

Apparently, Figure 6 also illustrates the dilemma of defining a “proper” finite size transition temperature for a system with a highly degenerate low energy domain structure. On the one hand, one could naively try to adjust to such a value that both minima of are of equal height. This choice precisely corresponds to the “equal height rule” for . However, at such a temperature, one observes a noticeable slope in the central “flat region” of , which signals that phase coexistence is not well established. In a plot of the quasi-Gaussian function this and other delicate features outside the peak regions do not give themselves away to the naked eye, since they are exponentially suppressed. On the other hand, choosing the inverse temperature to agree with the ratio-of-weights temperature with , which, as discussed above, converges exponentially fast to the exactly known inverse bulk transition temperature [33]

 β0=ln(1+√q)\lx@stackrel(q=30)=1.86829, (52)

the flat central region of is found within numerical precision to be horizontal, i.e. with vanishing slope, but now one notices a pronounced difference in height between the two minima at energy densities , which diminishes with growing system size. Numerically, this height difference is seen to approach

 λ(L)(ev,β0)−λ(L)(ec,β0)≈N−1lnq. (53)

Physically we can interpret these findings as follows. Suppose that precisely at the inverse transition temperature the system initially starts out in an ordered equilibrium state. Then, the probability to generate a fluctuation yielding a mixed state where half of the available volume is turned into a disordered state separated from the ordered part by a straight interfacial line, should differ from the “inverse” probability to produce the same state starting from the disordered state by a factor of , since in the latter case possible ordered states are available, while only one disordered configuration can be formed in the former one. Taking the logarithm of the fraction of these probabilities then produces (53). This completes the picture, since, as could have been anticipated from the practically perfect Gaussian nature of the two peaks in , the ratio-of-weights and ratio-of-heights temperatures are found to numerically agree for practical purposes. Without going into the details we also note that the equal height temperature can also be shown to coincide with the inverse temperature found by imposing a microcanonical Maxwell construction, i.e. chosing the value of for which suitable defined areas obtained from integrating between and coincide [46].

At this point, a few additional comments concerning the true nature of the above “Landau potential”, the non-monotonous behaviour of and the resulting appearance of branches of “negative specific heat” are in order. In fact, in the thermodynamic limit, indeed decreases with up to , stays constant at up to , and then decreases further, as it should be: no trace of any metastable states (or even “unstable states”) is left in the curve. Of course, this must be so: apart from statistical errors, Monte Carlo simulations yield the equilibrium statistical mechanics of any such model Hamiltonian exactly, and metastable or unstable curves cannot be the output of exact calculations in the framework of equilibrium statistical mechanics. So for the minimum position of moves towards (and its depth vanishes); similarly, the maximum position moves towards (and its height, relative to , vanishes as well). It is interesting to recall the physical significance of these extrema: for the finite system is still homogeneous, and the minimum is the signature of the first appearance of a “bubble” of the disordered phase within the otherwise homogeneously ordered phase (cf. Figure 3), while the maximum is the signature of the first appearance of a “droplet” of the ordered phase within the otherwise homogeneous disordered phase (Figure 2). As long as such “heterophase fluctuations” are absent, finite size effects are small in the curve shown in Figure 5; the strong finite size effects in between and are due to interfacial contributions to the “Landau potential” (Figure 6), which are of relative order in Figure 6. Thus, the branch of negative “specific heat” resulting from Figure 5 in the region where is an increasing function of is not at all an unphysical result, but simply reflects the importance of interfacial free energies in finite microcanonical systems. In addition, our use of the nomenclature “Landau potential” merely refers to the incomplete character of the Legendre transform (51), but should not mislead the reader to confuse this potential with “Landau potentials” of similar appearance as they are constructed in mean-field theory. In fact, our potential , whose information content is, after all, identical to that of the full microcanonical entropy density , describes the thermodynamics of two-phase coexistence in an inhomogeneous finite system without any approximation, and thus is conceptually quite different from a mean-field potential, which is constructed under the implicit constraint that the system is in a homogeneous phase throughout, whereas such states are thermodynamically unstable in reality.

We can gain confidence in the overall correctness and general quality of our data by comparing the exact value of the reduced planar interface tension as calculated from (27) to the one obtained by a finite size extrapolation of our data. In fact, as there are two minima of at energy densities whose values differ by as discussed above, there are two corresponding sets of data , corresponding to the difference between the central barrier taken at some energy density and the left and right minima , , respectively (see Figure 6). A standard finite size scaling extrapolation of these data to in the form

 β0σc,v(L)=β0σc,v(∞)+const/L, (54)

which is displayed in Figure 7, gives the two values and for the left and right difference, respectively, whose average differs by less than from the exact value computed from Formula (27).

From (4), (9) and (22)-(24) it is now straightforward to compute and for any prescribed total energy density . Note, however, that these formulae are only valid for a spherical phase separation geometry. The corresponding approximate density regions within which one may expect states which on average resemble spherical droplets and bubbles to dominate in the sampled microscopic system configurations may be found by visually inspecting the slopes of and in Figures 5 and 6, respectively, and cross-checking these ranges by examining corresponding snapshots taken during the course of the simulation. The total energy density regions for which we may expect the appearance of spherical droplets and bubbles are indicated in colour in Figures 5 and 6.

Our results for the interface tension at the surface of tension are gathered in Figures 8 and 9. In correctly interpreting these results, it is quite important to understand that they have been obtained from (4), (9) and (22)-(24) under the assumption of spherical geometry. The corresponding ranges of inverse radii for which one can expect this assumption to be valid have been marked in colour in these figures. Outside of these ranges, the data do not accurately describe a physical interface tension, but merely serve as a guide to the eye.

Looking at these results, one instantly notices the large finite size effects, manifesting themselves in the considerable offsets between the consecutive considered -values, which strikingly fail to collapse onto a common “master curve”. Currently, we find it difficult to understand the origin of this behaviour. For the planar interface tension, which is described by our data quite accurately as discussed above, strong but regular finite size effects are indeed expected. They may be heuristically understood in terms of the -dependent truncation of the wave vector spectrum of capillary waves running parallel to the interface. However, for a spherical bubble or droplet, identical radii should yield identical values of the interface tension, once the surrounding box has been chosen large enough to kill finite correlation length effects, which are however expected to be vanishingly small for a strong first order phase transition.

At the moment we do not have a clear explanation of these strong finite size effects which prevent us from a further meaningful analysis of the curvature dependence of . In similar approaches to study the interface tension of curved interfaces in a truncated 3d Lennard-Jones fluid [47] and a 3d fcc lattice gas model [48], we also find certain finite size effects, but they are much less pronounced than those of the present case. However, we have also observed finite size effects of comparable size in computing the interface tension from canonical simulations of a 2d Ising model. Thus, we believe that the large magnitude of the finite size effects has nothing to do with our microcanonical approach, but is rather related to the fact that both systems are two-dimensional. With growing linear system size , the gaps between consecutive -curves recorded in Figures 9 and 9 obviously diminish, so these curves are expected to eventually collapse onto a single “master curve” for large and . On the other hand, for we report that our Monte Carlo production runs failed to be sufficiently ergodic, indicating large residual entropic barriers with respect to “hidden” observables beyond our one-dimensional energy-based sampling, while at the same time the barriers observed in had already risen to a value of . To study much larger systems would thus require to overcome these additional hidden barriers, presumably by employing much more elaborate sampling techniques than the ones we are using here (see e.g. [49] or [50] for promising approaches). However, in our current work we are not interested in exceedingly large droplet sizes and their accompanying huge entropy and free energy barriers. Rather, out intention is to focus on the behavior of droplet and bubbles of moderate size, as this is the only regime that is of practical relevance for nucleation related questions. In any case, the origin and nature of the encountered finite size effects must currently be left to further study.

7 Conclusions

In this paper, we have addressed the investigation of phase coexistence of systems with a more-component parameter in the context of computer simulations, which necessarily involve systems of finite size. Such simulations of phase coexistence often are done with the motivation to extract information on the interfacial tension of flat and curved interfaces. While for systems with a scalar (i.e. one-component) order parameter this problem is normally considered in the grand-canonical and canonical ensemble of statistical mechanics, we have given a concise discussion of this approach to show that its extension to the multi-component case is formally possible but practically unfeasible. We then have presented, as an alternative, a microcanonical approach based on the number of states of energy for a system having a finite volume . In the entropy versus energy curve for the finite system there is a convex intruder (Figure 1), and the idea we follow in the present paper is to carefully analyze this intruder as a function of system size, in order to extract information on interfacial tensions. We exemplify our approach for the two-dimensional -state Potts model with a large number of states (), proposing also an extension of the heat bath algorithm from the canonical to the microcanonical ensemble. From these simulations we obtain very precise information on and also an effective potential ) per lattice site, being the energy density and the inverse temperature where in the thermodynamic limit () the first-order transition from the ordered phase to the disordered phase occurs (Figure 6). Also the derivative is obtained with meaningful accuracy (Figure 5). We have shown that the loop in such vs.  curves has nothing to do with the “van der Waals-like” loop of mean-field theories: in the latter, such loops describe a path of homogeneous states connecting the two phases between which the transition occurs; in reality, our loops (Figure 5) reflect two phase coexistence in finite systems, all parts of the loop describe full stable thermal equilibrium; any interpretation in terms of metastable or unstable states would be completely misleading. There is nothing mysterious about the “negative specific heat” that often is attributed to such loops - the whole loop just reflects interfacial effects, just as the ”hump” in between the two minima of the “Landau potential” in Figure 6; all these features disappear proportional to in the limit , and the correct horizontal parts in between and remain, as it should be. Thus, one should not be mislead by mean-field concepts when discussing first-order phase transitions in finite system in the microcanonical ensemble.

We have found that in the flat region in the center of Figure 6 the data allow an accurate estimation of the interfacial tension of flat interfaces between ordered and disordered phase (Figure 7), although also in this case finite size effects are clearly rather pronounced, and an extrapolation to is mandatory. However, the analysis of the ascending parts of in Figure 6 in terms of the radius-dependent interface tension of droplets (Figure 8) and bubbles (Figure 9) is more subtle: again huge finite size effects occur, and it is not possible at fixed radius to extrapolate to , because due to the droplet (bubble) evaporation/condensation transitions droplets at fixed radius are only stable in a rather restricted range of . While naively one could expect that different choices of yield mutually compatible results for , as approximately happens for one-component systems in dimensions, this is not the case here. Of course, our analysis does not explicitly consider the fact that at a given value of and the corresponding average value of (Figure 5) in the two-phase coexistence region at a given value of the droplet (or bubble) is strongly fluctuating both with respect to its size and its shape (Figures 2 and 3). We assume the shape of the droplet or bubble to be spherical, otherwise the information recorded does not suffice to extract . Future work along such lines must analyze this problem of droplet (bubble) fluctuations more closely, possibly by recording additional observables related to the droplet or bubble, to allow estimation of within reasonable error limits. In fact, in the fluctuations are found to be indeed much less pronounced, and – at last in the one-component case – meaningful results for are accessible [47, 48].

A. Tröster acknowledges support by the Austrian Science Fund (FWF): P22087-N16 and is grateful to Prof. J. Henderson for enlightening emails and to D. Reith for valuable computer assistance. K. Binder received support from the Deutsche Forschungsgemeinschaft (DFG) under grant No SPP 1296/BI 314/19-2.

8 References

References

1. D. W. Oxtoby, Journal of Physics: Condensed Matter 4, 7627 (1992).
2. J. L. Katz, Pure Appl. Chem. 64, 1661 (1992).
3. P. G. Debenedetti, Metastable Liquids (Princeton Univ. Press, Princeton, New Jersey, 1996).
4. D. Kashchiev, Nucleation: Basic Theory with Applications (Butterworth Heinemann, Oxford, UK, 2000).
5. M. P. Anisimov, Russ. Chem. Rev. 72, 591 (2003).
6. K. Binder and D. Stauffer, Advances in Physics 25, 343 (1976).
7. K. Binder, Rep. Prog. Phys. 50, 783 (1987).
8. O. Penrose and J. Lebowitz, J. Stat. Phys. 3, 211 (1971).
9. O. Penrose and J. Lebowitz, in Towards a Rigorous Molecular Theory of Metastability, edited by E. Montroll and J. Lebowitz (North Holland, Amsterdam, 1987), Chap. 5.
10. K. Binder and M. H. Kalos, Journal of Statistical Physics 22, 363 (1980), 10.1007/BF01014648.
11. H. Furukawa and K. Binder, Phys. Rev. A 26, 556 (1982).
12. K. Binder, Physica A: Statistical Mechanics and its Applications 319, 99 (2003).
13. D. Winter, P. Virnau, and K. Binder, Phys. Rev. Lett. 103, 225703 (2009).
14. B. J. Block, S. K. Das, M. Oettel, P. Virnau, and K. Binder, The Journal of Chemical Physics 133, 154702 (2010).
15. M. Schrader, P. Virnau, and K. Binder, Phys. Rev. E 79, 061104 (2009).
16. S. M. Thompson, K. E. Gubbins, J. P. R. B. Walton, R. A. R. Chantry, and J. S. Rowlinson, The Journal of Chemical Physics 81, 530 (1984).
17. P. R. ten Wolde and D. Frenkel, The Journal of Chemical Physics 109, 9901 (1998).
18. G. Kharlamov, A. Onischuk, S. Vosel, and P. Purtov, Colloids and Surfaces A: Physicochemical and Engineering Aspects 379, 10 (2011).
19. S. Ono and S. Kondo, in Molecular Theory of Surface Tension in Liquids, Vol. 10 of Handbuch der Physik, edited by S. Flügge (Springer, Berlin, 1960), pp. 134–280.
20. J. Rowlinson and B. Widom, Molecular Theory of Capillarity (Dover Publications, Inc, Mineola, New York, USA, 1982).
21. V. I. Kalikmanov, Statistical Physics of Fluids (Springer, Berlin, 2001).
22. H. Vehkamaki, Classical Nucleation Theory in Multicomponent Systems (Springer, Berlin, Heidelberg, New York, 2006).
23. G. S. Boltachev, V. G. Baidakov, and J. W. P. Schmelzer, Journal of Colloid and Interface Science 264, 228 (2003).
24. C. Borgs and R. Kotecky, J. Stat. Phys. 61, 79 (1990).
25. A. M. Ferrenberg and R. H. Swendsen, Phys. Rev. Lett. 61, 2635 (1988).
26. W. Thirring, Zeitschrift f. Physik A 235, 339 (1970).
27. D. H. E. Gross, Microcanonical Thermodynamics: Phase Transitions in ’Small’ Systems (World Scientific, Singapore, 2001).
28. H. Callen, Thermodynamics and an Introduction to Thermostatistics (Wiley & Sons, New York, 1985).
29. M. S. S. Challa, D. P. Landau, and K. Binder, Phys. Rev. B 34, 1841 (1986).
30. C. Borgs and R. Kotecky, Phys. Rev. Lett. 68, 1734 (1992).
31. F. Y. Wu, Rev. Mod. Phys. 54, 235 (1982).
32. T. Kihara, Y. Midzuno, and T. Shizume, J. Phys. Soc. Japan 9, 681 (1954).
33. R. J. Baxter, J. Phys. C 6, L445 (1973).
34. W. Janke, Phys. Rev. B 47, 14757 (1993).
35. C. Borgs and W. Janke, J. de Physique I 2, 2011 (1992).
36. C. Borgs and W. Janke, Phys. Rev. Lett. 68, 1738 (1992).
37. M. Fujimoto, J. Phys. A: Math. Gen. 30, 3779 (1997).
38. R. Zia and D. Wallace, J. Phys. A: Math. Gen. 8, 1495 (1975).
39. M. Weigel, Physics Procedia 3, 1499 (2010), proceedings of the 22th Workshop on Computer Simulation Studies in Condensed Matter Physics (CSP 2009).
40. T. Neuhaus and J. S. Hager, J. Stat. Phys. 113, 47 (2003).
41. D. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics (Cambridge University Press, Cambridge, 2000).
42. M. Newman and G. Barkema, Monte Carlo Methods in Statistical Physics (Claredon Press, Oxford, 1999).
43. M. P. Taylor, W. Paul, and K. Binder, Phys. Rev. E 79, 050801 (2009).
44. A. Tröster, C. Dellago, and W. Schranz, Phys. Rev. B 72, 094103 (2005).
45. A. Tröster, Phys. Rev. B 76, 012402 (2007).
46. W. Janke, Nucl. Phys. B (Proc. Suppl.) 63A-C, 631 (1998).
47. A. Tröster, M. Oettel, B. Block, P. Virnau, and K. Binder, (in preparation) (2011).
48. A. Tröster and K. Binder, Phys. Rev. Lett. accepted for publication (2011).
49. V. Martin-Mayor, Phys. Rev. Lett. 98, 137207 (2007).
50. B. Bauer, E. Gull, S. Trebst, M. Troyer, and D. Huse, J. Stat. Mech. 1, 01020 (2010).
100818