text:
# 1 twist

UCLA/00/TEP/06

INLO-PUB 02/00

Computation of the Vortex Free Energy in SU(2) Gauge Theory
Tamás G. Kovács

Department of Physics, Instituut-Lorentz for Theoretical Physics, P.O.Box 9506, 2300 RA, Leiden, The Netherlands
e-mail: kovacs@lorentz.leidenuniv.nl

and

E. T. Tomboulis

Department of Physics, UCLA, Los Angeles, CA 90095-1547
e-mail: tombouli@physics.ucla.edu

Abstract

We present the first measurement of the vortex free-energy order parameter at weak coupling for SU(2) in simulations employing multihistogram methods. The result shows that the excitation probability for a sufficiently thick vortex in the vacuum tends to unity. This is rigorously known to provide a necessary and sufficient condition for maintaining confinement at weak coupling in $`SU(N)`$ gauge theories.

The vortex free energy (also known as magnetic-flux free energy) order parameter in gauge theories is defined as the ratio of the partition function in the presence of a topologically trapped vortex excitation (introduced by a singular gauge transformation) to that without it. Its Fourier transform w.r.t. to the center ($`Z(N)`$) of the gauge group ($`SU(N)`$) defines the so-called electric-flux free energy which is rigorously known to provide an upper bound on the Wilson loop. These flux order parameters can characterize all possible phases of a (pure) gauge theory, and furthermore do this in terms of the behavior of the excitation expectation for a vortex. They were first considered in the study of gauge theories in , though the use of the analogous quantities in statistical mechanics goes back much further . The idea that vortex configurations underlie confinement at weak coupling has a long history, and has been the subject of intense recent activity. (We refer to Ref.  for a review of recent developments and references to early and recent work.)

In view of the physical significance of the magnetic-flux free energy, it may appear surprising that it has not been measured in simulations over the last twenty years. Accurate determination of (differences of) free energies in gauge theories, however, is well-known to be difficult. In fact, it is at first not quite clear how one should go about computing such totally nonlocal (lattice-length) quantities. We present here a computation for the group $`SU(2)`$ based on multihistogram methods . Such a method was recently used in Ref.  to compute the free energy of a pair of $`Z(N)`$ monopoles, a quantity related to the ’t Hooft loop operator. Our result demonstrates that the excitation expectation for a sufficiently extended ‘thick’ vortex at large $`\beta `$ is essentially unity. This is the feature responsible for maintaining the confining phase in $`SU(N)`$ gauge theories even at weak coupling.

We work on a $`d`$-dimensional hypercubic lattice $`\mathrm{\Lambda }`$ of size $`L_1\times \mathrm{}\times L_d`$ with periodic boundary conditions in all directions. We generally denote bonds by $`b`$, plaquettes by $`p`$, cubes by $`c`$, etc. The plaquette action is denoted by $`A_p(U_p)`$, where, as usual, $`U_p=_{bp}U_b`$, the product of the bond variables around the plaquette; for the minimal (Wilson) action $`A_p(U_p)=\beta \mathrm{Re}\mathrm{tr}U_p`$. The trace ”tr” is defined to include a $`1/N`$ normalization.

A coclosed set of plaquettes (2-cells) is a closed set of $`(d2)`$-cells on the dual lattice. Thus, in $`d=3`$, it is a closed loop of dual bonds; in $`d=4`$, a closed two-dimensional surface of dual plaquettes. For fixed $`\mu `$, $`\nu `$, let $`𝒱_{\mu \nu }`$ denote a coclosed set of plaquettes that winds through every 2-dim $`[\mu \nu ]`$-plane of $`\mathrm{\Lambda }`$, i.e. a topologically nontrivial plaquette set wrapped around the periodic lattice ($`d`$-torus) in the $`(d2)`$ directions $`\lambda \mu ,\nu `$ perpendicular to $`\mu ,\nu `$ This is depicted in figure 1(a), where the short lines represent the plaquettes in $`𝒱`$, with the horizontal axis representing the $`x^\mu `$, $`x^\nu `$ directions, and the vertical axis the remaining $`(d2)`$ perpendicular directions.

Define the partition function

$$Z_\mathrm{\Lambda }(\tau _{\mu \nu })=\underset{b}{}dU_b\mathrm{exp}\left(\underset{p𝒱_{\mu \nu }}{}A_p(U_p)\underset{p𝒱_{\mu \nu }}{}A_p(\tau _{\mu \nu }U_p)\right),$$
(1)
where the plaquette action $`A_p(U_p)`$ is replaced by the ‘twisted’ action $`A_p(\tau _{\mu \nu }U_p)`$ for each plaquette of $`𝒱_{\mu \nu }`$. Here the ‘twist’ $`\tau _{\mu \nu }Z(N)`$ is an element of the center. There are thus $`(N1)`$ different nontrivial choices for $`\tau _{\mu \nu }`$. The trivial element $`\tau _{\mu \nu }=1`$ is the ordinary partition function $`Z_\mathrm{\Lambda }(1)Z_\mathrm{\Lambda }`$.

As indicated by the notation on the l.h.s. of (1), the exact position or shape of $`𝒱_{\mu \nu }`$ is irrelevant; the only dependence is on the presence of the $`Z(N)`$ flux winding through each $`[\mu \nu ]`$-plane. It is indeed easily seen that $`𝒱_{\mu \nu }`$ can be moved around and distorted by a shift of integration variables, but not removed; it is rendered topologically stable by winding completely around the lattice (figure 1(b)). By the same token introducing two twists, $`\tau _{\mu \nu }`$ on $`𝒱_{\mu \nu }`$ and $`\tau _{\mu \nu }^{}`$ on $`𝒱_{\mu \nu }^{}`$ in (1), is equivalent to introducing one twist $`\tau _{\mu \nu }^{\prime \prime }=\tau _{\mu \nu }\tau _{\mu \nu }^{}`$ since $`𝒱_{\mu \nu }`$ and $`𝒱_{\mu \nu }^{}`$ can be brought together by a shift of integration variables (figure 2). This expresses the mod $`N`$ conservation of the $`Z(N)`$ flux introduced by the twist. Thus, for $`N=2`$, any odd number of such (nontrivial) twists is equivalent to one, and any even number to none.

The magnetic-flux free energy order parameter is now defined as

$`\mathrm{exp}(F_{\mathrm{mg}}(\tau _{\mu \nu }))`$ $`=`$ $`{\displaystyle \frac{Z_\mathrm{\Lambda }(\tau _{\mu \nu })}{Z_\mathrm{\Lambda }}}`$ (2)
$`=`$ $`\mathrm{exp}\left({\displaystyle \underset{p𝒱_{\mu \nu }}{}}\left(A_p(\tau _{\mu \nu }U_p)A_p(U_p)\right)\right).`$
Generalizations of (1)-(2) may be considered by introducing sets $`𝒱_{\kappa \lambda }`$ for several or all of the $`\frac{1}{2}d(d1)`$ possible distinct choices of planes $`[\kappa \lambda ]`$.

The twist amounts to a discontinuous (singular) $`SU(N)`$ gauge transformation on the configurations in (1) with multivaluedness in $`Z(N)`$ (so it is single-valued in $`SU(N)/Z(N)`$)), i.e. the introduction of a $`\pi _1(SU(N)/Z(N))=Z(N)`$ vortex. The set $`𝒱_{\mu \nu }`$ represents the topological obstruction to having singlevaluedness everywhere. (1) is then the partition sum for the system with a topologically stable vortex completely winding around the lattice; and (2) is the normalized expectation for the excitation of such a vortex. Hence, it is also referred to as the vortex free energy.

Choosing, say, $`[\mu \nu ]=[12]`$ in (1), we now drop the $`\mu \nu `$ subscript. One is interested in the behavior of (2) in the large volume limit (in the van Hove sense), i.e. as the size of the lattice increases in any power law fashion, e.g. $`L_\mu =2^{la_\mu }`$ for some fixed choice of positive exponents $`a_\mu `$, integer $`l\mathrm{}`$. Let $`A=L_1L_2`$ be the area of each $`[12]`$-plane, and $`L=L_3\mathrm{}L_d`$ the lattice volume in the perpendicular directions. One is interested, in particular, in $`LA`$. The twist introduces a cost in action localized on the plaquettes in $`𝒱`$. This cost, proportional to $`L`$, may be lowered if there are configurations that contribute with finite measure in the integral (2), and allow the flux introduced by the twist to spread in the two directions perpendicular to $`𝒱`$, so that the action is closer to its minimum; in other words, if there is finite probability for exciting a ‘thick’ vortex.

For sufficiently large lattices, there are then three possibilities ($`\tau 1`$):

1. $`\mathrm{exp}(F_{\mathrm{mg}}(\tau ))\mathrm{exp}(\alpha (\beta ,\tau )L)`$
2. $`\mathrm{exp}(F_{\mathrm{mg}}(\tau ))\mathrm{exp}(\beta c(\tau ){\displaystyle \frac{L}{A}})`$
3. $`\mathrm{exp}(F_{\mathrm{mg}}(\tau ))\mathrm{exp}(cLe^{\rho (\beta ,\tau )A})`$

In case (a) the magnetic flux stays focused in a thin vortex; this describes a Higgs phase. In (b) the flux can spread in a Coulomb-like fashion lowering the free-energy cost; this describes a massless Coulomb phase, where the long distance behavior is accurately given by weak coupling pertubative expansion. In (c) the gain in thickening the vortex is exponential; this characterizes the confinement phase. It is important to note that, in contrast to (a)-(b), only (c) gives a value which survives and in fact tends (exponentially) to unity for all ways of taking the thermodynamic limit as described above; this is the signature of the confinement phase.

Since our computation below is for $`N=2`$, we now write explicit formulae only for this case. The Fourier transform of (2) w.r.t. $`Z(N)`$ is known as the electric-flux free energy. For $`N=2`$ this is simply:

$$\mathrm{exp}(F_{\mathrm{el}})=\underset{\tau =1,1}{}\tau \mathrm{exp}(F_{\mathrm{mg}}(\tau ))=1\mathrm{exp}(F_{\mathrm{mg}}(1)).$$
(3)
Consider now a rectangular loop $`C`$ in a $`[12]`$-plane. Then, for any reflection positive plaquette action, the Wilson loop obeys the bound :

$$\mathrm{tr}(U[C])\left(\mathrm{exp}(F_{\mathrm{el}})\right)^{{\scriptscriptstyle \frac{A_C}{A}}},$$
(4)
where $`A_C`$ is the minimal area bounded by $`C`$. (4) shows that confining behavior (c) for the vortex free energy implies area-law for the Wilson loop with string tension bounded from below by the excitation expectation for a vortex. So confining behavior for the vortex free energy is a sufficient condition for linear asymptotic quark confinement.

Placing suitable constraints in the functional measure in (1) which forbid the spreading of flux across $`[12]`$-planes, thus eliminating the occurance of thick vortices, results in nonconfining behavior of type (a) above . In this case (4) cannot tell us anything about the Wilson loop. To show loss of confining behavior for the Wilson loop itself in the presence of such constraints, one needs a lower bound on it which exhibits perimeter-law. This was recently proven for large $`\beta `$ in . Thus the occurance of thick vortices is also a necessary condition for confinement at weak coupling.

Our measurement of (2) for $`SU(2)`$ was done by an application of a multihistogram method . From now on we restrict the form of the action to the Wilson action

$$A_p(U_p)=\beta \mathrm{tr}U_p,$$
(5)
which was used in the measurement. The basic quantity in our procedure is the density of states $`w(S)`$ as a function of the total action $`S`$ along the twisted plaquettes. This is defined as

$$w(S)=𝑑U_b\mathrm{exp}\left(\beta \underset{p𝒱}{}\mathrm{tr}U_p\right)\delta (S+\underset{p𝒱}{}\mathrm{tr}U_p).$$
(6)
If $`w(S)`$ is known, the partition function can be easily computed for any coupling $`\beta _𝒱`$ along $`𝒱`$ as

$$Z(\beta _𝒱)=𝑑Sw(S)e^{\beta _𝒱S}.$$
(7)
In particular, we are interested in $`Z(\beta )`$, the untwisted, and $`Z(\beta )`$, the twisted partition function. The problem is that the dominant contribution for $`Z(\beta _𝒱)`$ comes from different regions of $`S`$, depending on $`\beta _𝒱`$. Therefore one needs to know $`w(S)`$ to a good accuracy in a wide range of $`S`$. A simulation done at a certain value of $`\beta _𝒱`$, however will give accurate information on $`w(S)`$ only in a narrow neighbourhood of $`S_{\beta _𝒱}`$. The main idea of the Ferrenberg-Swendsen multihistogram method is to combine information on $`w(S)`$ coming from simulations at different $`\beta _𝒱`$’s to obtain $`w(S)`$ in a wide range of $`S`$ accurately. This can be done by noting that for a given $`\beta _𝒱`$ the probability distribution of $`S`$, $`P(S,\beta _𝒱)`$, goes as

$$P(S,\beta _𝒱)\frac{1}{Z(\beta _𝒱)}e^{\beta _𝒱S},$$
(8)
and that $`P(S,\beta _𝒱)`$ can be directly measured by making a histogram of the action along $`𝒱`$. In this way, any simulation at a certain $`\beta _𝒱`$ gives an estimate for $`w(S)`$,

$$w(S)=P(S,\beta _𝒱)e^{\beta _𝒱S}Z(\beta _𝒱).$$
(9)
These estimates coming from simulations with different $`\beta _𝒱`$’s (say $`\beta _1`$, $`\beta _2`$,…$`\beta _K`$) can then be averaged with suitable ($`S`$ dependent) weights to minimise the error in $`w(S)`$ over a given range of $`S`$. This results in the following set of coupled equations:

$`w(S)={\displaystyle \frac{_{n=1}^KP(S,\beta _n)}{_{n=1}^K\frac{\mathrm{exp}\left(\beta _nS\right)}{Z\left(\beta _n\right)}}}`$ (10)
$`Z(\beta _n)={\displaystyle 𝑑Se^{\beta _nS}w(S)},`$ (11)
which can be solved by iteration starting from $`Z(\beta _n)=1`$. To optimize the procedure, one needs a sufficient overlap between the $`P(S,\beta _n)`$ distributions corresponding to successive $`\beta _n`$’s. Since the distributions quickly become narrower with increasing lattice size, the number of simulations, $`K`$, also needs to be increased accordingly. This makes our measurement very expensive on large lattices. For the largest lattices we typically used $`K=4080`$ simulations with the $`\beta _n`$’s equally spaced in the $`\beta `$ $`+\beta `$ range.

The result of the computation for (2) is shown in figure 3. We have performed the computation on lattices of equal linear size in all directions for three different values of $`\beta `$. The lattice spacings are $`a=0.165`$ fm, $`a=0.119`$ fm and $`a=0.085`$ fm for $`\beta =2.3`$, $`\beta =2.4`$, and $`\beta =2.5`$, resp.

Notice that, with the lattice size expressed in physical units, the measurements for different $`\beta `$’s fall on the same curve, as they should. This indicates that the universal curve has been reached, and will not change at larger beta. Also, the onset of the sharp rise around $`0.7`$ fm is in the region of the finite temperature deconfining phase transition providing another indirect consistency check.

The approach to unity for sufficiently large lattice size in figure 3 is striking. In comparison, for Coulomb-like massless behavior, an upper bound obtained by action minimizing within the spin-wave approximation gives $`\mathrm{exp}(\beta (\pi /2)^2)0.085`$ at $`\beta =2.3`$. The points forming the upper part of the plot are well within the confinement region. The string tension values extracted from the vortex free energy in the confining region are consistent with the values from heavy-quark potential calculations (see eg. ), though still better precision in the measurement of the vortex free energy is required for precise quantitative comparisons.

In conclusion, the result of our computation clearly demonstrates that the weighted expectation for the excitation of a sufficiently thick vortex in the vacuum tends to one. In this sense the vacuum can indeed be viewed as having a ‘condensate’ of thick long vortices. This is sufficient for maintaining confinement at large $`\beta `$ in $`SU(N)`$ gauge theories. As mentioned, rigorous results also show it to be necessary: were the behavior for (2) exhibited in figure 3 not to occur, confinement at large beta would be lost.

We are very grateful to P. de Forcrand for correspondence. This research was supported by FOM (T.G.K.) and by NSF, Grant No. NSF-PHY-9819686 (E.T.T.).


format: text

text:
# The influence of critical behavior on the spin glass phase

## I Introduction

Despite over two decades of work, the controversy concerning the nature of the ordered phase of short range Ising spin glasses continues. For a few years, Monte Carlo simulations appeared to be providing evidence for replica symmetry breaking (RSB) in these systems. However, recent developments have cast doubt on this interpretation of the Monte Carlo data. In a series of papers on the Ising spin glass within the Migdal-Kadanoff approximation (MKA), we showed that the equilibrium Monte Carlo data in three dimensions that had been interpreted in the past as giving evidence for RSB can actually be interpreted quite easily within the droplet picture, with apparent RSB effects being attributed to a crossover between critical behavior and the asymptotic droplet-like behavior for small system sizes. We also showed that system sizes well beyond the reach of current simulations would probably be required in order to unambiguously see droplet-like behavior. The finding that the critical-point effects can still be felt at temperatures lower than those accessible by Monte Carlo simulations is supported by the Monte Carlo simulations of Berg and Janke who found critical scaling working reasonably well down to $`T=0.8T_c`$ for system sizes upto $`L=8`$ in three dimensions. The zero temperature study of Pallasini and Young also suggests that the ground- state structure of three-dimensional Edwards-Anderson model is well described by droplet theory, though the existence of low energy excitations not included in the conventional droplet theory remains an open question. Thus, while puzzles do remain, the weight of the evidence seems to be shifting towards a droplet-like description of the ordered phase in short range Ising spin glasses.

However, it is expected that critical point effects are less dominant in four dimensions than in three dimensions. Our aim in this paper is to quantify the extent of critical point effects in the low temperature phase of the four-dimensional Edwards-Anderson spin glass. We do this by providing results for the four-dimensional Ising spin glass in the MKA and compare these with existing Monte Carlo work. In particular, we study the Parisi overlap function and the link overlap function for system sizes up to $`L=16`$ and temperatures as low as $`T=0.16T_c`$. We find that for system sizes and temperatures comparable to those of the Monte Carlo simulations, the Parisi overlap distribution shows also in MKA the sample-to-sample fluctations and the stationary behavior at small overlap values, that are normally attributed to RSB. It is only for larger system sizes (or for lower temperatures), that the asymptotic droplet-like behavior becomes apparent. For the link overlap, we find similar double-peaked curves as those found in Monte-Carlo simulations. This double peak structure is expected on quite general grounds independent of the nature of the low temperature phase. However, we show that two peaks in the link overlap in MKA occur because of a difference between domain-wall excitations (which cross the entire system) and droplet excitations (which do not cross the entire system). We argue that for small system sizes, the effect of domain walls increases with increasing dimension, making it necessary to go very far below $`T_c`$ to see the asymptotic droplet behavior.

This paper is organized as follows: in section II, we define the quantities discussed in this paper, and the droplet-model predictions for their behavior. In section III, we describe the MKA, and our numerical methods of evaluating the overlap distribution. In section IV, we present our numerical results for the Parisi overlap distribution, and compare to Monte-Carlo data. The following section studies the link overlap distribution. Finally, section VI contains the concluding remarks, including some on the effects of critical behavior on the dynamics in the spin glass phase. Again we suspect that arguments which have been advanced against the droplet picture on the basis of dynamical studies have failed to take into account the effects arising from proximity to the critical point.

## II Definitions and Scaling Laws

The Edwards-Anderson spin glass in the absence of an external magnetic field is defined by the Hamiltonian

$$H=\underset{i,j}{}J_{ij}S_iS_j,$$

where the Ising spins can take the values $`\pm 1`$, and the nearest-neighbor couplings $`J_{ij}`$ are independent from each other and Gaussian distributed with a standard deviation $`J`$.

It has proven useful to consider two identical copies (replicas) of the system, and to measure overlaps between them. This gives information about the structure of the low-temperature phase, in particular about the number of pure states. The quantities considered in this paper are the Parisi overlap function $`P(q,L)`$ and the link overlap function $`P(q_l,L)`$. They are defined by

$$P(q,L)=\left[\delta \left(\underset{ij}{}\frac{S_i^{(1)}S_i^{(2)}+S_j^{(1)}S_j^{(2)}}{2N_L}q\right)\right],$$
(1)
and

$$P(q_l,L)=\left[\delta \left(\underset{ij}{}\frac{S_i^{(1)}S_i^{(2)}S_j^{(1)}S_j^{(2)}}{N_L}q_l\right)\right].$$
(2)
Here, the superscripts $`(1)`$ and $`(2)`$ denote the two replicas of the system, $`N_L`$ is the number of bonds, and $`\mathrm{}`$ and $`\left[\mathrm{}\right]`$ denote the thermodynamic and disorder average respectively. We use $`P(q,L)`$ and $`P(q_l,L)`$ to denote the overlap functions for a finite system of size $`L`$, reserving the more standard notation $`P(q)`$ and $`P(q_l)`$ for the limit $`lim_L\mathrm{}P(q,L)`$ and $`lim_L\mathrm{}P(q_l,L)`$.

In the mean-field RSB picture, $`P(q)`$ is nonzero in the spin glass phase in the entire interval $`[q_{EA},q_{EA}]`$, while it is composed only of two delta functions at $`\pm q_{EA}`$ in the droplet picture. Similarly, $`P(q_l)`$ is nonzero over a finite interval $`[q_l^{min},q_l^{max}]`$ in mean-field theory, while it is a delta-function within the droplet picture.

Much of the evidence for RSB for three- and four-dimensional systems comes from observing a stationary $`P(q=0,L)`$ for system sizes that are generally smaller than 20 in 3D and smaller than 10 in 4D, and at temperatures of the order of $`0.7T_c`$. However, even within the droplet picture one expects to see a stationary $`P(q=0,L)`$ for a certain range of system sizes and temperatures. The reason is that at $`T_c`$ the overlap distribution $`P(q,L)`$ obeys the scaling law

$$P(q,L)=L^{\beta /\nu }\stackrel{~}{P}(qL^{\beta /\nu }),$$
(3)
$`\beta `$ being the order parameter critical exponent, and $`\nu `$ the correlation length exponent. Above the lower critical dimension (which is smaller than 3), $`\beta /\nu `$ is positive, leading to an increase $`P(q=0,L)`$ as a function of $`L`$ (at $`T=T_c`$). On the other hand, for $`TT_c`$, the droplet model predicts a decay

$$P(q=0,L)1/L^\theta $$

on length scales larger than the (temperature–dependent) correlation length $`\xi `$, $`\theta `$ being the scaling exponent of the coupling strength $`J`$. A few words are in order here on what we mean by the correlation length. In the spin glass phase, all correlation functions fall off as a power law at large distances. However, within the droplet model, this is true only asymptotically, and the general form of the correlation function for two spins a distance $`r`$ apart, at a temperature $`TT_c`$, is $`r^\theta f(r/\xi )`$ where $`k_B`$ is the Boltzmann constant and $`f`$ is a scaling function. Thus, for $`r\xi `$ there are corrections to the algebraic long-distance behavior and the above expression defines the temperature-dependent correlation length. Note that for $`TT_c`$ this correlation length is expected to diverge with the exponent $`\nu `$.

Thus, for temperatures not too far below $`T_c`$, one can expect an almost stationary $`P(q=0,L)`$ for a certain range of system sizes. In three dimensions both $`\beta /\nu 0.3`$ and $`\theta 0.17`$ are rather small, this apparent stationarity may persist over a considerable range of system sizes $`L`$. However, in four dimensions, $`\beta /\nu 0.85`$  and $`\theta 0.65`$  and one would expect the crossover region to be smaller. In the present paper we shall investigate these crossover effects in four dimensions by studying $`P(q,L)`$ for the Edwards-Anderson spin glass within the MKA. It turns out that they are surprisingly persistent even at low temperatures, due to the presence of domain walls.

Monte-Carlo simulations of the link overlap distribution show a nontrivial shape with shoulders or even a double peak, which seems to be incompatible with the droplet picture, where the distribution should tend towards a delta-function. For sufficiently low temperatures and large length scales, the droplet picture predicts that the width of the link overlap distribution scales as

$$\mathrm{\Delta }q_l\sqrt{kT}L^{d_sd\theta /2},$$

where $`d_s`$ is the fractal dimension of a domain wall. Below, we will show that the nontrivial shape and the double peak reported from Monte-Carlo simulations are also found in MKA in four dimensions, and we will present strong evidence that it is due to the different nature of droplet and domain wall excitations. As the weight of domain walls becomes negligible in the thermodynamic limit, the droplet picture is regained on large scales.

## III Migdal-Kadanoff approximation

The Migdal-Kadanoff approximation (MKA) is a real-space renormalization group the gives approximate recursion relations for the various coupling constants. Evaluating a thermodynamic quantity in MKA in $`d`$ dimensions is equivalent to evaluating it on an hierarchical lattice that is constructed iteratively by replacing each bond by $`2^d`$ bonds, as indicated in Fig. 1. The total number of bonds after $`I`$ iterations is $`2^{dI}`$. $`I=1`$, the smallest non-trivial system that can be studied, corresponds to a system linear dimension $`L=2`$, $`I=2`$ corresponds to $`L=4`$, $`I=3`$ corresponds to $`L=8`$ and so on. Note that the number of bonds on hierarchical lattice after $`I`$ iterations is the same as the number of sites of a $`d`$-dimensional lattice of size $`L=2^I`$. Thermodynamic quantities are then evaluated iteratively by tracing over the spins on the highest level of the hierarchy, until the lowest level is reached and the trace over the remaining two spins is calculated . This procedure generates new effective couplings, which have to be included in the recursion relations.

In , it was proved that in the limit of infinitely many dimensions (and in an expansion away from infinite dimensions) the MKA reproduces the results of the droplet picture.

As was discussed in , the calculation of $`P(q,L)`$ is made easier by first calculating its Fourier transform $`F(y,L)`$, which is given by

$$F(y,L)=\left[\mathrm{exp}[iy\underset{ij}{}\frac{(S_i^{(1)}S_i^{(2)}+S_j^{(1)}S_j^{(2)})}{2N_L}]\right].$$
(4)
The recursion relations for $`F(y,L)`$ involve two- and four-spin terms, and can easily be evaluated numerically because all terms are now in an exponential. Having calculated $`F(y)`$ one can then invert the Fourier transform to get $`P(q,L)`$.

Similarly, $`P(q_l,L)`$ is calculated by first evaluating

$$F(y_l,L)=\left[\mathrm{exp}[iy_l\underset{ij}{}\frac{(S_i^{(1)}S_i^{(2)}S_j^{(1)}S_j^{(2)})}{N_L}]\right].$$
(5)
Before presenting our numerical results for the Parisi overlap and the link overlap, let us discuss the flow of the coupling constant $`J`$ in the low-temperature phase, as obtained in MKA. In order to obtain this flow, we iterated the MKA recursion relation on a set of $`10^6`$ bonds. At each iteration, each of the new set of $`10^6`$ bonds was generated by randomly choosing 16 bonds from the old set and taking the trace over the inner spins (with a bond arrangement as in Fig. 1). Figure 2 shows $`J/T`$ as function of $`L`$ for different initial values of the coupling strength. The critical point is at $`T_c2.1J`$. The first curve begins at $`J/T=0.5`$, which is close to the critical point, and it reaches the low-temperature behavior only at lengths around 1000. For an initial $`J/T=0.7`$, the asymptotic slope is already reached at $`L`$ around 40, and for $`J/T=3.0`$, which corresponds to $`T0.16T_c`$ the entire curve shows the asymptotic slope. The asymptotic slope is identical to the above-mentioned exponent $`\theta `$ and has the value $`\theta 0.75`$. In contrast to $`d=3`$ , we did not succeed in fitting the crossover regime by doing an expansion around the zero-temperature fixed point. The reason is that dimension 4 is too far above the lower critical dimension, such that the critical temperature is not small.

Note that for each temperature the length scale beyond which the flows of the coupling constants show the asymptotic behavior yields one estimate for the correlation length mentioned above. We have considered the flow to be in the asymptotic regime when its slope was within 90% of its asymptotic value. However, this estimate is specific to the flows of the coupling constant, and other quantities may show their asymptotic behavior later. In fact, as we shall see below, the convergence of the overlap distributions is much slower than that of the couplings, and we will have to give reasons for this.

## IV The Parisi overlap

We now discuss our results for the Parisi overlap. First, let us briefly describe the critical behavior. Fig. 3 shows a scaling plot for $`P(q,L)`$ for $`L=4,8,16`$ at $`T=T_c2.1J`$. We find a good data collapse if we use the value $`\beta /\nu =0.64`$, thus confirming the finite-size scaling ansatz Eq. 3.

We next move on to the low-temperature phase. In Fig. 4 we show $`P(q,L)`$ at $`T=0.5T_c`$ and $`L=8`$ for three different samples.

As one can see there are substantial differences between the samples. This sensitivity to samples for system sizes around 10 is in  interpreted as evidence for RSB. In our case, where we know that the droplet model is exact, it has to be considered a finite size effect. Note that we have not chosen the three samples in any particular manner. By comparing to the curves obtained for $`L=16`$ (not shown), we can even see the trend to an increasing number of peaks, just as in . Thus, one feature commonly associated with RSB is certainly present in within the MKA for temperatures and system sizes comparable to those studied in simulations.

Let us now focus on the behavior of $`P(q=0,L)`$ for different system sizes and temperatures. But before exhibiting our own data, we discuss the Monte Carlo data of Reger, Bhatt and Young who were the first to study $`P(q=0,L)`$ for the Edwards-Anderson spin glass. They studied system sizes $`L=2,3,4,5,6`$ at temperatures down to $`T=0.68T_c`$. At $`T=T_c`$ they found the expected critical scaling, $`P(q=0)L^{\beta /\nu }`$ with $`\beta /\nu 0.75`$. Then, as the temperature was lowered, the curves for $`P(q=0)`$ as a function of $`L`$ showed a downward curvature for the largest system sizes, which they interpreted as the beginning of the crossover between critical behavior and the low temperature behavior. At $`T=0.8T_c`$, $`P(q=0)`$ seemed to be roughly constant or decreasing slowly. However, the striking part of their data was that at $`T=0.68T_c`$ they found that $`P(q=0,L)`$ initially decreased as a function of system size for $`L=2,3,4`$ and then saturated for $`L=4,5,6`$. They interpreted this as suggestive of RSB. They admitted however that other explanations are possible.

The most recent Monte-Carlo simulation data for the 4d Ising spin glass are those in . These authors focus on $`T0.6T_c`$, and they find an essentially stationary $`P(q=0,L)`$ for system sizes up to the largest simulated size $`L=10`$. They argue, that stationarity over such a large range of $`L`$ values is most naturally interpreted as evidence for RSB. However, as can be seen from Fig. 3, the correlation length is of the order of 16 for these temperatures and therefore comparable to the system size.

In Fig.5, we show the MKA data for $`P(q=0,L)`$. We have calculated $`P(q=0,L)`$ for system sizes $`L=4`$,$`8`$,$`16`$ at temperatures $`T=T_c`$, $`0.68T_c`$, $`0.33T_c`$, and $`0.16T_c`$. At $`T=T_c`$, $`P(q=0,L)`$ grows as $`L^{\beta /\nu }`$ with $`\beta /\nu 0.64`$, in agreement with Fig.3. At $`T=0.68T_c`$ (the lowest temperature studied in , and not far from the lowest temperature studied by ), we do not see a clear decrease even for $`L=16`$. The curve for $`P(q=0)`$ looks more or less flat, though one could say that there is slight increase between $`L=4`$ and $`L=8`$ and a slight decrease between $`L=8`$ and $`L=16`$. This flat behavior is similar to what was found in  and . The deviation of the $`L=2`$ and $`L=3`$ data from the flat curve in  can probably be ascribed to artifacts at very small system sizes, which are also found elsewhere . For lower temperatures, where the correlation length is smaller than the system size, there is a clear decrease of $`P(q=0)`$ although the decrease is not asympotic even at a temperatures as low as $`T_c/6`$.

We conclude that the observed stationarity of $`P(q=0,L)`$ in Monte-Carlo data is due to the effects of a finite system size and finite temperature. Similarly, Monte-Carlo simulations at $`T0.5T_c`$ and at system sizes around 10, should be able to show the negative slope in $`P(q=0,L)`$. In the not too far future, it should become possible to perform these simulations.

The fact that $`P(q=0,L)`$ does not show asymptotic behavior even at $`T=T_c/6`$ for the system sizes that we have studied, is surprising, and is different from our findings in $`d=3`$ . That $`P(q=0,L)`$ converges slower towards the asymptotic behavior than the flow of the coupling constant (see Fig. 2), can be understood in the following way: A Parisi overlap value close to zero can be generated by a domain wall excitation. For large system sizes and low temperatures, such an excitation occurs with significant weight only in those samples where a domain wall excitation costs little energy. These are exactly the samples with a small renormalized coupling constant at system size $`L`$. As the width of the probability distribution function of the couplings increases with $`L^\theta `$, the probability for obtaining a small renormalized coupling decreases as $`L^\theta `$. This is the argument that predicts that $`P(q=0,L)L^\theta `$. However, for smaller system sizes and higher temperatures, there are corrections to this argument. Thus, even samples with a renormalized coupling that is not small can contribute to $`P(q=0,L)`$ by means of large or multiple droplet excitations, or of thermally activated domain walls. For this reason, $`P(q=0,L)`$ can be expected to converge towards asymptopia slower that the coupling constant itself. Furthermore, as we shall see in the next section, the superposition of domain wall excitations and droplet excitations leads to deviations from simple scaling, which may further slow down the convergence towards asymptotic scaling behavior.

## V The Link Overlap

The link overlap gives additional information about the spin glass phase that is not readily seen in the Parisi overlap. The main qualitative differences between the Parisi overlap and the link overlap are (i) that flipping all spins in one of the two replicas changes the sign of $`q`$ but leaves $`q_l`$ invariant, and (ii) that flipping a droplet of finite size in one of the two replicas changes $`q`$ by an amount proportional to the volume of the droplet, and $`q_l`$ by an amount proportional to the surface of the droplet. Thus, the link overlap contains information about the surface area of excitations.

First, let us study $`P(q_l,L)`$ as function of temperature, for a given system size $`L=4`$. Fig. 6 shows our curves for $`T=0.8T_c`$, $`0.67T_c`$, $`0.56T_c`$, $`0.48T_c`$, and $`0.33T_c`$. They appear to result from the superposition of two different peaks, with their distance increasing with decreasing temperature, and the weight shifting from the left peak to the right peak.

Fig. 7 shows $`P(q_l,L)`$ for fixed $`T=0.33T_c`$ and for different $`L`$. One can see that with increasing system size the peaks move closer together, and the weight of the left-hand peak decreases.

These results are similar to what we found in MKA in three dimensions , however, in four dimensions the peaks are more pronounced. Monte-Carlo simulations of the four-dimensional Ising spin glass also show two peaks for certain system sizes and temperatures . This feature is attributed by the authors to RSB. However, as it is also present in MKA, there must be a different explanation. The width of the curves shrinks with increasing system size in , just as it does in MKA and as is expected from the droplet picture. If the RSB scenario were correct, the width would go to a finite value in the limit $`L\mathrm{}`$.

In the following we present evidence that the left peak corresponds to configurations where one of the two replicas has a domain wall excitation, and the right peak to configurations where one of the two replicas has a droplet excitation. In MKA, domain wall excitations involve flipping of one side of the system, including one of the two boundary spins of the hierarchical lattice, while droplet excitations involve flipping of a group of spins in the interior. If the sign of the renormalized coupling is positive (negative), the two boundary spins are parallel (antiparallel) in the ground state. By plotting separately the contributions from configurations with and without flipped boundary spins, we can separate domain wall excitations from droplet excitations. Fig. 8 shows the three contributions from configurations where none, one, or both replicas have a domain wall. Clearly, the left peak is due to domain wall excitations, and the right peak to droplet excitations.

Similar curves are obtained for other values of the parameters. We thus have shown that the qualitative differences between droplet and domain wall excitations are sufficient to explain the structure of the link overlap distribution, and no other low-lying excitations like those invoked by RSB are needed.

The weight with which domain-wall excitations occur is in agreement with predictions from the droplet model. The probability of having a domain wall in a system of size $`L`$ is according to the droplet picture of the order of

$$(T/J)L^\theta ,$$

which is $`0.25`$ at $`T=0.33T_c`$ and $`L=4`$, and $`0.15`$ at $`T=0.33T_c`$ and $`L=8`$. From our simulations, we find that the relative weights of domain walls for these two situations are $`0.12`$ and $`0.076`$, which fits the droplet picture very well if we include a factor 1/2 in the above expression. Domain walls become negligible only when the product $`(T/J)L^\theta `$ becomes small. In higher dimensions, the critical value of $`T/J`$ becomes larger, and for a given relative distance from the critical point, the weight of domain walls therefore also becomes larger. This explains why the effect of domain walls is more visible in 4 dimensions than in 3 dimensions. However, with increasing system size, domain walls should become negligible more rapidly in higher dimensions, due to the larger value of the exponent $`\theta `$.

## VI Conclusions

Our results for the Parisi overlap distribution in four dimensions show that there are rather large finite size effects in four dimensions which give rise to phenomena normally attributed to RSB. The system sizes needed to see the beginning of droplet like behavior within the MKA are larger, and the temperatures are lower, than those studied by Monte Carlo simulations. However, at temperatures not too far below those studied in Monte Carlo simulations ($`T=0.5T_c`$), the weight of the Parisi overlap distribution function $`P(q=0,L)`$ within the MKA appears to decrease, albeit with an effective exponent different from the asymptotic value. Thus, simulations at these temperatures for the Ising spin glass on a cubic lattice might resolve the controversy regarding the nature of the ordered state in short range spin glasses. However, the MKA is a low dimensional approximation and it is possible that the system sizes needed to see asymptotic behavior for a hypercubic lattice in four dimensions are different from what is indicated by the MKA. So, any comparison of the MKA with the Monte Carlo data should be taken with a pinch of salt.

Recently, a modified droplet picture was suggested by Houdayer and Martin , and by Bouchaud . Within this picture, excitations on length scales much smaller than the system size are droplet-like, however, there exist large-scale excitations that extend over the entire system and that have a small energy that does not diverge with increasing system size. As we have demonstrated within MKA, the double-peaked curves for the link overlap distribution, can be fully explained in terms of two types of excitations that contribute to the low-temperature behavior, namely domain-wall excitations and droplet excitations. We therefore believe that there is no need to invoke system-wide low-energy excitations that are more relevant than domain walls.

Finally, the whole field of dynamical studies of spin glasses is thought by many  to provide a strong reason for believing the RSB picture. A very recent study of spin glass dynamics on the hierarchical lattice , on which the MKA is exact, indicates that no ageing occurs at low temperatures in the response function whereas in Monte-Carlo simulations on the Edwards-Anderson model and in spin glass experiments ageing is seen in the response function. We suggest that the ageing behavior found in Monte-Carlo simulations and experiment are in fact often dominated by critical point effects, and not by droplet effects. Indeed we would expect that if the simulations of Ref.  were performed at temperatures closer to the critical temperature then ageing effects would be seen in the response function, since near the critical point of even a ferromagnet such ageing effects occur . The reason why experiments and simulations on the Edwards-Anderson model see ageing in the response function is that they are probing time scales that may be less than the critical time scale, which is given by

$$\tau =\tau _0(\xi /a)^z,$$

with $`a`$ the lattice constant and $`\tau _0`$ the characteristic spin-flip time. The dynamical critical exponent $`z6`$ in 3 dimensions . Only for droplet reversals which take place on time scales larger than $`\tau `$ (i.e for reversals of droplets whose linear dimensions exceed $`\xi `$) will droplet results for the dynamics be appropriate. However, because of the large values of $`\xi `$ down to temperatures of at least $`0.5T_c`$ and the large value of $`z`$, $`\tau `$ may be very large in the Monte-Carlo simulations and experiments. Thus if $`\xi /a`$ is 100, then $`\tau /\tau _0`$ is $`10^{12}`$, which would make droplet like dynamics beyond the reach of a Monte-Carlo simulation. In practice, most data will be in a crossover regime leading to an apparently temperature dependent exponent $`z(T)`$ (see for example Ref ).

###### Acknowledgements.

We thank A. P. Young for discussions and for encouraging us to write this paper. Part of this work was performed when HB and BD were at the Department of Physics, University of Manchester, supported by EPSRC Grants GR/K79307 and GR/L38578. BD also acknowledges support from the Minerva foundation.


text:
# Density functional theory of phase coexistence in weakly polydisperse fluids

## Abstract

The recently proposed universal relations between the moments of the polydispersity distributions of a phase-separated weakly polydisperse system are analyzed in detail using the numerical results obtained by solving a simple density functional theory of a polydisperse fluid. It is shown that universal properties are the exception rather than the rule.

Département de Physique des Matériaux (UMR 5586 du CNRS),
Université Claude Bernard-Lyon1, 69622 Villeurbanne Cedex, France
<sup>††</sup> Physique des Polymères, Université Libre de Bruxelles,
Campus Plaine, CP 223, B-1050 Brussels, Belgium

PACS numbers: 05.70.-a, 64.75.+g, 82.60.Lf

Many, natural or man-made, systems are mixtures of similar instead of identical objects. For example, in a colloidal dispersion the size and surface charge of the colloidal particles are usually distributed in an almost continuous fashion around some mean value. When this distribution is very narrow the system can often be assimilated  to a one-component system of identical objects. Such a system is usually called monodisperse whereas otherwise it is termed polydisperse. Since polydispersity is a direct consequence of the physico-chemical production process it is an intrinsic property of many industrial systems. Therefore, many author have included polydispersity into the description of a given phase of such systems. More recently, a renewed interest can be witnessed for the study of phase transitions occuring in weakly polydisperse systems. The phase behavior of polydisperse systems is of course much richer than that of its monodisperse counterpart. It is also more difficult to study theoretically, essentially because one has to cope with an infinity of thermodynamic coexistence conditions. Therefore, several authors have proposed approximation schemes which try to bypass this difficulty. In the present study we take to opposite point of view by solving numerically the infinitely many thermodynamic coexistence conditions for a simple model polydisperse system. On this basis we have studied the radius of convergence of the weak polydispersity expansion used in ref.4 and found that their “universal law of fractionation” and some of their conclusions have to be modified in several cases.

The statistical mechanical description of a polydisperse equilibrium system is equivalent to a density functional theory for a system whose number density, $`\rho (𝐫,\sigma )`$, depends besides the position variable $`𝐫`$ (assuming spherical particles) also on at least one polydispersity variable $`\sigma `$ (which we consider to be dimensionless). Such a theory is completely determined once the intrinsic Helmholtz free-energy per unit volume, $`f[\rho ]`$, has been specified as a functional of $`\rho (𝐫,\sigma )`$ (for notational convenience the dependence on the temperature $`T`$ will not be indicated explicitly). For the spatially uniform fluid phases considered here (and also implicitly in ref.4) we have, $`\rho (𝐫,\sigma )\rho (\sigma )`$, and the pressure can be written as, $`p[\rho ]=𝑑\sigma \rho (\sigma )\mu (\sigma ;[\rho ])f[\rho ]`$, where $`\mu (\sigma ;[\rho ])=\delta f[\rho ]/\delta \rho (\sigma )`$, is the chemical potential of “species” $`\sigma `$. When a parent phase of density $`\rho _0(\sigma )`$ phase separates into $`n`$ daughter phases of density $`\rho _i(\sigma )`$ ($`i=1,\mathrm{},n`$) the phase coexistence conditions imply that, $`p[\rho _1]=p[\rho _2]=\mathrm{}=p[\rho _n]`$, and $`\mu (\sigma ;[\rho _1])=\mu (\sigma ;[\rho _2])=\mathrm{}=\mu (\sigma ;[\rho _n])`$. For simplicity we consider here only the case of two daughter phases ($`n=2`$) and rewrite moreover $`\rho _i(\sigma )=\rho _ih_i(\sigma )(i=0,1,2)`$ in terms of the average density $`\rho _i`$ and a polydispersity distribution $`h_i(\sigma )`$ such that $`𝑑\sigma h_i(\sigma )=1`$. Since the ideal gas contribution to $`f[\rho ]`$ is exactly known one has, $`\mu (\sigma ;[\rho ])=k_BT\mathrm{ln}\{\mathrm{\Lambda }^3(\sigma )\rho (\sigma )\}+\mu _{ex}(\sigma ;[\rho ])`$, where $`k_B`$ is Boltzmann’s constant, $`\mathrm{\Lambda }(\sigma )`$ is the thermal de Broglie wavelength of species $`\sigma `$ and $`\mu _{ex}`$ the excess (ex) contribution to $`\mu `$. This allows us to rewrite the equality of the chemical potentials of the two daughter phases as, $`h_1(\sigma )=h_2(\sigma )A(\sigma )`$, where $`A(\sigma )`$ is a shorthand notation for:

$$A(\sigma )=\frac{\rho _2}{\rho _1}\mathrm{exp}\beta \left\{\mu _{ex}(\sigma ;[\rho _2])\mu _{ex}(\sigma ;[\rho _1])\right\}$$
(1)
with $`\beta =1/k_BT`$. The polydispersity distributions are further constrained by the relation, $`x_1h_1(\sigma )+x_2h_2(\sigma )=h_0(\sigma )`$, which expresses particle number conservation. The number concentration of phase 1, $`x_1=1x_2`$, is given by the lever rule: $`x_1=\frac{\rho _1}{\rho _1\rho _2}\frac{\rho _0\rho _2}{\rho _0}`$. Combining these two relations one finds:

$$h_2(\sigma )h_1(\sigma )=h_0(\sigma )H(\sigma )$$
(2)
where $`H(\sigma )(1A(\sigma ))/(x_2+x_1A(\sigma ))`$. Eq.(2) is the starting point to relate the difference between the moments of the daughter phases, $`\mathrm{\Delta }_k=𝑑\sigma \sigma ^k(h_2(\sigma )h_1(\sigma ))`$, to the moments, $`\xi _k=𝑑\sigma \sigma ^kh_0(\sigma )`$ ($`k=0,1,2,\mathrm{}`$), of the parent phase distribution $`h_0(\sigma )`$. Indeed, when $`\sigma `$ is chosen such that $`h_0(\sigma )`$ tends to the Dirac delta function $`\delta (\sigma )`$ in the monodisperse limit, $`\mathrm{\Delta }_k`$ can be obtained from (2) by expanding $`H(\sigma )`$ around $`\sigma =0`$, $`H(\sigma )=_{l=0}^{\mathrm{}}a_l\sigma ^l`$, yielding for a weakly polydisperse system, $`\mathrm{\Delta }_k=_{l=0}^{\mathrm{}}a_l\xi _{l+k}`$. The normalization of the $`h_i(\sigma )`$ ($`i=0,1,2`$) implies $`\mathrm{\Delta }_0=0`$, $`\xi _0=1`$ or $`a_0=_{l=1}^{\mathrm{}}a_l\xi _l`$, and eliminating $`a_0`$ from $`\mathrm{\Delta }_k`$ yields the general moment relation:

$$\mathrm{\Delta }_k=a_1\xi _{k+1}+\underset{l=2}{\overset{\mathrm{}}{}}a_l(\xi _{k+l}\xi _l\xi _k).$$
(3)
where we took moreover into account that $`\sigma `$ can always be chosen such that $`\xi _1=0`$. When only the first term in the r.h.s. of (3) is retained we recover the universal law $`\mathrm{\Delta }_k/\mathrm{\Delta }_l=\xi _{k+1}/\xi _{l+1}`$, put forward in ref.4. The question left unanswered by the study of ref.4 concerns the radius of convergence of the weak polydispersity expansion (3). In order to study this problem in more detail we now consider a simple model system for which we can determine the $`h_i(\sigma )(i=1,2)`$ numerically and compare the results with (3). The free energy density functional chosen here corresponds to a simple van der Waals (vdW) model for the liquid-vapor transition in polydisperse systems of spherical particles of variable size:

$`f[\rho ]=`$ $`k_BT{\displaystyle 𝑑\sigma \rho (\sigma )\left\{\mathrm{ln}(\frac{\mathrm{\Lambda }^3(\sigma )\rho (\sigma )}{E[\rho ]})1\right\}}`$ (4)
$`+{\displaystyle \frac{1}{2}}{\displaystyle 𝑑\sigma 𝑑\sigma ^{}V(\sigma ,\sigma ^{})\rho (\sigma )\rho (\sigma ^{})}`$
where, $`E[\rho ]=1𝑑\sigma v(\sigma )\rho (\sigma )`$, describes the average excluded volume correction for particles of radius $`R_\sigma `$ and volume $`v(\sigma )=\frac{4\pi }{3}R_\sigma ^3`$, while $`V(\sigma ,\sigma ^{})=𝑑𝐫V(r;\sigma ,\sigma ^{})`$ is the integrated attraction between two particles of species $`\sigma `$ and $`\sigma ^{}`$, for which we took the usual vdW form, $`V(r;\sigma ,\sigma ^{})=ϵ_0(R_\sigma +R_\sigma ^{})^6/r^6`$ for $`rR_\sigma +R_\sigma ^{}`$ and zero otherwise, $`ϵ_0`$ being the amplitude of the attraction at the contact of the two particles. The size-polydispersity can be described in terms of the dimensionless variable, $`\sigma =R_\sigma /R1`$, with $`R`$ the mean value of $`R_\sigma `$ in the parent phase, hence $`\xi _1=𝑑\sigma \sigma h_0(\sigma )=0`$. The thermodynamics is given in terms of $`h_0(\sigma )`$, the dimensionless temperature $`t=k_BT/ϵ_0`$ and the dimensionless density $`\eta =v_0\rho `$, with $`v_0=\frac{4\pi }{3}R_0^3`$ and $`R_0`$ the value of $`R_\sigma `$ in the monodisperse limit. The coexistence conditions are integral equations which can be solved numerically using, for instance, an iterative algorithm for any $`t`$, $`\eta _0=v_0\rho _0`$ and $`h_0(\sigma )`$. For $`h_0(\sigma )`$ we took a Schulz distribution with zero mean. The normalized distribution is given, for $`1\sigma <\mathrm{}`$, by $`h_0(\sigma )=\alpha ^\alpha (1+\sigma )^{\alpha 1}e^{\alpha (1+\sigma )}/\mathrm{\Gamma }(\alpha )`$, with $`\mathrm{\Gamma }(\alpha )`$ the gamma function and $`1/\alpha `$ a width parameter which measures the distance to the monodisperse limit, $`h_0(\sigma )\delta (\sigma )`$ when $`\alpha \mathrm{}`$. We then have: $`\xi _0=1`$, $`\xi _1=0`$, $`\xi _2=1/\alpha `$, $`\xi _3=2/\alpha ^2`$, $`\xi _4=\frac{3}{\alpha ^2}+\frac{6}{\alpha ^3}`$, $`\xi _5=\frac{20}{\alpha ^3}+\frac{24}{\alpha ^4}`$, etc. For a weakly polydisperse system we retain only the dominant terms of (3) in a $`1/\alpha `$ expansion. From (3) we obtain then: $`\mathrm{\Delta }_1=a_1(\mathrm{})\xi _2+O(1/\alpha ^2)`$, $`\mathrm{\Delta }_2=a_1(\mathrm{})\xi _3+a_2(\mathrm{})(\xi _4\xi _2^2)+O(1/\alpha ^3)=\{a_1(\mathrm{})+a_2(\mathrm{})\}\xi _3+O(1/\alpha ^3)`$, $`\mathrm{\Delta }_3=a_1(\mathrm{})\xi _4+O(1/\alpha ^3)`$, etc, where $`a_l(\mathrm{})`$ are the values of $`a_l`$ for $`\alpha \mathrm{}`$. Using the vdW expression (4) to evaluate (1) one finds, for ex. for $`t=1.0`$ and $`\eta _0=0.5`$, $`a_1(\mathrm{})=1.75`$ and $`a_2(\mathrm{})=2.68`$. Using the corresponding numerical solutions found for $`h_1(\sigma )`$ and $`h_2(\sigma )`$ (see Fig.1) it can be seen from Fig.2 that $`\mathrm{\Delta }_1/\xi _21.75`$, $`\mathrm{\Delta }_2/\xi _30.93`$ and $`\mathrm{\Delta }_3/\xi _41.75`$ are obeyed to within ten percent for $`\alpha `$ larger than, respectively, 40, 80 and 150. We can conclude thus that the weak polydispersity expansion (3) is valid (to dominant order) for Schulz distributions $`h_0(\sigma )`$ with a dispersion $`\left((\xi _2\xi _1^2)^{1/2}\right)`$ smaller than, say, 0.1 ($`\alpha 100`$). These values do of course depend on the thermodynamic state but the case considered here ($`t=1`$, $`\eta _0=0.5`$) is representative of other $`t,\eta _0`$ values. Note also that we have verified numerically that the radius of convergence of (3) with respect to $`1/\alpha `$ is fairly sensitive to the total amount of polydispersity present. Allowing, for instance, the amplitude $`ϵ_0`$ of the pair potential $`V(r;\sigma ,\sigma ^{})`$ to depend on $`\sigma `$ and $`\sigma ^{}`$ does reduce the radius of convergence of (3) considerably. From the above it follows that, $`\frac{\mathrm{\Delta }_3}{\mathrm{\Delta }_1}`$ follows the universal law, $`\frac{\mathrm{\Delta }_3}{\mathrm{\Delta }_1}=\frac{\xi _4}{\xi _2}`$, put forward in ref.4 whereas $`\frac{\mathrm{\Delta }_2}{\mathrm{\Delta }_1}`$ follows the non-universal law, $`\frac{\mathrm{\Delta }_2}{\mathrm{\Delta }_1}=\{1+\frac{a_2(\mathrm{})}{a_1(\mathrm{})}\}\frac{\xi _3}{\xi _2}`$. We have verified that similar results can be obtained for different $`h_0(\sigma )`$ distributions. Taking, for instance, a Gaussian for $`h_0(\sigma )`$ similar results are found although $`\xi _3=0`$ for this case. This invalidates the conclusion of ref.4 that a particular importance should be attached to the skewness of $`h_0(\sigma )`$. In conclusion, the general moment relation (3) can yield useful information about the phase behavior of weakly polydisperse systems but this information is in general not universal.

Figure Captions

FIG. 1. The polydispersity distributions $`h_n(\sigma )`$ of the parent phase ($`n=0`$: full curve) (a Schulz distribution with the width parameter $`\alpha =50`$), the low-density (n=1: dotted curve) and the high-density (n=2: circles) daughter phases, as obtained by numerically solving the coexistence conditions of the van der Waals model of eq.(4) for $`t=1`$, $`\eta _0=0.5`$. The corresponding dimensionless densities of the coexisting daughter phases are $`\eta _1=0.106`$, $`\eta _2=0.521`$ whereas for the monodisperse system one has, $`\eta _1=0.103`$, $`\eta _2=0.608`$. Also shown are $`h_1(\sigma )h_0(\sigma )`$ (dashed curve) and $`[h_2(\sigma )h_0(\sigma )]50`$ (triangles).

FIG. 2. The ratio $`\mathrm{\Delta }_k/\xi _{k+1}`$ ($`k=1,2,3`$) versus $`1/\alpha `$ as obtained from the numerical solution of the van der Waals model of eq.(4) for $`t=1`$, $`\eta _0=0.5`$ and a Schulz distribution for $`h_0(\sigma )`$. The symbols are as follows: circles(k=1), squares(k=2) and triangles(k=3). The dotted lines indicate their asymptotic ($`\alpha \mathrm{}`$) values. The arrows indicate for each case the radius of convergence of the weak polydispersity expansion of eq.(3).