# Fractionation of polydisperse systems: multi-phase coexistence

## Abstract

The width of the distribution of species in a polydisperse system is employed in a small-variable expansion, to obtain a well-controlled and compact scheme by which to calculate phase equilibria in multi-phase systems. General and universal relations are derived, which determine the partitioning of the fluid components among the phases. The analysis applies to mixtures of arbitrarily many slightly-polydisperse components. An explicit solution is approximated for hard spheres.

It is vital to gain an understanding of polydispersity, due to its ubiquity in both synthetic and biological complex fluids. A polydisperse substance is a mixture of infinitely many components, and can, in general, separate into arbitrary numbers of coexisting phases. These properties typically engender great mathematical complexity, and have been a stumbling block to the concise formulation of polydisperse thermodynamics. Experimental  and simulational  studies of polydisperse polymeric fluids and colloidal suspensions have catalogued diverse behavior and intricate phase diagrams. Until recently, theoretical treatments of polydispersity relied on uncontrolled approximations  or idealised models , while generic schemes and fundamental understanding remained elusive.

The phase behavior of pure (i.e. monodisperse) systems is (in principle at least) relatively straightforward to analyse. The standard method, formulated last century , involves integrating (by various approximate methods) the Boltzmann factor over all configurations to construct the Helmholtz free energy as a function of temperature, density and volume. From this, the densities of coexisting phases can be calculated and the phase diagram deduced. One source of difficulty in analysing the phase equilibria of polydisperse systems is that the density alone does not fully characterize a phase. Instead, we wish to know the entire composition of each phase. That requires the evaluation of an infinite set of variables.

Two systematic schemes were developed recently to solve the polydisperse phase equilibria problem. The powerful ‘annealed moments’ method of Sollich and Cates  and Warren  applies to a large subset of model systems, and will not be further discussed here. The second scheme, which applies to real systems, was developed by the present author . It uses the width of the distribution of species as a small expansion parameter, and is therefore valid for slightly polydisperse systems. In other words, the scheme is applicable whenever the polydisperse property (e.g. the size or charge of a particle) varies only a little throughout the system. (N.B. Here, ‘particle’ is used to denote the polydisperse fluid elements, which could be polymer molecules, colloidal latices, etc.) The method was used  to find the complete distributions of species in two coexisting phases of any slightly polydisperse system, resulting in a universal law of fractionation. In this paper, the method is applied to coexistence between arbitrary numbers of phases — a situation of importance to many polydisperse substances .

A slightly polydisperse system (one with a narrow distribution of species) is in principle very different from a truly pure one, which has no mixing entropy and whose distribution is a Dirac delta function. Nevertheless, one would expect the physical properties of the two systems to be very similar. That similarity motivates this study, since the pure system is vastly simpler to analyse than its polydisperse counterpart. To exploit that simplicity, a formalism is required which treats mono- and poly-disperse systems on an equal footing. Such a formalism is now derived (extending the method for two-phase equilibrium in Ref. ).

Let us first define a number $`\epsilon _i`$ to characterize each of the $`N`$ particles in the system (with $`i=1,\mathrm{},N`$) . For size polydispersity, this is the fractional difference $`\epsilon _i(R_iR_0)/R_0`$ of the particle’s radius $`R_i`$ from some reference length $`R_0`$ (with the obvious generalization to charge polydispersity etc.). Henceforth, $`\epsilon `$ shall be referred to as the size parameter, for definiteness. The population of species in the system is characterized by a continuous distribution $`f(\epsilon )`$, which is unnormalized so

$$_{\mathrm{}}^{\mathrm{}}f(\epsilon )𝑑\epsilon =N.$$

In general, the free energy $`F`$ of a polydisperse system is a complicated functional of $`f(\epsilon )`$. It will be expressed in units of $`k_BT`$ where $`k_B`$ is Boltzmann’s constant and $`T`$ is temperature. For a polydisperse ideal gas, the free energy $`F^{\mathrm{id}}`$ is easily shown  to be, per unit volume,

$$\frac{F^{\mathrm{id}}}{V}=𝑑\epsilon \frac{f(\epsilon )}{V}\left[\mathrm{ln}\frac{f(\epsilon )}{V}1\right]$$
(1)   
which is the usual ideal gas free energy, summed over all species. As this expression contains the mixing entropy, it is useful to write the free energy of a non-ideal system as

$$FF^{\mathrm{id}}+F^{\mathrm{ex}}.$$
(2)   
Here $`F^{\mathrm{ex}}`$ is the ‘excess’ part of the free energy (over and above the ideal part), deriving from interactions.

Let us consider a system whose ‘initial’ population (before phase separation) is known. This will be called the ‘parent’ distribution $`f_P(\epsilon )`$. In a system where this parent is partitioned into $``$ coexisting phases, we wish to determine the distribution $`f(\epsilon )_𝒜`$ in each phase $`𝒜=1,\mathrm{},`$. By conservation of matter,

$$\underset{𝒜=1}{\overset{}{}}f(\epsilon )_𝒜=f_P(\epsilon ).$$
(3)   
At equilibrium, the chemical potential is equal in all coexisting phases. This statement applies for each species of particles, so the equation

$$\mu (\epsilon )_𝒜=\mu (\epsilon )_{}\epsilon $$
(4)   
represents an uncountable infinity of thermodynamic constraints for any pair of phases $`𝒜`$ and $``$. Since there is a continuum of species, the chemical potential is a functional derivative of the free energy

$$\mu (\epsilon )\frac{\delta F[f(\epsilon )]}{\delta f(\epsilon )}.$$
(5)   
From Eq. 2, $`\mu (\epsilon )`$ can be written in two parts

$$\mu (\epsilon )=\mu ^{\mathrm{id}}(\epsilon )+\mu ^{\mathrm{ex}}(\epsilon ).$$
(6)   
Functional differentiation of Eq. 1 yields

$$\mu ^{\mathrm{id}}(\epsilon )=\mathrm{ln}\left[\frac{f(\epsilon )}{V}\right].$$
(7)   
Collecting together Eqs. 4, 6 and 7 gives the ratios of densities in any two of the $``$ coexisting phases

$$\frac{f(\epsilon )_{}/V_{}}{f(\epsilon )_𝒜/V_𝒜}=\mathrm{exp}\left(\mu ^{\mathrm{ex}}(\epsilon )_𝒜\mu ^{\mathrm{ex}}(\epsilon )_{}\right)$$
(8)   
in terms of the excess parts of the chemical potentials. Thus, all but one distribution can be eliminated from Eq. 3, yielding the solution for any given phase

$$f(\epsilon )_𝒜=f_P(\epsilon )/\underset{=1}{\overset{}{}}\frac{V_{}}{V_𝒜}\mathrm{exp}\left(\mu ^{\mathrm{ex}}(\epsilon )_𝒜\mu ^{\mathrm{ex}}(\epsilon )_{}\right)$$
(10)   
$$\text{where }\mu ^{\mathrm{ex}}(\epsilon )\frac{\delta F^{\mathrm{ex}}[f(\epsilon )]}{\delta f(\epsilon )}.$$
(11)   
Given a knowledge of $`F^{\mathrm{ex}}`$, which specifies the interactions in the system, (and of the phase volumes), Eqs. 10 and 11 represent a complete solution to the problem. However, they constitute an uncountable infinity of non-linear simultaneous equations. This is the source of the mathematical complexity mentioned earlier.

Some simplification is achieved by making a change of variables. Rather than expressing a thermodynamic state in terms of the densities of the individual species of particles, $`f(\epsilon )/V`$, let us use moments of this distribution (as in ). The thermodynamic variables

$$\rho _\alpha _{\mathrm{}}^{\mathrm{}}\epsilon ^\alpha \frac{f(\epsilon )}{V}𝑑\epsilon ;\alpha =0,1,\mathrm{},\mathrm{}$$
(12)   
will be called ‘moment densities’. Note that $`\rho _\alpha =\overline{\epsilon ^\alpha }\rho `$, so that $`\rho _0`$ is the overall number density $`\rho `$. \[Mean powers of the size parameter, $`\overline{\epsilon ^\alpha }`$, are moments of the normalized distribution $`p(\epsilon )f(\epsilon )/N`$.\] Each moment density, being a linear combination of conserved species densities, is itself conserved and, accordingly, respects the usual equilibrium conditions. For instance, each ‘moment chemical potential’, defined by $`\mu _\alpha (F/V)/\rho _\alpha `$, is equal in coexisting phases. This is clear from expanding the species chemical potential in partial derivatives

$$\mu (\epsilon )\frac{\delta F}{\delta f(\epsilon )}=\underset{\alpha =0}{\overset{\mathrm{}}{}}\frac{F}{\rho _\alpha }\frac{\delta \rho _\alpha }{\delta f(\epsilon )}=\underset{\alpha =0}{\overset{\mathrm{}}{}}\mu _\alpha \epsilon ^\alpha .$$
(13)   
Thus, equality of $`\mu (\epsilon )`$ in coexisting phases requires equality of each $`\mu _\alpha `$.

We now have a discrete set of thermodynamic variables, and can substitute the power series expression (Eq. 13) for $`\mu (\epsilon )`$ into Eq. 10, yielding

$$f(\epsilon )_𝒜=f_P(\epsilon )/\underset{=1}{\overset{}{}}\frac{V_{}}{V_𝒜}\mathrm{exp}\left(\underset{\alpha =0}{\overset{\mathrm{}}{}}(\mu _{\alpha 𝒜}^{\mathrm{ex}}\mu _\alpha ^{\mathrm{ex}})\epsilon ^\alpha \right)$$
(15)   
$$\text{with }\mu _\alpha ^{\mathrm{ex}}\frac{F^{\mathrm{ex}}/V}{\rho _\alpha }$$
(16)   
which, with Eq. 12, form a countable infinity of simultaneous equations. The excess free energy is now a function $`F^{\mathrm{ex}}(\rho _0,\rho _1,\mathrm{})`$ of the moment densities.

The equations thus far are perfectly general, but the advantages of this formalism become apparent when we consider a narrow distribution of sizes, i.e. a system which is close to monodisperse. If the origin for the parameter $`\epsilon `$ is chosen (by fixing the reference $`R_0`$) to be close to the centre of the narrow distribution, then $`\epsilon `$ is a small number for most if not all particles. Hence in Eq. 15, $`f_P(\epsilon )`$ vanishes for large $`\epsilon `$, so the power series in the denominator becomes a well-controlled expansion. The results have a particularly simple form if the origin is chosen to be the mean of the parent distribution, so that $`\overline{\epsilon }_P0`$. Henceforth this choice is assumed.

The solution to Eqs. 12, 15 and 16 is now calculated to first order in $`\epsilon `$. This will yield the exact phase equilibria in the limit of a narrow parent, $`\overline{\epsilon ^2}_P0`$. Expanding Eq. 15 to first order and integrating over $`\epsilon `$ gives

$$\underset{=1}{\overset{}{}}\frac{V_{}}{V_𝒜}\mathrm{exp}(\mu _{0𝒜}^{\mathrm{ex}}\mu _0^{\mathrm{ex}})=\frac{N}{N_𝒜}\left[1+O(\epsilon ^2)\right].$$

To zeroth order, Eq. 8 gives

$$\frac{V_{}}{V_𝒜}\mathrm{exp}(\mu _{0𝒜}^{\mathrm{ex}}\mu _0^{\mathrm{ex}})=\frac{N_{}}{N_𝒜}\left[1+O(\epsilon )\right].$$

Note the different orders of expansion. Substituting these expressions back into 15 yields

$$\frac{f(\epsilon )_𝒜}{N_𝒜}=\frac{f_P(\epsilon )}{N}\left[1\epsilon \mu _{1𝒜}^{\mathrm{ex}}+\frac{\epsilon }{N}\underset{=1}{\overset{}{}}N_{}\mu _1^{\mathrm{ex}}+O(\epsilon ^2)\right].$$
(17)   
To obtain Eq. 17, the prefactor of $`f_P(\epsilon )`$ in Eq. 15 is expanded to first order in $`\epsilon `$, but the distribution $`f_P(\epsilon )`$ itself remains exact. Thus, other than narrowness, no limitations are put on the form of $`f_P(\epsilon )`$. Any distribution can be treated, however asymmetric or discontinuous. The population may even contain finite amounts of some components, contributing delta functions to $`f_P(\epsilon )`$.

In Eq. 17 we see that the distribution in any given phase $`𝒜`$ depends, as one would expect, on the properties of all the other $``$ phases with which it coexists. However, taking the difference (denoted $`\mathrm{\Delta }`$) between the normalized distributions in any two of the $``$ coexisting phases, we find the strikingly simple expression

$$\mathrm{\Delta }p(\epsilon )\epsilon p_P(\epsilon )\mathrm{\Delta }\mu _1^{\mathrm{ex}}$$
(18)   
in the limit as $`\overline{\epsilon ^2}_P0`$, where $`p_P(\epsilon )`$ is the normalized parent distribution. \[Note that the solution for each phase (Eq. 17) is recoverable from the neater sum (Eq. 3) and difference (Eq. 18) equations.\] Surprisingly we have found that, in the multi-phase system, the difference in compositions of any pair of phases is identical to the expression found earlier  for two-phase coexistence. Thus the same universal laws follow , relating any pair of phases. This is not an obvious result, since the parent appearing in Eq. 18 is the combined population of the whole system, not just of the two phases in question as it is in the two-phase coexistence problem.

Equation 18 is very generally applicable. It is valid for any system with a narrow distribution (that is, narrower than the range of linearizability of the fugacity), whatever particles or interactions it comprises. Furthermore, recall that $`\epsilon `$ need not parameterize size deviations, but could represent charge, mass or any other sole polydisperse quantity. By analysing multi-phase coexistence, we have found that Eq. 18 does not even depend on $``$, the number of phases present.

We have considered a system in which a slightly polydisperse fluid component is partitioned among several phases. The coexistence of more than two phases may be the result of tuning the temperature to the triple point of the monodisperse reference system. Alternatively, the slightly polydisperse particles may be in the presence of other, dissimilar components which, by the Gibbs phase rule (which states that an $`n`$-component mixture can exhibit up to $`n+1`$ coexisting phases at arbitrary temperature), can induce multi-phase coexistence . Within such a multi-component system, a particular, slightly polydisperse component will respect the above relations, which may be tested by an experimental probe which is ‘blind’ to the other components. For instance, near-monodisperse colloidal particles in the presence of ‘depletant’ species  exhibit multiple phases. Light scattered only from the near-monodisperse colloid contains information on its fractionation , which should obey Eq. 18. As an illustration, a multi-phase colloidal sample with the composition shown in the figure, will obey the above relations, applied only to those particles in range $`X`$, with the origin of $`\epsilon `$ defined at its centre. The relations are equally applicable to particles in range $`Y`$, if we are blind to all other particles (e.g. they may be made invisible by matching their refractive index to that of the solvent), and redefine $`\epsilon =0`$ appropriately.

The form of the solution in Eq. 18 is of interest in itself, not least for the non-appearance of $``$. However, one quantity remains unknown: the constant of proportionality $`\mathrm{\Delta }\mu _1^{\mathrm{ex}}`$. That constant is system-dependent. For some substances, $`\mu _1^{\mathrm{ex}}`$ can be calculated using thermodynamic perturbation theory . Unfortunately, this is not possible for a system of hard-spheres, as it’s Hamiltonian is non-differentiable. Since the hard-sphere system is of great practical interest for modelling systems with repulsive interactions, the constant of proportionality is now calculated for that case.

The excess part of the free energy of the polydisperse hard-sphere fluid can be Taylor expanded in the small size parameter of each of the $`N`$ particles of interest thus

$$F^{\mathrm{ex}}=F_{\mathrm{mono}}^{\mathrm{ex}}+\underset{i=1}{\overset{N}{}}\epsilon _i\frac{F^{\mathrm{ex}}}{\epsilon _i}|_{\epsilon _i=0}+O(\epsilon ^2)$$

where $`F_{\mathrm{mono}}^{\mathrm{ex}}`$ is the excess free energy of the reference component of monodisperse hard spheres (in the presence of the rest of the system — see Fig. 1). In the reference component, all particles are alike, so the differentiation may be performed on particle number 1 only, without loss of generality, giving

$$F^{\mathrm{ex}}=F_{\mathrm{mono}}^{\mathrm{ex}}+N\overline{\epsilon }\frac{F^{\mathrm{ex}}}{\epsilon _1}|_{\epsilon _1=0}+O(\epsilon ^2).$$

The change in the identity (the species) of particle 1 when its size is varied affects only $`F^{\mathrm{id}}`$. The excess free energy contains the physical effect of the particle’s size on the rest of the system. By its presence in the container, particle 1 simply excludes other particles from a volume $`V_{\mathrm{excl}}`$, given that its interactions are purely hard and repulsive. Thus, increasing its size reduces the effective system volume, so

$$\frac{F^{\mathrm{ex}}}{\epsilon _1}=\frac{F^{\mathrm{ex}}}{V}\frac{dV_{\mathrm{excl}}}{d\epsilon _1}.$$
(19)   
In fact the volume from which particle 1 excludes other particles, $`V_{\mathrm{excl}}`$, depends on their species, so the quantity in Eq. 19 is a net effective value, defined by the equation. For the special case of an almost pure hard sphere system (not in the presence of other, dissimilar components), $`V_{\mathrm{excl}}=\frac{4}{3}\pi \overline{R}_P^3(2+\epsilon _1)^3`$ at low density (correct up to second virial coefficient). At high density, the geometry of high-order inter-particle interactions modifies this. In any case, $`dV_{\mathrm{excl}}/d\epsilon _1`$ is of order a particle volume. The resulting excess free energy density of a polydisperse hard sphere fluid is

$$\frac{F^{\mathrm{ex}}}{V}=\frac{F_{\mathrm{mono}}^{\mathrm{ex}}}{V}+12\rho _1P^{\mathrm{ex}}V^{\mathrm{eff}}+O(\epsilon ^2)$$
(20)   
where $`V^{\mathrm{eff}}`$ is some (unknown) effective volume, of order the volume of an average sphere, and exactly that for a near-pure, low-density system. Applying Eq. 16 yields

$$\mu _1^{\mathrm{ex}}=12P^{\mathrm{ex}}V^{\mathrm{eff}}$$
(21)   
in terms of the system’s excess pressure $`P^{\mathrm{ex}}`$ over an ideal gas. Since coexisting phases have the same total pressure, it follows that $`\mathrm{\Delta }\mu _1^{\mathrm{ex}}=12V^{\mathrm{eff}}\mathrm{\Delta }P^{\mathrm{id}}`$. So the constant of proportionality in Eq. 18 is

$$\mathrm{\Delta }\mu _1^{\mathrm{ex}}=12V^{\mathrm{eff}}\mathrm{\Delta }\rho $$
(22)   
for hard spheres in ergodic (fluid) phases. This calculation contains the lowest-order effects of polydispersity. Once the polydispersity is sufficient to alter the mode of packing (e.g. small particles preferentially filling the gaps between big ones), higher-order analysis is needed.

It is apparent that combining a moment description with a small-variable expansion in the distribution’s width is a productive way to analyse polydisperse systems. While the applications of this study are clearly wide-ranging, it is intended to extend its scope by analysing correlation functions and multiply-polydisperse systems . In addition, some work is required, using higher-order analysis, to establish the radius of convergence of the expansion, and quantify more precisely the method’s regime of validity.

Acknowledgements Many thanks for informative discussions go to Michael Cates, Peter Sollich, David Fairhurst, Patrick Warren and Wilson Poon. The work was funded by the EPSRC (GR/K56025) and a Royal Society of Edinburgh SOEID Research Fellowship.