{'fraction_non_alphanumeric': 0.04162544375254612, 'fraction_numerical': 0.010097189082232439, 'mean_word_length': 4.467149220489977, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Perturbative polydispersity: Phase equilibria of near-monodisperse systems

The conditions of multi-phase equilibrium are solved for generic polydisperse systems. The case of multiple polydispersity is treated, where several properties (e.g. size, charge, shape) simultaneously vary from one particle to another. By developing a perturbative expansion in the width of the distribution of constituent species, it is possible to calculate the effects of polydispersity alone, avoiding difficulties associated with the underlying many-body problem. Explicit formulae are derived in detail, for the partitioning of species at coexistence and for the shift of phase boundaries due to polydispersity. `Convective fractionation' is quantified, whereby one property (e.g. charge) is partitioned between phases due to a driving force on another. To demonstrate the ease of use and versatility of the formulae, they are applied to models of a chemically-polydisperse polymer blend, and of fluid-fluid coexistence in polydisperse colloid-polymer mixtures. In each case, the regime of coexistence is shown to be enlarged by polydispersity.


I. INTRODUCTION
Materials whose constituent elements are much larger than atoms or simple molecules are now ubiquitous. Such substances include polymer blends, colloidal suspensions and emulsions. In the manufacture of such large components, it is impossible to produce a pure strain of truly identical particles. Instead, each colloidal latex is minutely different in size from its colleagues, and the number of monomers on each long polymer molecule inevitably varies across some distribution. These systems are said to be 'polydisperse', though the word is applied only loosely to the polymeric case, since not every molecule is unique. Some intriguing experimental [1][2][3] and simulational [4][5][6][7][8] data have been collected from polydisperse systems. For instance, unusual textures have been found in polydisperse polymer blends [9], and coexisting phases of hard spheres have been shown to contain unequal proportions of the various particle sizes, in both simulation [6,7] and experiment [1,10]. It is important, then, to understand the statistical thermodynamics of polydisperse mixtures.
It is well known how to calculate the equilibrium state of a thermodynamic system [11]. At fixed temperature, an expression must be found for the Helmholtz free energy. It is then straightforward to determine any property of the system at equilibrium: its value is that which minimises the free energy. To establish the concentrations of the system's various constituents, the minimisation must be performed subject to constraints of global conservation. The constitution of a uniform system cannot vary, because the amount of each component is conserved, but the free energy can sometimes be lowered by partitioning into coexisting phases of differing compositions.
The non-trivial part of this procedure is to establish a formula for the Helmholtz free energy. In principle, the subsequent minimisation to derive phase equilibria is simple. For a polydisperse system the story is quite different. Here, the number of components is thermodynamically large. In general, this does not greatly complicate the (already difficult) many-body problem of formulating the free energy. However, since that free energy is now a function of infinitely many conserved variables, the minimisation procedure, though formally understood [12,13], becomes intractable.
In the past, many calculations have been performed to study the effects of polydispersity, using pedagogical model systems with a convenient form of Hamiltonian [14,15], species distribution [16][17][18], or approximate free energy [13,[19][20][21][22][23]. Also, approximate solutions can be found by making ad hoc guesses about the compositions of coexisting phases, and then minimising with respect to just a few variables such as the overall density of each phase [14,24]. Alternatively, the minimisation can be performed numerically, by replacing the continuous distribution of concentrations with an arbitrary, finite set of variables [25].
The above methods invariably involve some way of approximating the many-body problem inherent in the thermodynamics of the particular system. The present study differs in that we shall calculate only the effects of polydispersity. The analytical procedure was outlined recently [26][27][28], and is presented in more detail in the present article, and applied in new ways to solve a range of problems in the physics of polydispersity. The method is a perturbative one, whereby polydisperse phase equilibria are derived from a reference state that is monodisperse and also in multi-phase equilibrium. We make no assumptions as to the nature of the system, and require no special properties or approximations to ensure tractability of the statistical mechanics, since we shall not address the underlying many-body problem. Instead, we assume that the many-body properties of the monodisperse reference system are known (whether from experiment or theory), and derive some exact thermodynamic relations, applicable to most systems, for the changes in its properties due to the introduction of polydispersity.
The aim then is to treat a polydisperse system as a perturbation to a monodisperse reference state, using the variables that characterise the polydisperse particles (for instance, size and charge) as expansion parameters. There are two major obstacles to such a theory. Firstly, the reference state is singular. To calculate thermodynamic equilibria, a knowledge of the chemical potential of each species is required. In the reference system, the population is zero for all species but one. Unfortunately, the chemical potential of a non-populated species is (as a rule) negative infinity. The second obstacle is that a narrow distribution of species is not necessarily a smooth one, and therefore may not be conducive to small-variable expansion. To proceed, then, we must isolate the badly behaved functions for exact treatment, and expand only smoothly varying quantities. The distribution of species is not expanded (contrast [16]), so that any such distribution may be treated, including a set of delta peaks representing a finite number of components. The canonical ensemble is used throughout, as it is of greatest practical relevance. Thus, the overall mix of species (the 'parent distribution') is specified a priori, and the pressure and chemical potentials are derived quantities, used to calculate the resulting distribution adopted by each coexisting phase after partitioning.
Perturbative methods suggested previously [29,30] have been dogged by a proliferation of variables, leading to unwieldy expressions. That problem is avoided here by relating all properties of the polydisperse system to just two or three functions that parameterise a generic free energy. An alternative approach for perturbing about a monodisperse reference is to assume that all species but one are dilute, and to use their concentrations as small expansion parameters [13], in the manner of a virial series. In the present study, by expanding in the width, rather than height, of the distribution of polydisperse components, we are using a reference state that correctly treats the many-particle interactions of concentrated, phase-separated systems. For narrow distributions therefore, this series is expected to converge more quickly than a virial-like expansion. An exception is at critical points, where both methods fail.
The formalism and notation are set out in the next section, where a generic free energy is perturbatively expanded to low order, to derive the pressure and chemical potentials of a single polydisperse phase. These expressions are used in section III to analyse polydisperse phase equilibria. The normalised distributions, giving the relative amount of each species in each phase, are calculated in section III A. Their differences (due to fractionation) are found to obey a very simple law. For systems that are simultaneously polydisperse in two properties (e.g. size and charge), this law is shown to describe 'convective fractionation', whereby one property is partitioned between phases due to a driving force on the other. In section III B, the effect of polydispersity on the total number densities at coexistence is calculated. The resulting shift of a phase boundary along the density-axis of a phase diagram is shown to be proportional to the variance (standard deviation squared) of the parent distribution, in the limit of a narrow distribution. Explicit expressions are given for the cloud-and shadow-point densities, as well as more general binodals. All the derived formulae quantifying polydisperse phase equilibria are expressed in terms of a couple of basic parameters of the free energy. In case the free energy is not known for the polydisperse system in question, a method for evaluating the relevant parameters is given in section IV A for soft interaction potentials, and in section IV B for hard spheres and related systems. The method for soft potentials is extended to determine the lowest-order effect of polydispersity on correlation functions in appendix A. To demonstrate the ease of use of the formulae for phase equilibria, they are applied in section V to a Flory-Huggins model of chemically polydisperse polymers, and in section VI to fluid-fluid coexistence in colloid-polymer mixtures, where slight polydispersity is shown to widen the coexistence region.

II. PROPERTIES OF A SINGLE PHASE
Consider a polydisperse phase of overall number density ρ, whose constituent particles are distinguished by a number ν of properties (e.g. mass, charge, diameter, oblateness. . . ). Each particle is characterised by an ν-component vector ǫ, drawn from a narrow distribution p(ǫ) normalised over the whole of ǫ-space. Let each component of ǫ represent the (dimensionless) fractional deviation of each property from some reference value. For instance, a hard sphere of radius r in a polydisperse sample would be characterised by a scalar ǫ = (r − r 0 )/r 0 , in terms of the reference size r 0 . Since the distribution p(ǫ) is narrow, the components of ǫ are always small numbers.
To isolate the badly behaved part of the chemical potential, we shall write the phase's free energy density (which is a functional of the distribution of species densities ρp(ǫ)) in two parts This is an exact statement, as it serves only to define F ex as the excess free energy density, over and above that of a polydisperse ideal gas of densities ρp(ǫ). The ideal part is simply the sum, over the continuum of species, of the free energy density of an ideal gas of each species, The chemical potentials for the continuum of species are given by a functional derivative [12] of F with respect to the density of each species, which, from Eq. (1), splits into two parts, By functional differentiation of Eq. (2), the ideal part of the chemical potential is which is singular in the monodisperse limit p(ǫ) → δ(ǫ), because it derives from the entropy of mixing. The other part, µ ex (ǫ), remains well behaved, since it describes the physics of interactions which vary little from one species to the next. The obstacles to a small-variable expansion, discussed above, have now been isolated in the badly-behaved but well-defined function µ id (ǫ) in Eq. (5), and the narrow, but possibly rapidly varying distribution p(ǫ), also appearing in Eq. (5). These quantities will be treated exactly, while the effects of interactions, embodied in the excess free energy density F ex , are expanded in the small variable ǫ.

A. Expansion of the excess free energy
The excess free energy density in a phase at equilibrium is a function(al) F ex ([p(ǫ)], ρ) of the intensive variables p(ǫ) and ρ. The dependence of F ex on temperature and other fields is suppressed in the notation, and F ex is assumed to be absolutely minimised with respect to any non-conserved order parameters.
Any distribution p(ǫ) is uniquely defined by its moments { ǫ , ǫǫ , ǫǫǫ , . . .} which, for vectorial ǫ, are averages over p(ǫ) of outer products of ǫ. For instance, Since p(ǫ) is narrow, higher moments are increasingly small. So a controlled small-variable expansion of F ex is where F ex (ρ, 0, 0, . . .) is the free energy of a phase of monodisperse particles all characterised by ǫ = 0. We shall write this as F ex m (ρ), where the subscript m denotes a quantity in the monodisperse reference system. The coefficients A, B and C are vector and tensor functions of the overall density ρ. The notation O(ǫ 3 ) indicates terms of third order in ǫ and higher. Equation (6) is a generic expansion, not specific to any particular system. The use of a non-singular expansion, involving only integer powers of small quantities, has the status of a conjecture. An alternative motivation for the form of Eq. (6) is given in section IV.
B. Expansion of the excess chemical potential Following Eq. (3), Eq. (6) is differentiated, using the identity where δ ǫ m /δ[ρp(ǫ)] = (ǫ m − ǫ m )/ρ. In this notation, ∂F ex /∂ ǫ is a vector whose components are the derivatives of F ex with respect to each component of ǫ . Similarly, ∂F ex /∂ ǫǫ is a tensor. Thus we obtain a general expression, truncated at second order in ǫ, for the excess chemical potential µ ex (ǫ) of the species with property ǫ, in a phase characterised by {ρ, p(ǫ)}, where µ ex m (ρ) is the excess chemical potential, dF ex m /dρ of a monodisperse reference phase at density ρ. Equation (7) will be central to the calculations of phase equilibria in section III.

C. Expansion of the pressure
The pressure of a polydisperse phase with density ρ and distribution p(ǫ) can be found from which is the continuum analogue of the standard thermodynamic relation for the pressure of a mixture. Using Eqs. (1) and (4) it follows that where P id is the ideal pressure of the polydisperse mixture. It is easily confirmed from Eqs. (2) and (5) that this is the usual ideal gas pressure P id = ρ. Substituting the series expressions for the excess free energy and chemical potential (Eqs. (6) and (7)) into Eq. (8), and using the pressure of the monodisperse reference system at density ρ, P m (ρ) = ρ − F ex m (ρ) + ρ µ ex m (ρ), yields the expansion for the pressure of the polydisperse phase, III. PHASE EQUILIBRIA

A. Fractionation
At coexistence between two or more phases, the distribution p(ǫ) of properties is typically different for each phase. We shall assume that the overall distribution of all particles in all phases is known (since the particles might have been synthesised en masse, before phase separation was induced). This will be called the parent distribution p p (ǫ). Without loss of generality, the origin of ǫ will henceforth be set at the mean of this parent distribution, so that ǫ p ≡ 0. The aim of this section is to calculate the distribution p α (ǫ) in any given phase α, of M coexisting phases.
The number of particles in phase α is some fraction n α of the total number in the system. So conservation of material is expressed as M β n β = 1. In fact, a stronger criterion than this holds. The number of particles belonging to each species is conserved. That is, M β n β p β (ǫ) = p p (ǫ) (10) for all ǫ. This is a convenient point at which to introduce a special notation which will be useful later. If Q is any property of a phase then, for that phase, let ∆[Q] be its deviation from the mean of Q over all phases, weighted by the number of particles in each. Hence, in phase α, The notation [. . .] α β indicates the difference between the quantity evaluated in phase α and in phase β. Returning to the problem of coexistence, note that, in any two of the M coexisting phases, α and β say, the chemical potential is equal, i.e. µ α (ǫ) = µ β (ǫ) for all ǫ. Combining Eqs. (4) and (5), this implies which can be used to eliminate each p β (ǫ) in Eq. (10), yielding We shall now use Eq. (7) for the excess chemical potential, truncated at first order, µ ex (ǫ) = µ ex m (ρ) + ǫ · (A ′ − A/ρ) + ǫ · A/ρ + O(ǫ 2 ) and expand the exponential in Eq. (12) as So Eq. (12) can be written where (13) is integrated with respect to ǫ, then the term linear in ǫ disappears, due to the condition ǫ p ≡ 0. We are left with 1 = 1/S α + O(ǫ 2 ), or S α = 1 + O(ǫ 2 ), which can be substituted back into Eq. (13). Finally, integrating Eq. (11) with µ ex (ǫ) expanded to zeroth order yields which we also substitute into Eq. (13). The result is This is the result we wanted: an expression for the distribution of species in each phase. Not surprisingly, it is almost the same as the parent distribution (which remains unapproximated), since all species have similar properties (ǫ is small for each). The lowest-order correction to p p (ǫ) is expressed in terms of A, a (vectorial) parameter of the excess free energy (Eq. (6)), and ρ β , the density of each phase. In principle, the quantities on the right hand side of Eq. (14) are known. The densities of the phases ρ α will differ a little (by an amount δ α , say) from the binodals of the monodisperse reference system ρ m α , so that ρ α = ρ m α + δ α , and the numbers of particles n α will have similar small changes. However, in Eq. (14) one may substitute either the monodisperse or the polydisperse values of ρ and n, whichever are available from theory or experiment, since the small differences affect only higher-order terms.
A more elegant expression is obtained if we find the difference between the normalised distributions p(ǫ) of two of the M coexisting phases. From Eq. (14), which, together with Eq. (10), contains the same information as Eq. (14). The arrow(→) indicates the strict asymptotic limit as ǫ 2 p → 0. The parameter A/ρ in Eq. (15), is the lowest-order coefficient of ǫ in Eq. (7), so we can write A gradient in ǫ-space is denoted ∇ ǫ , and it is evaluated for phases α, β in Eq. (16) at the mean species ǫ = 0. Equation (16) says that, at coexistence, the distributions separate along the steepest gradient of µ ex α − µ ex β in ǫ-space. It is informative to take moments of Eq. (15), since these will tell us how the various properties of the particles are correlated within the phases. One might wish to examine averages of the individual components ǫ i of ǫ, such as ǫ 1 the mean value of some property such as the size of particles, or ǫ 1 ǫ 2 the cross-correlation between e.g. sizes and charges of the particles. The most general moment, from which this information can be extracted, is the mean of an outer product of m vectors ǫ. That is the mth rank tensor m whose components are all possible correlators of order m. From Eq. (15), the difference of such a moment between two of the coexisting phases is which involves the next moment of the parent distribution, and one summation (the scalar product) over the ν polydisperse properties. The scalar version of Eq. (18) was derived previously for substances that are polydisperse in a single property, at two-phase [26] and multi-phase [27] coexistence.
It has been noted [31] that for some distributions, odd moments ǫ 2k+1 p are not of order (2k+1) in the width of the distribution but, due to cancellation in the integral (Eq. (17)), are of higher order, so that the RHS of Eq. (18) is zero to O(ǫ m+1 ) and the unevaluated O(ǫ m+2 ) terms are the dominant order 1 . Specifically, this is the case for distributions which become symmetric in the narrow limit (e.g. Gaussian, Schultz). For strongly asymmetric distributions though, excepting unlikely cancellations, the mth moment is of mth order in the width. Equations (15) and (16) remain valid for all narrow distributions.
As the notation of Eq. (18) is rather abstract, let us consider a specific example: a system of particles with a range of radii r and charges q. At equilibrium the particles are partitioned into two coexisting phases denoted α and β, and we are interested in the difference between the average radii in these phases. The overall average radius throughout the system is called r p and the average charge q p . Applying Eq. (18) with m = 1, ǫ 1 = r/ r p − 1 and ǫ 2 = q/ q p − 1, and taking the first component gives the answer The first term says that, in the absence of charges, the amount of size fractionation is proportional to the overall size variance (the square of the standard deviation), with a coefficient a indicating the difference between each phase's affinity for large particles. With charge polydispersity, the second term indicates that size fractionation can occur even if a = 0 so that neither phase favours larger particles. If size and charge are cross-correlated in the parent, i.e. bigger particles tend to have bigger (or smaller) charges, then a 'convective' fractionation of sizes will be driven by b, the 'driving force' separating charges. These two contributions are additive for a system with a narrow distribution.

B. Shift of the phase boundaries
So far we have established the lowest order change (with width of the distribution) of the distribution of properties in each phase, relative to the original, 'parent' distribution (Eq. (14)). However, a knowledge of the normalised distributions is not a complete characterisation of the coexisting phases. We would also like to know the overall densities of the phases, ρ α . For a monodisperse system, the binodal densities are independent of the number of particles in each phase. This is not so when the system is polydisperse, since the size of one phase's share of the available particles influences the shape of its distribution of species and, hence, its thermodynamic properties. One coexistence of particular interest is the 'cloud point'. A phase is at its cloud point when its density lies on the 'cloud curve' in the phase diagram. It then coexists with an infinitesimal amount of another phase, whose density is on the corresponding 'shadow' curve. The cloud point is important from a practical point of view, since it marks the threshold at which a separation of phases is first observable, and the resulting morphology and fractionation can have useful applications. The cloud point also inspires theoretical interest since it is relatively easy to treat in a two-phase system, where the distribution of species in the majority phase is simply equal to the parent.
In the present analysis, as in the previous section, we shall calculate the conditions of general phase equilibria, not just the cloud and shadow curves. Having found the general solution, we shall return to the special case of the cloud point.
To establish the phase equilibria, we shall demand that the pressure and chemical potentials are equal in all phases, using the expressions for pressure and chemical potential derived in section II. Before proceeding further, let us take moments of Eq. (14), as this will lead to some substantial simplifications. We see that the second moment, ǫǫ , in any given phase is, to lowest order, equal to that in the parent. However, the mean, ǫ , in that phase is also of second order in the parent's width, since the mean of the parent vanishes by definition. In orders of the width of the parent, and So the terms involving ǫ ǫ in Eqs. (7) and (9) can be dropped, since they are of order σ 4 p . Let us clarify this point. The excess chemical potentials in Eq. (7) and the pressure in Eq. (9) were calculated consistently to second order in the deviations ǫ of the properties of a phase's constituent particles from some reference. In general, then, to evaluate these thermodynamic potentials for an arbitrary phase, with respect to an arbitrary reference, terms involving ǫ ǫ must be preserved. However, the choice ǫ p = 0, to fix the arbitrary reference, allows us to drop the terms in Eqs. (7) and 9 with coefficient C, for coexisting phases. This is because we have found that when several phases coexist, their mean properties differ very little from each other, so that ǫ α − ǫ p is second-order small.
To calculate the densities of the coexisting phases, we first demand that the pressure P (ρ, ǫ , ǫǫ ) is equal in any two phases α, β of the M coexisting phases. We shall write this as and use the expression for P (ρ, ǫ , ǫǫ ) in Eq. (9). Equation (20) is a condition on ρ α , the binodal densities which we want to find. As mentioned above, in the limit of a narrow parent, these densities are close to the monodisperse binodals ρ m α , differing by an amount δ α , so that ρ α = ρ m α + δ α . We shall determine the small shift δ α . In Eq. (9) we can Taylor-expand the monodisperse pressure P m (ρ α ) = P m (ρ m α ) + P ′ m (ρ m α ) δ α . It will transpire that δ α is of order σ 2 p , so that it is sufficient to truncate this Taylor expansion at first order in δ α . The condition of pressure balance (Eq. (20)) becomes We see that, due to the low order of expansion, the change in density δ affects only P m , the 'monodisperse contribution' to the pressure. Effectively, each of the polydisperse phases can be replaced by a monodisperse phase that is subject to an 'external field' ( ǫ · (ρA ′ − A) + ǫǫ : (ρB ′ − B)). In this interpretation, each (fictitious) monodisperse phase responds to the field with a density change δ, governed by the response function P ′ m . The same interpretation can be applied to the chemical potential balance below.
We can eliminate the first term because the monodisperse reference phases at densities ρ m α are in pressure balance, that is [P m (ρ m )] α β = 0. The density derivative of the pressure in the monodisperse system, P ′ m (ρ m ), can be expressed in terms of the phase's isothermal compressibility, Hence pressure balance dictates for any pair (α, β) of the M phases. Another constraint on the shifts in density δ α,β can be derived from the equality of all chemical potentials between any pair of phases, which leads to Eq. (11) for the excess parts. We substitute from Eq. (7) for µ ex (ǫ), and further substitute Expanding the exponential of small quantities and integrating gives Now, applying Eq. (21) and the condition [µ m (ρ m )] α β = 0 yields which, with Eq. (22), forms a pair of simultaneous equations in two unknowns, δ α and δ β . As an aside, it is interesting that, as well as fixing the M unknowns δ α , the 2(M − 1) constraints in Eqs. (22) and (23) impose extra conditions on the reference system for M > 2. The extra constraints are analogous to the Gibbs phase rule: a single component, whether monodisperse or near-monodisperse, can exist as a pair of coexisting phases under a range of conditions, whereas triple or higher coexistence requires a special choice of external fields [32]. It is intriguing that the extra conditions imposed on the monodisperse reference system constrain the values of A(ρ) and B(ρ) (defined in Eq. (6)), which, though defined in the monodisperse system, are not manifestly relevant to its phase equilibria.
It is straightforward to solve the linear equations (22 and 23) to yield the value of δ α , expressed as , where expressions for ǫ , ǫǫ can be substituted from Eqs. (19). However, the expression is inelegant, since there are M phases, of which β is in no way special. The shift δ α in the α-binodal should have a symmetric dependence on all other phases. Since β is arbitrary, we may symmetrize Eq. (24) by averaging it over all β, with respect to any weight factor z αβ : . The final answer for the shift δ of the density of one 2 of the M coexisting phases due to polydispersity is plus terms of order σ 3 p . As this is a lowest-order expression, the densities ρ appearing on the R.H.S. may be taken as either the monodisperse or the polydisperse binodal values, whichever is most convenient.
As yet, we do not know the values of the parameters A, B, A ′ and B ′ in Eq. (25). They will be discussed in Section IV. Nevertheless, the structure of Eq. (25) is informative. As stated above, δ is indeed of order σ 2 p . If, as is often the case [6,7,23,24,33,34], a phase diagram is drawn for some system, showing the densities of coexisting phases against the width σ p of some parent distribution, then Eq. (25) gives the leading-order Taylor expansion of the shape of the phase boundaries 3 . By definition, they meet the (σ p = 0)-axis at the monodisperse values ρ m , and we have shown that they do so perpendicularly, since there is no term linear in σ p . According to Eq. (25), the curvature of each phase boundary is proportional to the isothermal compressibility κ of the phase. Hence the small-variable expansion breaks down at a critical point, where the compressibility diverges. This is not unexpected, since the position of the critical point can itself be shifted by the presence of polydispersity, for instance to a different temperature. Criticality in the monodisperse system is a coexistence between two (or more) phases of equal densities. Under the same conditions, slightly-polydisperse phases, being non-critical, will either coexist at finitely-different densities, or not demix at all. In either case, the thermodynamic state is far from the monodisperse reference, so the divergence of Eq. (25) at criticality is correct.
Equation (25) exhibits a property typical of polydisperse systems. For a given parent distribution, the positions of the phase boundaries depend on the number of particles in each phase (since this affects the averages in the definition of ∆[. . .]). This contrasts with monodisperse systems, for which tie-lines can be drawn, connecting a set of points in the phase diagram's forbidden region to the same coexisting end-points. It is interesting that this signature of polydispersity [43] is evident at the lowest order in σ p . Hence, as soon as any shift in the phase boundaries due to polydispersity is observable in a system, that shift should vary as a tie-line is traversed.
As mentioned above, a coexisting state of particular interest is that between an infinitesimal amount of one phase (the 'shadow' phase) and a majority phase defined to be at its 'cloud' point. Such a state defines the extreme boundary of the coexistence region. The movement of the cloud point 4 , δ c is of great importance, since this tells an experimenter where, in the phase diagram, to first expect phase separation to occur. Let the portion of particles in the shadow phase be the infinitesimally small number n s = s. Then the fraction belonging to the majority, cloud phase is n c = 1 − s. For any property Q, ∆[Q] vanishes in the majority phase as s → 0, since the cloud phase defines the number-weighted mean. In the cloud phase, s denotes the difference between the values of Q in the cloud and shadow phases, Q c − Q s . In the shadow phase, Substituting these expressions into Eq. (25) gives the shift of the cloud point as The coexisting shadow phase has a density ρ m s + δ s where It is not always appropriate to express the position of a phase boundary in terms of the number density of particles, ρ = ρ m + δ. For instance, the phase diagram is often represented in terms of the volume fraction of particles, φ. This is the product of the number density and the volume of a particle. Since fractionation occurs, so that the mean volume of a particle differs from phase to phase, φ and ρ are not equivalent, and the coexistence region has a different shape in the two representations of the phase diagram. Quantities such as φ are easily derived from the above results as follows. For spherical particles, for instance, using the scalar ǫ to denote fractional deviations in particle radii: with the quantities on the RHS given in Eqs. (19) and (25).
In describing the phase diagram, one final quantity of interest is the width of the coexistence region, sometimes called the miscibility gap. In terms of ρ, this is the range of densities of the system as a whole, for which it separates into more than one phase. Say the miscibility gap is bounded by two phases, called α and β in the reference system. Quenching to one edge of this gap, the cloud point of phase α, will cause an infinitesimal amount of phase β to form. At the other extreme, a vanishing amount of phase α coexists with majority phase β. So the miscibility gap is given by application of Eq. (26) alone, to establish the two cloud points of phases α and β. The miscibility gap [ρ c ] α β is then given, in terms of the gap in the monodisperse reference system, by Note that, with the gap defined to be positive, [−1/ρ] α β is positive. In principal, the lowest-order change in the miscibility gap due to polydispersity, given in Eq. (29), may be positive or negative, depending on the system. We shall return to this topic in section V.

IV. EVALUATING THE PARAMETERS
We have determined the state of thermodynamic equilibrium for polydisperse phases, in terms of certain properties of the monodisperse reference phases. Some of these properties, the density ρ and compressibility κ, are well established for most common systems. Others, the parameters A(ρ), A ′ (ρ), B(ρ) and B ′ (ρ), require further discussion. These parameters appear in the excess free energy, Eq. (6), and we require their values in the monodisperse system (where ǫ = ǫǫ = 0) at some density ρ. One might already have an expression to hand for the free energy of a particular system, in the form of Eq. (6), in which case it is trivial to identify the appropriate parameters and substitute them into the various expressions above. The results for two commonly-studied systems are given in sections V and VI. Alternatively, the required numbers may be established by experiment, or may require derivation from first principles.
We shall now find some generic expressions for A(ρ), B(ρ) and, for completeness, also C(ρ), which will inform our physical interpretation of these quantities. Then, in section IV A, we shall apply those expressions, to construct a method for first-principles calculation of the parameters from thermodynamic perturbation theory. That theory is applicable whenever the Hamiltonian of the system can be differentiated with respect to the properties ǫ of the particles. Though many systems meet this criterion, one system of particular theoretical and pedagogical interest, a set of hard spheres, does not. Its Hamiltonian is discontinuous, being zero for all physical configurations and rising discontinuously to infinity if any spheres overlap. However, hard spheres belong to a special class of systems (scalable systems) for which the parameter A, which controls the emergence of fractionation (see Eqs. (14,15,18)) is calculable by an alternative method, detailed in section IV B.
To find general expressions for the quantities A, B and C, let us Taylor expand a generic excess free energy with respect to the properties ǫ i of each constituent particle separately, where i labels the N particles. We expand to second order in these N vectorial variables thus where the subscript m indicates evaluation in the system where ǫ i = 0 for all i. In that system the particles are all the same, so the derivative ∂F ex /∂ǫ i takes the same vectorial value for any particle i. Hence we may evaluate it for particle number 1, without loss of generality. Similarly, the second derivate ∂ 2 F ex /∂ǫ i ∂ǫ j takes just two tensorial values, depending on whether i and j label the same particle (number 1, say) or different particles (1 and 2, say). Hence, dividing the above equation by the system volume V gives the excess free energy density Note that, although the final term contains a factor N , it remains intensive, because ∂ 2 F ex /∂ǫ 1 ∂ǫ 2 ∼ 1/V since particles 1 and 2 interact less as the system grows. Comparing Eq. (30) with Eq. (6), we can identify expressions for the coefficients A, B and C, which we also express in terms of derivatives of Eqs. (6) and (7), with all formulae evaluated in the monodisperse reference system. To clarify the meanings of the various derivatives, let us read Eq. (31a) from left to right. It states that the parameter A, which is a property of a monodisperse phase, is the rate of change of excess free energy density as the mean property ǫ of all particles in the phase is changed. The next equality asserts that A is the overall number density ρ times the rate at which the extensive excess free energy changes when the property ǫ of particle number 1 alone is changed. Finally from Eq. (31a), this is the same as ρ times the variation in excess chemical potential for different species. (Though different species are absent from the monodisperse system, their excess chemical potentials remain well defined.) In sections IV A and IV B two methods are developed, using the various relations in Eqs. (31) to ascertain the values of A, B and C in any given system.

A. Thermodynamic perturbation theory
In this section, polydispersity is treated as a perturbation to the Hamiltonian of a system. The Hamiltonian of a polydisperse system (neglecting the trivial kinetic part) is H(Γ), a function of the set Γ of particle coordinates (both positional and internal) which defines a configuration of the system. Given that a monodisperse reference system in the same configuration Γ has Hamiltonian H m (Γ), we define the perturbationĤ(Γ) by By Boltzmann and Gibbs, the Helmholtz free energy of the polydisperse system is (with units k B T ≡ 1) and that of the monodisperse system is In these two equations we may set H and H m to zero to find expressions for the ideal part, which we subtract, leaving the excess part of the free energies, By applying Eq. (32), we see that the difference between these expressions is where the argument of the logarithm is the thermal average of exp(−Ĥ(Γ)) in the monodisperse system. Let us denote the thermal average of a stochastic variable ζ by ζ . This is a Boltzmann-weighted average over configurations, denoted by double brackets to distinguish it from averages over the distribution of species p(ǫ), for which we have been using single brackets . . . . Hence, the excess free energy of a polydisperse system, in terms of a perturbationĤ to a monodisperse system, is where the subscript m indicates that the thermal average is performed in the monodisperse system. Equation (33) has the familiar form of thermodynamic perturbation theory but, for polydispersity, it applies only to the excess part of the free energy. The ideal parts of the mono-and polydisperse free energies differ by a non-perturbative amount. Note that we are not making any approximation of the interactions, unlike many perturbation theories which use an ideal or harmonic reference state. Hence, as with our previous analysis, Eq. (33) treats the effects arising purely due to polydispersity, in a fully interacting system.
We now apply the second equality in Eq. (31a) to the perturbative expression for the excess free energy, Eq. (33), to find A by varying the properties of a single particle in an otherwise monodisperse system. The first term in Eq. (33), F ex m is a constant and therefore does not contribute. We find where the 'hat' (ˆ) has been dropped from the Hamiltonian H since, by definition, the perturbationĤ is the only part that varies with ǫ 1 . Similarly, for Eqs. (31b) and (31c), In principle, the derivatives of the Hamiltonian appearing in Eqs. (34a), (34b) and (34c) are known from the microscopic physics of a given system. For instance, in a system of particles with central, symmetric, pairwise additive interactions, whose inter-particle potential is U (r, ǫ i , ǫ j ), the Hamiltonian is In terms of the matrices of derivatives of the potential, Here, g m (r) = ρ(0)ρ(r) m /ρ 2 is the radial distribution function in the monodisperse reference phase. To obtain Eq. (37a), we have used the fact that the density of particle centres in the system at any instant is ρ(r) = N i=1 δ (3) (r− r i ). Equations (34b, 34c) can similarly be evaluated for central, symmetric, pairwise additive interactions, giving m (r, r ′ , r ′′ ) ≡ ρ(0)ρ(r)ρ(r ′ )ρ(r ′′ ) m /ρ 4 . Notice that Eqs. (37a, b, c) require a knowledge of spatial correlations in the monodisperse system only. In the appendix, it is shown how general thermal averages are perturbed by polydispersity, for such systems with soft, pairwise interactions.
Equations (34) and (37) allow the parameters affecting phase equilibria (to lowest order in the width of the distribution) to be calculated for any system with soft (i.e. differentiable) interactions. We shall put these expressions to use on some real systems in sections V and VI.

B. Scalable systems: The special case of hard spheres
In this section we shall consider a system of prevalent theoretical interest, polydisperse hard spheres. Hard spheres interact via a potential which is zero, except for configurations where the spheres overlap, which are forbidden and hence have infinite potential energy. Hard-sphere systems attract interest because they exemplify substances with short-range repulsive interactions, while lacking a characteristic energy scale. This leads to temperature-independent behaviour which is purely entropy-driven, and engenders simplicity due to the small number of tunable parameters. Unfortunately, hard spheres are an exceptional case to which the thermodynamic perturbation theory in section IV A cannot be applied, since the internal energy is zero and the interaction potential is discontinuous. Hence the values of the parameters A, B and C cannot be found by that method, though their formulations in Eqs. (31) still hold.
In this section we shall exploit a special property of the hard-sphere system to establish the exact value of the parameter A (previously approximated [35]). In fact the method can be applied to any system of particles whose interactions are 'scalable' -a term which will be elucidated below. With a knowledge of A, the partitioning of hard spheres between slightly-polydisperse phases is fully determined in Eq. (14), (15) or (18). Unfortunately it is not clear how such a method might be used to determine B, the other parameter on which the phase boundaries depend. In section VI its value is taken from an approximate expression for the free energy of polydisperse hard spheres, of which there are many examples in the literature [36,14,23].
To calculate A, the coefficient of ǫ in the free energy expansion (Eq. (6)), for polydisperse hard spheres, we use the first equality of Eq. (31a). Let the radius r of each sphere be measured relative to a reference length r 0 thus: r = (1 + ǫ) r 0 . As this is the only property that varies from sphere to sphere, the value ǫ which characterises the particles is a scalar, so Eq. (31a) reduces to a scalar equation. That equation states that A is the rate at which the free energy density changes when the first moment ǫ of size deviations is increased, while the other moments remain stationary 5 . When this hypothetical change is made, the particles all grow by the same amount, so the system remains monodisperse. In terms of the unique radius r of the monodisperse particles, the derivative is which will be evaluated at r = r 0 . One can imagine making this change in two stages: firstly the whole system is scaled up, so that the particles, and the space between them, and the volume of the system all increase, while the concentration φ (the fraction of space occupied by particles) is held constant. Then the system is compressed back to its original volume while the particles remain at their new large size, so the concentration increases. The resulting change in the extensive excess free energy (F ex = V F ex ) is given by The last term here is due to the reversible work done against excess pressure when the system is compressed, with (∂V /∂r) φ,N = 3V /r being a conversion factor between length and volume. So we have The first term on the right hand side is trivial to calculate because there exist no length scales (such as the range of an interaction) that remain fixed as the particles grow. Hence the growth is simply an overall scaling [17]. This is what was meant by the term "scalable system" above. By writing the free energy in terms of the usual configurational integral over particle positions measured in units of the thermal de Broglie wavelength, it is easy to establish that this term vanishes, since the excess part does not vary with the overall scale factor. Hence, for hard spheres [37], where, as usual, k B T has been set to unity. Since, at equilibrium, P is a constant for all phases, the coefficient in Eqs. (15,18) may be written [A/ρ] α β = 3P [ρ −1 ] α β for size-polydispersity in scalable systems. Though it is not apparent that an analogous method exists for calculating B, the parameter C is calculable from the first equality of Eq. (31c), by taking the second derivative of excess free energy density with respect to particle size in a monodisperse scalable system. We shall not perform that calculation here.
For monodisperse hard spheres, there is a transition [38,48] from a crystal of volume fraction 0.545 to a fluid phase of volume fraction 0.494 which, according to the Carnahan-Starling equation of state [39], exists at a pressure 6.17 ± 0.02, in units scaled by k B T and by the volume of a particle. Hence, from Eq. (38), the coefficient [A/ρ] crystal fluid , which appears in Eq. (18), evaluates to 3.51 ± 0.04. This is consistent with simulation results in which the average sizes of polydisperse hard spheres were measured in coexisting fluid and crystal phases, and the difference plotted against the variance of the overall distribution [40]. The gradient was found to be 3.55 ± 0.01, and changed by only a few percent up to polydispersities of 6% and more.

V. EXAMPLE: CHEMICALLY POLYDISPERSE POLYMERS
It has been necessary to plough through a fair amount of mathematics in order to arrive at some simple and practical formulae. The convenience of this approach is that the lengthy derivations have been performed once and for all. They will not require repetition for each application. Lest the reader has lost sight of the wood for the trees, let us demonstrate how easily the formulae may be applied to a model of polymeric phase equilibria.
We consider a Flory-Huggins-style [41] description of a polymer blend. The reference system is a binary mixture of polymers (species 'A' and 'B') of equal size L, but two different chemical types. The chemical nature of each species is parameterised by a number in the interval (−1, 1), which may be interpreted as the hydrophobicity of the molecule. We make species B chemically 'neutral' (hydrophobicity zero), while the value for species A is a. If the concentration of A-polymers is φ then that of B-polymers is 1 − φ, and the mean-field free energy density is given by where φ 1 = aφ is the 'hydrophobic concentration'. The Flory-Huggins interaction parameter χ determines the strength of attraction between chemically similar species. This model system separates into coexisting A-rich and A-poor phases for χ ≡ χa 2 > 2, the concentrations of which are easily found e.g. from a double-tangent construction [11] on the free energy (Eq. (39)). The resulting phase diagram in the (φ, χ)-plane is shown by the solid line in Fig. 1.
If component A is now made chemically polydisperse, by varying slightly the constituent monomers on each polymer molecule, while component B remains monodisperse and chemically neutral, acting only as a polymeric solvent, then the overall hydrophobic concentration of polymer A is the mean of a distribution: φ 1 = φ a p(a) da = φ a and the free energy density becomes The model free energy of Eq. (40) was studied previously to demonstrate a sophisticated approximation scheme for systematically reducing the dimensionality of polydisperse phase equilibria problems [43]. Accurate cloud (o) and shadow (*) points were calculated by that method, and are reproduced in Fig. 1 (where χ ≡ χa 2 0 ) for a Gaussian parent with a p ≡ a 0 = 0.5 and polydispersity σ = 8% [42].
For comparison, our lowest-order perturbative formulae for the cloud and shadow points (Eqs. (26,27)) can be evaluated with ease. Note that the last term of Eq. (40) has the form of an ideal free energy density if we identify φ with ρ in Eq. (2) (up to an irrelevant term linear in ρ). Hence, by writing a = a 0 (1 + ǫ) to expand about the monodisperse reference system with a = a 0 = 0.5, Eq. (40) can be cast in the form of Eq. (6), with all of which are scalar quantities since only one property (hydrophobicity) is polydisperse. The cloud and shadow points (φ c = φ m c + δ c , φ s = φ m s + δ s ) are obtained by substituting these values into Eqs. (26,27) together with the isothermal compressibility derived from Eq. (39) where φ m c , φ m s are the respective points on the monodisperse binodal (solid line in Fig. 1). Eqs. (41,42) are plotted in Fig. 1 as dashed and dash-dotted lines respectively. As expected, the perturbative expansion scheme breaks down near the critical point, where the compressibility diverges and the separate reference phases disappear. This invalid regime extends over only a small range of χ, and the phase boundaries rapidly converge, with increasing χ, to a very accurate answer.
Notice that the miscibility gap is broadened by polydispersity. Equation (29) shows that the polydisperse perturbation to the gap must be positive whenever B = B ′ = 0, i.e. when there is no explicit dependence of the excess free energy on the variance of the distribution of species. We also have formulae prescribing the chemical properties of the two phases. By definition, the distribution of polymers in the cloud phase is just the same as the parent distribution that was input to the system. The shadow phase, on the other hand, has chosen its own preferred mix. The mean hydrophobicity of type-A polymers in this phase is given by Eq. (19a) (or equivalently Eq. (14)) as a s = a 0 (1 + ǫ s ) = a 0 (1 + 2σ 2 χ(φ s − φ c ) + O(σ 3 )), which is also in good agreement with exact numerical data.

VI. EXAMPLE: STERICALLY-STABILISED COLLOID-POLYMER MIXTURES
To further illustrate the utility of the formulae derived in this article, let us make a case study of a familiar complex fluid: a colloid-polymer mixture [44]. The phase behaviour of the monodisperse system is first briefly reviewed. A theoretical understanding of its behaviour is well-established [45]. However, the effects of colloidal polydispersity on the system have not previously been calculated, despite the fact that the results of colloidal synthesis are inevitably slightly polydisperse. The effects of polymeric polydispersity on the system, which have been modelled previously [46], will not be analysed.

A. Established monodisperse behaviour
In this system, small balls (typically of diameter 10 nm to 1 µm) of PMMA ('Perspex') [44] or silica [47] are suspended in an organic liquid. To avoid the balls (or 'colloidal latices') clumping due to van der Waals attraction, they are each given a short brush-like coat of 'steric' polymer chains, making the latices behave as hard spheres, with no long-range interactions but a very short-range strong repulsion. Hard spheres of this kind are known [48] to have a simple phase diagram, and carry no latent heat. The only controlling parameter is the fraction of space φ occupied by the spheres. A single phase transition exists, between an amorphous fluid state and an FCC crystal, which coexist at φ ≈ 0.494 and φ ≈ 0.545 respectively [38]. The phase behaviour is enriched by the addition of free polymer chains to the mixture. Each of these chains, being very long, bunches itself into the form of a self-avoiding random walk. Such a 'ball of string' cannot penetrate a colloidal latex, and so interacts with the colloid particles as if it too were a hard sphere. On the other hand, when two randomly coiled polymer molecules meet, they can interpenetrate and, with a suitable choice of suspending solvent, can be made almost entirely non-interacting 6 , like an ideal gas [44,45]. The ideal gas pressure that they exert on the surface of an isolated colloidal hard sphere is isotropic. However, if a second hard sphere is so close to the first that the non-penetrating polymer coils cannot fit between them, then the polymeric pressure is not felt on the adjacent parts of the colloids, so they are pushed together.
Asakura and Oosawa [50] showed that this mechanism could be described in terms of an effective attraction between colloidal spheres, known as the depletion interaction. Thus, rather than considering a mixture of two components (colloid and polymer) in the solvent, one can imagine a single-component system, consisting of spheres with a hard core repulsion and a longer-range pairwise 7 attraction. This effective attractive potential U (r) between two latices of radius a and centre-to-centre separation r is simply the product of the (ideal) polymeric osmotic pressure Π p and the volume of space from which the centres of polymer coils are excluded by the hard spheres [48]. At r = 2a (contact of the two hard spheres) the potential rises discontinuously to infinity, and for r > 2(a + r g ), where r g is the radius of gyration of a polymer chain, there is no interaction, since the diameter of a polymer coil (2r g ) sets the range of the effective potential.
The phase diagram of this system in the (φ, Π p )-plane has been calculated for the monodisperse case [45], using results from the Percus-Yevick integral equation theory for hard spheres [52]. It was found (in agreement with experiments [44]) to resemble the phase diagram of simple atomic substances such as argon, with Π −1 p playing the rôle of temperature. As well as the colloidal crystal, there are low-and high-concentration amorphous phases, dubbed 'colloidal gas and liquid', which meet at a critical point. We shall study the 'fluid-fluid' coexistence of these amorphous phases, thus avoiding problems of non-ergodicity arising in the polydisperse crystal [51].

B. Calculating the polydisperse phase equilibria
The strength of the treatment of polydispersity in this article is that it uses, as input, properties of the monodisperse system. The phase behaviour of monodisperse colloid-polymer mixtures is well understood, so application of the new formulae for polydisperse phase equilibria will be straightforward. In addition to the monodisperse phase diagram, we require the functions A(ρ) and B(ρ), defined by Eq. (6), which are scalars in this case since the particle radii alone are polydisperse.
Unlike the previous example in section V, we are not given an expression for the free energy of the polydisperse system. According to Eq. (31a), A(ρ) can be found from the monodisperse free energy, for which we use the expression given in Ref. [45]. It is simply the derivative of the free energy with respect to the radius of monodisperse colloidal particles (scaled by the radius), at fixed number density. Hence an expression for A(ρ) is straightforwardly found, so long as care is taken not to overlook the dependence of concentration on particle size, φ = 4 3 πa 3 ρ, and to keep the polymer size r g fixed while taking the derivative.
It is not so trivial to evaluate B(ρ) without a prescribed polydisperse free energy. Firstly, we must specify the difference between the Hamiltonian of the system with slightly polydisperse colloid, and that of a monodisperse reference system. Let us write the Hamiltonian of the colloid-polymer mixture with polydisperse colloid as a sum of two parts: (43) the Hamiltonian H HS of the hard-sphere colloid, plus effective interactions H eff dep due to the depletion of polymer. Following the analysis [45] of the monodisperse case, we assume that the resulting free energy is also separable into two parts: where F HS is approximately the free energy of polydisperse hard spheres in the absence of polymer. This is consistent with the Weeks-Chandler-Andersen (WCA) approximation [53] which states that the spatial distribution of repulsive particles (in this case hard spheres) is not significantly altered by the introduction of additional interactions which are purely attractive (here contained in H eff dep ). Accordingly, the polymeric contribution to the free energy F dep will be approximated by averaging the depletion Hamiltonian H eff dep over the pure hard sphere distribution. Note that the structure of Eq. (44) implies that the desired function B(ρ) also splits into two terms, The depletion interaction is soft, so the thermodynamic perturbation theory of section IV A is suitable for calculating its contribution in B dep . The hard sphere potential, on the other hand, is non-differentiable and therefore inappropriate for description in terms of an energetic perturbation, as noted above. Other methods will be used to deal with the hard sphere part B HS . Let us tackle the polymeric (depletion) contribution first. As stated above, the effective Hamiltonian (i.e. the polymeric free energy for a fixed configuration of colloidal spheres) H eff dep is simply the product of the ideal polymeric osmotic pressure and the volume of space from which centres of polymer coils are excluded. Let us define the polymercolloid size ratio as ξ ≡ r g /a, and v s ≡ 4 3 πa 3 to be the volume of a reference sphere with the mean radius. An isolated colloidal sphere of radius a(1 + ǫ) excludes polymer coils from a volume v s (1 + ξ + ǫ) 3 . Hence, even for well-separated colloids beyond the interaction range, there is a bulk contribution to the depletion Hamiltonian, where the polymeric osmotic pressure has been re-scaled:Π p ≡ v s ξ 3 Π p , so that the (φ,Π p )-plane will correspond to phase diagrams presented in Ref. [45] for colloid-polymer mixtures. For monodisperse colloid, H bulk is usually ignored, since it is a constant, independent of volume fraction φ. For a polydisperse system, on the other hand, it cannot be neglected, as it contains ǫ-dependence. From Eq. (34b), its contribution to B(ρ) is When two colloidal latices, of radii (1 + ǫ 1 )a and (1 + ǫ 2 )a, have a centre-centre separation (2 + ǫ 1 + ǫ 2 )a < r < (2 + 2ξ + ǫ 1 + ǫ 2 )a, their individual regions of polymeric exclusion overlap, so that the total volume available to polymer increases, giving rise to the depletion attraction. It can be shown by some elementary geometry that the overlap volume of the two spheres of depletion; radii r 1 = (1 + ξ + ǫ 1 )a and r 2 = (1 + ξ + ǫ 2 )a; is The resulting interaction potential U (r, ǫ 1 , ǫ 2 ) = −Π p V overlap gives a pairwise contribution to the effective Hamiltonian, as in Eq. (35). Its derivatives, as defined in Eqs. (36), are which can be substituted into Eq. (37b) for B dep . Additionally, that equation requires the pair and three-point distribution functions g m (r), g m (r, r ′ ) for the monodisperse system. At the level of the WCA approximation, these are replaced by the pure hard-sphere (HS) distribution functions, g ≈ g HS , for which we use the Percus-Yevick expression [54]. As a simplifying assumption, let us also use the Kirkwood superposition approximation [55] for the three-body correlations, It remains only to find the hard-sphere contribution to B. Unfortunately, the methods developed here do not enable a first-principles derivation of this quantity, so we must look to other analyses of polydisperse hard spheres for an evaluation of B HS . In particular, we refer to the free energy expression for the polydisperse hard-sphere fluid due to Boublik, Mansoori, Carnahan, Starling and Lealand [36] (BMCSL), which is known to reduce to the Carnahan-Starling free energy [39] in the monodisperse limit. The application of this expression perhaps requires some clarification. The use of results from other studies of polydispersity by no means depreciates the present analysis. In cases where a polydisperse free energy is already known, one is spared the application of first-principles methods such as section IV. However, the results of section III (or equivalent methods) are still required, in order to derive (to lowest-order) the phase equilibria from that free energy. Indeed, the gap between a knowledge of the free energy and a knowledge of the phase equilibria is evidenced by the BMCSL free energy, whose regime of thermodynamic instability was only recently established [20].
To extract the required function, the BMCSL free energy is cast in the form of Eq. (6) by expanding to order ǫ 2 , at a density ρ = φ 0 /v s , giving The zeroth-order term is the excess free energy of a monodisperse hard-sphere fluid, and recovers the Carnahan-Starling equation of state. By comparison with Eq. (6), the coefficient of ǫ 2 is the required function v s B HS . For definiteness, let us consider a polymer-colloid size ratio ξ = 0.4. For a given polymeric osmotic pressurê Π p , the coexisting concentrations φ of monodisperse colloidal fluid phases are given in Ref. [45]. These values for the concentrations of coexisting monodisperse gas (φ g m ) and liquid (φ l m ) phases are substituted into the expressions calculated above to obtain values of A and B for the coexisting phases. The monodisperse gas-liquid phase boundaries φ g m , φ l m are shown as a solid line in Fig. 2 from the critical point atΠ p ≈ 0.41 to the triple point atΠ p ≈ 0.54. The calculated values of A and B and their derivatives are used in Eqs. (26,27) to find the change in each phase's cloud-and shadow-point densities due to polydispersity, which are translated into concentrations via Eq. (28). These shifts in concentration are added to the monodisperse phase boundary in Fig. 2 to give the resulting cloud and shadow points φ c,s = φ c,s m + δφ c,s for polydispersity σ = 8%, shown as dashed and dash-dotted curves respectively. We see from the figure that the expansion has remained well controlled (has not blown up) even for values ofΠ p very close to criticality of the reference system. The theory predicts a widening of the coexistence region (the miscibility gap) so that liquid will condense from a polydisperse gas of lower concentration than in the monodisperse case. We also observe that when the gas phase exists in an infinitesimal quantity (at the shadow point), it is more concentrated than at its cloud point. The same applies to the liquid phases. Similar results are also found for different size ratios ξ (studied for 0.3 < ξ < 1). The above calculation of B(ρ) for the colloid-polymer mixture involves a degree of inexactitude, since the radial distribution function is approximated by the Percus-Yevick hard-sphere distribution function, while three-point correlations are even more crudely estimated. To assess the importance of precision in the correlations, the cloud and shadow points were re-calculated using a Heaviside step function g(r) = Θ(r − 2a) in place of the Percus-Yevick radial distribution. The polydispersity-induced shift in the high-density cloud point changes sign when the correlations are neglected, indicating, not surprisingly, that the correct correlations are important in determining the high-density phase boundary. Though the positions of the other phase boundaries are also affected by the change in g(r), certain qualitative features remain unaltered, and are therefore expected to be reliable predictions of our approximate model. In particular, each shadow point remains more concentrated than the equivalent cloud point and, far from the critical point, the gases are shifted to lower concentrations by polydispersity.
According to Eqs. (15) and (18), the difference between the normalised populations of the coexisting phases is proportional to −[A/ρ] l g and the fractional difference in mean particle size is [ ǫ ] l g = −[A/ρ] l g σ 2 p where σ p is the overall polydispersity. So the coefficient −[A/ρ] l g determines the amount of partitioning or fractionation. It is plotted for the gas-liquid phase boundary (fromΠ p = 0.41 to 0.54) in Fig. 3. Note that it is positive, indicating that, on average, the liquid phase contains larger particles than the gas. The hard-sphere contribution is negative, but is outweighed by the depletion part.

VII. SUMMARY
The results derived here are exact in the limit of low polydispersity, and can therefore provide accurate information on a great many systems, whose manufacture attempts to approximate the ideal of monodisperse constituents. Additionally, the results will be of pragmatic use even in systems with a wide scatter of particle properties, in providing qualitative predictions and estimates for the effects of polydispersity. In any case, minimal effort is required to substitute values into the formulae.
Equations (14) and, equivalently, (15) show that, in the narrow limit, the difference (due to fractionation) between normalised distributions in coexisting phases does not depend on the volumes of the phases, and is in fact remarkably simple, requiring only one system-dependent parameter per polydisperse property. From Eq. (18), the mth moment about the centre of the distribution differs, between coexisting phases, by an amount proportional to the (m + 1)th moment of the parent 8 . For hard spheres, the constant of proportionality has been determined (Eq. (38)), and agrees with data from simulation [40]. The calculation is indistinguishable from the data at 2% polydispersity, and deviates by only 5% at 4% polydispersity.
Additionally, at the nematic cloud point of the polydisperse Zwanzig model of hard rods [33], the form of Eq. (18) has been shown [56] to hold to within 5% up to the remarkably large polydispersity of 50%. In section III A, it was shown that this simple rule can lead to 'convective fractionation' for multiply-polydisperse systems, whereby one property is partitioned between phases due to a driving force on another.
The expressions found for the shift in the phase boundaries due to polydispersity, at a general coexistence (Eq. (25)) and at the cloud (Eq. (26)) and shadow points (Eq. (27)), show that the shift in the density of a phase is proportional to its isothermal compressibility and to the square of the polydispersity. This explains why these phase boundaries are parabolic near σ = 0 when plotted as a function of polydispersity σ in a range of models; see e.g. [6,7,23,24,33,34]. The other parameters on which these shifts depend can be extracted directly from the form of the polydisperse free energy if it is known or, if not, from the Hamiltonian via the thermodynamic perturbation theory in section IV A. In the case of pairwise interactions, the relevant parameters are obtained from a knowledge of only two-and threebody correlations in the reference system, even when larger clusters are responsible for the existence of the phase transition in question. This illustrates a strength of the small-polydispersity expansion: that the underlying reference system takes care of complicated many-body interactions, so they do not need to be calculated in the analysis of the polydispersity.
The application of the methods to the fluid-fluid coexistence of colloid-polymer mixtures in section VI revealed that colloidal polydispersity leads to a widening of the coexistence region, and favours larger particles living in the liquid phase. Similar methods are used in appendix A to find how polydispersity alters correlation functions.
Other methods exist for calculating polydisperse phase equilibria. Whatever the theoretical formalism by which one chooses to analyse a polydisperse system, the results, if correct, will tend in the limit of low polydispersity to Eq. (15) for the degree of fractionation, and Eq. (25) for the movement of the binodals. In polydisperse systems, correlation functions are easier to calculate than phase boundaries, as there is no partitioning of the sample to consider. The thermal average of a stochastic quantity ζ in a single polydisperse phase is given by the usual perturbative expression with respect to a monodisperse reference system We shall require the second-order expansion of this expression, in the small quantityĤ, which must therefore be small (in units of k B T ), Let us restrict the discussion to a system of particles interacting via a pairwise-additive, symmetric, isotropic, central potential U (r, ǫ 1 , ǫ 2 ), for which the Hamiltonian is given in Eq. (35). With the Hamiltonian Taylor-expanded to second order in ǫ, the perturbation can be written in terms of the derivatives of the interaction potential defined in Eqs. (36). Substituting this expression into Eq. (A2), we keep terms to second order in ǫ, and measure the deviations ǫ with respect to the mean properties of particles in the single-phase system, so that ǫ = 0. This yields ζ = ζ m − ǫǫ : where Ψ(r) ≡ U 11 (r) − 1 2 U 1 (r)U 1 (r).