text:
# Moment free energies for polydisperse systems

## I Introduction

The thermodynamics of mixtures of several chemical species is, since Gibbs, a well established subject (see e.g. ). But many systems arising in nature and in industry contain, for practical purposes, an infinite number of distinct, though similar, chemical species. Often these can be classified by a parameter, $`\sigma `$, say, which could be the chain length in a polymeric system, or the particle size in a colloid; both are routinely treated as continuous variables. In other cases (see e.g. ) $`\sigma `$ is instead a parameter distinguishing species of continuously varying chemical properties. In the most general case, several attributes may be required to distinguish the various particle species in the system (such as length and chemical composition in length-polydisperse random copolymers) and $`\sigma `$ is then a collection of parameters . The thermodynamics of polydisperse systems (as defined above) is of crucial interest to wide areas of science and technology; it is sometimes also referred to as “continuous thermodynamics” (see e.g. ).

Standard thermodynamic procedures  for constructing phase equilibria in a system of volume $`V`$ containing $`M`$ different species can be understood geometrically in terms of a free energy surface $`f(\rho _j)`$ (with $`f=F/V`$, and $`F`$ the Helmoltz free energy) in the $`M`$-dimensional space of number densities $`\rho _j`$. Tangent planes to $`f`$ define regions of coexistence, within which the free energy of the system is lowered by phase separation. The volumes of coexisting phases follow from the well-known “lever rule”. Here “surface” and “plane” are used loosely, to denote manifolds of appropriate dimension. This procedure becomes unmanageable, both conceptually and numerically, in the limit $`M\mathrm{}`$ which formally defines a polydisperse system. There is now a separate conserved density $`\rho (\sigma )`$ for each value of $`\sigma `$; $`\rho (\sigma )`$ is in fact a density distribution and the overall number density of particles is written as $`\rho =𝑑\sigma \rho (\sigma )`$. The free energy surface is $`f=f[\rho (\sigma )]`$ (a functional) which resides in an infinite-dimensional space. Gibbs’ rule allows the coexistence of arbitrarily many thermodynamic phases.

To make phase equilibria in polydisperse systems more accessible to both computation and physical intuition, it is clearly helpful to reduce the dimensionality of the problem. Theoretical work in this area has made significant progress in finding simplified forms of the conditions for phase coexistence, thus making these more numerically tractable . Our aim is to achieve a similar reduction in dimensionality on a higher level, that of the free energy itself. We show that, for a large class of models, it is possible to construct a reduced free energy depending on only a small number of variables. From this, meaningful information on phase equilibria can be extracted by the usual tangent plane procedure, with obvious benefits for a more intuitive understanding of the phase behavior of polydisperse systems. As an important side effect, this procedure also leads to robust algorithms for calculating polydisperse phase equilibria numerically. In particular, such algorithms can handle coexistence of more than two phases with relative ease compared to those used previously.

A clue to the choice of independent variables for such a reduced free energy comes from the representation of experimental data. We recall first the definition of a “cloud point” (see e.g. ): This is the point at which, for a system with a given density distribution $`\rho (\sigma )`$, phase separation first occurs as the temperature $`T`$ or another external control parameter is varied. The corresponding incipient phase is called the “shadow”. Now consider diluting or concentrating the system, i.e., varying its overall density $`\rho `$ while maintaining a fixed “shape” of polydispersity $`n(\sigma )=\rho (\sigma )/\rho `$. We will find it useful later to refer to the collection of all systems $`\rho (\sigma )=\rho n(\sigma )`$ which can be obtained by this process as a “dilution line” (in the space of all density distributions $`\rho (\sigma )`$). Plotting the cloud point temperature $`T`$ versus density $`\rho `$ then defines the cloud point curve (CPC), while plotting $`T`$ versus the density of the shadow gives the shadow curve . The differences in the shapes $`n(\sigma )`$ of the density distributions in the different phases are hidden in this representation; only the overall densities $`\rho `$ in each phase are tracked. The density $`\rho `$ is a particular moment of $`\rho (\sigma )`$ (the zeroth one); higher order moments would be given by $`𝑑\sigma \sigma ^i\rho (\sigma )`$. Generalizing slightly, this suggests that our reduced free energy should depend on several (generalized) moment densities $`\rho _i=𝑑\sigma w_i(\sigma )\rho (\sigma )`$, defined by appropriate (linearly independent) weight functions $`w_i(\sigma )`$. Ordinary power-law moment densities are included as the special case $`w_i(\sigma )=\sigma ^i`$. Whatever the choice of moments, we insist on $`w_0=1`$ so that in all cases the zeroth moment density $`\rho _0`$ coincides with the overall number density $`\rho `$.

As an illustration, consider the simplest imaginable case where the true free energy $`f`$ already has the required form, i.e., it depends only on a finite set of $`K`$ moment densities of the density distribution:

$$f=f(\rho _i),i=1\mathrm{}K\mathrm{o}ri=0\mathrm{}K1$$
(1)
where the range of $`i`$ depends on whether $`\rho _0\rho `$ is among the $`K`$ moment densities on which $`f`$ depends. (From now on, we use the notation $`i=1\mathrm{}K`$ inclusively, to cover both possibilities.) In coexisting phases one demands equality of particle chemical potentials, defined as $`\mu (\sigma )=\delta f/\delta \rho (\sigma )`$, for all $`\sigma `$. Because

$`\mu (\sigma ){\displaystyle \frac{\delta f}{\delta \rho (\sigma )}}={\displaystyle \underset{i}{}}{\displaystyle \frac{f}{\rho _i}}w_i(\sigma )={\displaystyle \underset{i}{}}\mu _iw_i(\sigma )`$

this implies that all “moment” chemical potentials, $`\mu _if/\rho _i`$, are likewise equal among phases. The second requirement for phase coexistence is that the pressures or osmotic pressures , $`\mathrm{\Pi }`$, of all phases must also be equal. But from

$`\mathrm{\Pi }=f{\displaystyle 𝑑\sigma \mu (\sigma )\rho (\sigma )}=f{\displaystyle \underset{i}{}}\mu _i\rho _i`$

one sees that this again involves only the moment densities $`\rho _i`$ and their chemical potentials $`\mu _i`$. Finally, the moment densities also obey the “lever rule”, as follows. Let the overall density distribution of a system of volume $`V`$ be $`\rho ^{(0)}(\sigma )`$; we call this the “parent” distribution. If (after a lowering of temperature, for example) this parent has split into $`p`$ coexisting “daughter” phases with $`\sigma `$-distributions $`\rho ^{(\alpha )}(\sigma )`$, each occupying a fraction $`v^{(\alpha )}`$ of the total volume ($`\alpha =1\mathrm{}p`$), then particle conservation implies the usual lever rule (or material balance) among species: $`_\alpha v^{(\alpha )}\rho ^{(\alpha )}(\sigma )=\rho ^{(0)}(\sigma ),\sigma `$. Multiplying this by a weight function $`w_i(\sigma )`$ and integrating over $`\sigma `$ shows that the lever rule also holds for the moment densities:

$$\underset{\alpha =1}{\overset{p}{}}v^{(\alpha )}\rho _i^{(\alpha )}=\rho _i^{(0)}$$
(2)
These results express the fact that any linear combination of conserved densities (a generalized moment density) is itself a conserved density in thermodynamics. We have shown, therefore, that if the free energy of the system depends only on $`K`$ moment densities $`\rho _1\mathrm{}\rho _K`$ we can view these as the densities of $`K`$ “quasi-species” of particles, and construct the phase diagram via the usual construction of tangencies and the lever rule. Formally this has reduced the problem to finite dimensionality, although this is trivial here because $`f`$, by assumption, has no dependence on any variables other than the $`\rho _i`$ ($`i=1\mathrm{}K`$).

Of course, it is uncommon for the free energy $`f`$ to obey (1). In particular, the entropy of an ideal mixture (or, for polymers, the Flory-Huggins entropy term) is definitely not of this form. On the other hand, in very many thermodynamic (especially mean field) models the excess (i.e., non-ideal) part of the free energy does have the simple form (1). In other words, if we decompose the free energy as (setting $`k_\mathrm{B}=1`$)

$$f=Ts_{\mathrm{id}}+\stackrel{~}{f},s_{\mathrm{id}}=𝑑\sigma \rho (\sigma )\left[\mathrm{ln}\rho (\sigma )1\right]$$
(3)
then the excess free energy $`\stackrel{~}{f}`$ is a function of some moment densities only:

$$\stackrel{~}{f}=\stackrel{~}{f}(\rho _i),i=1\mathrm{}K$$
(4)
Examples of models of this kind include polydisperse hard spheres (within the generalization by Salacuse and Stell  of the BMCSL equation of state ), polydisperse homo- and copolymers , and van der Waals fluids with factorized interaction parameters . With the exception of a brief discussion in Sec. VI, this paper concerns models with free energies of the form (3,4), which we call “truncatable”. (This terminology emphasizes that the number $`K`$ of moment densities appearing in the excess free energy of truncatable models is finite, while for a non-truncatable model the excess free energy depends on all details of $`\rho (\sigma )`$, corresponding to an infinite number of moment densities.) In what follows, we regard each model free energy as given, and do not discuss the issue of how good a description of the real system it offers, nor how or whether it can be derived from an underlying microscopic Hamiltonian. Whenever we refer to “exact” results, we mean the exact thermodynamics of such a model as specified by its free energy.

Different authors refer in different terms to $`s_{\mathrm{id}}`$ in (3). Throughout this paper will call the quantity

$`s_{\mathrm{id}}={\displaystyle 𝑑\sigma \rho (\sigma )\left[\mathrm{ln}\rho (\sigma )1\right]}`$

(which is, up to a factor of $`T`$, the ideal part of the free energy density) the “entropy of an ideal mixture”. By writing $`\rho (\sigma )=\rho n(\sigma )`$, where $`n(\sigma )`$ is the normalized distribution of $`\sigma `$, this can be decomposed as

$$s_{\mathrm{id}}=\rho (\mathrm{ln}\rho 1)\rho 𝑑\sigma n(\sigma )\mathrm{ln}n(\sigma )$$
(5)
The first term on the right hand side is the “entropy of an ideal gas”, while the second term gives the “entropy of mixing”. The prefactor $`\rho `$ reflects the fact that both are expressed per unit volume rather than per particle.

In principle, the entire phase equilibria for any truncatable system (obeying (4)) can be computed exactly by a finite algorithm. Specifically, the spinodal stability criterion involves a $`K`$-dimensional square matrix  whereas calculation of $`p`$-phase equilibrium involves solution of $`(p1)(K+1)`$ coupled nonlinear equations (see Sec. III E). This method has certainly proved useful , but is cumbersome, particularly if one is interested mainly in cloud- and shadow-curves, rather than coexisting compositions deep within multi-phase regions . Various ways of simplifying the procedure exist , but there has previously not been a systematic alternative to the full computation. Note also that the nonlinear phase equilibrium equations permit no simple geometrical interpretation or qualitative insight akin to the familiar rules for constructing phase diagrams from the free energy surface of a finite mixture. This further motivates the search for a description of polydisperse phase equilibria in terms of a reduced free energy which depends on only a small number of density variables.

In previous work, the authors originally arrived independently at two definitions of a reduced free energy in terms of moment densities (moment free energy for short) . Though based on distinctly different principles, the two approaches led to very similar results. We explain this somewhat surprising fact in the present work; at the same time, we describe the two methods in more detail and explore the relationship between them. We also discuss issues related to practical applications and give a number of example results for simple model systems.

The first route to a moment free energy takes as its starting point the conventional form of the ideal part of the free energy, $`f_{\mathrm{id}}=T𝑑\sigma \rho (\sigma )[\mathrm{ln}\rho (\sigma )1]`$, as in (3). This can be thought as defining a (hyper-)surface, with the “horizontal” co-ordinate being $`\rho (\sigma )`$ and the “height” of the surface $`f_{\mathrm{id}}[\rho (\sigma )]`$. As explained in Sec. II A, this surface can then be projected geometrically onto one with only a finite set of horizontal co-ordinates: these are chosen to be the moment densities appearing in the excess free energy. Physically, this corresponds to minimizing the free energy with respect to all degrees of freedom in $`\rho (\sigma )`$ that do not affect these moment densities. We call this first route the “projection” method.

The second approach rederives the entropy of mixing in the ideal part of the free energy in a form that depends explicitly only on the chosen moment densities. As described in Sec. II B, the expression that results is intractable in general because it still contains the full complexity of the problem. However, in situations where there are only infinitesimal amounts of all but one of the phases in the systems, the entropy of mixing can be evaluated in a closed form. Applying this functional form in regimes where it is not strictly valid (i.e., when phases of comparable volumes coexist) generates a moment free energy by this “combinatorial” method.

We show the equivalence of the two approaches in Sec. II C. There, we first demonstrate that the general form of the entropy of mixing obtained by the combinatorial method can be transformed to the standard expression $`𝑑\sigma n(\sigma )\mathrm{ln}n(\sigma )`$. Second, we show that the moment free energies arrived at by the two methods are in fact equal, with the projection method being slightly more generally applicable.

In Sec. III, we then discuss the properties of the moment free energy (as obtained from either the projection or the combinatorial method). By construction, it depends only on a (finite) number of moment densities, and so achieves the desired reduction in dimensionality. The moment densities can be treated as densities of “quasi-species” of particles, and the standard procedures of the thermodynamics of finite mixtures can be applied to the moment free energy to calculate phase equilibria. This is only useful, however, if the results are faithful to the actual phase behavior of the underlying polydisperse system (as modeled by the given truncatable free energy). We show that this is so: In fact, the construction of our moment free energy is such that exact binodals (cloud-point and shadow curves), critical (and multi-critical) points and spinodals are obtained. Beyond the onset of phase separation, where coexisting phases occupy comparable volumes, the results are not exact, but can be refined arbitrarily by adding extra moment densities. This procedure is necessary also to ensure that, in regions of multi-phase coexistence, the correct number of phases is found. In Sec. IV, we discuss the practical implementation of our method, followed by a number of examples (Sec. V). Our results are summarized in Sec. VI, where we also outline perspectives for future work.

Note that throughout this work, we focus on the case of phase coexistence at fixed volume. However, as described in App. A, the formalism can easily be applied to scenarios where the (mechanical or osmotic) pressure is fixed instead.

## II Derivation of moment free energy

### A Projection method

The starting point for this method is the decomposition (3) for the free energy of truncatable systems (4):

$$f=T𝑑\sigma \rho (\sigma )\left[\mathrm{ln}\frac{\rho (\sigma )}{R(\sigma )}1\right]+\stackrel{~}{f}(\rho _i)$$
(6)
As explained in the previous section, truncatable here means that the excess free energy $`\stackrel{~}{f}`$ depends only on $`K`$ moment densities $`\rho _i`$. Note that, in the first (ideal) term of (6), we have included a dimensional factor $`R(\sigma )`$ inside the logarithm. This is equivalent to subtracting $`T𝑑\sigma \rho (\sigma )\mathrm{ln}R(\sigma )`$ from the free energy. Since this term is linear in densities, it has no effect on the exact thermodynamics; it contributes harmless additive constants to the chemical potentials $`\mu (\sigma )`$. However, in the projection route to a moment free energy, it will play a central role.

We now argue that the most important moment densities to treat correctly are those that actually appear in the excess free energy $`\stackrel{~}{f}(\rho _i)`$. Accordingly we divide the infinite-dimensional space of density distributions into two complementary subspaces: a “moment subspace”, which contains all the degrees of freedom of $`\rho (\sigma )`$ that contribute to the moment densities $`\rho _i`$, and a “transverse subspace” which contains all remaining degrees of freedom (those that can be varied without affecting the chosen moment densities $`\rho _i`$. Physically, it is reasonable to expect that these “leftover” degrees of freedom play a subsidiary role in the phase equilibria of the system, a view justified a posteriori below. Accordingly, we now allow violations of the lever rule, so long as these occur solely in the transverse space. This means that the phase splits that we calculate using this approach obey the lever rule for the moment densities, but are allowed to violate it in other details of the density distribution $`\rho (\sigma )`$ (see Fig. 10 in Sec. V B 1 for an example). These “transverse” degrees of freedom, instead of obeying the strict particle conservation laws, are chosen so as to minimize the free energy: they are treated as “annealed”. If, as assumed above, $`\stackrel{~}{f}=\stackrel{~}{f}(\rho _i)`$ only depends on the moment densities retained, this amounts to maximizing the ideal mixture entropy ($`s_{\mathrm{id}}`$) in (6), while holding fixed the values of the moment densities $`\rho _i`$. Note that we are allowed, if we wish, to include among the retained densities “redundant” moments on which $`\stackrel{~}{f}`$ has a null dependence. We will have occasion to do this later on.

At this point, the factor $`R(\sigma )`$ in (6), which is immaterial if all conservation laws are strictly obeyed, becomes central. Indeed, maximizing the entropy over all distributions $`\rho (\sigma )`$ at fixed moment densities $`\rho _i`$ yields

$$\rho (\sigma )=R(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right)$$
(7)
where the Lagrange multipliers $`\lambda _i`$ are chosen to give the desired moment densities

$$\rho _i=𝑑\sigma w_i(\sigma )R(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right)$$
(8)
The corresponding minimum value of $`f`$ then defines our projected free energy as a function of the moment densities $`\rho _i`$:

$$f_{\mathrm{pr}}(\rho _i)=Ts_{\mathrm{pr}}(\rho _i)+\stackrel{~}{f}(\rho _i),s_{\mathrm{pr}}=\rho _0\underset{i}{}\lambda _i\rho _i$$
(9)
Here $`s_{\mathrm{pr}}`$ is the projected entropy of an ideal mixture. The first term appearing in it, $`\rho _0=𝑑\sigma \rho (\sigma )`$, is the “zeroth moment”, which is identical to the overall particle density $`\rho `$ defined previously. If this is among the moment densities used for the projection (or more generally, if it is a linear combination of them), then the term $`T\rho _0`$ is simply a linear contribution to the projected free energy $`f_{\mathrm{pr}}(\rho _i)`$ and can be dropped because it does not affect phase equilibrium calculations. Otherwise, $`\rho _0`$ needs to be expressed – via the $`\lambda _i`$ – as a function of the $`\rho _i`$ and its contribution cannot be ignored. We will see an example of this in Sec. V.

The projection method yields a free energy $`f_{\mathrm{pr}}(\rho _i)`$ which only depends on the chosen set of moment densities. These can now be viewed as densities of ‘quasi-species’ of particles, and a finite-dimensional phase diagram can be constructed from $`f_{\mathrm{pr}}`$ according to the usual rules, ignoring the underlying polydisperse nature of the system. Obviously, though, the results now depend on $`R(\sigma )`$ which is formally a “prior distribution” for the entropy maximization. To understand its role, we recall that the projected free energy $`f_{\mathrm{pr}}(\rho _i)`$ was constructed as the minimum of $`f[\rho (\sigma )]`$ at fixed $`\rho _i`$; that is, $`f_{\mathrm{pr}}`$ is the lower envelope of the projection of $`f`$ onto the moment subspace. Crucially, the shape of this envelope depends on how, by choosing a particular prior distribution $`R(\sigma )`$, we “tilt” the infinite-dimensional free energy surface before projecting it. This geometrical point of view is illustrated, for a mixture of only two species, in Fig. 1.

To understand the effect of the prior $`R(\sigma )`$ physically, we note that the projected free energy is simply the free energy of phases of a system in which the density distributions $`\rho (\sigma )`$ are of the form (7). The prior $`R(\sigma )`$ determines which distributions lie within this “maximum entropy family” (or “family” for short), and it is the properties of phases with these distributions that the projected free energy represents. Typically, one is interested in a system where a fixed overall “parent” (or “feed”) distribution $`\rho ^{(0)}(\sigma )`$ becomes subject to separation into various phases. In such circumstances, we should generally choose this parent distribution as our prior, $`R(\sigma )=\rho ^{(0)}(\sigma )`$, thereby guaranteeing that it is contained within the family (7). Having done this, we note that the projection procedure will be exactly valid, to whatever extent the density distributions actually arising in the various coexisting phases of the system under study are members of the corresponding family

$$\rho (\sigma )=\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right).$$
(10)
In fact, the condition just described holds whenever all but one of a set of coexisting phases are of infinitesimal volume compared to the majority phase. This is because the density distribution, $`\rho ^{(0)}(\sigma )`$, of the majority phase is negligibly perturbed, whereas that in each minority phase differs from this by a Gibbs-Boltzmann factor, of exactly the form required for (10); we show this formally in Sec. III. Accordingly, our projection method yields exact cloud point and shadow curves. By the same argument, critical points (which in fact lie at the intersection of these two curves) are exactly determined; the same is true for tricritical and higher order critical points. Finally, spinodals are also found exactly. We defer explicit proofs of these statements to Sec. III.

The projection method does, however, give only approximate results for coexistences involving finite amounts of different phases. This is because linear combinations of different density distributions from the family (10), corresponding to two (or more) phases arising from the same parent $`\rho ^{(0)}(\sigma )`$, do not necessarily add to recover the parent distribution itself:

$$\underset{\alpha }{}v^{(\alpha )}\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _i^{(\alpha )}w_i(\sigma )\right)\rho ^{(0)}(\sigma )$$
(11)
unless all except one of the $`v^{(\alpha )}`$ are infinitesimal. Moreover, according to Gibbs’ phase rule, a projected free energy depending on $`n`$ moment densities will not normally predict more than $`n+1`$ coexisting phases, whereas a polydisperse system can in principle separate into an arbitrary number of phases. We explain in Sec. III how both of these shortcomings can be overcome by systematically including extra moment densities within the projection procedure. How quickly convergence to the exact results occurs depends on the choice of weight functions for the extra moment densities; we discuss this point further in the context of the examples of Sec. V.

### B Combinatorial method

Now we turn to the second method of constructing a moment free energy. As before, we recognize from the outset that the physics of the problem is contained in the excess free energy and therefore the most important variables are the moment densities which feature in this. On the other hand, the ideal part of the free energy (essentially the entropy of mixing, up to the trivial ideal gas term) is a book-keeping device that accounts for the number of ways of partitioning a priori a given size distribution between two or more coexisting phases. (To focus the discussion, we identify $`\sigma `$ as particle size in this section; the arguments do of course remain valid for a general polydispersity parameter $`\sigma `$.) The entropy of mixing has its origin in $`1/N!`$ factors in the partition function (where $`N`$ is a particle number) which are usually derived from the classical limit of quantum statistics. The origin of the $`1/N!`$ factors in classical statistical mechanics has been the subject of debate ever since Gibbs . The problem is also connected with the question of the extensivity of entropy, and the Gibbs paradox. Below we present a completely classical derivation of the entropy of mixing which has the additional benefit of indicating the appropriate generalisation for moment densities. To start with, follow Gibbs and define a non-extensive free energy $`F^{}`$ as a configuration space ($`\mathrm{\Gamma }`$) integral

$$e^{F^{}/T}=𝑑\mathrm{\Gamma }e^{H/T}.$$
(12)
No $`1/N!`$ appears in this as all particles are distinguishable in classical statistical mechanics.

Now consider two phases in coexistence as one joint system. Assume they occupy volumes $`V^{(1)}`$ and $`V^{(2)}`$, respectively, and contain $`N^{(1)}`$ and $`N^{(2)}`$ particles. Following literally the prescription in eq. (12), the free energy of the joint system of $`N^{(1)}+N^{(2)}=N`$ particles is found from

$$e^{F^{}/T}=\underset{\mathrm{prtns}}{}e^{(F_{}^{}{}_{}{}^{(1)}+F_{}^{}{}_{}{}^{(2)})/T}$$
(13)
where the configuration space integral has been done in two parts. Firstly, for each way of partitioning the particles between the two phases, the individual configuration space integrals give a product of the individual partition functions. Secondly, and crucially, one must sum over the $`N!/N^{(1)}!N^{(2)}!`$ ways of partitioning the particles between phases. Now define an extensive free energy by reinserting the $`1/N!`$ as though the particles were indistinguishable:

$$e^{F_{\mathrm{in}}/T}=\frac{1}{N!}𝑑\mathrm{\Gamma }e^{H/T},$$
(14)
Eq. (13) can then be written as

$$e^{F_{\mathrm{in}}/T}=e^{(F_{\mathrm{in}}^{(1)}+F_{\mathrm{in}}^{(2)})/T}_{\mathrm{prtns}}$$
(15)
where the average is taken over all partitions, with equal a priori probabilities. This result is the cornerstone of the combinatorial method. If we separate $`F_{\mathrm{in}}`$ into the entropy of an ideal gas and the remaining non-ideal (excess) parts, it is written as

$$F_{\mathrm{in}}=NT[\mathrm{ln}(N/V)1]+\stackrel{~}{F}$$
(16)
with $`\stackrel{~}{F}`$ being the excess free energy. Note that $`F_{\mathrm{in}}`$ does not yet contain an entropy of mixing term because it treats the particles as though they are indistinguishable. The entropy of mixing will reappear when we come to do the average over partitions in (15). Note also that conventional thermodynamics follows from (15) if all particles are identical as far as their mutual interactions are concerned. There is then no entropy of mixing to consider, so that $`F_{\mathrm{in}}`$ is just the conventional free energy $`F`$; and since the average over partitions is trivial one has $`F=F^{(1)}+F^{(2)}`$. Similarly, as we show in Sec. II C, one can recover from (15) the conventional form of the entropy of mixing as given in (5).

We now show how the average over partitions in (15) results in an expression for the entropy of mixing which depends explicitly on moment variables. The key is to note that, by our assumption that we are dealing with a truncatable model, the excess free energy $`\stackrel{~}{F}=V\stackrel{~}{f}(\rho _i)`$ depends only on a limited number of moment densities $`\rho _i`$ (as well as volume $`V`$, temperature $`T`$, and other state variables which we suppress below). If the density $`\rho _0`$ is among the moment densities $`\rho _i`$ – which we assume throughout this section – then specifying $`V`$ and the $`\rho _i`$ is equivalent to specifying $`N`$, $`V`$ and the normalized moments $`m_i=\rho _i/\rho _0`$ ($`i>0`$). So we can think of $`\stackrel{~}{F}`$ as a function of $`N`$, $`V`$ and the $`m_i`$. From (16), $`F_{\mathrm{in}}`$ then depends on the same set of variables. The choice of the $`m_i`$, which are moments of the normalized particle size distribution $`n(\sigma )`$, as independent variables is natural in the present context because we are considering the particle number $`N`$ and volume $`V`$ of the coexisting phases to be fixed.

To avoid notational complexity, and without loss of generality, we now specialize to the case where there is only one normalized moment $`m_i`$, a generalized mean size

$`m={\displaystyle 𝑑\sigma w(\sigma )n(\sigma )}`$

We shall also write $`\rho _m=\rho m=𝑑\sigma w(\sigma )\rho (\sigma )`$ for the corresponding moment density, and $`m^{(1)}`$ and $`m^{(2)}`$ for the values of the generalized mean size in the first and second phase. The overall (“parent”) size distribution is $`n^{(0)}(\sigma )`$ with mean $`m^{(0)}`$. In each of the two phases, in which $`N`$ and $`V`$ are held fixed, $`F_{\mathrm{in}}`$ only depends on $`m`$, and so the average over partitions becomes

$$e^{(F_{\mathrm{in}}^{(1)}+F_{\mathrm{in}}^{(2)})/T}_{\mathrm{prtns}}𝑑m^{(1)}𝒫(m^{(1)})e^{(F_{\mathrm{in}}^{(1)}+F_{\mathrm{in}}^{(2)})/T}$$
(17)
Here $`𝒫(m^{(1)})`$ is the probability distribution for the generalized mean size in the first phase, taken over partitions with fixed $`N^{(1)}`$ and $`N^{(2)}`$, with equal a priori probabilities. Note that given $`m^{(1)}`$, $`m^{(2)}`$ is fixed in the second phase by the moment equivalent of particle conservation: $`N^{(1)}m^{(1)}+N^{(2)}m^{(2)}=Nm^{(0)}`$. The integral in (17) can be replaced by the maximum of the integrand in the thermodynamic limit, because $`\mathrm{ln}𝒫(m^{(1)})`$ is an extensive quantity. Introducing a Lagrange multiplier $`\mu _m`$ for the above moment constraint then shows that the quantity $`\rho _m`$ has the same status as the density $`\rho \rho _0`$ itself: both are thermodynamic density variables. This reinforces the discussion in the introduction, where we showed that moment densities can be regarded as densities of “quasi-species” of particles.

In what follows, we only need to refer to $`m^{(1)}`$ (not $`m^{(2)}`$) and therefore drop the superscript for brevity. Since $`\mathrm{ln}𝒫(m)`$ is extensive, we can write $`\mathrm{ln}𝒫(m)=Ns(m)`$ where $`s`$ is the entropy of mixing per particle expressed as a function of the moment $`m`$ (or, more generally, of the full set of normalized moments). The total free energy of the system then takes the form

$$F_{\mathrm{in}}=\underset{m}{\mathrm{min}}F_{\mathrm{in}}^{(1)}+F_{\mathrm{in}}^{(2)}NTs(m),$$
(18)
and we need to calculate the entropy of mixing $`s(m)`$. We recognize that this quantity depends not only on the generalized mean size, but also on $`x=N^{(1)}/N`$, the fraction of particles in the first phase. A formal result for $`s(m)`$ can be obtained by first calculating the joint probability distribution $`𝒫(xm,x)`$ of $`xm`$ and $`x`$. From this, we can find $`𝒫(m)`$ according to

$$𝒫(m)=\frac{𝒫(xm,x)}{𝑑m^{}𝒫(xm^{},x)}=\frac{𝒫(xm,x)}{𝒫(x)/x}$$
(19)
Writing the log probabilities in (19) as $`\mathrm{ln}𝒫(xm,x)=Ns(xm,x)`$ and $`\mathrm{ln}𝒫(x)=Ns(x)`$, we then have in the thermodynamic limit

$$s(m)=s(xm,x)s(x)$$
(20)
The method for calculating $`𝒫(xm,x)`$ proceeds by introducing an indicator function $`ϵ_i`$ for each particle, deemed to be $`1`$ if the particle is in the first phase, and $`0`$ if the particle is in the second . Then $`Nx=_{i=1}^Nϵ_i`$ and $`Nxm=_{i=1}^Nϵ_iw(\sigma _i)`$. We now write the moment generating function for $`𝒫(xm,m)`$ as follows

$`{\displaystyle d(xm)𝑑x𝒫(xm,x)\mathrm{exp}[N(\theta x+\lambda xm)]}`$ $`=\mathrm{exp}\left[\theta _{i=1}^Nϵ_i+\lambda _{i=1}^Nϵ_iw(\sigma _i)\right]_{\mathrm{prtns}}`$ (23)
$`={\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{exp}\{ϵ_i[\theta +\lambda w(\sigma _i)]\}_{\mathrm{prtns}}=2^N{\displaystyle \underset{i=1}{\overset{N}{}}}\{1+\mathrm{exp}[\theta +\lambda w(\sigma _i)]\}`$
$`=2^N\mathrm{exp}\left({\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{ln}\{1+\mathrm{exp}[\theta +\lambda w(\sigma _i)]\}\right).`$
In the second equality we have used the fact that the $`ϵ_i`$ are independent discrete random variables taking the values $`0`$ and $`1`$ with a priori equal probabilities. Taking logarithms of (23) and dividing by $`N`$, we obtain on the right hand side an average over $`i=1\mathrm{}N`$ of a function of $`\sigma _i`$. Since the $`\sigma _i`$ are drawn independently from $`n^{(0)}(\sigma )`$, the law of large numbers guarantees that this average tends to the corresponding average over $`n^{(0)}(\sigma )`$ in the limit $`N\mathrm{}`$ (compare the discussion in Ref. ). We therefore have the final result

$$\frac{1}{N}\mathrm{ln}d(xm)𝑑x\mathrm{exp}\{N[s(xm,x)+\theta x+\lambda xm]\}=𝑑\sigma n^{(0)}(\sigma )\mathrm{ln}\{1+\mathrm{exp}[\theta +\lambda w(\sigma )]\}$$
(24)
where we have dropped an unimportant additive constant ($`\mathrm{ln}2`$) and have replaced $`𝒫(xm,x)`$ by $`\mathrm{exp}[Ns(xm,x)]`$.

Now, the left hand side can be evaluated, by the steepest descent method, in the limit $`N\mathrm{}`$. Defining

$`\psi (\theta ,\lambda )={\displaystyle 𝑑\sigma n^{(0)}(\sigma )\mathrm{ln}\{1+\mathrm{exp}[\theta +\lambda w(\sigma )]\}}`$

we get $`s(xm,x)+\theta x+\lambda xm=\psi (\theta ,\lambda )`$ at the point where $`s/x+\theta =s/(xm)+\lambda =0`$. We recognize that the relationship between $`s(x,xm)`$ and $`\psi (\theta ,\lambda )`$ is a double Legendre transform. Inverting this shows that

$$s(xm,x)=\psi (\theta ,\lambda )\theta x\lambda xm$$
(25)
where

$$x=\psi /\theta xm=\psi /\lambda $$
(26)
Eq. (25) is essentially the desired result; to obtain $`s(m)`$ from (20), only $`s(x)`$ remains to be determined. It is obtained from $`𝒫(x)`$ via $`\mathrm{exp}[Ns(x)]=𝒫(x)=d(xm)𝒫(xm,x)`$ which by inspection of the moment generating function is seen to correspond to the point $`\lambda =0`$. At this point, $`\psi =\mathrm{ln}(1+e^\theta )`$ and $`x=\psi /\theta =e^\theta /(1+e^\theta )`$. We find that $`\theta =\mathrm{ln}[x/(1x)]`$ and after some manipulation,

$$s(x)=\psi x\theta =x\mathrm{ln}x(1x)\mathrm{ln}(1x)$$
(27)
This is recognized as the standard entropy of mixing that is lost when a total of $`N`$ particles is partitioned into $`N^{(1)}=xN`$ and $`N^{(2)}=(1x)N`$ particles in two phases. The simple form of the result shows that the calculation leading to (25) is essentially a generalization of the Stirling approximation.

Inserting (25) and (27) into (20), we finally have the desired result for $`s(m)`$,

$$s(m)=\psi (\theta ,\lambda )\theta x\lambda xm+x\mathrm{ln}x+(1x)\mathrm{ln}(1x)$$
(28)
Unfortunately this is mainly formal because neither the integral defining $`\psi `$ nor the Legendre transform are likely to be tractable. However, we show in Sec. II C that (28) is equivalent to the more conventional form of the entropy of mixing as given by the second term in in eq. (5). From a conceptual point of view, it should be noted that the conventional form is normally derived by “binning” the distribution of particle sizes $`\sigma `$ and taking the number of bins to infinity after the thermodynamic limit $`N\mathrm{}`$ has been taken (see e.g. ). In our above derivation, on the other hand, we have assumed that even for finite $`N`$ all particles have different sizes $`\sigma _i`$, drawn randomly from $`n^{(0)}(\sigma )`$; the “polydisperse limit” is thus taken simultaneously with the thermodynamic limit. The relation between these two approaches – which lead to the same results – has been discussed in detail by Salacuse . Note that the first limit is physically more plausible for many homopolymer systems (where there may only be thousands or millions of species, with many particles of each) whereas the second limit is more natural for colloidal materials (and also some random copolymers) in which no two particles present are exactly alike, even in a sample of macroscopic size.

Although the full result (28) is intractable, progress can be made for $`x0`$, when the number of particles in one phase is much smaller than in the other. The limit $`x0`$ implies $`\theta \mathrm{}`$ and hence $`\psi (\theta ,\lambda )\mathrm{exp}[\theta +h(\lambda )]`$ where $`h=\mathrm{ln}𝑑\sigma n^{(0)}(\sigma )e^{\lambda w(\sigma )}`$ is a generalized cumulant generating function for $`n^{(0)}(\sigma )`$. From this we derive $`x=\psi /\theta =\psi `$, giving $`\theta =\mathrm{ln}xh`$. The generalized mean size is given by $`xm=\psi /\lambda =\psi h/\lambda `$ and hence $`m=h/\lambda `$. Inserting these results into (25) then shows that

$`s(xm,x)=xx\mathrm{ln}x+x(h\lambda m),m={\displaystyle \frac{h}{\lambda }}`$

and hence from (28)

$$s(m)=x(h\lambda m)+𝒪(x^2).$$
(29)
Noting that $`xN=N^{(1)}`$ is the number of particles in the small phase, the total free energy (18) can then be written as

$$F_{\mathrm{in}}=F_{\mathrm{in}}^{(2)}+F_{\mathrm{in}}^{(1)}N^{(1)}T(h\lambda m)$$
(30)
The term $`(h\lambda m)`$, being multiplied by $`N^{(1)}T`$, can now be interpreted as the entropy of mixing *per particle in the small phase*. It arises from the deviation of the (generalized) mean size $`mm^{(1)}`$ in the small phase from the mean size $`m^{(0)}`$ of the parent, and is given by the Legendre transform of the generalized cumulant generating function:

$$s_{\mathrm{mix}}(m)=h\lambda m,\text{where}m=\frac{h}{\lambda },h=\mathrm{ln}𝑑\sigma n^{(0)}(\sigma )e^{\lambda w(\sigma )}.$$
(31)
Let us now examine the relation between $`s_{\mathrm{mix}}(m)`$, thereby defined, and $`s(m)`$ introduced previously. By construction, $`Ns(m)`$ is the entropy of mixing of the system as a whole; both phases contribute to this. The result (29), which can also be written as

$$Ns(m)=N^{(1)}s_{\mathrm{mix}}(m)$$
(32)
appears to contradict this, seeing as it does not contain a term proportional to $`N^{(2)}`$. The resolution of this paradox comes from the neglected $`𝒪(x^2)`$ terms in (29). In fact, using the single-phase entropy of mixing defined in eq. (31) – and reinstating on the r.h.s. the superscript on $`mm^{(1)}`$ – we can write (32) in the more symmetric form

$`Ns(m)=N^{(1)}s_{\mathrm{mix}}(m^{(1)})+N^{(2)}s_{\mathrm{mix}}(m^{(2)})`$

This is still correct to leading (linear) order in $`x`$ because the added term is $`𝒪(x^2)`$. \[To see this, use the fact that $`m^{(2)}=m^{(0)}xm^{(1)}/(1x)`$ and hence $`m^{(2)}m^{(0)}=𝒪(x)`$. From (31) it can then be deduced that $`s_{\mathrm{mix}}(m^{(2)})(m^{(2)}m^{(0)})^2=𝒪(x^2)`$.\] Similarly, we can rewrite the total free energy (30) as

$$F_{\mathrm{in}}=[F_{\mathrm{in}}^{(1)}N^{(1)}Ts_{\mathrm{mix}}(m^{(2)})]+[F_{\mathrm{in}}^{(2)}N^{(2)}Ts_{\mathrm{mix}}(m^{(2)})]$$
(33)
This shows that in the limit where one of the two phases is much larger than the other one, the free energy of each of the phases – now including the entropy of mixing – is given by

$$F=F_{\mathrm{in}}NTs_{\mathrm{mix}}(m)=NT\left[\mathrm{ln}(N/V)+1+s_{\mathrm{mix}}(m)\right]+\stackrel{~}{F}$$
(34)
This expression, which only depends on the particle size distribution through the generalized mean size $`m`$, is our desired moment free energy.

The Legendre transform result (31) for the (single-phase) entropy of mixing is appealingly simple; we will illustrate its application in Sec. V A, using polydisperse Flory-Huggins theory as an example. As discussed briefly in App. B, eq. (31) also establishes an interesting connection to large deviation theory. However, the most important aspect of the result (31,34) is that – as we have shown – it gives exact results for the limiting case $`N^{(1)}N^{(2)}`$ (that is, $`x0`$). This includes two important classes of problems, that are also handled exactly by the projection method (for essentially the same reason). The first is the determination of spinodal curves and critical points. Intuitively these are exact because they are related to the stability of the system with respect to small variations in its density or composition, which can be probed by allowing fluctuations to take place in a vanishingly small subregion. The second concerns the cloud point and shadow curves. Again intuitively, these are exact because by definition only an infinitesimal amount of a second phase has appeared. Formal proofs of these statements are given in Sec. III.

If the result (31,33,34) is applied in the regime where it is no longer strictly valid, i.e., where the two coexisting phases contain comparable number of particles (or, equivalently, occupy comparable volumes), then approximate two-phase coexistences can be calculated. It is also straightforward to show that the analog of eqs. (31,33,34) holds in the case where two or more small phases coexist with a much larger phase; application in the regime of comparable phase volumes then provides approximate results for multi-phase coexistence.

To end this section, let us state in full the analogue of (31,34) for the case of several moment densities, restoring the notation used in the previous sections. The square bracket on the r.h.s. of (34) is the moment expression for the entropy of an ideal mixture. If, as in Sec. II A, we measure this entropy per unit volume (rather than per particle, as previously in the current Section) and generalize to several moment densities, we find by the combinatorial approach the following moment free energy:

$$f_{\mathrm{comb}}=Ts_{\mathrm{comb}}+\stackrel{~}{f},s_{\mathrm{comb}}=\rho _0(\mathrm{ln}\rho _01)+\rho _0\left(h\underset{i0}{}\lambda _im_i\right),$$
(35)
with

$$h=\mathrm{ln}𝑑\sigma n^{(0)}(\sigma )\mathrm{exp}\left(\underset{i0}{}\lambda _iw_i(\sigma )\right)$$
(36)
Here $`\rho _0=N/V`$ is the particle density as before.

### C Relation between the two methods

Although eq. (35) looks somewhat different from the corresponding result (9) for $`f_{\mathrm{pr}}`$ obtained by the projection method, the two methods are, mathematically, almost equivalent, as we now show. (For brevity we return to the case of a single moment $`m`$ and continue to refer to the polydisperse feature $`\sigma `$ as ‘size’.) First, consider the exact expression for the entropy of mixing (28) derived within the combinatorial method. This should be equivalent to the conventional result used as the starting point (5) in the projection method. To see this, write the Legendre transform conditions (26) out explicitly:

$$x=𝑑\sigma n^{(0)}(\sigma )\frac{\mathrm{exp}[\theta +\lambda w(\sigma )]}{1+\mathrm{exp}[\theta +\lambda w(\sigma )]},xm^{(1)}=𝑑\sigma w(\sigma )n^{(0)}(\sigma )\frac{\mathrm{exp}[\theta +\lambda w(\sigma )]}{1+\mathrm{exp}[\theta +\lambda w(\sigma )]}$$
(37)
Here we have reinstated the superscript to show that $`m^{(1)}`$ is the value in phase one of the (generalized) moment $`m`$:

$`m^{(1)}={\displaystyle 𝑑\sigma w(\sigma )n^{(1)}(\sigma )}`$

Comparing with (37), we can identify the (normalized) particle size distribution in that phase as

$`n^{(1)}(\sigma )={\displaystyle \frac{n^{(0)}(\sigma )}{x}}{\displaystyle \frac{\mathrm{exp}[\theta +\lambda w(\sigma )]}{1+\mathrm{exp}[\theta +\lambda w(\sigma )]}}`$

Particle conservation $`xn^{(1)}(\sigma )+(1x)n^{(2)}(\sigma )=n^{(0)}(\sigma )`$ implies a similar form for the size distribution in phase two:

$`n^{(2)}(\sigma )={\displaystyle \frac{n^{(0)}(\sigma )}{1x}}{\displaystyle \frac{1}{1+\mathrm{exp}[\theta +\lambda w(\sigma )]}}.`$

It then takes only a few lines of algebra to show that the entropy of mixing (28) can be written (for any $`x`$) as

$`s(m)`$ $`=`$ $`x{\displaystyle 𝑑\sigma n^{(1)}(\sigma )\mathrm{ln}\frac{n^{(1)}(\sigma )}{n^{(0)}(\sigma )}}(1x){\displaystyle 𝑑\sigma n^{(2)}(\sigma )\mathrm{ln}\frac{n^{(2)}(\sigma )}{n^{(0)}(\sigma )}}`$ (38)
$`=`$ $`x{\displaystyle 𝑑\sigma n^{(1)}(\sigma )\mathrm{ln}n^{(1)}(\sigma )}(1x){\displaystyle 𝑑\sigma n^{(2)}(\sigma )\mathrm{ln}n^{(2)}(\sigma )}+{\displaystyle 𝑑\sigma n^{(0)}(\sigma )\mathrm{ln}n^{(0)}(\sigma )}`$ (39)
The second equation, with the third (constant) term discarded, is the standard result (compare eqs. (5,6,A1)). However the first expression appears more natural, and is retained, within the combinatorial derivation; it embodies the intuitively reasonable prescription that entropy is best measured relative to the parent distribution $`n^{(0)}(\sigma )`$. This was also an essential ingredient of the projection method, as described in Sec. II A, where the “prior” in the entropy expression was likewise identified with the parent. Retaining the parent as prior also avoids subtleties with the definition of the integrals in (39) in the case where the phases contain monodisperse components, corresponding to $`\delta `$-peaks in the density distributions .

Having established that the rigorous starting points of the combinatorial and projection method are closely related, we now show that the subsequent approximations also lead to essentially the same results. The relevant approximations are use of (35) for $`f_{\mathrm{comb}}`$ (which adopts the small $`x`$ form for $`s_{\mathrm{mix}}`$) in the combinatorial case and use of (9) for $`f_{\mathrm{pr}}`$ (which minimizes over transverse degrees of freedom) in the projection case. First note that the density $`\rho _0`$ of phases within the family (10) is given by

$`\rho _0={\displaystyle 𝑑\sigma \rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right)}`$

Comparing this with (36), and using the fact that $`n^{(0)}(\sigma )=\rho ^{(0)}(\sigma )/\rho _0^{(0)}`$, one sees that

$`\rho _0={\displaystyle 𝑑\sigma \rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right)}=\rho _0^{(0)}e^{\lambda _0+h}`$

where the sum over $`i`$ now includes $`i=0`$. Solving for $`h`$, one has $`h=\lambda _0+\mathrm{ln}(\rho _0/\rho _0^{(0)})`$. The ideal mixture entropy derived by the combinatorial route, eq. (35), can thus be rewritten as

$`s_{\mathrm{comb}}=\rho _0(\mathrm{ln}\rho _01)+\rho _0\mathrm{ln}(\rho _0/\rho _0^{(0)})\lambda _0\rho _0{\displaystyle \underset{i0}{}}\lambda _i\rho _0m_i=\rho _0{\displaystyle \underset{i}{}}\lambda _i\rho _i\rho _0\mathrm{ln}\rho _0^{(0)}.`$

This is identical to the projected entropy $`s_{\mathrm{pr}}`$, eq. (9), except for the last term. But by construction, the combinatorial entropy assumes that $`\rho _0`$ – the overall density – is among the moment densities retained in the moment free energy. The difference $`s_{\mathrm{comb}}s_{\mathrm{pr}}=\rho _0\mathrm{ln}\rho _0^{(0)}`$ is then linear in this density, and the combinatorial and projection methods therefore predict exactly the same phase behavior.

In summary, we have shown that the projection and combinatorial methods for obtaining moment free energies give equivalent results. The only difference between the two approaches is that within the projection approach, one need not necessarily retain the zeroth moment, which is the overall density $`\rho _0=\rho `$, as one of the moment densities on which the moment free energy depends. If $`\rho _0`$ does not appear in the excess free energy, this reduces the minimum number of independent variables of the moment free energy by one (see Sec. V for an example).

## III Properties of the moment free energy

In the previous sections, we have derived by two different routes (namely (9) and (35)) our moment free energy for truncatable polydisperse systems. We now investigate the properties of this moment free energy, which we henceforward denote $`f_\mathrm{m}`$ and write in the form (9)

$$f_\mathrm{m}(\rho _i)=Ts_\mathrm{m}(\rho _i)+\stackrel{~}{f}(\rho _i),s_\mathrm{m}=\rho _0\underset{i}{}\lambda _i\rho _i;$$
(40)
here $`s_\mathrm{m}`$ is the “moment entropy” of an ideal mixture. In particular, we formally compare the phase behavior predicted from this free energy (by treating the $`\rho _i`$ as densities of “quasi-species” of particles, and applying the usual tangency construction) with that obtained from the exact free energy of the underlying (truncatable) model

$$f[\rho (\sigma )]=T𝑑\sigma \rho (\sigma )\left[\mathrm{ln}\rho (\sigma )1\right]+\stackrel{~}{f}(\rho _i).$$
(41)
Let us first collect a few simple properties of the moment free energy (40) which will be useful later. Recall that (40) faithfully represents the free energy density of any phase with density distribution in the family

$$\rho (\sigma )=\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _iw_i(\sigma )\right)$$
(42)
where $`\rho ^{(0)}(\sigma )`$ is the density distribution of the parent. The moment densities $`\rho _i`$ are then related to the Lagrange multipliers $`\lambda _i`$ by

$$\rho _i=𝑑\sigma w_i(\sigma )\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{j}{}\lambda _jw_j(\sigma )\right)$$
(43)
If we regard $`\rho _0`$ as a function of the $`\lambda _i`$, then from (43) we have $`\rho _0/\lambda _i=\rho _i`$. Together with eq. (40), this implies that the moment entropy $`s_\mathrm{m}`$ has the structure of a Legendre transform . For the first derivatives of $`s_\mathrm{m}`$ with respect to the moment densities, this yields

$`{\displaystyle \frac{s_\mathrm{m}}{\rho _i}}=\lambda _i`$

while the matrix of second derivatives is the negative inverse of the matrix of “second-order moment densities” $`\rho _{ij}`$ :

$$\frac{^2s_\mathrm{m}}{\rho _i\rho _j}=(𝑴^1)_{ij},(𝑴)_{ij}=\frac{^2\rho _0}{\lambda _i\lambda _j}=𝑑\sigma w_i(\sigma )w_j(\sigma )\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{k}{}\lambda _kw_k(\sigma )\right)\rho _{ij}.$$
(44)
The chemical potentials $`\mu _i`$ conjugate to the moment densities follow as

$$\mu _i=\frac{f_\mathrm{m}}{\rho _i}=T\lambda _i+\frac{\stackrel{~}{f}}{\rho _i}=T\lambda _i+\stackrel{~}{\mu }_i$$
(45)
and their derivatives w.r.t. the $`\rho _i`$, which give the curvature of the moment free energy, are

$$\frac{\mu _i}{\rho _j}=\frac{^2f_\mathrm{m}}{\rho _i\rho _j}=T(𝑴^1)_{ij}+\frac{^2\stackrel{~}{f}}{\rho _i\rho _j}.$$
(46)
The pressure, finally, is given by

$$\mathrm{\Pi }_\mathrm{m}=f_\mathrm{m}\underset{i}{}\mu _i\rho _i=T\rho _0+\stackrel{~}{f}\underset{i}{}\stackrel{~}{\mu }_i\rho _i.$$
(47)
Eqs. (45-47) are all calculated via the moment free energy. The three corresponding quantities obtained from the exact free energy (41) are, first, the chemical potentials conjugate to $`\rho (\sigma )`$:

$$\mu (\sigma )=\frac{\delta f}{\delta \rho (\sigma )}=T\mathrm{ln}\rho (\sigma )+\underset{i}{}\stackrel{~}{\mu }_iw_i(\sigma )$$
(48)
second, their derivatives w.r.t. $`\rho (\sigma )`$:

$$\frac{\delta \mu (\sigma )}{\delta \rho (\sigma ^{})}=\frac{\delta ^2f}{\delta \rho (\sigma )\delta \rho (\sigma ^{})}=\frac{T\delta (\sigma \sigma ^{})}{\rho (\sigma )}+\underset{i,j}{}\frac{^2\stackrel{~}{f}}{\rho _i\rho _j}w_i(\sigma )w_j(\sigma ^{})$$
(49)
and, third, the resulting expression for the pressure:

$$\mathrm{\Pi }=f𝑑\sigma \mu (\sigma )\rho (\sigma )=T\rho _0+\stackrel{~}{f}\underset{i}{}\stackrel{~}{\mu }_i\rho _i.$$
(50)
The last of these is identical to the result (47) derived from the moment free energy.

We note one important consequence of the form of the exact chemical potentials (48) for truncatable models: If two phases $`\rho ^{(1)}(\sigma )`$ and $`\rho ^{(2)}(\sigma )`$ have the same chemical potentials $`\mu (\sigma )`$, then the ratio of their density distributions can be written as a Gibbs-Boltzmann factor:

$$\frac{\rho ^{(1)}(\sigma )}{\rho ^{(2)}(\sigma )}=\mathrm{exp}\left[\underset{i}{}\beta \left(\stackrel{~}{\mu }_i^{(2)}\stackrel{~}{\mu }_i^{(1)}\right)w_i(\sigma )\right]$$
(51)
This implies that if one of the density distributions is in the family (42), then so is the other. The same argument obviously applies if there are several phases with equal chemical potentials. Conversely, we have for the chemical potential difference between any two phases in the family (42)

$$\mathrm{\Delta }\mu (\sigma )=\underset{i}{}(T\mathrm{\Delta }\lambda _i+\mathrm{\Delta }\stackrel{~}{\mu }_i)w_i(\sigma )=\underset{i}{}\mathrm{\Delta }\mu _iw_i(\sigma ).$$
(52)
The chemical potentials $`\mu (\sigma )`$ of two such phases are therefore equal if and only if their moment chemical potentials $`\mu _i`$ are equal. Combining this with the fact that the pressure (47) derived from the moment free energy is exact (compare (50)), we conclude: Any set of (two or more) coexisting phases calculated from the moment free energy obeys the exact phase equilibrium conditions. That is, if they were brought into contact, they would genuinely coexist. (However, they will not necessarily obey the lever rule with respect to the parent $`\rho ^{(0)}(\sigma )`$ of the given family, unless all but one of them – which then coincides with the parent – are of vanishing volume; see (11).)

### A General criteria for spinodals and (multi-) critical points

We now demonstrate that, for any truncatable model, the moment free energy gives exact spinodals and (multi-) critical points. By this we mean that, using the moment free energy, the values of external control parameters – such as temperature – at which the parent phase $`\rho ^{(0)}(\sigma )`$ becomes unstable or critical can be exactly determined. Our argument treats spinodals and (multi-) critical points in a unified fashion, using the fact that all of them occur when the difference between phases with equal chemical potentials becomes infinitesimal. Because the truncatable models that we are considering are generally derived from mean-field approaches, we do not concern ourselves with critical point singularities, assuming instead that the free energy is a smooth function of all order parameters (the densities in the system). Likewise, we do not discuss the subtle question of how, beyond mean-field theory, free energies can actually be defined in spinodal and unstable regions .

Let us first recap the general criteria for spinodals and critical points in multi-component systems. In order to treat the critera derived from the exact and moment free energies simultaneously, we use the common notation $`𝝆`$ for the vector of densities specifying the system: for the exact free energy, the components of $`𝝆`$ are the values $`\rho (\sigma )`$; for the moment free energy, they are the reduced set $`\rho _i`$. We write the corresponding vector of chemical potentials as

$`𝝁(𝝆)=f(𝝆).`$

This notation emphasizes that the chemical potentials are functions of the densities. The criterion for a spinodal at the parent phase $`𝝆^{(0)}`$ is then that there is an incipient instability direction $`\delta 𝝆`$ along which the chemical potentials do not change:

$$(\delta 𝝆)𝝁(𝝆^{(0)})=(\delta 𝝆)f(𝝆^{(0)})=\mathrm{𝟎}.$$
(53)
As the second form of the criterion shows, an equivalent statement is that the curvature of the free energy along the direction $`\delta 𝝆`$ vanishes. Eq. (53) can also be written as

$$𝝁(𝝆^{(0)}+ϵ\delta 𝝆)𝝁(𝝆^{(0)})=𝒪(ϵ^2)$$
(54)
which is closely related to the critical point criterion that we describe next.

Near a critical point, the parent $`𝝆^{(0)}`$ coexists with another phase that is only slightly different; if, as we assume here, the free energy function is smooth, these two phases are separated – in $`𝝆`$-space – by a “hypothetical phase” which has the same chemical potentials but is (locally) thermodynamically unstable. (This is geometrically obvious even in high dimensions; between any two minima of $`f(𝝆)𝝁𝝆`$, at given $`𝝁`$, there must lie a maximum or a saddle point, which is the required unstable “phase”.) Now imagine connecting these three phases by a smooth curve in density space $`𝝆(ϵ)`$. At the critical point, all three phases collapse, and the variation of the chemical potential around $`𝝆(ϵ=0)=𝝆^{(0)}`$ must therefore obey

$`𝝁(𝝆(ϵ))𝝁(𝝆^{(0)})=𝒪(ϵ^3)`$

Similarly, if $`n`$ phases coexist, we can connect them and the $`n1`$ unstable phases in between them by a curve $`𝝆(ϵ)`$ . If we define an $`n`$-critical point as one where all these phases become simultaneously critical, we obtain the criterion

$$𝝁(𝝆(ϵ))𝝁(𝝆^{(0)})=𝒪(ϵ^{2n1}).$$
(55)
This formulation was proposed by Brannock ; the cases $`n=2`$ and $`n=3`$ correspond to ordinary critical and tricritical points, respectively. The spinodal criterion (54) is of the same form as (55) if one chooses the curve $`𝝆(ϵ)=𝝆^{(0)}+ϵ\delta 𝝆`$.

In summary, we have that the phase $`𝝆^{(0)}`$ is a spinodal or $`n`$-critical point if there is a curve $`𝝆(ϵ)`$ with $`𝝆(ϵ=0)=𝝆^{(0)}`$ such that

$$\mathrm{\Delta }𝝁𝝁(𝝆(ϵ))𝝁(𝝆^{(0)})=𝒪(ϵ^l).$$
(56)
where $`l=2`$ for a spinodal, and $`l=2n1`$ for an $`n`$-critical point. Brannock  has shown that these criteria are equivalent to the determinant criteria introduced by Gibbs . The above forms are more useful for us because they avoid having to define infinite-dimensional determinants (as otherwise required to handle the exact free energy $`f(𝝆)`$ of a polydisperse system, whose argument $`𝝆`$ denotes an infinite number of components, even in the truncatable case). They also show the analogy with the standard criteria for single-species systems (where $`𝝆`$ has only a single component) more clearly.

We can now apply (56) to the exact free energy (41) and show that the resulting criteria are identical to those obtained from the moment free energy. First, note that the curve $`𝝆(ϵ)\rho (\sigma ;ϵ)`$ can always be chosen to lie within the family (42). This follows from the derivation of the criterion (56): The curve $`𝝆(ϵ)`$ is defined as passing through $`l`$ phases with equal chemical potentials, one of them being the parent $`𝝆^{(0)}`$. As shown in (51), all these phases are therefore within the family (42). But then (52) implies that the condition (56) that $`\mathrm{\Delta }\mu (\sigma )`$ must be zero to $`𝒪(ϵ^l)`$ is equivalent to the same requirement for the differences $`\mathrm{\Delta }\mu _i`$ in the moment chemical potential. This proves that the conditions for spinodals and (multi-) critical points derived from the exact and moment free energies are equivalent. Intuitively, one can understand this as follows: Having shown that the spinodal/critical point conditions can be formulated solely in terms of density distributions $`\rho (\sigma )`$ within the family (42), it is sufficient to know the free energy of those density distributions. This is exactly the moment free energy.

Finally, we note that the exactness of these stability and critical point conditions holds not just for the parent $`\rho ^{(0)}(\sigma )`$, but for all phases within the family (42). This is clear, because any such phase could itself be chosen as the parent without changing the family; the moment free energy would change only by irrelevant terms linear in the moment densities. In general, one will not necessarily be interested in the properties of such “substitute” parents. As an exception, substitute parents that differ from the parent only by a change of the overall density are of interest: they lie on the dilution line of distributions from which cloud point and shadow curves are calculated. The dilution line is included in the family (42) if the overall density $`\rho _0`$ (with corresponding weight function $`w_0(\sigma )=1`$) is retained in the moment free energy; in the combinatorial derivation, this is automatically the case.

### B Spinodals

For completeness, we now give the explicit form  of the spinodal criterion (53) for truncatable systems; see also  for an equivalent derivation using the combinatorial approach. Using (46) and abbreviating the matrix of second derivatives of the excess free energy as $`\stackrel{\mathbf{~}}{𝑭}`$, the spinodal condition (53) becomes

$$\left|\stackrel{\mathbf{~}}{𝑭}+T𝑴^1\right|=0,(\stackrel{\mathbf{~}}{𝑭}+T𝑴^1)𝜹𝝆=\mathrm{𝟎}.$$
(57)
As before, the (nonzero) vector $`𝜹𝝆`$ with components $`\delta \rho _i`$ gives the direction of the spinodal instability; the moment densities $`\rho _i`$ and $`\rho _{ij}`$ that appear in $`\stackrel{\mathbf{~}}{𝑭}`$ and $`𝑴`$ are to be evaluated for the parent distribution $`\rho ^{(0)}(\sigma )`$ being studied. Note that (57) is valid for any $`\rho ^{(0)}(\sigma )`$; no specific assumptions about the parent were made in the derivation. We can thus simply drop the “(0)” superscript: for any phase with density distribution $`\rho (\sigma )`$, the point where (57) first becomes zero, as external control parameters are varied, locates a spinodal instability.

More convenient forms of (57) that avoid matrix inversions are obtained after multiplication by the second order moment matrix $`𝑴`$ (which is positive definite for linearly independent weight functions $`w_i(\sigma )`$ and therefore has nonzero determinant):

$$Y=\left|\mathrm{𝟏}+\beta 𝑴\stackrel{\mathbf{~}}{𝑭}\right|=0,(\mathrm{𝟏}+\beta 𝑴\stackrel{\mathbf{~}}{𝑭})𝜹𝝆=\mathrm{𝟎}.$$
(58)
Here $`\beta =1/T`$ in the standard notation. From our general statements in Sec. III A, the spinodal criterion derived from the exact free energy (41) must be identical to this; this is shown explicitly in App. C. Note that the spinodal condition depends only on the (first-order) moment densities $`\rho _i`$ and the second-order moment densities $`\rho _{ij}`$ of the distribution $`\rho (\sigma )`$ (given by (43) and (44)); it is independent of any other of its properties. This simplification, which has been pointed out by a number of authors , is particularly useful for the case of power-law moments (defined by weight functions $`w_i(\sigma )=\sigma ^i`$): If the excess free energy only depends on the moments of order 0, 1$`\mathrm{}`$ $`K1`$ of the density distribution, the spinodal condition involves only $`2K1`$ moments \[up to order $`2(K1)`$\].

The general discussion in Sec. III A as well as the explicit calculation in App. C show that the spinodal instability direction lies within (or more precisely, is tangential to) the family (42) of density distributions. This fact has a simple geometrical interpretation: By construction, the family lies along a “valley” of the free energy surface (compare Fig. 1). Away from the valley floor, any change in $`\rho (\sigma )`$ increases the free energy, corresponding to a positive curvature; the spinodal direction, for which the curvature vanishes, must therefore be along the valley floor.

### C Critical points

Next, we show the explicit form of the critical point criterion for truncatable systems. The general condition (56) for an ordinary critical point ($`n=2`$) was shown by Brannock  to be equivalent to

$$(\delta 𝝆)f(𝝆^{(0)})=\mathrm{𝟎},(\delta 𝝆)^3f(𝝆^{(0)})=0$$
(59)
The first part of this is simply the spinodal criterion, as expected. To evaluate the second part for the moment free energy (40), we need the third derivative of $`s_\mathrm{m}`$ w.r.t. the moment densities $`\rho _i`$. Writing (44) as

$`{\displaystyle \underset{j}{}}\rho _{ij}{\displaystyle \frac{^2s_\mathrm{m}}{\rho _j\rho _k}}=\delta _{ik}`$

and differentiating w.r.t. one of the moment densities, one finds after a little algebra

$`{\displaystyle \frac{^3s_\mathrm{m}}{\rho _i\rho _j\rho _k}}={\displaystyle \underset{lmn}{}}(𝑴^1)_{il}(𝑴^1)_{jm}(𝑴^1)_{kn}\rho _{lmn}`$

with third-order moment densities $`\rho _{lmn}`$ defined in the obvious way. Thus, evaluating (59) for the moment free energy (40), we find that critical points have to obey

$$T\underset{ijk}{}v_iv_jv_k\rho _{ijk}+\delta \rho _i\delta \rho _j\delta \rho _k\frac{^3\stackrel{~}{f}}{\rho _i\rho _j\rho _k}=0,v_i=\underset{j}{}(𝑴^1)_{ij}\delta \rho _j$$
(60)
in addition to the spinodal condition (58). As before, the criterion has to be evaluated for the parent phase $`\rho ^{(0)}(\sigma )`$ under consideration; but because $`\rho ^{(0)}(\sigma )`$ can be chosen arbitrarily, it applies to all density distributions $`\rho (\sigma )`$.

As expected from the general discussion in Sec. III A, the criterion (60) can also be derived from the exact free energy; an alternative form involving the spinodal determinant $`Y`$ is given in App. D. Eq. (60) shows that the location of critical points depend only on the moment densities $`\rho _i`$, $`\rho _{ij}`$ and $`\rho _{ijk}`$ . For a system with an excess free energy depending only on power-law moments up to order $`K1`$, the critical point condition thus involves power-law moments of the parent only up to order $`3(K1)`$.

### D Onset of phase coexistence: Cloud point and shadow

So far in this section, we have shown that the moment free energy gives exact results for spinodals and (multi-) critical points. Now we consider the onset of phase coexistence, where (on varying the temperature, for example) a parent phase with density distribution $`\rho ^{(0)}(\sigma )`$ first starts to phase separate. As explained in Sec. I, this temperature together with the overall density $`\rho ^{(0)}\rho _0^{(0)}`$ of the parent defines a “cloud point”; the density of the incipient daughter phase gives the “shadow”. If the parent begins to coexist with $`p`$ phases simultaneously ($`p=2`$ at a triple point, for example), there will be $`p`$ such shadows.

At the onset of phase coexistence, one of the coexisting phases is by definition the parent $`\rho ^{(0)}(\sigma )`$; the lever rule does not yet play any role because the daughters $`\rho ^{(\alpha )}(\sigma )`$ ($`\alpha =1\mathrm{}p`$) occupy an infinitesimal fraction of the total volume. It then follows from (51) that all daughters lie within the family (42) of density distributions. As shown in (52), the condition for equality of chemical potentials $`\mu (\sigma )`$ between any two of the coexisting phases (parent and daughters) then becomes

$`\mathrm{\Delta }\mu (\sigma )={\displaystyle \underset{i}{}}\mathrm{\Delta }\mu _iw_i(\sigma )`$

and is satisfied if and only if the moment chemical potentials $`\mu _i`$ are equal in all phases. Likewise, the exact and moment free energies give the same condition for equality of pressure $`\mathrm{\Pi }`$ in all phases, because they yield identical expressions (47,50) for $`\mathrm{\Pi }`$. In summary, we see that the conditions for the onset of phase coexistence are identical for the exact and moment free energies; the moment free energy therefore gives exact cloud points and shadows.

### E Phase coexistence beyond onset

As stated in Sec. II A, the moment free energy does not give exact results beyond the onset of phase coexistence, i.e., in the regime where the coexisting phases occupy comparable fractions of the total system volume. As shown in Sec. III A, the calculated phases will still be in exact thermal equilibrium; but the lever rule will now be violated for the “transverse” degrees of freedom of the density distributions. This is clear from (11): In general, no linear combination of distributions from this family can match the parent $`\rho ^{(0)}(\sigma )`$ exactly.

A more detailed understanding of the failure of the moment free energy beyond phase coexistence can be gained by comparing with the formal solution of the exact phase coexistence problem. Assume that the parent $`\rho ^{(0)}(\sigma )`$ has separated into $`p`$ phases numbered by $`\alpha =1\mathrm{}p`$. The condition (51), which follows from equality of the chemical potentials $`\mu (\sigma )`$ in all phases, implies that we can write their density distributions $`\rho ^{(\alpha )}(\sigma )`$ as

$$\rho ^{(\alpha )}(\sigma )=\stackrel{~}{R}(\sigma )\mathrm{exp}\left(\underset{i}{}\lambda _i^{(\alpha )}w_i(\sigma )\right)$$
(61)
for some function $`\stackrel{~}{R}(\sigma )`$. If phase $`\alpha `$ occupies a fraction $`v^{(\alpha )}`$ of the system volume, particle conservation $`_\alpha v^{(\alpha )}\rho ^{(\alpha )}(\sigma )=\rho ^{(0)}(\sigma )`$ then gives

$$\stackrel{~}{R}(\sigma )=\frac{\rho ^{(0)}(\sigma )}{_\alpha v^{(\alpha )}\mathrm{exp}\left(_i\lambda _i^{(\alpha )}w_i(\sigma )\right)},\rho ^{(\alpha )}(\sigma )=\rho ^{(0)}(\sigma )\frac{\mathrm{exp}\left(_i\lambda _i^{(\alpha )}w_i(\sigma )\right)}{_\beta v^{(\beta )}\mathrm{exp}\left(_i\lambda _i^{(\beta )}w_i(\sigma )\right)}$$
(62)
If there are $`K`$ moment densities $`\rho _i`$, $`i=1\mathrm{}K`$, then this exact solution is parameterized by $`(p1)(K+1)`$ independent parameters: $`(p1)K`$ parameters $`\lambda _i^{(\alpha )}`$ (noting that the $`\lambda _i`$ of one phase can be fixed arbitrarily), and $`p1`$ parameters $`v^{(\alpha )}`$ (noting that one phase volume is fixed by the constraint $`_\alpha v^{(\alpha )}=1`$). Comparing (61) with (51), one sees that the $`K`$ quantities

$`T\lambda _i^{(\alpha )}+\stackrel{~}{\mu }_i^{(\alpha )}(i=1\mathrm{}K)`$

must be the same in all phases $`\alpha `$; the same is true for the pressures $`\mathrm{\Pi }^{(\alpha )}`$, and this gives the required total number of $`(p1)(K+1)`$ constraints. This formally defines the exact solution of the phase coexistence problem for truncatable systems. Its practical value is limited by the difficulties of finding the solution numerically, as pointed out in Sec. I; see also Sec. IV. There is also no interpretation of the result in terms of a free energy depending on a small number of densities.

Nevertheless, eq. (62) is useful for a comparison with the solution provided by the moment free energy method. From (61), the exact coexisting phases are all members of a family of density distributions of the general form (7); this is simply a consequence of the requirement of equal chemical potentials (51). This “exact coexistence family” has an “effective” prior  $`\stackrel{~}{R}(\sigma )\rho ^{(0)}(\sigma )`$ and is therefore different from the “original” family (42). The moment free energy only gives us access to distributions from the original family, not the exact coexistence family, and therefore cannot yield exact solutions for phase coexistence (beyond its onset). Fig. 2 illustrates this point explicitly in the simplified context of a bidisperse system.

It is now easy to see, however, that – as stated in Sec. II A – the exact solution can be approached to arbitrary precision by including extra moment densities in the moment free energy. (This leaves the exactness of spinodals, critical points, cloud-points and shadows unaffected, because none of our arguments excluded a null dependence of $`\stackrel{~}{f}`$ on certain of the $`\rho _i`$.) Indeed, by adding further moment densities one can indefinitely extend the family (42) of density distributions, thereby approaching with increasing precision the actual distributions in all phases present; this yields phase diagrams of ever-refined accuracy.

The roles of the “original” moment densities (those appearing in the excess free energy) and the extra ones are quite different, however. To see this, note from (45) that equality among phases of the moment chemical potentials implies that, for the extra moments only, the corresponding Lagrange multipliers must themselves be equal in all coexisting phases. (This is because there is, by construction, no excess part to the “extra” chemical potentials.) We can therefore drop the phase index on these Lagrange multipliers and write the density distribution in phase $`\alpha `$ as

$$\rho ^{(\alpha )}(\sigma )=\rho ^{(0)}(\sigma )\mathrm{exp}\left(\underset{\mathrm{extra}i}{}\lambda _iw_i(\sigma )\right)\mathrm{exp}\left(\underset{\mathrm{original}i}{}\lambda _i^{(\alpha )}w_i(\sigma )\right)$$
(63)
Comparing with (7) and (61), we see that the extra Lagrange multipliers can be thought of as providing a “flexible prior” $`R(\sigma )`$, allowing a better approximation to the effective prior $`\stackrel{~}{R}(\sigma )`$ required by particle conservation. The “tuning” of the prior with these extra Lagrange multipliers also has an important effect on the number of coexisting phases that can be found by the moment method (see Sec. III G).

### F Global and local stability

Most numerical algorithms for phase coexistence calculations, including ours, initially proceed by finding a solution to the equilibrium conditions of equal chemical potentials and pressures in all phases, rather than by a direct minimization of the total free energy of the system. It is then crucial to verify whether this solution is stable, both locally (i.e., with respect to small fluctuations in the compositions of the phases) and globally with respect to splitting into a larger number (or different) phases (see Fig. 3). Global stability is a particularly important issue in our context because there is in principle no limit on the number of coexisting phases in a polydisperse system.

A useful tool for stability calculations is the “tangent plane distance” . Let us first define this generically for a system with (a vector of) densities $`𝝆`$, free energy density $`f(𝝆)`$ and chemical potentials $`𝝁(𝝆)=f(𝝆)`$; this notation is the same as in Sec. III A. Assume we have found a candidate “phase split”, that is a collection of $`p`$ phases $`𝝆^{(\alpha )}`$ ($`\alpha =1\mathrm{}p`$) which satisfy the phase equilibrium conditions of equal chemical potentials $`𝝁=𝝁^{(\alpha )}`$ and pressures $`\mathrm{\Pi }=\mathrm{\Pi }^{(\alpha )}`$. They occupy fractions $`v^{(\alpha )}`$ of the total volume, and the overall density distribution is thus $`𝝆_{\mathrm{tot}}=_\alpha v^{(\alpha )}𝝆^{(\alpha )}`$. Note that we do not yet assume that the lever rule is satisfied, i.e., that $`𝝆_{\mathrm{tot}}=𝝆^{(0)}`$; as explained above, this equality will generally not hold for phase splits calculated from the moment free energy. We now define “global tangent plane (TP) stability” for such a phase split as the property that there is no other phase split (i.e., no other tangent plane) that gives a lower total free energy for the same overall density distribution $`𝝆_{\mathrm{tot}}`$. In more intuitive language, this means that if we were to put the phases $`𝝆^{(\alpha )}`$ into contact with each other, the resulting system would be thermodynamically stable; neither the composition nor the number of phases would change over time. Note, however, that since in general $`𝝆_{\mathrm{tot}}=_\alpha v^{(\alpha )}𝝆^{(\alpha )}`$ is not equal to $`𝝆^{(0)}`$, the phase split that is globally TP-stable need not accurately reflect the number and composition of phases into which the parent $`𝝆^{(0)}`$ would actually split under the chosen thermodynamic conditions.

To define the tangent plane distance (TPD), note first that by virtue of the coexistence conditions, all phases $`𝝆^{(\alpha )}`$ lie on a tangent plane to the free energy surface. Points $`(𝝆,f)`$ on this tangent plane obey the equation $`f𝝁𝝆+\mathrm{\Pi }=0`$, with $`𝝁`$ and $`\mathrm{\Pi }`$ the chemical potentials and pressure common to all phases. For a generic phase with density distribution $`𝝆`$ and free energy $`f(𝝆)`$, the same expression will have a nonzero value which measures how much “below” or “above” the tangent plane it lies. This defines the TPD

$$t(𝝆)=f(𝝆)𝝁^{(\alpha )}𝝆+\mathrm{\Pi }^{(\alpha )}.$$
(64)
Here we have added the superscript $`\alpha `$ to emphasize that the chemical potential and osmotic pressure used in the calculation of the TPD are those of the calculated phase equilibrium (and hence of any of the participating phases $`\alpha `$), rather than those of the test phase $`𝝆`$. It is then clear intuitively – and can be shown more formally  – that the calculated phase coexistence is globally TP-stable if the TPD is non-negative everywhere. Geometrically, this simply means that no part of the free energy surface must protrude beneath the tangent plane; otherwise the total free energy of the system could be lowered by constructing a new tangent plane that touches the protruding piece. Global TP-stability of course encompasses local stability; the latter simply corresponds to the requirement that the TPD be a local minimum (with the value $`t=0`$) at each of the phases in the candidate solution.

When verifying global TP-stability, it is obviously sufficient to check the value of the TPD at all of its stationary points. By differentiating (64) w.r.t. $`𝝆`$, one sees that at these points, the chemical potentials $`𝝁(𝝆)`$ are the same as in the calculated coexisting phases. For our truncatable polydisperse systems, it then follows from eq. (51) that we only need to consider the TPD for test phases $`\rho (\sigma )`$ which are in the same family (7) as the calculated coexisting phases. This is a crucial point: Even though global TP-stability is a statement about stability in the infinite-dimensional space of density distributions $`\rho (\sigma )`$, it can be checked by only considering density distributions from a $`K`$-dimensional family. In the (generally hypothetical) case of an exactly calculated phase split, this family would be (62), with the effective prior $`\stackrel{~}{R}(\sigma )`$ . For a phase split calculated from the moment free energy, it is the family (42) with the parent as prior, with any Lagrange multipliers for extra moments fixed to their values in the calculated coexisting phases. The TPD of such a test phase is (using (41) for the exact free energy and (48) for the chemical potentials $`\mu ^{(\alpha )}(\sigma )`$ of the coexisting phases)

$`t[\rho (\sigma )]`$ $`=`$ $`f[\rho (\sigma )]{\displaystyle 𝑑\sigma \mu ^{(\alpha )}(\sigma )\rho (\sigma )}+\mathrm{\Pi }^{(\alpha )}`$ (65)
$`=`$ $`\stackrel{~}{f}(\rho _i)+{\displaystyle 𝑑\sigma \rho (\sigma )\left\{T\left[\mathrm{ln}\rho (\sigma )1\right]T\mathrm{ln}\rho ^{(\alpha )}(\sigma )\underset{i}{}\stackrel{~}{\mu }_i^{(\alpha )}w_i(\sigma )\right\}}+\mathrm{\Pi }^{(\alpha )}`$ (66)
$`=`$ $`\stackrel{~}{f}(\rho _i)+T{\displaystyle 𝑑\sigma \rho (\sigma )\left[\mathrm{ln}\frac{\rho (\sigma )}{\rho ^{(0)}(\sigma )}1\right]}{\displaystyle \underset{i}{}}\left[\lambda _i^{(\alpha )}+\stackrel{~}{\mu }_i^{(\alpha )}\right]\rho _i+\mathrm{\Pi }^{(\alpha )}`$ (67)
$`=`$ $`f_\mathrm{m}(\rho _i){\displaystyle \underset{i}{}}\mu _i^{(\alpha )}\rho _i+\mathrm{\Pi }^{(\alpha )}.`$ (68)
A comparison with (64) shows that this is identical to the TPD that one would derive from the moment free energy alone. Here again, the underlying polydisperse nature of the problem can therefore be disregarded once the moment free energy has been obtained. In summary, one can determine whether a phase split calculated from the moment free energy is globally TP-stable (which means that phases of the predicted volumes and compositions, would, if placed in contact, indeed coexist) using only the TPD derived from the moment free energy; the same is trivially true of the weaker requirement of local stability .

Recall, however, that the overall density distribution $`\rho _{\mathrm{tot}}(\sigma )`$ for a phase split found from the moment free energy only has the same moment densities $`\rho _i^{(0)}`$ as the parent $`\rho ^{(0)}(\sigma )`$, but differs in other details (the transverse degrees of freedom). Global TP-stability thus guarantees that such an approximate phase split is thermodynamically stable, but does not imply that it is identical (in either number or composition of phases) to the exact one. Nevertheless, progressively increasing the number of extra moment densities in the moment free energy will make $`\rho _{\mathrm{tot}}(\sigma )`$ a progressively better approximation to $`\rho ^{(0)}(\sigma )`$, and so the exact phase split, with the correct number of phases, must eventually be recovered to arbitrary accuracy; see Fig. 4 for an illustration.

### G Geometry in density distribution space

In the above discussion of the properties of the moment free energy, we have focused on obtaining the phase behavior of a system with a given parent distribution $`\rho ^{(0)}(\sigma )`$. This is the point of view most relevant for practical applications of the method, and the rest of this Section is not essential for understanding such applications. Nevertheless, from a theoretical angle, it is also interesting to consider the global geometry of the space of all density distributions $`\rho (\sigma )`$, without reference to a specific parent: the exact free energy (41) of the underlying (truncatable) model induces two-phase tie lines and multi-phase coexistence regions in this space, and one is led to ask how the moment free energy encodes these properties.

As pointed out in Sec. II A, the definition of the moment free energy depends on a prior $`R(\sigma )`$ and represents the properties of systems with density distributions $`\rho (\sigma )`$ in the corresponding maximum entropy family (7). Instead of identifying $`R(\sigma )=\rho ^{(0)}(\sigma )`$, we now allow a general prior $`R(\sigma )`$. Conceptually, it then makes sense to associate the moment free energy with the family (7) rather than the specific prior. This is because any density distribution from (7) can be chosen as prior, without changing the moment free energy (apart from irrelevant linear terms in the moments $`\rho _i`$), or the identity of the remaining family members. The construction of the moment free energy thus partitions the space of all $`\rho (\sigma )`$ into (an infinite number of) families (7); different families give different moment free energies that describe the thermodynamics “within the family”. This procedure gives meaningful results because, as shown at the end of Sec. III A, coexisting phases are always members of the same family. By considering the entire ensemble of families (7) and their corresponding moment free energies, one can thus in principle recover the exact geometry of the density distribution space and the phase coexistences within it. Note that this includes regions with more than $`K+1`$ phases, even when the excess (and thus moment) free energy depends on only $`K`$ moment densities; see Fig. 4. Consistent with Gibbs’ phase rule, however, the families for which this occurs are exceptional: they occupy submanifolds (of measure zero) in the space of all families. Accordingly, to find the families on which such “super-Gibbs” multi-phase coexistences occur, the corresponding prior must be very carefully tuned. Indeed, the probability of finding such priors “accidentally”, without either adding extra moments to the moment free energy description or solving the exact phase coexistence problem, is zero. Generically, one needs to retain at least $`n`$ moment densities in the moment free energy to find $`n+1`$ phases in coexistence. From (63), this corresponds to having $`nK`$ parameters available for tuning the location of the family (or equivalently, its prior) in density distribution space. But note that, once the extra moments have been introduced, the required tuning need not be done “by hand”: it is achieved implicitly by requiring the lever rule to be obeyed not only for the $`K`$ original moment densities, but also for the $`nK`$ extra ones. (The solution of the phase equilibrium conditions thus effectively proceeds in an $`n`$-dimensional space; but once a solution has been found, its (global TP-)stability can still be checked by computing the TPD only within a $`K`$-dimensional family of distributions. See Sec. III F.)

As pointed out above, phase splits calculated from the moment free energy allow violations of the lever rule. In the global view, this fact also has a simple geometric interpretation: the families (7) are (generically) curved. In other words, for two distributions $`\rho ^{(1)}(\sigma )`$ and $`\rho ^{(2)}(\sigma )`$ from the same family, the straight line $`\rho (\sigma )=ϵ\rho ^{(1)}(\sigma )+(1ϵ)\rho ^{(2)}(\sigma )`$ connecting them lies outside the family. More generally, if $`p`$ coexisting phases $`\rho ^{(\alpha )}(\sigma )`$ have been identified from the moment free energy, then the overall density distribution $`\rho _{\mathrm{tot}}(\sigma )=_\alpha v^{(\alpha )}\rho ^{(\alpha )}(\sigma )`$ is different from the parent $`\rho ^{(0)}(\sigma )`$; geometrically, it lies on the hyperplane that passes through the phases $`\rho ^{(\alpha )}(\sigma )`$.

## IV Practical implementation of the moment method

The application of the moment free energy method to the calculation of spinodal and critical points is straighforward using conditions (58), and (60) or (D3), respectively, and is further illustrated in Sec. V below. We therefore focus in this Section on phase coexistence calculations.

Recall that in the moment approach, each phase $`\alpha `$ is parameterized by Lagrange multipliers $`\lambda _i^{(\alpha )}`$ for the original moments (the ones appearing in the excess free energy of the system) and the fraction $`v^{(\alpha )}`$ of system volume that it occupies. If extra moments are used, there is one additional Lagrange multiplier $`\lambda _i`$ for each of them; these are common to all phases. These parameters have to be chosen such that the pressure (47) and the moment chemical potentials $`\mu _i`$ given by (45) are equal in all phases. Furthermore, the (fractional) phase volumes $`v^{(\alpha )}`$ have to sum to one, and the lever rule has to be satisfied for all moments (both original and extra):

$`{\displaystyle \underset{\alpha }{}}v^{(\alpha )}\rho _i^{(\alpha )}=\rho _i^{(0)}.`$

In a system with $`K`$ original moments that is being studied using an $`n`$-moment free energy (i.e., with $`nK`$ extra moments), and for $`p`$ coexisting phases, one has $`p(K+1)+nK=(p1)(K+1)+n+1`$ parameters and as many equations. Starting from a suitable initial guess, these can, in principle, be solved by a standard algorithm such as Newton-Raphson . Generating an initial point from which such an algorithm will converge, however, is a nontrivial problem, especially when more than two phases coexist.

To simplify this task, we work with a continuous control parameter such as temperature, density $`\rho _0`$ of the parent phase, or interaction parameter $`\chi `$ for polymers. Taking the latter case as an example, we start the calculation at a small value of $`\chi `$ where we are sure to be in a single-phase region; thus the parent phase under consideration is stable. The basic strategy is then to increment $`\chi `$ and detect potential new phases as we go along. At each step, the phase equilibrium conditions are solved by Newton-Raphson for the current number of phases. Then we check for local stability of the solution by calculating the Hessian of the TPD around each of the coexisting phases and verifying that it is positive definite. If an instability (i.e., a negative eigenvalue of the Hessian) is found, we search for local minima starting from points displaced either way along the instability direction. If two new local minima are found in this way, we add them to list of phases and delete the unstable phase; if only one new local minimum is uncovered, we add this but retain the old phase. Finally, we check for global (TP-)stability. As explained in Sec. III F, this involves scanning all “test” phases from the same maximum entropy family for possible negative values of the TPD. However, the test phase will have the same extra Lagrange multipliers, and is therefore parameterized in terms of the $`\lambda _i`$ for the original moments only. The extra Lagrange multipliers are held fixed, so their number is irrelevant for the computational cost of the stability check. Put differently, within our algorithm the precision of the coexistence curves should depend on the total number of moments retained, $`n`$, rather than the number of moments in the excess free energy, $`K`$. Computational effort, on the other hand, is dominated by the global stability check, and thus mainly sensitive to $`K`$ (the dimension of the space to be searched for new phases), rather than $`n`$.

This is a substantial efficiency gain, but despite it, an exhaustive search for local TPD minima over the $`K`$-dimensional space of original Lagrange multipliers is unrealistic except for $`K=1`$. We use instead a Monte-Carlo-type algorithm  to sample the TPD at representative points. If any negative values of the TPD are encountered, we choose the smallest such value, find the nearest local minimum of the TPD and add this point to the list of phases. If any new phases have been found during the (local and global) stability checks, we assign each of them a default phase volume ($`v=0.1`$, say), reduce the phase volumes of the old phases accordingly, update the number of phases and loop back to the Newton-Raphson solution of the phase equilibrium conditions, using the current list of phases as an initial guess. This process is repeated until the calculated solution is found to be stable. Our implementation of this basic scheme also contains some additional elements (such as an adaptive choice of the stepsize for the control parameter $`\chi `$, and checks for very small phase volumes) which are useful near points where new phases appear or old ones vanish.

Comparing our own approach, which is based on the moment free energy with extra moment densities as just outlined, to (62), which is the exact solution for truncatable systems, we see that the former actually uses more parameters (essentially one phase-independent Lagrange multiplier per extra moment) to represent the solution. At first sight, this may appear counter-productive. However, the Lagrange multipliers $`\lambda _i^{(\alpha )}`$ and the (fractional) phase volumes $`v^{(\alpha )}`$ are much less strongly coupled in the moment free energy approach. The phase volumes are only determined by the lever rule, while in the exact solution they “feed back” into the effective prior $`\stackrel{~}{R}(\sigma )`$ and therefore into the equilibrium conditions of equal chemical potentials and pressures. In the examples studied below, we have found that the numerical advantages of this decoupling (in terms of stability, robustness and convergence of our algorithm) can easily outweigh the larger number of parameters that it requires.

As with any numerical algorithm, it is important to develop robust criteria by which the convergence of the solution can be judged. This has a novel aspect, in the moment method, since violations of the lever rule are allowed: alongside normal numerical convergence criteria one needs a method for deciding whether the effect of these on the predicted phase behaviour is significant. We develop appropriate criteria in Sec. V B in the context of a specific example. The basic ideas is that since any state of phase coexistence predicted by the method represents the exact behaviour of some parent, then so long as this parent is close enough to the true one, the predicted behaviour will lie within the range of uncertainty that arises anyway, from not knowing the true parent to arbitrary experimental precision.

## V Examples

We now illustrate how the moment method is applied and demonstrate its usefulness for several examples. The first two (Flory-Huggins theory for length-polydisperse homopolymers and dense chemically polydisperse copolymers, respectively) contain only a single moment density in the excess free energy and are therefore particularly simple to analyse and visualize. In the third example (chemically polydisperse copolymers in a polymeric solvent), the excess free energy depends on two moment densities and this will gives us the opportunity to discuss the appearance of more complex phenomena such as tricritical points.

### A Homopolymers with length polydispersity

Let us start with the simple but well studied example of polydisperse Flory-Huggins theory . One considers a system of homopolymers with a distribution of chain lengths; the polydisperse feature $`\sigma `$ is simply the chain length $`L`$, i.e., the number of monomers in each chain. (We treat this as a continuous variable.) The density distribution $`\rho (L)`$ then gives the number density of chains as a function of $`L`$. We choose the segment volume $`a^3`$ as our unit of volume, making $`\rho (L)`$ dimensionless. The volume fraction occupied by the polymer is then simply the first moment of this distribution, $`\varphi \rho _1=𝑑LL\rho (L)`$. Within Flory-Huggins theory, the free energy density is (in units such that $`k_\text{B}T=1`$)

$$f=𝑑L\rho (L)[\mathrm{ln}\rho (L)1]+(1\rho _1)\mathrm{ln}(1\rho _1)+\chi \rho _1(1\rho _1)$$
(69)
where the Flory $`\chi `$-parameter plays essentially the role of an inverse temperature. Before analysing this further, we note that for a bidisperse system with chain lengths $`L_1`$ and $`L_2`$ and number densities $`r_1`$ and $`r_2`$, the corresponding expression would be

$$f=r_1(\mathrm{ln}r_11)+r_2(\mathrm{ln}r_21)+(1\rho _1)\mathrm{ln}(1\rho _1)+\chi \rho _1(1\rho _1),\rho _1=L_1r_1+L_2r_2$$
(70)
This free energy, with $`L_1=10`$ and $`L_2=20`$, was used to generate the examples shown in Figs. 1 and 2.

Returning now to (69), we note first that the excess free energy

$$\stackrel{~}{f}=(1\rho _1)\mathrm{ln}(1\rho _1)+\chi \rho _1(1\rho _1)$$
(71)
depends only on the moment density $`\rho _1`$. If no extra moments are used, the moment free energy is therefore a function of a single density variable, $`\rho _1`$. Let us work out its construction explicitly for the case of a parent phase with a Schulz-distribution of lengths, given by

$$\rho ^{(0)}(L)=\rho _0^{(0)}\frac{1}{\mathrm{\Gamma }(a)b^a}L^{a1}e^{L/b}$$
(72)
Here $`a`$ is a parameter that determines how broad or peaked the distribution is; it is conventionally denoted by $`\alpha `$, but we choose a different notation here to prevent confusion with the phase index used in the general discussion so far. The chain number density of the parent is $`\rho _0^{(0)}`$, and the normalized first moment $`m_1^{(0)}=\rho _1^{(0)}/\rho _0^{(0)}`$ is simply the (number) average chain length

$$L_Nm_1^{(0)}\frac{\rho _1^{(0)}}{\rho _0^{(0)}}=ab.$$
(73)
The parent is thus parameterized in terms of $`a`$, $`\rho _0^{(0)}`$ and $`\rho _1^{(0)}`$ (or $`b`$). The density distributions in the family (42) are given by

$$\rho (L)=\rho ^{(0)}(L)e^{\lambda _1L}=\rho _0^{(0)}\frac{1}{\mathrm{\Gamma }(a)b^a}L^{a1}e^{(b^1\lambda _1)L}$$
(74)
and the zeroth and first moment densities are easily worked out, for members of this family, to be

$`\rho _0=\rho _0^{(0)}\left({\displaystyle \frac{1}{1b\lambda _1}}\right)^a,\rho _1=\rho _1^{(0)}\left({\displaystyle \frac{1}{1b\lambda _1}}\right)^{a+1}.`$

Because the family only has a single parameter $`\lambda _1`$, these two moment densities are of course related to one another:

$$\rho _0=c\rho _1^{a/(a+1)},c=\rho _0^{(0)}\left(\rho _1^{(0)}\right)^{a/(a+1)}.$$
(75)
Now we turn to the moment entropy, given by (40): $`s_\mathrm{m}=\rho _0\lambda _1\rho _1`$, and to be considered as a function of $`\rho _1`$. After a little algebra (and eliminating $`b`$ in favor of $`\rho _1^{(0)}`$) one finds

$$s_\mathrm{m}=\rho _0^{(0)}\left[(1+a)\left(\frac{\rho _1}{\rho _1^{(0)}}\right)^{a/(a+1)}a\frac{\rho _1}{\rho _1^{(0)}}\right]$$
(76)
The last term is linear in $`\rho _1`$ and can be disregarded for the calculation of phase equilibria. This then gives the following simple result for the moment free energy

$$f_\mathrm{m}=(a+1)c\rho _1^{a/(a+1)}+\stackrel{~}{f}.$$
(77)
(In conventional polymer notation, this result would read $`f_\mathrm{m}=(\alpha +1)c\varphi ^{\alpha /(\alpha +1)}+\stackrel{~}{f}`$.) From this we can now obtain the spinodal condition, for example, which identifies the value of $`\chi `$ where the parent becomes unstable. The general criterion (53) simplifies in our case of a single moment density to

$$\frac{d^2f_\mathrm{m}}{d\rho _1^2}|_{\rho _1=\rho _1^{(0)}}=0\frac{1}{1\rho _1^{(0)}}2\chi +\frac{a}{a+1}\frac{\rho _0^{(0)}}{(\rho _1^{(0)})^2}=0$$
(78)
The same condition for a phase with general density distribution $`\rho (L)`$ – rather than our specific Schulz parent $`\rho ^{(0)}(L)`$ – follows from (58) as

$$\frac{1}{1\rho _1}2\chi +\frac{1}{\rho _2}$$
(79)
As expected, this becomes equivalent to (78) for parents of the Schulz form (72), which obey $`\rho _2^{(0)}=\rho _0^{(0)}a(a+1)b^2`$ and eq. (73). Note that in standard polymer notation, eq. (79) would be written as

$$\frac{1}{1\varphi }2\chi +\frac{1}{L_W\varphi }$$
(80)
where $`L_W=\rho _2/\rho _1`$ is the weight average chain length.

The critical point condition is obtained similarly. For a single density variable, eq. (59) reduces to

$$\frac{d^3f_\mathrm{m}}{d\rho _1^3}|_{\rho _1=\rho _1^{(0)}}=0\frac{1}{(1\rho _1^{(0)})^2}\frac{a(a+2)}{(a+1)^2}\frac{\rho _0^{(0)}}{(\rho _1^{(0)})^3}=0$$
(81)
For a general distribution, the critical point criterion (60) becomes instead

$`{\displaystyle \frac{\rho _3}{\rho _2^3}}+{\displaystyle \frac{d^3\stackrel{~}{f}}{d\rho _1^3}}=0`$

and inserting the Flory-Huggins excess free energy (71) gives

$$\frac{1}{(1\rho _1)^2}\frac{\rho _3}{\rho _2^3}=0$$
(82)
Again, this simplifies to (81) for Schulz parents, for which $`\rho _3^{(0)}=\rho _0^{(0)}a(a+1)(a+2)b^3`$. The traditional statement of (82), in terms of $`\varphi \rho _1`$, $`L_W`$ and the $`Z`$-average chain length $`L_Z=\rho _3/\rho _2`$ reads

$$\frac{1}{(1\varphi )^2}\frac{L_Z}{L_W^2}\frac{1}{\varphi ^2}=0.$$
(83)
Fig. 5 shows the moment free energy (77) for a given parent and different values of $`\chi `$. The value of $`\chi `$ for which there is a double tangent at the point $`\rho _1=\rho _1^{(0)}`$ representing the parent gives the cloud point; the second tangency point gives the polymer volume fraction $`\rho _1`$ in the coexisting phase, hence the shadow point. Repeating this procedure for different values of $`\rho _1^{(0)}`$ (while maintaining $`L_N=\rho _1^{(0)}/\rho _0^{(0)}`$ and $`a`$ constant), one obtains the full CPC and shadow curve for a given normalized parent length distribution $`n(L)=\rho ^{(0)}(L)/\rho _0^{(0)}`$.

We now consider the properties of the moment free energy with the (chain number) density $`\rho _0`$ retained as an extra moment. This provides additional geometrical insight into the properties of polydisperse chains (while for the numerical determination of the CPC and shadow curve, the above one-moment free energy is preferable). To construct the two-moment free energy, we proceed as before. The family (42) is now

$$\rho (L)=\rho ^{(0)}(L)e^{\lambda _0+\lambda _1L}=\rho _0^{(0)}\frac{1}{\mathrm{\Gamma }(a)b^a}L^{a1}e^{\lambda _0(b^1\lambda _1)L}$$
(84)
with zeroth and first moment densities

$`\rho _0=\rho _0^{(0)}e^{\lambda _0}\left({\displaystyle \frac{1}{1b\lambda _1}}\right)^a,\rho _1=\rho _1^{(0)}e^{\lambda _0}\left({\displaystyle \frac{1}{1b\lambda _1}}\right)^{a+1}.`$

The relations can be inverted to express the Lagrange multipliers $`\lambda _0`$ and $`\lambda _1`$ in terms of the moment densities as

$`\lambda _0`$ $`=`$ $`(a+1)\mathrm{ln}{\displaystyle \frac{\rho _0}{\rho _0^{(0)}}}a\mathrm{ln}{\displaystyle \frac{\rho _1}{\rho _1^{(0)}}}`$
$`\lambda _1`$ $`=`$ $`{\displaystyle \frac{1}{b}}\left(1{\displaystyle \frac{\rho _0/\rho _0^{(0)}}{\rho _1/\rho _1^{(0)}}}\right)`$
and one obtains the moment entropy

$$s_\mathrm{m}=\rho _0\lambda _0\rho _0\lambda _1\rho _1=(a+1)\rho _0\mathrm{ln}\frac{\rho _0}{\rho _0^{(0)}}+a\rho _0\mathrm{ln}\frac{\rho _1}{\rho _1^{(0)}}+(a+1)\rho _0a\frac{\rho _0^{(0)}}{\rho _1^{(0)}}\rho _1.$$
(85)
As expected, this reduces to (76) for systems with density distributions from our earlier one-parameter family (74), where $`\rho _0`$ can be expressed as a function of $`\rho _1`$ according to (75). With both $`\rho _0`$ and $`\rho _1`$ retained as moment densities in the moment free energy, a number of linear terms in (85) can be dropped, giving the final result

$$f_\mathrm{m}=(a+1)\rho _0\mathrm{ln}\rho _0a\rho _0\mathrm{ln}\rho _1+\stackrel{~}{f}.$$
(86)
Note that the dependence on the parent distribution is now only through $`a`$. This can be understood from the general discussion in Sec. III G: The family (84) of density distributions now contains all Schulz-distributions with the given $`a`$, and the moment free energy is insensitive to which member of this family (specified by $`\rho _0^{(0)}`$ and $`\rho _1^{(0)}`$) is used as the parent.

The result (86) can of course also be obtained via the combinatorial method, as follows. The normalized Schulz parent distribution is given by $`n^{(0)}(L)=[\mathrm{\Gamma }(a)b^a]^1L^{a1}e^{bL}`$. The cumulant generating function is thus $`h(\lambda _1)=a\mathrm{ln}(1b\lambda _1)`$, giving $`m_1=h/\lambda _1=ab/(1b\lambda _1)`$. Dropping constants and linear terms in the average chain length $`m_1`$ (which do not affect the phase behavior), the generalized entropy of mixing per particle, $`s_{\mathrm{mix}}=h\lambda _1m_1`$ then becomes $`s_{\mathrm{mix}}=a\mathrm{ln}m_1`$. Adding the ideal gas term and the excess part of the free energy, we thus find the moment free energy density

$$f_\mathrm{m}=\rho _0\mathrm{ln}\rho _0a\rho _0\mathrm{ln}(\rho _1/\rho _0)+\stackrel{~}{f}$$
(87)
in agreement with (86); it is a simple two component free energy. In the conventional notation, eq. (87) would read

$`f_\mathrm{m}=\rho \mathrm{ln}\rho \alpha \rho \mathrm{ln}(\varphi /\rho )+(1\varphi )\mathrm{ln}(1\varphi )+\chi \varphi (1\varphi ).`$

The spinodal curve (SC) and critical point (CP) condition may now be calculated from eq. (87), using either the methods outlined in Sec. III A or the more traditional determinant conditions  (see also App. D). One finds

$`{\displaystyle \frac{1}{1\rho _1}}2\chi +{\displaystyle \frac{a}{a+1}}{\displaystyle \frac{\rho _0}{\rho _1^2}}=0,`$ $`(\mathrm{SC})`$ (88)
$`{\displaystyle \frac{1}{(1\rho _1)^2}}{\displaystyle \frac{a(a+2)}{(a+1)^2}}{\displaystyle \frac{\rho _0}{\rho _1^3}}=0.`$ $`(\mathrm{CP})`$ (89)
Evaluated at the parent, these conditions are identical to eqs. (78) and (81) – which were derived from the one-moment free energy – as they must be. But, as explained at the end of Sec. III A, they also hold more generally for all systems with Schulz density distributions with the given value of $`a`$; see the remarks after eq. (86). Of particular interest among these systems are those that differ from the parent only in their overall density while having the same normalized length distribution $`n(L)=n^{(0)}(L)`$. They form the “dilution line”; their number averaged chain length satisfies $`m_1=m_1^{(0)}=𝑑LLn^{(0)}(L)`$, which translates to $`\rho _1=L_N\rho _0`$ (where $`L_Nm_1^{(0)}=ab`$ is, as before, the number average chain length of the parent). Inserting this “dilution line constraint” into the above equations, we obtain as the spinodal and critical point conditions for systems on the dilution line

$`{\displaystyle \frac{1}{1\rho _1}}2\chi +{\displaystyle \frac{a}{a+1}}{\displaystyle \frac{1}{L_N\rho _1}}=0,`$ $`(\mathrm{SC})`$ (90)
$`{\displaystyle \frac{1}{(1\rho _1)^2}}{\displaystyle \frac{a(a+2)}{(a+1)^2}}{\displaystyle \frac{1}{L_N\rho _1^2}}=0.`$ $`(\mathrm{CP})`$ (91)
Again, agreement with the general results (80) and (83) is easily verified by noting that for Schulz length distributions, $`L_W=\rho _2/\rho _1=(a+1)b=L_N(a+1)/a`$ and $`L_Z=\rho _3/\rho _2=(a+2)b=L_N(a+2)/a`$ for the weight and $`Z`$-average chain lengths, respectively.

The above discussion focussed on spinodals, critical points, and cloud and shadow curves, all of which are found exactly by the moment approach. (As usual, “exact” is subject to the assumed validity of the original, truncatable free energy expression (71).) Within the moment approach, we have also calculated the full phase behavior in the $`(\rho _1,\rho _0)`$ plane. This was done numerically from the two-moment free energy (87), facilitated by a partial Legendre transformation from $`\rho _0`$ to the associated chemical potential $`\mu _0`$, which can be performed analytically in this case (and effectively brings one back to the one-moment free energy (77)). Fig. 6(a)–(c) shows binodal curves, tielines, spinodal curves and critical points in the $`(\rho _1,\rho _0)`$ plane, for three values of $`\chi `$. These plots give an appealing geometric insight into the problem of phase separation in this system. For instance, the slope of the tielines indicates a size partitioning effect: the more dense phases are enriched in long chains.

The heavy line in Fig. 6(a)–(c) is the dilution line constraint $`m_1=m_1^{(0)}`$ or $`\rho _1=L_N\rho _0`$. As discussed already, we are interested mainly in systems whose mean composition lies on this line. (With $`\chi `$ as an extra variable, this line becomes a plane ($`\rho _1=L_N\rho _0,\chi `$) in the space $`(\rho _1,\rho _0,\chi )`$.) Not all of the phase behavior shown in Fig. 6 is then accessible. The extremities of phase separation on the $`\rho _1=L_N\rho _0`$ line are points where phase separation just starts to occur, and the locus of these points in the $`(\rho _1=L_N\rho _0,\chi )`$ plane defines the cloud curve, which bounds the phase separation region. This is shown in Fig. 6(d). This plot also includes the spinodal stability curve and the critical point from eqs. (90) and (91). The spinodal curve and the cloud curve touch at the critical point, which no longer lies at the minimum of either. This distorted behavior is a well known feature of polydisperse systems. Here it is seen to be due to the way that *regular* phase behavior in $`(\rho _1,\rho _0,\chi )`$ space is cut through by the dilution line constraint. Note that only two moments are needed to understand this qualitative effect, and that Fig. 6 gives direct geometrical insight into it, which would be very hard to extract from any exact solution based on (61) and (62).

The compositions of the phases that just start to appear as the phase separation region is entered do not in general lie on the dilution line $`\rho _1=L_N\rho _0`$. These phases are the “shadows”; they lie at the other ends of the cloud point tielines in Fig. 6(a)–(c). This is illustrated for $`\chi =0.62`$ in Fig. 7. The shadow phase compositions may be projected onto the ($`\rho _1=L_N\rho _0`$, $`\chi `$) plane to give shadow curves, but the projection is not unique. The shadow curve obtained by ignoring the value of $`\rho _0`$ and projecting onto the $`\rho _1`$-axis, for example, is shown in Fig. 8(a). That obtained by ignoring the value of $`\rho _1`$ and projecting onto the $`\rho _0`$-axis is shown in Fig. 8(b). The different shapes are due to the fact that the locus of shadow phase compositions in general only intersects the $`(\rho _1=L_N\rho _0,\chi )`$ plane at a single point—the critical point.

In conclusion, we emphasize again that all the curves shown in Fig. 6(d) and Fig. 8 are exact, even though they have been constructed from the phase behavior shown in Fig. 6(a)–(c) which is itself a projection of the true phase behavior of the fully polydisperse system.

### B Copolymer with chemical polydispersity

Next, we apply the moment free energy formalism to random AB-copolymers in the melt or in solution . For simplicity, we assume that the copolymer chain lengths are monodisperse, so that each chain contains an identical number $`L`$ of monomers. Even with this simplification, the problem is quite challenging because, as we shall see, a given parent can split into an arbitrarily large number of phases as the interaction strength is increased.

For this problem we choose as our polydisperse feature, $`\sigma `$, the difference between the fractions of A and B-monomers on a chain. This takes values in the range $`\sigma =1\mathrm{}1`$, with the extreme values $`\pm 1`$ corresponding to pure A and B chains, respectively. We now measure volumes in units of the chain volume (which is $`L`$ times larger than the monomeric volume used in the previous example). The total density $`\rho _0=𝑑\sigma \rho (\sigma )`$ is then just the copolymer volume fraction $`\varphi `$. (Note that for the case of length polydispersity, $`\varphi `$ was related to the first moment of the length distribution; here it is the zeroth moment of the distribution of chain compositions $`\sigma `$.) Within Flory-Huggins theory, the free energy density (in units of $`k_\text{B}T`$ per chain volume) is

$$f=𝑑\sigma \rho (\sigma )[\mathrm{ln}\rho (\sigma )1]+\frac{L}{L_\mathrm{s}}(1\rho _0)\mathrm{ln}(1\rho _0)+L\stackrel{~}{\chi }\rho _0(1\rho _0)L\stackrel{~}{\chi }^{}\rho _1^2$$
(92)
The interaction between A and B monomers enters through $`\rho _1=𝑑\sigma \sigma \rho (\sigma )`$. In (92), we have specialized to the case where the interactions of the A and B monomers with each other and with the solvent are symmetric in A and B; otherwise, there would be an additional term proportional to $`\rho _0\rho _1`$. We have also allowed the solvent to be a polymer of length $`L_\mathrm{s}`$; the case of a monomeric solvent is recovered for $`L_\mathrm{s}=1`$. In the following, we absorb factors of $`L`$ into the interaction parameters by setting

$`\chi =L\stackrel{~}{\chi },\chi ^{}=L\stackrel{~}{\chi }^{}`$

Defining also $`r=L/L_\mathrm{s}`$, the excess free energy takes the form

$$\stackrel{~}{f}=r(1\rho _0)\mathrm{ln}(1\rho _0)+\chi \rho _0(1\rho _0)\chi ^{}\rho _1^2.$$
(93)
The moment free energy (40), with only $`\rho _0`$ and $`\rho _1`$ retained, is then

$$f_\mathrm{m}=\lambda _0\rho _0+\lambda _1\rho _1\rho _0+r(1\rho _0)\mathrm{ln}(1\rho _0)+\chi \rho _0(1\rho _0)\chi ^{}\rho _1^2.$$
(94)
If extra moment densities with weight functions $`w_i(\sigma )`$ are included, these simply add a term $`\lambda _i\rho _i`$ each.

#### 1 Copolymer without solvent

We first consider the special case of a copolymer melt, for which the overall density (i.e., copolymer volume fraction) is constrained to take everywhere its maximum value $`\rho _0=1`$. The free energy (94) then simplifies to

$$f_\mathrm{m}=\lambda _0+\lambda _1\rho _1\chi ^{}\rho _1^2$$
(95)
up to irrelevant constants . Note that the same expression for the excess free energy of the dense system is obtained for a model in which the interaction energy between two species $`\sigma ,\sigma ^{}`$ varies as $`(\sigma \sigma ^{})^2`$:

$$\stackrel{~}{f}=\frac{1}{2}\chi ^{}𝑑\sigma 𝑑\sigma ^{}\rho (\sigma )\rho (\sigma ^{})(\sigma \sigma ^{})^2=\chi ^{}\rho _1^2+𝑑\sigma \rho (\sigma )\sigma ^2$$
(96)
The second term on the r.h.s. is linear in $`\rho (\sigma )`$ and can be discarded for the purposes of phase equilibrium calculations. The model (95) can therefore also be viewed as a simplified model of chemical fractionation, with $`\sigma `$ being related to aromaticity, for example, or another smoothly varying property of the different polymers. Eq. (96) suggests that the model should show fractionation into an ever-increasing number of phases as $`\chi ^{}`$ is increased: only the entropy of mixing prevents each value of $`\sigma `$ from forming a phase by itself, which would minimize this form of $`\stackrel{~}{f}`$. It is therefore an interesting test case for the moment free energy approach (and the method of adding further moment densities), yet simple enough for exact phase equilibrium calculations to remain feasible , allowing detailed comparisons to be made.

As a concrete example, we now consider phase separation from parent distributions of the form $`\rho ^{(0)}(\sigma )\mathrm{exp}(\gamma \sigma )`$ (for $`1\sigma 1`$, otherwise zero). The shape parameter $`\gamma `$ is then a fixed function of the parental first moment density $`\rho _1^{(0)}=𝑑\sigma \sigma \rho ^{(0)}(\sigma )`$. Fig 9 shows the exact coexistence curve for $`\rho _1^{(0)}=0.2`$, along with the predictions from our moment free energy with $`n`$ moment densities ($`\rho _i=𝑑\sigma \sigma ^i\rho (\sigma )`$, $`i=1\mathrm{}n`$) retained. Comparable results are found for other $`\rho _1^{(0)}`$. Even for the minimal set of moment densities ($`n=1`$) the point where phase separation first occurs on increasing $`\chi ^{}`$ is predicted correctly (this is the cloud point for the given parent). As more moments are added, the calculated coexistence curves approach the exact one to higher and higher precision. As expected, the precision decreases at high $`\chi ^{}`$, where fractionated phases proliferate; in this region, the number of coexisting phases predicted by the moment method increases with $`n`$. However, it is not always equal to $`n+1`$, as one might expect from a naive use of Gibbs’ phase rule; three-phase coexistence, for example, is first predicted for $`n=4`$. In fact, for fractionation problems such as this (but not more generally), study of the low temperature limit (large $`\chi ^{}`$) suggests that to obtain $`n+1`$ phases one may have to include up to $`2n`$ moments. We show below that using localized weight functions (rather than powers of $`\sigma `$) for the extra moments can reduce this number, but also gives less accurate coexistence curves.

In Fig. 10, we show how, for a given value of $`\chi ^{}`$, the lever rule is satisfied more and more accurately as the number of moment densities $`n`$ in the moment free energy is increased: the total density distribution $`\rho _{\mathrm{tot}}(\sigma )=_\alpha v^{(\alpha )}\rho ^{(\alpha )}(\sigma )`$ approaches the parent $`\rho ^{(0)}(\sigma )`$. This is by design, because the moment method forces the moment densities $`\rho _i=𝑑\sigma w_i(\sigma )\rho _{\mathrm{tot}}(\sigma )`$ of $`\rho _{\mathrm{tot}}(\sigma )`$ to be equal to those of the parent, for $`i=1\mathrm{}n`$; as $`n`$ is increased, therefore, $`\rho _{\mathrm{tot}}(\sigma )`$ approaches $`\rho ^{(0)}(\sigma )`$.

The above results for the dense copolymer model raise two general questions. First, how large does $`n`$ (the number of moment densities retained) have to be for the moment method to give reliable results for the coexisting phases (and, in particular, the correct number of phases)? And second, how should the weight functions for the extra moment densities be chosen? Regarding the first question, note first that the theoretical framework developed so far only says that the exact results (and therefore the correct number of phases) will be approached as $`n`$ gets large; it does not say for which value of $`n`$ a reasonable approximation is obtained. In fact it is clear that no “universal” (problem-independent) value of $`n`$ can exist. The model studied in this section already provides a counterexample: As $`\chi ^{}`$ is increased, the (exact) number $`p`$ of coexisting phases increases without limit; therefore the minimum value of $`n`$ ($`n=p1`$) required to obtain this correct number of phases from a moment free energy calculation can also be arbitrarily large.

Nevertheless, we see already from Fig. 9 that, for the current example, the robustness of the results to the addition of extra moment densities provides a reasonable qualitative check of the accuracy of the coexistence curves. To study the convergence of the moment free energy results to the exact ones more quantitatively, especially in cases where the latter are not available for comparison, we need a measure of the deviation between $`\rho _{\mathrm{tot}}(\sigma )`$ and $`\rho ^{(0)}(\sigma )`$. To obtain a dimensionless quantity that does not scale with the overall density of the parent, we consider the normalized distributions $`n_{\mathrm{tot}}(\sigma )=\rho _{\mathrm{tot}}(\sigma )/\rho _{\mathrm{tot}}`$ and $`n^{(0)}(\sigma )=\rho ^{(0)}(\sigma )/\rho ^{(0)}`$ and define as our measure of deviation the “log-error”

$$\delta =𝑑\sigma n^{(0)}(\sigma )\left(\mathrm{ln}\frac{n^{(0)}(\sigma )}{n_{\mathrm{tot}}(\sigma )}\right)^2$$
(97)
This becomes zero only when $`n_{\mathrm{tot}}(\sigma )=n^{(0)}(\sigma )`$; otherwise it is positive. When the deviations between the two distributions are small, the logarithm becomes to leading order the relative deviation $`n_{\mathrm{tot}}(\sigma )/n^{(0)}(\sigma )1`$ and we can identify $`\delta ^{1/2}`$ as the root-mean-square relative deviation. We work with the logarithm rather than directly with the relative deviation because the former gives more sensible behavior for larger deviations; in particular, isolated points where $`n_{\mathrm{tot}}(\sigma )`$ is nonzero while $`n^{(0)}(\sigma )`$ is close to zero do not lead to divergences of $`\delta `$ .

Fig. 11 shows the $`n`$-dependence of the log-error $`\delta `$ for three fixed values of $`\chi ^{}`$, with the same parent as in Fig 9 . The correct number of phases (3, 4, 5, respectively) is reached at values of $`\delta 4\times 10^3`$, $`5\times 10^4`$, $`2\times 10^3`$ for $`\chi =5`$, 10, 15, respectively. This suggests that the log-error $`\delta `$ might generally be a useful heuristic criterion for guiding the choice of $`n`$, with values of $`\delta `$ around $`10^4`$ (corresponding to an average deviation between $`n_{\mathrm{tot}}(\sigma )`$ and $`n^{(0)}(\sigma )`$ of around 1%) ensuring that the correct number of phases has been detected. Note that, although for the current example we did know the correct number of phases (because the simplicity of the problem made an exact calculation feasible), the definition of $`\delta `$ does not require us to know in advance any properties of the exact solution of the coexistence problem. Also, note that in real physical materials it is extremely unlikely for the parent or “feedstock” distribution $`\rho ^{(0)}(\sigma )`$ to be known to better than 1% accuracy. Hence even for systems, such as polymers, where a mean-field (and thus truncatable) free energy model is thought to be highly accurate, deviations from the lever rule, at about the 1% level, will almost never be the main source of uncertainty in phase diagram prediction.

We now turn to the second general question, regarding the choice of weight functions for the extra moment densities. Comparing (63) with the formally exact solution (62) of the coexistence problem tells us at least in principle what is required: the log-ratio $`\mathrm{ln}\stackrel{~}{R}(\sigma )/\rho ^{(0)}(\sigma )`$ of the “effective prior” and the parent needs to be well approximated by a linear combination of the weight functions of the extra moment densities. However, the effective prior is unknown (otherwise we would already have the exact solution of the phase coexistence problem), and so this criterion is of little use .

To make progress, let us try to gain some intuition from the example at hand (dense random copolymers). So far, we have simply taken increasing powers of $`\sigma `$ for the extra weight functions. On the other hand, looking at Fig. 10 one might suspect that localized weight functions might be more suitable for capturing the fractionating of the parent into the various daughter phases. We therefore tried the simple “triangle” and “bell” weight functions shown in Fig. 12(a). For a given $`n`$, there are $`n1`$ extra moment densities (the first moment density, with weight function $`w_1(\sigma )=\sigma `$, is prescribed by the excess free energy) and we take these to be evenly spaced across the range $`\sigma =1\mathrm{}1`$, with overlaps as indicated in Fig 12(a). (For the case of triangular weight functions, the log ratio $`\mathrm{ln}\rho (\sigma )/\rho ^{(0)}(\sigma )=_i\lambda _iw_i(\sigma )`$ between the density distribution of any phase and that of the parent is thus approximated by a “trapezoidal” function across $`n`$ equal intervals spanning $`1\mathrm{}1`$.) Fig. 12(b) shows the corresponding coexistence curves for the triangle case with $`n=1,2,3,4`$. Note that multi-phase coexistence is generally detected for smaller $`n`$ than for the case of power-law weight functions (Fig. 9), in agreement with the expectations outlined above. However, the predicted number of phases does not vary monotonically with $`n`$: in the range of $`\chi ^{}`$ shown in Fig. 12(b), up to four phases are predicted for $`n=3`$, whereas $`n=4`$ gives no more than three. This effect illustrates a significant difference between power-law and triangle weight functions. The former, but not the latter, make up an incremental set: increasing $`n`$ by one corresponds to adding a new weight function while leaving the existing set of weight functions unchanged. The family (42) of accessible density distributions $`\rho (\sigma )`$ is thus enlarged without losing any of its former members. Intuitively, one then expects that the number of phases cannot decrease as $`n`$ is increased. But for nonincremental sets, such as triangle weight functions (and bells, which give qualitatively similar results to triangles), each increase in $`n`$ creates a larger family but completely discards the previous one, so that there is then no reason for the results to be monotonic in $`n`$.

We regard the fact that the number of phases can go up and down with increasing $`n`$ as a disadvantage of non-incremental weight functions. In the current example, we also find that they generally give higher values of the log-error than power-law weight functions; see Fig. 12(c). Finally, incremental sets of weight functions also have the benefit that the results for $`n`$ can be used to initialize the calculation for $`n+1`$, thus helping to speed up the numerics (though we have not yet exploited this property). On the basis of these observations, we would generally recommend the use of incremental sets of weight functions for extra moment densities. Beyond this, the optimal choice of extra weight functions remains largely an open problem. One avenue worth exploring in future work might be the possibility of choosing extra weight functions adaptively. One could imagine monitoring the log-error in a phase coexistence calculation, and increasing the number of extra moments by one every time it becomes larger than a certain threshold ($`10^4`$, say). The new weight function could then be chosen as a best fit to the current lever-rule violation $`\mathrm{ln}\rho _{\mathrm{tot}}(\sigma )/\rho ^{(0)}(\sigma )`$ (for example as a linear combination of weight functions from some predefined “pool”). We are planning to investigate this method in future work.

To conclude this section on the dense random copolymer model, we briefly discuss the spinodal criterion, and ask whether critical or multi-critical points can exist for a general parent distribution $`\rho ^{(0)}(\sigma )`$ (with $`\rho _0^{(0)}=1`$). The criterion (56), applied to our one-moment free energy $`f_\mathrm{m}(\rho _1)`$ becomes simply

$$f_\mathrm{m}(\rho _1)=𝒪\left((\rho _1\rho _1^{(0)})^{l+1}\right)$$
(98)
where $`l=2`$ for a spinodal and $`l=2n1`$ for an $`n`$-critical point. From eq. (95), we have that

$$f_\mathrm{m}=\lambda _0+\lambda _1\rho _1\chi ^{}\rho _1^2=s_\mathrm{m}\chi ^{}\rho _1^2$$
(99)
The excess part $`\chi ^{}\rho _1^2`$ is quadratic in $`\rho _1`$ and therefore only enters the spinodal criterion; the additional conditions for critical points only depend on the moment entropy $`s_\mathrm{m}=\lambda _0\lambda _1\rho _1`$. Because of the dense limit constraint $`\rho _1=1`$, $`\lambda _0`$ is related to $`\lambda _1`$ by

$`\lambda _0=\mathrm{ln}{\displaystyle 𝑑\sigma \rho ^{(0)}(\sigma )e^{\lambda _1\sigma }}h(\lambda _1)`$

As expected from the general discussion in Sec. II C, $`\lambda _0`$ is thus the cumulant generating function $`h(\lambda _1)`$ of $`\rho ^{(0)}(\sigma )`$ ($`n^{(0)}(\sigma )`$ in the dense limit we are considering), introduced in eq. (36) within the combinatorial approach. The moment entropy $`s_\mathrm{m}=h\lambda _1\rho _1`$ is its Legendre transform; its behavior for $`\rho _1\rho _1^{(0)}`$ reflects that of $`h`$ for small $`\lambda _1`$. We can thus use the expansion $`h=_{j=1}^{\mathrm{}}(c_j/j!)\lambda _1^j`$, where the $`c_j`$ are the cumulants of $`\rho ^{(0)}(\sigma )`$; $`c_1\rho _1^{(0)}`$ is its mean, $`c_2`$ its variance etc.. Using the properties of the Legendre transform, it then follows that the leading term of the moment entropy is $`s_\mathrm{m}=(\rho _1\rho _1^{(0)})^2/(2c_2)`$. Inserting this into (99), we see immediately that the spinodal condition becomes

$$2\chi ^{}c_2=2\chi ^{}\left(\rho _2^{(0)}(\rho _1^{(0)})^2\right)=1$$
(100)
To find the additional conditions for critical points, consider what happens if the cumulants $`c_3`$, …, $`c_l`$ are zero, hence $`h=\rho _1^{(0)}\lambda _1+c_2\lambda _1^2/2+𝒪(\lambda _1^{l+1})`$. It can then easily be shown that the moment entropy behaves as $`s_\mathrm{m}=(\rho _1\rho _1^{(0)})^2/(2c_2)+𝒪((\rho _1\rho _1^{(0)})^{l+1})`$; conversely, one can show that this behavior of $`s_\mathrm{m}`$ implies that the cumulants $`c_3`$, …, $`c_l`$ are zero. We thus have the simple condition that for an $`n`$-critical point to be observed, the cumulants $`c_3`$, $`\mathrm{}`$, $`c_{2n1}`$ of the parent distribution have to be zero; the value of $`\chi ^{}`$ at this critical point is given by the spinodal condition (100). (For an ordinary critical point, $`n=2`$, the condition $`c_3=0`$ is well known in this context .) We therefore come to the somewhat surprising conclusion that even in a very simple polydisperse system such as this – with a single moment density occurring in the excess free energy – multi-critical points of arbitrary order can occur, at least in principle. In practice, the required fine-tuning of the parent will of course make experimental study of these points difficult.

#### 2 Copolymer with solvent

We now consider the random copolymer model in the presence of solvent, i.e., for a copolymer volume fraction $`\rho _1<1`$. We are not aware of previous work on this model in the literature, but will briefly discuss below the link to models of homopolymer/copolymer mixtures . The excess free energy (93) then depends on two moment densities, rather than just one as in all previous examples. For simplicity, we restrict ourselves to the case of a neutral solvent which does not in itself induce phase separation; this corresponds to $`\chi =0`$, making the excess free energy

$$\stackrel{~}{f}=r(1\rho _0)\mathrm{ln}(1\rho _0)\chi ^{}\rho _1^2$$
(101)
This is still sufficiently simple that the exact spinodal condition (58) – for a system with a general density distribution $`\rho (\sigma )`$ – can be worked out analytically; one finds

$$2\chi ^{}=\frac{r\rho _0+1\rho _0}{r(\rho _2\rho _0\rho _1^2)+\rho _2(1\rho _0)}$$
(102)
In the dense limit $`\rho _01`$, the earlier result (100) is recovered. To study critical and multi-critical points, we use the moment free energy corresponding to (101),

$$f_\mathrm{m}=\lambda _0\rho _0+\lambda _1\rho _1\rho _0+r(1\rho _0)\mathrm{ln}(1\rho _0)\chi ^{}\rho _1^2$$
(103)
and start from the criterion (56). The moment chemical potentials are given by

$`\mu _0`$ $`=`$ $`\lambda _0r[\mathrm{ln}(1\rho _0)1]`$ (104)
$`\mu _1`$ $`=`$ $`\lambda _12\chi ^{}\rho _1`$ (105)
and the curve $`𝝆(ϵ)`$ referred to in (56) – along which the critical phases merge – is simply parameterized in terms of $`\lambda _0(ϵ)`$ and $`\lambda _1(ϵ)`$. The restriction that the curve must pass through the parent at $`ϵ=0`$ means that $`\lambda _0(0)=\lambda _1(0)=0`$, so to get the small $`ϵ`$ behavior of the chemical potentials (105) we can expand for small $`\lambda _0`$ and $`\lambda _1`$. To simplify the algebra, we restrict ourselves to parent distributions $`\rho ^{(0)}(\sigma )`$ which are symmetric about $`\sigma =0`$. Using this symmetry and the fact that $`ϵ`$ is simply a dummy parameter which can be reparameterized arbitrarily, we can then set $`\lambda _1=ϵ`$ and $`\lambda _0=a_2ϵ^2+a_4ϵ^4+𝒪(ϵ^6)`$.

To have a spinodal, the moment chemical potentials (105) must differ from those of the parent by no more than $`𝒪(ϵ^2)`$. We find that this gives the condition $`2\chi ^{}=1/\rho _2^{(0)}`$, in agreement with (102): symmetric parents have $`\rho _1^{(0)}=0`$. For a critical point, the chemical potential difference must be decreased to $`𝒪(ϵ^3)`$. This gives no new condition (but determines $`a_2`$), implying that parents with symmetric distributions are automatically critical at their spinodal. This result has its origin in the symmetry $`\sigma \sigma `$ of the excess free energy, and can be understood by considering the $`\rho _0`$, $`\rho _1`$, $`\chi ^{}`$ phase diagram of the moment free energy (103). In this phase diagram there is a two-dimensional spinodal surface, and within this surface lies a line of critical points. The symmetry $`\sigma \sigma `$ (corresponding to $`\rho _1\rho _1`$, $`\lambda _1\lambda _1`$) forces this critical line to lie within the plane $`\rho _1=0`$ of symmetrical distributions, implying that for such distributions spinodal and critical points coincide.

For a tricritical point, finally, where the chemical potential difference must be further reduced to $`𝒪(ϵ^3)`$, we find after some algebra the condition

$`\rho _4^{(0)}={\displaystyle \frac{3r(\rho _2^{(0)})^2}{r\rho _0^{(0)}+1\rho _0^{(0)}}}`$

As an example, consider a uniform parent distribution $`\rho ^{(0)}(\sigma )=\text{const}`$, which can be written as $`\rho ^{(0)}(\sigma )=\rho _0^{(0)}/2`$, using the fact that $`\rho _0^{(0)}=_1^1𝑑\sigma \rho ^{(0)}(\sigma )`$. Then $`\rho _2^{(0)}=\rho _0^{(0)}/3`$ and $`\rho _4^{(0)}=\rho _0^{(0)}/5`$ and so a tricritical point occurs if the overall density (i.e., the copolymer volume fraction) is $`\rho _0^{(0)}=3/(2r+3)`$. Fig. 13 shows the coexistence curve calculated for this parent (with $`r=1`$), which clearly shows the tricritical point at the predicted value $`\chi ^{}=1/(2\rho _2^{(0)})=r+3/2=2.5`$. Our numerical implementation manages to locate the tricritical point and follow the three coexisting phases without problems; we take that as a signature of its robustness . Note that the tricritical point that we found is closely analogous to that studied by Leibler  for a symmetric blend of two homopolymers and a symmetric random copolymer that is, nonetheless, chemically monodisperse (in the sense that $`\sigma =0`$ for all copolymers present). In fact, in our notation, the scenario of Ref.  simply corresponds to a parent density of the form $`\rho ^{(0)}(\sigma )\delta (\sigma 1)+\delta (\sigma +1)`$, with the copolymer ($`\sigma =0`$) now playing the role of the neutral solvent.

We conclude this section with another illustration of the geometrical intuition which the moment free energy can provide. In Fig. 14, we plot the moment free energy $`f_\mathrm{m}(\rho _0,\rho _1)`$, eq. (103), for the uniform parent $`\rho ^{(0)}(\sigma )`$ studied above. Two values of $`\chi ^{}`$ are shown: the tricritical value ($`\chi ^{}=2.5`$), and a value just in the three-phase region ($`\chi ^{}=2.65`$); the developing three minima of the free energy can clearly be seen. Even though the underlying polydisperse system has an infinite number of degrees of freedom $`\rho (\sigma )`$, the moment method thus gives us a tool for understanding and visualizing the occurrence of the tricritical point in terms of a simple free energy surface with two independent variables.

## VI Conclusion

Polydisperse systems contain particles with an attribute $`\sigma `$ which takes values in a continuous range (such as particle size in colloids or chain length in polymers). They therefore have an infinite number of conserved densities, corresponding to a density distribution $`\rho (\sigma )`$. The fact that the free energy can depend on all details of $`\rho (\sigma )`$ makes the analysis of the phase behavior of such systems highly non-trivial. However, in many (especially mean-field) models the excess free energy only depends on a finite number of moment densities, i.e., (generalized) moments of $`\rho (\sigma )`$; the only dependence on other details of $`\rho (\sigma )`$ is through the ideal part of the free energy. For such models, which we call truncatable, we showed how a moment free energy can be constructed, which only depends on the moment densities appearing in the excess free energy. Two possible approaches for the construction of the moment free energy were described in Sec. II. The first, which we call the projection method (Sec. II A), is primarily geometrical; the second is mainly combinatorial (Sec. II B) and starts from first principles to derive an expression for the ideal part of the free energy in terms of moment densities. In Sec. II C, we showed that these two approaches give essentially equivalent results.

In Sec. III, we explored the properties of the moment free energy in detail. The moment free energy depends only on a finite number of moment densities, and can be used to predict phase behavior using the conventional methods known from the thermodynamics of finite mixtures. We showed that, for any truncatable model, this procedure yields the same spinodals, critical and multi-critical points as the original free energy, even though the latter depends on all details of the density distribution $`\rho (\sigma )`$. It also correctly predicts the onset of phase coexistence, i.e., the cloud points and corresponding shadows. Beyond the onset of phase coexistence (i.e., when one or more coexisting phases occupy comparable volumes), the results from the moment free energy are approximate, but can be made arbitrarily accurate by retaining additional moment densities in the moment free energy. We also explained how to check the local and global stability of a phase coexistence calculated from the moment free energy, and how the moment free energy reflects overall phase behavior in the space of density distributions $`\rho (\sigma )`$.

We discussed aspects of the numerical implementation of the moment free energy method in Sec. IV, highlighting its many advantages in terms of robustness, efficiency and stability. The method allows violations of the lever rule, but then ensures that these be kept small by adding $`nK`$ extra moments (where $`K`$ is the number of moment densities actually present in the excess free energy of the underlying truncatable model). The solution of the phase equilibrium conditions thus proceeds in an $`n`$-dimensional space; but global (tangent plane) stability can always be verified by minimimization (of the tangent plane distance) in a space of dimension $`K`$. For $`K>1`$ exact minimization is not practicable, but standard search methods can instead be employed; we chose a Monte-Carlo type algorithm for this.

We then surveyed several sample applications in Sec. V. For length-polydisperse homopolymers (treated within Flory-Huggins theory), we showed explicitly how the moment free energy with only a single moment density (the volume fraction of polymer) gives the correct spinodals and critical points. By retaining the chain number density as an additional moment density, we could rationalize the behavior of cloud point and shadow curves as cuts through conventional two-species phase diagrams. As a second example, the case of copolymers with chemical polydispersity was considered, both without and with solvent. This gave us the opportunity to test the performance of the moment free energy method for multi-phase coexistence, and in particular to study how the accuracy of the predictions in the coexistence region increases as more and more moment densities are retained. In the presence of solvent, the same theory predicts a tricritical point under appropriate conditions, which was detected without problems, along with the associated three phase region, by our numerical algorithm. A plot of the two-moment free energy gave a simple geometrical picture of the origin of this tricritical point.

The results for the copolymer model showed that the log-error (97) can be used to decide when enough moment densities have been included to locate the correct number of phases and their compositions. For general purposes we suggest a log-error tolerance of order $`10^4`$; this corresponds to limiting the rms lever-rule violations at the predicted phase coexistence (found by global minimization of the moment free energy) to 1% or so of the parent distribution. Given that in practice the parent is itself not known to such high precision, these violations lie well within an acceptable range of error. Indeed, every state of coexistence actually predicted by our algorithm represents the exact phase equilibrium of some or other parent, indistinguishable from the real one to within experimental accuracy. Here, as elsewhere in this paper, “exact” means exactly as predicted by a truncatable model free energy; such models are, of course, themselves approximate, which is another reason not to expend resources reducing the log-error below the level suggested above.

In summary, with exact results for cloud point and shadow curves, critical points and spinodals, as well as refinably accurate coexistence curves and multi-phase regions, the moment free energy method allows rapid and accurate computation of the phase behavior of many polydisperse systems. Moreover, by establishing the link to a projected free energy $`f_\mathrm{m}(\rho _i)`$ which is a function of a finite set of conserved moment densities $`\rho _i`$, it restores to the problem much of the geometrical interpretation and insight (as well as the computational methodology) associated with the phase behavior of finite mixtures. This contrasts with many procedures previously in use for truncatable polydisperse systems . Some previous approximations to the problem have used (generalized) moments as thermodynamic coordinates; see e.g. . The moment free energy method provides a rational basis for these methods and, so long as the generalized moments are correctly chosen, guarantees that many properties of interest are found exactly. Note also that the commonplace method of “binning” the $`\sigma `$-distribution into discrete “pseudo-components” can be seen as a special case of the moment method in which each weight function is zero outside the corresponding bin; but since these moments do not, in general, coincide with those contained in the excess free energy the advantages described above are then lost.

We finish with some remarks about extensions to the moment method, within and outside the sphere of phase diagram prediction. First, methods of this type may extend to nontruncatable models, for which the excess free energy cannot be written directly in terms of a finite number of moments as in (6). For example, many mean-field theories correspond to a variational minimization of the free energy: $`FE_0TS_0`$, where the subscript $`0`$ refers to a trial Hamiltonian . In such a case, one might choose to first make a physically motivated decision about which (and how many) moment densities $`\rho _i`$ to keep, and then include among the variational parameters the “transverse” degrees of freedom of the density distribution $`\rho (\sigma )`$. We have not pursued this approach, though it may form a promising basis for future developments. Secondly, we note that physical insight based on moment free energies may increasingly play a part in understanding kinetic problems. For example, one of us  has argued that in many systems the zeroth moment (mean particle density) can relax more rapidly (by collective diffusion) than can the higher moments (by interdiffusion of different species), and that this may lead to novel kinetics in certain areas of the phase diagram. And finally, we point to recent work  in which the moment method has been extended, via density functional theory, to allow study of inhomogeneous states; this opens the way to studying the effect of polydispersity on interfacial tensions and other interfacial thermodynamic properties.

Acknowledgements: We thank P. Bladon, N. Clarke, R. M. L. Evans, T. McLeish, P. Olmsted, I. Pagonabarraga, W. C. K. Poon and R. Sear for helpful discussions. PS is grateful to the Royal Society for financial support. This work was also funded in part under EPSRC GR/M29696 and GR/L81185.

## A Moment (Gibbs) free energy for fixed pressure

In this appendix, we explain how the construction of the moment free energy is extended to scenarios where the pressure (rather than the overall system volume) is held constant. One then describes the system by the particle number distribution $`N(\sigma )=V\rho (\sigma )`$; this is normalized so that its integral gives the total number of particles, $`N=𝑑\sigma N(\sigma )`$. The relevant thermodynamic potential is the Gibbs free energy $`G=\mathrm{min}_VF+\mathrm{\Pi }V`$, which we now construct. For a truncatable system, whose free energy density obeys eqs. (3,4) by definition, the Helmholtz free energy can be written as

$`F=T{\displaystyle 𝑑\sigma N(\sigma )\left[\mathrm{ln}\frac{N(\sigma )}{V}1\right]}+\stackrel{~}{F}`$

where the excess free energy $`\stackrel{~}{F}(N_i,V)=V\stackrel{~}{f}(\rho _i)`$ depends on volume $`V`$ and on the moments $`N_i=V\rho _i=𝑑\sigma N(\sigma )w_i(\sigma )`$ of the particle number distribution. If $`\stackrel{~}{F}=0`$, a Legendre transform gives the Gibbs free energy of an ideal mixture:

$`G_{\mathrm{id}}=\underset{V}{\mathrm{min}}F_{\mathrm{id}}+\mathrm{\Pi }V=T{\displaystyle 𝑑\sigma N(\sigma )\mathrm{ln}\frac{\mathrm{\Pi }N(\sigma )}{NT}}=NT{\displaystyle 𝑑\sigma n(\sigma )\mathrm{ln}n(\sigma )}+NT\mathrm{ln}{\displaystyle \frac{\mathrm{\Pi }}{T}}`$

Here $`n(\sigma )=N(\sigma )/N=\rho (\sigma )/\rho `$ is the normalized particle distribution. In the general case, we split off the ideal part

$`G=G_{\mathrm{id}}+\stackrel{~}{G}=\underset{V}{\mathrm{min}}F+\mathrm{\Pi }V`$

and obtain for the excess Gibbs free energy

$`\stackrel{~}{G}=\underset{V}{\mathrm{min}}NT\mathrm{ln}{\displaystyle \frac{NT}{\mathrm{\Pi }V}}+\stackrel{~}{F}(N_i,V)+\mathrm{\Pi }V`$

This is a function of $`N`$ and the $`N_i`$ (as well as the fixed control parameters $`\mathrm{\Pi },T`$, which we suppress in our notation). We can therefore write the excess Gibbs free energy per particle, $`\stackrel{~}{G}(N_i,N)/N=\stackrel{~}{g}(m_i)`$, as a function of the moments $`m_i=N_i/N`$ ($`i=1\mathrm{}K`$) of the normalized particle distribution $`n(\sigma )`$. Altogether, the Gibbs free energy per particle of a truncatable system becomes

$$g[n(\sigma )]=G/N=T𝑑\sigma n(\sigma )\mathrm{ln}n(\sigma )+T\mathrm{ln}\frac{\mathrm{\Pi }}{T}+\stackrel{~}{g}(m_i)$$
(A1)
This has again a truncatable structure; the excess part $`\stackrel{~}{g}`$ “inherits” its moment structure from $`\stackrel{~}{f}`$. Note that because of the normalization of the $`m_i`$, $`\stackrel{~}{g}`$ does not depend on the density $`\rho \rho _0`$. It is therefore normally a function of one less variable than $`\stackrel{~}{f}`$, unless $`\stackrel{~}{f}`$ is already independent of $`\rho `$ (as is the case in Flory-Huggins theory, eq. (71), for example). The chemical potentials $`\mu (\sigma )`$ follow from (A1) as

$`\mu (\sigma )`$ $`=`$ $`{\displaystyle \frac{\delta }{\delta N(\sigma )}}Ng[N(\sigma )/N]`$ (A2)
$`=`$ $`{\displaystyle \frac{\delta g}{\delta n(\sigma )}}+g{\displaystyle 𝑑\sigma ^{}n(\sigma ^{})\frac{\delta g}{\delta n(\sigma ^{})}}`$ (A3)
$`=`$ $`T\mathrm{ln}n(\sigma )+{\displaystyle \underset{i}{}}w_i(\sigma ){\displaystyle \frac{\stackrel{~}{g}}{m_i}}+\left[T\mathrm{ln}{\displaystyle \frac{\mathrm{\Pi }}{T}}+\stackrel{~}{g}{\displaystyle \underset{i}{}}m_i{\displaystyle \frac{\stackrel{~}{g}}{m_i}}\right]`$ (A4)
Multiplying by $`n(\sigma )`$ and integrating over $`\sigma `$, one has $`g=𝑑\sigma \mu (\sigma )n(\sigma )`$ as it should be (compare eq. (A1)). Comparing (A4) with (48), one deduces that

$`{\displaystyle \frac{\stackrel{~}{g}}{m_i}}=\stackrel{~}{\mu }_i`$

and that the density $`\rho `$ of a phase with moments $`m_i`$ is given by

$`T\mathrm{ln}\left({\displaystyle \frac{T\rho }{\mathrm{\Pi }}}\right)=\stackrel{~}{g}{\displaystyle \underset{i}{}}m_i\stackrel{~}{\mu }_i.`$

To construct the moment Gibbs free energy, one can now proceed as in the constant volume case: in the ideal contribution, we replace $`\mathrm{ln}n(\sigma )\mathrm{ln}[n(\sigma )/n^{(0)}(\sigma )]`$, with a normalized parent distribution $`n^{(0)}(\sigma )`$, and then minimize $`g`$ at fixed values of the $`m_i`$. The minimum occurs for distributions $`n(\sigma )`$ from the family

$$n(\sigma )=n^{(0)}(\sigma )\mathrm{exp}\left(\lambda _0+\underset{i}{}\lambda _iw_i(\sigma )\right)$$
(A5)
The extra Lagrange multiplier $`\lambda _0`$ comes from the normalization condition $`m_0=𝑑\sigma n(\sigma )=1`$. Inserting into (A1), one obtains the moment Gibbs free energy

$$g_\mathrm{m}(m_i)=T\left(\lambda _0+\underset{i}{}\lambda _im_i\right)+T\mathrm{ln}\frac{\mathrm{\Pi }}{T}+\stackrel{~}{g}(m_i)$$
(A6)
where the $`m_i`$ are given by

$`m_i={\displaystyle 𝑑\sigma w_i(\sigma )n^{(0)}(\sigma )\mathrm{exp}\left(\lambda _0+\underset{i}{}\lambda _iw_i(\sigma )\right)}={\displaystyle \frac{𝑑\sigma w_i(\sigma )n^{(0)}(\sigma )\mathrm{exp}\left(_i\lambda _iw_i(\sigma )\right)}{𝑑\sigma n^{(0)}(\sigma )\mathrm{exp}\left(_i\lambda _iw_i(\sigma )\right)}}`$

The derivatives of $`g_\mathrm{m}`$ can again be identified as moment chemical potentials and obey

$$\mu _i\frac{g_\mathrm{m}}{m_i}=T\lambda _i+\frac{\stackrel{~}{g}}{m_i}=T\lambda _i+\stackrel{~}{\mu }_i.$$
(A7)
The constraint $`m_0=1`$ does not affect this result because the partial derivative is taken with the values of all moments other than $`m_i`$ fixed anyway. For the same reason, the second derivatives of the ideal part of $`g_\mathrm{m}`$ are (up to a factor $`T`$) given by the inverse of the matrix of second-order normalized moments $`m_{ij}`$, by analogy with (46). (Note that the first row and column of the second-order moment matrix, corresponding to $`i=0`$ or $`j=0`$, need to be retained; they can only be discarded after the inverse has been taken.)

From (A4), one deduces immediately that if one of a number of coexisting phases is in the family (A5), then so are all others. Equality of the exact chemical potentials $`\mu (\sigma )`$ in all phases is then equivalent to equality of the quantities

$$\stackrel{~}{\mu }_i+T\lambda _i=\mu _i(i=1\mathrm{}K),T\lambda _0+\stackrel{~}{g}\underset{i}{}m_i\stackrel{~}{\mu }_i.$$
(A8)
But from (A6,A7), these are simply the moment chemical potentials $`\mu _i=g_\mathrm{m}/m_i`$ and (up to an unimportant additive constant) the associated Legendre transform $`g_\mathrm{m}_im_i\mu _i`$. Therefore, phase equilibria can be constructed by applying the usual tangent plane construction to the moment Gibbs free energy, in analogy with the constant volume case. Note, however, that this tangent plane construction now takes place in the space of the normalized moments $`m_i`$, rather than that of the unnormalized moment densities $`\rho _i`$. As pointed out above, this generically reduces the dimension by one. For the trivial case of a monodisperse system, $`\stackrel{~}{f}`$ is a function of $`\rho `$ only and $`\stackrel{~}{g}`$ does not depend on any density variables at all, so that the tangent plane condition becomes void as expected.

The exactness statements in Sec. III can also be directly translated to the constant pressure case. The arguments above imply directly that the onset of phase coexistence is found exactly from the moment Gibbs free energy: all phases are in the family (A5), because one of them (the parent) is, and the requirement of equal chemical potentials is satisfied. Spinodals and (multi-) critical points are also found exactly. Arguing as in Sec. III A, and using the vector notation of (56), the criterion for such points is found as

$$\mathrm{\Delta }𝝁𝝁(𝒏(ϵ))𝝁(𝒏^{(0)})=𝒪(ϵ^l)$$
(A9)
where $`l=2`$ for a spinodal and $`l=2n1`$ for an $`n`$-critical point. Again, the curve $`𝒏(ϵ)`$ only has to pass through phases with equal chemical potentials and can therefore be chosen to lie within the family (A5). As discussed above (eq. (A8)), the chemical potential difference on the left hand side of (A9) is then zero to the required order in $`ϵ`$ if the same holds for the $`\mathrm{\Delta }\mu _i`$ and $`\mathrm{\Delta }h`$, where $`h=g_\mathrm{m}_im_i\mu _i`$. It may look as if the constraint on $`h`$ gives one more condition here than in the constant volume case; but in fact $`\mathrm{\Delta }\mu _i=𝒪(ϵ^l)`$ for all $`i`$ already implies that $`\mathrm{\Delta }g_\mathrm{m}=𝒪(ϵ^{l+1})`$ and hence $`\mathrm{\Delta }h=𝒪(ϵ^l)`$. The exact condition (A9) for spinodals and critical points is therefore again equivalent to the same condition obtained from the moment Gibbs free energy.

## B Moment entropy of mixing and large deviation theory

In this appendix, we discuss some interesting properties of the Legendre transform result (31) for the moment entropy of mixing, in particular its relation to large deviation theory (LDT).

We first note that (31), which we derived by taking the limit $`x0`$ of the result (28) for general $`x`$, can also be obtained by a more direct route. In the limit $`x0`$, the sizes of particles in the smaller system become independent random variables drawn from $`n^{(0)}(\sigma )`$; the second phase can be viewed as a reservoir to which the small phase is connected. One writes the moment generating function for $`𝒫(m)`$ in the small phase as a product of $`xN`$ independent moment generating functions of $`n^{(0)}(\sigma )`$, and evaluates the integral over $`𝒫(m)`$ by a saddle point method .

One can also view (31) as a generalisation of Cramér’s large deviation theory (LDT) . Recall that this theory treats the “wings” of a distribution like $`𝒫(m)`$ correctly, to which the central limit theorem (CLT) does not apply . However, the CLT approximation in this problem has quite a beautiful interpretation, and is worth describing separately. It can be derived by assuming that $`n^{(0)}(\sigma )`$ is itself a Gaussian in eq. (31), or just written down directly:

$`𝒫(m)\mathrm{exp}\left[N{\displaystyle \frac{(mm^{(0)})^2}{2v^{(0)}}}\right]`$

where $`v^{(0)}`$ is the variance of the parent distribution. We see that this gives a term in the free energy of the form

$`\mathrm{\Delta }f=T\rho s=\rho T{\displaystyle \frac{(mm^{(0)})^2}{2v^{(0)}}}.`$

This result states that there is an entropic spring term in the free energy that penalizes deviations of the mean size $`m`$ from the mean size in the parent distribution. The spring constant is inversely proportional to the variance in the parent distribution, thus if the parent is narrower, it is harder to move away from the parental mean size . However, as indicated above, the CLT is not sufficiently accurate for our purposes. For instance it gives a finite weight to $`m<0`$ even if $`w(\sigma )`$ is strictly positive. That said, the CLT may be attractive for narrow size distributions, and one might make a connection with the recent results of Evans et al  (even though the latter do not seem to rely on a specific shape of the distribution). We have not fully explored this avenue.

## C Spinodal criterion from exact free energy

In this appendix, we apply the spinodal criterion (53) to the exact free energy (41) and show that it can be expressed in a form identical to (58). This result has been given by a number of authors , but we include it here for the sake of completeness.

Choosing for $`𝝆`$ the vector $`\rho (\sigma )`$ (whose “components” are the values of $`\rho (\sigma )`$ for all $`\sigma `$), the spinodal criterion (53) applied to (41) becomes

$$𝑑\sigma ^{}\frac{\delta ^2f}{\delta \rho (\sigma )\delta \rho (\sigma ^{})}\delta \rho (\sigma ^{})=\underset{i,j}{}\frac{^2\stackrel{~}{f}}{\rho _i\rho _j}w_i(\sigma )\delta \rho _j+T\frac{\delta \rho (\sigma )}{\rho (\sigma )}=0$$
(C1)
The change of $`\mathrm{ln}\rho (\sigma )`$ along the instability direction (which is $`\delta \rho (\sigma )/\rho (\sigma )`$) is therefore a linear combination of the weight functions $`w_i(\sigma )`$. This means that the instability direction is within the family (42), consistent with our general discussion in Sec. III A. One can now rewrite (C1) in the form

$$\delta \rho (\sigma )=\beta \underset{i,j}{}w_i(\sigma )\rho (\sigma )\frac{^2\stackrel{~}{f}}{\rho _i\rho _j}\delta \rho _j$$
(C2)
and take the moment with the $`k`$-th weight function $`w_k(\sigma )`$ to get

$`\delta \rho _k+\beta {\displaystyle \underset{i,j}{}}\rho _{ki}{\displaystyle \frac{^2\stackrel{~}{f}}{\rho _i\rho _j}}\delta \rho _j=0.`$

As promised, this is identical to the spinodal condition (58) derived from the moment free energy, with the matrix multiplications written out explicitly.

## D Determinant form of critical point criterion

In this appendix, we give the form of the critical point criterion (60) that uses the spinodal determinant $`Y`$ from (58. At a critical point, the instability direction connects two neighboring points on the spinodal. The first order variation of the spinodal determinant $`Y`$ along this direction must therefore vanish. This gives

$$\delta Y=\underset{i}{}\frac{Y}{\rho _i}\delta \rho _i+\underset{ij}{}\frac{Y}{\rho _{ij}}\delta \rho _{ij}=0.$$
(D1)
Here the $`\delta \rho _{ij}`$ are the changes in the second order moment densities along the instability direction. We can express these as a function of the $`\delta \rho _i`$ via the change in the Lagrange multipliers $`\lambda _i`$: From the definition (43), a change in the $`\lambda _i`$ changes the first and second order moment densities according to

$$\frac{\rho _i}{\lambda _j}=\rho _{ij}=𝑴_{ij},\frac{\rho _{ij}}{\lambda _k}=\rho _{ijk}.$$
(D2)
Hence changes in second and first order moment densities are related by

$`\delta \rho _{kl}={\displaystyle \underset{i}{}}\rho _{kli}\delta \lambda _i={\displaystyle \underset{ij}{}}\rho _{kli}(𝑴^1)_{ij}\delta \rho _j`$

Inserting this into (D1) gives the determinant form of the critical point condition

$$\underset{i}{}\frac{Y}{\rho _i}\delta \rho _i+\underset{ij}{}\underset{kl}{}\frac{Y}{\rho _{kl}}\rho _{kli}(𝑴^1)_{ij}\delta \rho _j=0.$$
(D3)
Of course, if in an application of the moment free energy method one succeeds in calculating $`Y`$ explicitly as a function of the $`\rho _i`$ alone (rather than as a function of the $`\rho _i`$ and the $`\rho _{ij}`$), then one can determine critical points simply from the condition

$`\delta Y={\displaystyle \underset{i}{}}{\displaystyle \frac{dY}{d\rho _i}}\delta \rho _i=0.`$

We have written a total derivative sign here to indicate that the dependence of $`Y`$ on the $`\rho _{ij}`$ is accounted for implicitly.