Phase separation dynamics of polydisperse colloids: a mean-field lattice-gas theory

TL;DR

This paper introduces a mean-field lattice-gas theory for polydisperse colloids, providing analytical insights into phase separation dynamics, including spinodal behavior, fractionation effects, and long-lived heterogeneities.

Contribution

It develops a new dynamical mean-field model for polydisperse colloids that captures complex phase separation behaviors and predicts effects of fractionation and density fluctuations.

Findings

Analytical determination of critical points and spinodal curves.

Identification of qualitative changes in early-time dynamics across a quenched spinodal.

Observation of long-lived composition heterogeneities due to evaporation and fractionation.

Abstract

New insights into phase separation in colloidal suspensions are provided via a new dynamical theory based on the Polydisperse Lattice-Gas model. The model gives a simplified description of polydisperse colloids, incorporating a hard-core repulsion combined with polydispersity in the strength of the attraction between neighbouring particles. Our mean-field equations describe the local concentration evolution for each of an arbitrary number of species, and for an arbitrary overall composition of the system. We focus on the predictions for the dynamics of colloidal gas-liquid phase separation after a quench into the coexistence region. The critical point and the relevant spinodal curves are determined analytically, with the latter depending only on three moments of the overall composition. The results for the early-time spinodal dynamics show qualitative changes as one crosses a 'quenched' spinodal that excludes fractionation and so allows only density fluctuations at fixed composition. This effect occurs for dense systems, in agreement with a conjecture by Warren that, at high density, fractionation should be generically slow because it requires inter-diffusion of particles. We verify this conclusion by showing that the observed qualitative changes disappear when direct particle-particle swaps are allowed in the dynamics. Finally, the rich behaviour beyond the spinodal regime is examined, where we find that the evaporation of gas bubbles with strongly fractionated interfaces causes long-lived composition heterogeneities in the liquid phase; we introduce a two-dimensional density histogram method that allows such effects to be easily visualized for an arbitrary number of particle species.

Figures

Tables

Table 1. Table 1: Abbreviations: SD = Spinodal Decomposition; NG = Nucleation and Growth; fluct. = fluctuations.

Region	Stage 1	Stage 2
R1	SD in density	SD of ‘sizes’
R2	NG of density fluct.	SD of ‘sizes’
R3	—	SD of ‘sizes’
R4	—	NG of ‘size’ fluct.

Equations

$$H=-\sum _{⟨i,j⟩}\sum _{\alpha ,\beta }{\sigma }_{\alpha }{\sigma }_{\beta }{n}_i^{\alpha }{n}_j^{\beta }$$

$$\sum _{\gamma =0}^M{n}_i^{\gamma }=1\forall i$$

$$F=-\sum _{⟨i,j⟩}\sum _{\alpha ,\beta }{ϵ}_{\alpha \beta }{p}_i^{\alpha }{p}_j^{\beta }+T\sum _i\sum _{\gamma =0}^M{p}_i^{\gamma }\ln p_i^{\gamma }$$

$$f(\left\{p^{\alpha }\right\},T)=-\frac{z}{2}\sum _{\alpha ,\beta }{ϵ}_{\alpha \beta }{p}^{\alpha }{p}^{\beta }+T\sum _{\gamma =0}^M{p}^{\gamma }\ln p^{\gamma }$$

$$S_{\alpha \beta }=\frac{{\partial }^{2}f}{\partial p^{\alpha }\partial p^{\beta }}$$

$$T=z\left({\rho }_2-{\rho }_1^2\right).$$

$$T=z{\left(\frac{{\rho }_1}{\rho }\right)}^2\left(\rho -{\rho }^2\right).$$

$$\frac{\partial P(C,t)}{\partial t}=\sum _{C^{\prime }}\left[W(C^{\prime },C)P(C^{\prime },t)-W(C,C^{\prime })P(C,t)\right]$$

$$⟨Q⟩(t)=\sum _{C}Q(C)P(C,t)$$

$$\frac{d{p}_i^{\alpha }}{dt}=-\sum _{j\in \partial i}{J}_{ij}^{\alpha }$$

$$J_{ij}^{\alpha}=\sum\limits_{\gamma=0}^{M}\Bigl{\langle}n_{i}^{\alpha}n_{j}^{\gamma}w_{ij}^{\alpha\gamma}-n_{j}^{\alpha}n_{i}^{\gamma}w_{ji}^{\alpha\gamma}\Bigr{\rangle}$$

$$w_{ij}^{\alpha \gamma }=w^{\alpha \gamma }\left[\frac{1}{1+\exp \left(\Delta H_{ij}^{\alpha \gamma }/T\right)}\right]$$

$$H(C)=-{ϵ}_{\alpha \gamma }-\sum _{\beta }\sum _{k\in \partial i\backslash j}{ϵ}_{\alpha \beta }{n}_k^{\beta }-\sum _{\beta }\sum _{l\in \partial j\backslash i}{ϵ}_{\gamma \beta }{n}_l^{\beta }+\dots$$

$$H(C^{\prime })=-{ϵ}_{\alpha \gamma }-\sum _{\beta }\sum _{k\in \partial i\backslash j}{ϵ}_{\gamma \beta }{n}_k^{\beta }-\sum _{\beta }\sum _{l\in \partial j\backslash i}{ϵ}_{\alpha \beta }{n}_l^{\beta }+\dots$$

$$\Delta H_{ij}^{\alpha \gamma }=$$

$$-\sum _{\beta }\left(\sum _{k\in \partial i\backslash j}{ϵ}_{\gamma \beta }{n}_k^{\beta }-\sum _{l\in \partial j\backslash i}{ϵ}_{\gamma \beta }{n}_l^{\beta }\right).$$

$$\displaystyle\frac{dp_{i}^{\alpha}}{dt}=-\sum\limits_{j\in\partial i}\sum\limits_{\gamma=0}^{M}\left[\frac{p_{i}^{\alpha}p_{j}^{\gamma}w^{\alpha\gamma}}{1+\exp{\left(\bigl{\langle}\Delta H_{ij}^{\alpha\gamma}\bigr{\rangle}/T\right)}}-\frac{p_{j}^{\alpha}p_{i}^{\gamma}w^{\alpha\gamma}}{1+\exp{\left(\bigl{\langle}\Delta H_{ji}^{\alpha\gamma}\bigr{\rangle}/T\right)}}\right]$$

$$\displaystyle\bigl{\langle}\Delta H_{ij}^{\alpha\gamma}\bigr{\rangle}={}$$

$$-\sum _{\beta }\left(\sum _{k\in \partial i}{ϵ}_{\gamma \beta }{p}_k^{\beta }-\sum _{l\in \partial j}{ϵ}_{\gamma \beta }{p}_l^{\beta }\right).$$

$$p_i^{\alpha }=p^{\alpha }+{\delta }_i^{\alpha }$$

$$\frac{d\delta_{i}^{\alpha}}{dt}=\sum\limits_{\gamma=0}\mathscr{M}^{\alpha\gamma}\Delta_{\rm d}\Bigl{(}\mu_{i}^{\alpha}-\mu_{i}^{\gamma}\Bigr{)}$$

$${\mathcal{M}}^{\alpha \gamma }\equiv \frac{w^{\alpha \gamma }}{2T}{p}^{\alpha }{p}^{\gamma }$$

$${\mu }_i^{\alpha }=-\sum _{j\in \partial i}\sum _{\beta }{ϵ}_{\alpha \beta }{p}_j^{\beta }+T\ln \left(p_i^{\alpha }/p_i^0\right)$$

$${\Delta }_{\mathrm{d}}{g}_i=\sum _{j\in \partial i}\left(g_j-g_i\right)$$

$${\delta }_j^{\alpha }=\delta p^{\alpha }\exp \left[i\mathbf{k}\cdot {\mathbf{x}}_j+\omega t\right]$$

$$\omega =$$

$$\displaystyle\left.+w_{s}\left(\rho\rho_{2}-\rho_{1}^{2}\right)\biggr{]}\right\}-\frac{A(\mathbf{k})}{4T}\biggl{[}\left(T\left(w_{0}-w_{s}\right)\rho\right)^{2}$$

$$\displaystyle+\left(A(\mathbf{k})+z\right)^{2}\left(w_{0}\left(1-\rho\right)\rho_{2}+w_{s}\left(\rho\rho_{2}-\rho_{1}^{2}\right)\right)^{2}\biggr{.}$$

$$+2T\left(w_0-w_s\right)\left(A(\mathbf{k})+z\right)(w_0\left(1-\rho \right)\left(\rho {\rho }_2-2{\rho }_1^2\right)$$

$$\displaystyle\biggl{.}\left.+w_{s}\rho\left(\rho\rho_{2}-\rho_{1}^{2}\right)\right)\biggr{]}^{1/2}$$

$$f^{(0)}(\sigma )=\frac{1}{g!}{\left(\frac{g+1}{\bar{\sigma }}\right)}^{g+1}{\sigma }^g\exp \left[-\left(\frac{g+1}{\bar{\sigma }}\right)\sigma \right].$$

$$\bar{\sigma }{\rho }_0-{\rho }_1=(\bar{\sigma }-{\sigma }_A)p^A+(\bar{\sigma }-{\sigma }_B)p^B=d(p^B-p^A)$$

$$H=-\sum _{⟨i,j⟩}\sum _{\alpha ,\beta }{ϵ}_{\alpha \beta }{n}_i^{\alpha }{n}_j^{\beta }$$

Full text

\makeFNbottom

[TABLE]

††footnotetext: Disordered Systems Group, Department of Mathematics, King’s College London, London, United Kingdom. E-mail: pablo.decastro@kcl.ac.uk††footnotetext: † Electronic Supplementary Information (ESI) available: Animation shows evolution of phase-separation dynamics snapshots for a binary mixture (left) and their two-dimensional density histograms (right). Parameters: $p^{A}=p^{B}=0.435$, $T=0.5$, $d=0.25$, $w_{0}=1$, $w_{s}=0$, $L=150$, from $t=0$ to $t=4497.5$. See DOI: 10.1039/b000000x/

1 Introduction

Traditionally, thermodynamics and statistical mechanics deal with systems composed of identical particles, also called monodisperse systems. However, nature is often more complex: many soft matter systems, such as biological and industrial fluids (in the form of colloidal dispersions, liquid crystals, polymer solutions, etc.), are polydisperse in the sense that their constituent particles exhibit variation in terms of one (or several) attribute(s).1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 This polydisperse attribute, denoted here by $\sigma$, can be particle size, shape, charge, molecular weight, chemical nature, etc.16 In the thermodynamic limit, these systems are usually regarded as having effectively an infinite number of components (although the theory that will be presented here applies to systems with an arbitrary number of components).17, 18 Everyday and not-so-everyday examples include blood, paint, milk, clay, photonic crystals, shampoo, viruses, globular proteins, pharmaceuticals, and even sewage, among many others.19, 20, 21, 22, 23, 24, 25, 26, 27, 28 Therefore, understanding the impact of polydispersity on the phase behaviour of many-body systems is of fundamental, commercial, and practical interest.29 For instance, knowing under what conditions (and how) a multicomponent fluid will demix may be essential in determining the shelf life of a product.30, 31 On the fundamental end of the range, it is believed that (symmetry-breaking) phase transitions have occurred in the early universe via nucleation of vacuum-field bubbles, thus generating the different fundamental forces of nature; despite having a different sort of fluid dynamics and not being in the scope of the present work, it is interesting that the early-universe hot plasma can also be seen as a multicomponent fluid.32, 33 Our focus here is on colloidal fluids, for which the fluid particles are influenced by the thermal agitation from the solvent in which they are suspended, and quantum effects can be neglected. (See ref. 34 for an introduction to polydisperse colloidal systems and their phase behaviour.)

The state of a polydisperse system (or any of its phases) is described by a density distribution $\rho(\sigma)$, defined such that $\rho(\sigma)d\sigma$ is the number density of particles with $\sigma$-values in the range $\left[\sigma,\sigma+d\sigma\right]$. One can also consider a multicomponent system, composed of a finite number $M$ of particle species (each one labelled by its own value of $\sigma$). Each component has concentration $p^{\alpha}$, with $\alpha=1,2,\dots,M$. (Note that the $\alpha$ here is a superscript, not an exponent.) The continuous approach can be thought of as the limit $M\rightarrow\infty$, leading to $p^{\alpha}\rightarrow\rho(\sigma_{\alpha})d\sigma$, where $\sigma_{\alpha}$ is the value of $\sigma$ associated with species $\alpha$. In either case, the system has an overall composition, which specifies the ratios of the densities of different species. The simplest example is the one of a binary system, i.e. $M=2$: its composition could then be specified by the ratio $p^{A}/p^{B}$, for instance. We will often specialize to this binary case for ease of explanation.

Similarly to monodisperse fluids, polydisperse fluids can phase-separate into regions with higher and lower concentrations of particles. The typical experiment is to decrease the temperature of a homogeneous system to a value within the coexistence region, setting off a dynamical evolution towards separation into two (or more) equilibrium phases. (See Fig. 1 and the animation provided in the ESI.†) In polydisperse phase separation, however, these phases will not just differ in density but generally also in composition. This process is called fractionation (or partitioning) 17 because it implies that particles of different species distribute themselves unevenly into the new phases. Fractionation is responsible for much of the complexity in the phase behaviour of polydisperse systems; 17, 35, 36, 37 its effect on phase-separation kinetics will be our main focus.

Let us explain fractionation in more detail using the $M=2$ example. Fig. 2 shows a schematic example of a phase diagram of a binary fluid at a fixed temperature. The purple (middle) point indicates the species densities $(p^{A},p^{B})$ of the parent phase, i.e. of the initial homogeneous system. Any other phase of the same composition would lie on the dashed line through the parent and the origin. This is called the dilution line because in a colloidal fluid such phases can be prepared by adding or removing solvent. Phase separation can occur into pairs of coexisting phases – identified from the requirements of equal pressure and species chemical potentials – that are shown as end points of tielines. Due to particle conservation, the overall system composition must remain unchanged during phase separation, so the actual daughter phases generated must lie on the tieline passing through the parent. In the generic case these daughter phases both lie off the dilution line so have a composition that is different from the parent (and from each other). This is the phenomenon of fractionation.

In the case of general $M$ similar considerations apply. Here the parent is specified by a density distribution $\rho^{(0)}(\sigma)$. This can be decomposed as $\rho^{(0)}(\sigma)=\rho^{(0)}_{0}f^{(0)}(\sigma)$, where $\rho^{(0)}_{0}=N/V$ is the overall particle number density and $f^{(0)}(\sigma)$, the normalized parent shape function, specifies the composition. As $\rho^{(0)}_{0}$ is varied, $\rho^{(0)}(\sigma)$ traces out the dilution line in density distribution space. To obtain phase diagrams one needs to project from this $M$-dimensional space. Often only the density of coexisting phases is recorded, to recover the polydisperse analogue of a monodisperse density–temperature phase diagram (see Fig. 3a).

Fractionation causes the equilibrium phase behaviour of polydisperse systems to be much richer than that of their monodisperse counterparts.17 To see why, note that for large $M$, the classical tangent plane procedure used to find phase coexistence from a given free energy becomes awkward both conceptually and computationally. Also, the number of coexisting phases allowed by Gibbs’ phase rule grows with $M$, and becomes arbitrary in the limit $M\to\infty$. These challenges have been tackled using a number of theoretical methods, uncovering many qualitatively novel effects of polydispersity on phase behaviour. 17, 38, 39, 40 For instance, Evans has developed a perturbative approach for narrow distributions of the polydisperse attribute.16, 36, 41, 42, 43 This predicts, for example, that the difference in mean particle size between the two daughter phases is proportional to the variance of the parent distribution. It can also predict specific trends, e.g. that in the fluid-solid coexistence of size-disperse colloidal hard spheres, the solid (crystalline) phase contains on average larger particles than the fluid.36 Going beyond perturbative approaches, the full equilibrium phase diagram for polydisperse hard spheres has been found via accurate free energy expressions for the fluid and solid phases, showing very different equilibrium properties already for moderate particle size spreads.44 The physical intuition behind this is the fact that particles of different sizes cannot easily fit into a regular crystal lattice. These results were derived using a moment projection method.17, 45, 46, 47 This exploits a mapping of the full density distribution space onto moment densities, defined as $\rho_{n}\equiv\int\sigma^{n}\rho(\sigma)d\sigma$, where one integrates over all possible values of $\sigma$. Note that the zeroth moment density $\rho_{0}$ is just the overall density $N/V$. Another important parameter, the standard deviation of the polydisperse attribute $\sigma$ normalized by its mean, can be expressed in terms of moment densities as $[(\rho_{2}\rho_{0}/\rho_{1}^{2})-1]^{1/2}$. It is often referred to as degree of polydispersity or simply polydispersity.

The present work focuses on the dynamics towards equilibrium. Warren proposed that the phase-separation kinetics in polydisperse systems should proceed in two stages.48 To understand this two-stage scenario, he starts by considering an initially homogeneous system (suddenly placed within a two-phase region). The classical picture predicts two types of phase-separation kinetics. If placed not too far from the phase boundary encompassing the two-phase region in the phase diagram, the system would phase-separate via nucleation and growth. However, if placed further away from the boundary phase separation would proceed via spinodal decomposition. Lying entirely below the coexistence curve, the classical spinodal curve separates the two types of kinetics. In the case of phase separation via spinodal decomposition, arbitrarily small fluctuations in the homogeneous system decrease the free energy, and thus grow exponentially before eventually being limited by nonlinear effects. Warren highlights the fact that, although the spinodal curve is a convention that does not rigorously survive in a correct statistical mechanical treatment, it is a rough guide to the transition region between the two types of kinetic behaviour. He proceeds by asking what the effect of polydispersity is on this classical picture. Assuming that in dense systems, fractionation is potentially a very slow process, he suggests that two stages of kinetics might be expected. In a first stage, the system might phase-separate by relaxing the overall density to a phase equilibrium dictated by ‘quenching’ (i.e. ‘fixing’) the polydisperse distribution in any phase to the one in the parent phase. Over longer timescales, the system might then gradually redistribute particles between coexisting phases to reflect the possibility of further lowering the free energy by fractionating. Therefore, a second stage would begin when the system started to significantly fractionate. Referring to the polydisperse density distribution as the size distribution, Warren says: ‘During this slow process, the particle number density would be able to adjust itself continually to follow the current stage of the size partitioning process. The argument for this [two-stage] scenario is based on the observation that, in a dense system, the overall density can be much more easily relaxed by collective particle motion (collective diffusion), than can the size distribution be adjusted by individual particle rearrangements (self diffusion)’ (p. 2199). Experimental measurements of self- and collective diffusion constants support this picture.49

Warren’s two-stage scenario is certainly physically reasonable: density fluctuations can be created by moving groups of particles ‘in sync’ to form regions of higher density. Fractionation, on the other hand, requires particles from different species to ‘push past each other’ in opposite directions. The scenario also makes an interesting connection to moment densities: the zeroth moment $\rho_{0}$, which gives the total number density of all species, should relax much more rapidly at a local level than the higher moments, whose equilibration requires interdiffusion of different particle species. This suggests that moment densities can remain useful in understanding the kinetics of phase separation, and our results in this paper provide some support for this.

Existing theoretical approaches to polydisperse phase-separation dynamics have focussed primarily on polymeric systems. One could make some progress by binning the range of $\sigma$, which reduces the problem to the dynamics of a finite mixture. This approach of course becomes numerically challenging as the number of bins grows. Clarke has shown that the method is nonetheless useful for investigating the early-time phase-separation dynamics of polydisperse polymers.17, 50 Again in the context of polymeric materials, Pagonabarraga and Cates developed an analysis of the dynamics based on time evolution equations for polymer densities driven by chemical potential gradients. The form of these equations had been proposed phenomenologically in Clarke’s work, but the approach of ref. 51 allowed the associated mobility coefficients to be derived explicitly. Apart from the case of length polydispersity, where the results were more subtle, Pagonabarraga and Cates studied the case of chemical polydispersity, where different polymer chains have different monomer compositions and hence different interaction strengths. For this scenario the coupled dynamical equations could be fully solved in certain cases. Pagonabarraga and Cates51 also studied the mode spectrum of the various density fluctuations in a system undergoing spinodal decomposition; from this analysis, they concluded that depending on where in the phase diagram the system is placed the kinetics will proceed in accordance with Warren’s two-stage scenario.

Compared to the polymeric case, there is a lack of theoretical work designed to model fractionation effects in the phase-separation dynamics of spherical colloids (but see ref. 41 below). This paper is designed to fill this gap. The approach we use to investigate the phase-separation dynamics of mixtures is the theory described in refs. 52 and 53 by Plapp and Gouyet. These studies were, however, concerned with binary metallic alloys. In this context they addressed rather different questions from ours, based on assumptions about the dynamics and the particle interactions that are quite distinct from the colloidal case. Nonetheless, our development of the mean-field dynamical equations has close similarities with the methods of refs. 52 and 53.

There are a few other theoretical investigations of polydisperse colloidal dynamics in the literature, where simulations are implemented and some aspects of fractionation investigated.54, 55, 42 As a result, the dynamics of phase separation in polydisperse colloids remains a challenging (and mostly unsolved) problem. One of the difficulties is that the kinetics of phase separation could be so slow as to make the actual equilibrium phase compositions unobservable in experiments.17 It has been argued that this is the case for polydisperse hard-sphere crystals: once particles join to a crystal nucleus growing from the hard-sphere fluid, they essentially no longer diffuse on experimental timescales.41 The size distribution of particles in the crystal will thus ‘freeze in’, and will be determined by the mechanism of crystal growth rather than the conditions of thermodynamic equilibrium. Although recent advances have been made,56 such non-equilibrium effects on the experimentally observed phase behaviour of colloidal systems are definitely not fully understood.

In this paper we present a mean-field theory for the Polydisperse Lattice-Gas (PLG) model, which has been proposed as a simple description of polydisperse colloids.35 In Section 2 we present the model and its mean-field phase diagram. In Section 3 we endow the model with an appropriate dynamics and derive the mean-field evolution equations for this. Section 4 gives an early-time regime analysis of our equations. Then, in Section 5, we go beyond the spinodal regime and study the full late-time dynamics of our model. Section 6 summarizes our results and outlines some directions for future research.

2 PLG mean-field phase diagram

2.1 PLG model

Hereafter we use the PLG model to investigate the phase-separation dynamics in polydisperse colloidal systems.35 The model is described by the Hamiltonian

[TABLE]

where $i$ runs over the sites of a periodic lattice $i=1,\dots,L^{D}$, assumed simple cubic and $D$-dimensional in this work, with lattice spacing $a=1$, unless otherwise stated; the sum runs over all pairs $\langle i,j\rangle$ of nearest-neighbour sites; $\sigma_{\alpha}$ is the value of the polydisperse attribute associated with particle species $\alpha$, which controls the strength of interparticle interactions; it is a positive number for attractive interactions as considered here. We consider a mixture with $M$ species, with the summations over $\alpha$ and $\beta$ therefore running from $1$ to $M$. The (occupation) variable $n_{i}^{\alpha}$ simply counts the number of particles of species $\alpha$ at site $i$, for which a hard-core constraint is imposed:

[TABLE]

where $n_{i}^{0}$ refers to vacancy, i.e. $n_{i}^{0}=1$ indicates the presence of a vacancy at site $i$, or, equivalently, a solvent particle. Note that the solvent particles are ‘passive’51 in this framework, in such a way that any non-hydrodynamic effect caused by them has already been effectively included in the interaction between the particles as described by the model Hamiltonian. We also neglect the fact that colloidal particles may interact with one another via hydrodynamic interactions mediated by the solvent. We therefore use the term vacancy in the following, rather than solvent particle. Observe that $n_{i}^{0}$ can be expressed in terms of the other $n_{i}^{\alpha}$, from eqn (2). In summary, each lattice site may be either vacant or occupied by a single colloidal particle of polydisperse attribute $\sigma_{\alpha}$. The instantaneous density distribution follows as $p^{\alpha}=L^{-D}\sum_{i}n_{i}^{\alpha}$. We will use the letter $\gamma$ as the species index for summations running from [math] to $M$, while for summations running from $1$ to $M$, we use $\alpha$ (or $\beta$), unless otherwise specified.

Note that in the Hamiltonian (1) the interaction strength between any two neighbouring particles is assumed to be $\sigma_{\alpha}\sigma_{\beta}$, i.e. the product of their polydisperse attributes, though to preserve generality we will often write this in the form $\epsilon_{\alpha\beta}=\sigma_{\alpha}\sigma_{\beta}$. Thus the role of polydispersity in this model is to engender a spread of possible interaction strengths between particles, a situation which contrasts with the single interaction strength characterizing the simple Ising lattice-gas model. As observed in ref. 35, this allows the PLG model to capture the essential qualitative features that distinguish polydisperse fluids from their monodisperse counterparts. Nonetheless one has to bear in mind that, in real colloids, polydispersity often occurs in the size of the particles. This will have additional consequences, e.g. on the local packing of particles in dense regions, that a lattice model cannot account for. In principle one could extend the approach by allowing particles to occupy several contiguous lattice sites, thus explicitly representing their size. This is not a trivial extension as particle moves would then correspond to simultaneous changes of potentially many $n_{i}^{\alpha}$, but may be an interesting avenue for future work.

2.2 Mean-field phase diagram

We next present the mean-field phase diagram for the PLG model. This will serve as a useful reference point for our later discussion of the phase-separation dynamics. We first explain how to obtain the relevant curves of the diagram, starting from the spinodals and then moving on to the cloud (‘binodal’) and shadow curves.

2.2.1 Free energy and spinodals

Phase separation via spinodal decomposition occurs when the system is placed in a region of the phase diagram that is unstable to fluctuations; fluctuations of any size will then lead to phase separation. This contrasts with the case of nucleation and growth, where finite fluctuations – corresponding to a nucleus of a critical size – are required for the system to escape from a metastable state.

For the phase diagram of monodisperse systems, the spinodal curve (i.e. the curve below which spinodal decomposition occurs) can be calculated by joining up the inflection points of the free energy curve at each temperature. For the spinodal curve of a polydisperse system, one likewise needs to identify points where the free energy develops negative curvature.

In order to obtain an expression for the Helmholtz (inhomogeneous) free energy of the polydisperse lattice gas, $F=\left\langle H\right\rangle-TS$, we use a variational mean-field approximation. This is obtained from the Gibbs–Bogoliubov–Feynman variational bound,57, 58 using a variational approximation to the equilibrium distribution that is factorized over lattice sites. The latter is fully characterized by local densities $p_{i}^{\alpha}=\left\langle n_{i}^{\alpha}\right\rangle$. For the PLG Hamiltonian (1), this leads to

[TABLE]

where the Boltzmann constant has been set equal to $1$. Recalling that the volume of the system is equal to the total number of sites ($V=L^{D}$) and applying eqn (3) to a homogeneous configuration, the free energy density $f=F/V$ can be written as

[TABLE]

where $z$ is the lattice coordination number [i.e. $z=4$ for the square lattice (in $D=2$)] and $p^{\alpha}$ is the overall density of particles of species $\alpha$. Here the hard-core constraint (2) means that in the second, entropic term, $p^{0}=1-\sum_{\alpha}p^{\alpha}$. For a monodisperse system (i.e. $M=1$), the free energy density reduces to a function $f=f(\rho,T)$, where $\rho=N/V\equiv p^{1}$ is the overall density of particles. Here and in the following we drop the 0-subscript from our previous notation $\rho_{0}$ where the meaning is clear from the context. (Confusion with the density distribution $\rho(\sigma)$ should not arise as this is always written with its polydisperse attribute argument $\sigma$.) In this monodisperse case, the spinodal density at any temperature $T$ would be found by from the equation $\frac{\partial^{2}f}{\partial\rho^{2}}=0$. In the more general, polydisperse, case one needs to consider the Hessian matrix $S$ with elements

[TABLE]

and then solve the equation $\det(S)=0$ (or the equivalent spinodal criterion as written in ref. 47) to obtain the spinodal curve.

For polydisperse systems, one can, in the spirit of Warren’s two-stage scenario, define two types of spinodals.48 The first is the annealed spinodal curve, which is precisely the one given by $\det(S)=0$. For the PLG model, this leads to a simple expression in terms of the moment densities:

[TABLE]

The quenched spinodal curve, on the other hand, is defined as the spinodal curve that is calculated by assuming that the system is not allowed to fractionate. In other words, one treats the composition as fixed, and the overall density $\rho$ as the only variable. To calculate this, note that the free energy given by eqn (4) is a function of all densities $p^{\alpha}$. At fixed composition these are proportional to the overall density $\rho$, e.g. for $M=2$ one would have $p^{A}=p^{B}=\rho/2$ if the dilution line $p^{A}/p^{B}=1$ is considered. Inserting these $\rho$-dependencies into the free energy gives the quenched free energy, as a function of $\rho$; call this $f_{Q}(\rho,T)$, for example. To find the quenched spinodal one then only has to solve $\frac{\partial^{2}f_{Q}}{\partial\rho^{2}}=0$, which in our case yields

[TABLE]

One can easily check that this coincides with eqn (6) for monodisperse systems, where $\rho_{2}/\rho=(\rho_{1}/\rho)^{2}$.

The above distinction between two types of spinodal will prove useful later because, if Warren’s two-stage scenario holds, one would expect that the dynamics of phase separation would proceed, at least initially, as if the system ‘did not know’ it could further lower its free energy by fractionating. Therefore, one could expect that the early-stage dynamics of a system operating in the two-stage scenario should proceed in accordance with the quenched spinodal instead of the annealed one.

For completeness we note that the critical point lies on the (annealed) spinodal curve and is identified by an additional condition. This can be obtained using the methods of ref. 47 and reads $2\rho_{1}^{3}-3\rho_{1}\rho_{2}+\rho_{3}=0$.

2.2.2 Cloud and shadow curves

Because of fractionation, the conventional vapour-liquid binodal (or coexistence curve) of a monodisperse system splits into a cloud curve, marking the onset of (polydisperse) phase coexistence, and a shadow curve, giving the density of the incipient phase. The critical point appears at the intersection of these curves, rather than at the maximum of either. 17 This splitting is seen in experiments on polydisperse fluids (see e.g. ref. 59). Similarly to the spinodal curve, we can define annealed and quenched cloud curves. These will be distinct because the onset of phase coexistence will be delayed to lower temperature if the system is not allowed to lower its free energy by fractionation.

Our numerical data for cloud and shadow curves are determined by solving the equations for two-phase (bulk) phase equilibrium. (Three-phase coexistence can also occur in the PLG, but at lower temperatures than we consider here.) In the binary case ($M=2$), the phase equlibrium conditions at given $T$ are equality of the chemical potentials $\mu^{A}$, $\mu^{B}$ and the pressure between the two phases $(p^{A}_{\thinspace\thinspace\!{}_{\textup{I}}},p^{B}_{\thinspace\thinspace\!{}_{\textup{I}}})$ and $(p^{A}_{\thinspace\!{}_{\textup{II}}},p^{B}_{\thinspace\!{}_{\textup{II}}})$. Here the chemical potentials $\mu^{A}$ and $\mu^{B}$ are given by the partial derivatives of the free energy density $f(p^{A},p^{B},T)$ with respect to $p^{A}$ and $p^{B}$, respectively, and the pressure is $P=-f+\mu^{A}p^{A}+\mu^{B}p^{B}$. For $M>2$, the generalization is straightforward, giving $M+1$ coexistence equations to solve. (Alternatively, one could use the moment free energy method, since our free energy density (4) is truncatable.47) The cloud temperature for a given parent density is found by lowering $T$ and checking when phase coexistence with the parent as one of the coexisting phases is first found. The shadow curve identifies the density of the second coexisting phase, which at this point is present only in an infinitesimal fraction of the system volume. In a monodisperse system, cloud and shadow curve coincide and are then identical to the conventional binodal curve.

2.2.3 Phase diagram

Fig. 3a shows our results for the phase diagram of a binary mixture. We chose the dilution line $p^{A}=p^{B}$ and $\sigma$-values given by $\sigma_{A}=1+d$ and $\sigma_{B}=1-d$, where $d$ is a number between [math] and $1$. (Note that the greater the $d$, the greater the polydispersity in this bidisperse case.) The plot shows the annealed cloud and shadow curves, with the critical point at their intersection; the quenched binodal curve; and the annealed and quenched spinodal curves. Qualitatively these curves look as one would expect them to on general grounds.17

For our discussion of phase-separation dynamics most relevant are the cloud curve as it signals the onset of phase coexistence, and the spinodal curve. The annealed and quenched versions of each of these two curves divide the phase diagram into distinct regions where Warren’s hypothesis predicts different sequences of phase-separation dynamics. Inside the quenched spinodal, the system should initially phase-separate by spinodal decomposition in density only; between the quenched spinodal and quenched binodal the first stage should be nucleation and growth of density fluctuations. The second stage would then involve fractionation, again by spinodal decomposition or nucleation and growth depending on the position relative to the annealed spinodal.

As an example, Fig. 3b shows a zoomed-in portion of the phase diagram. Starting at $\rho=0.5$ and moving towards higher densities at constant temperature $T=0.94$, one visits four distinct regions with respect to the annealed and quenched cloud and spinodal curves. Labelling these by R1, R2, R3, R4 with increasing density gives the predictions in Table 1, where we have used the term ‘size’ generically to refer to the polydisperse attribute $\sigma$. In this paper we develop a mean-field theory that excludes fluctuations, hence does not capture nucleation and growth dynamics. Hence we will focus on the distinction between different spinodal regions, i.e. region R1 versus R2/R3. Note that we will use the term ‘early-time dynamics’ throughout to refer to the onset of phase-separation dynamics, irrespective of the relevant time scale. In particular, in region R2/R3 the early-time dynamics should be stage 2 spinodal decomposition, which happens on a slow time scale because it involves fractionation.

3 Kinetic PLG model

In this section we present the kinetics we assume for the PLG model with an arbitrary number of species $M$, and derive our mean-field dynamical equations.

We model the dynamics as resulting from jumps of particles to nearest-neighbour sites. This can be described in the general form of a master equation

[TABLE]

where $P(C,t)$ is the probability of finding the system in a configuration $C=\{n_{i}^{1},n_{i}^{2},\dots,n_{i}^{M};i=1,\dots,L^{D}\}$ at time $t$, and $W(C,C^{\prime})$ is the transition rate from a configuration $C$ to another one $C^{\prime}$. One can then define a time-dependent average for any observable ${Q}(C)$ as

[TABLE]

In particular we define a time-dependent local density for each species as $p_{i}^{\alpha}(t)=\langle n_{i}^{\alpha}\rangle(t)$. From the master equation it can be shown that each species follows a conservation law of the form

[TABLE]

where we have introduced the notation $j\in\partial i$, meaning that the summation has to be performed over all nearest-neighbour sites $j$ of site $i$; the current of $\alpha$-particles across the link $i\rightarrow j$ is given by

[TABLE]

where $w_{ij}^{\alpha\gamma}$ is the jump rate for an $\alpha$-particle at site $i$ to exchange positions with a $\gamma$-particle at site $j$ (or with a vacancy, if $\gamma=0$). Observe that, effectively, we are considering a kinetic model equipped with two elementary processes: (i) the jump of a particle from an occupied lattice site to an empty one, and also (ii) the direct interchange between two particles from arbitrary species. Throughout this process, the configuration of the system changes, but the overall number of particles of each species remains constant, which is the physical origin of the conservation law (10). It is worth noting that the physical elementary process are the jumps to empty sites, i.e. particle-vacancy exchanges. However, for the subsequent analysis it is useful to include also direct particle-particle swaps, so that we develop the theory initially for generic $w^{\alpha\gamma}_{ij}$.

Note from that the prefactors of $w^{\alpha\gamma}_{ij}$ in eqn (11) ensure the target site is occupied by a $\gamma$-particle (or empty, if $\gamma=0$) and the start site is occupied by an $\alpha$-particle. This is similar to the theory described by Plapp and Gouyet in ref. 52, but there they did not consider direct interchange of particles, and the number of species $M$ was set equal to $2$ (as previously mentioned); moreover, they considered Arrhenius law jump rates, which was argued to be more suitable in their context of phase separation in binary metallic alloys. Here, we use Glauber-like jump rates:

[TABLE]

where $\Delta H_{ij}^{\alpha\gamma}=H(C^{\prime})-H(C)$ is the energy difference associated with the jump, i.e. it is the energy difference between the configuration before the jump ($C$), in which there is an $\alpha$-particle at site $i$ and a $\gamma$-particle at site $j$, and after the exchange ($C^{\prime}$); note that $C$ and $C^{\prime}$ are identical except for the exchange between $\alpha$ and $\gamma$. The prefactor $w^{\alpha\gamma}$ is an ‘attempt rate’. This gives the actual jump rate $w_{ij}^{\alpha\gamma}$ when the energy change $\Delta H_{ij}^{\alpha\gamma}$ is large and negative, while it is reduced exponentially for large positive energy changes.

Let us write the energy of the system for the configuration before the jump, isolating only the contributions involving the two particles $\alpha$ and $\gamma$ that are swapping:

[TABLE]

Here the first term on the right-hand side is the contribution only from the interaction between the two particles that are swapping. The second term is the contribution from the interactions between the $\alpha$-particle at site $i$ and all of its neighbours except the one at site $j$. The third term similarly accounts for the interactions of the $\gamma$-particle. Similarly, we have that

[TABLE]

is the energy of the system after the exchange.

The energy change is then

[TABLE]

We now invoke a mean-field approximation to evaluate the currents (11). In this we neglect fluctuations of $\Delta H_{ij}^{\alpha\gamma}$, i.e. we replace this quantity by its average. This can be justified by thinking about a high-dimensional limit, where the number of nearest neighbour sites $z$ of any lattice site is large enough for the local ‘fields’ appearing in $\Delta H_{ij}^{\alpha\gamma}$ to average out fluctuations. In the same spirit we also drop the restrictions on the sums in $\Delta H_{ij}^{\alpha\gamma}$, which only changes the local fields by a relative amount $1/z$. It then only remains to perform the average of the kinetic prefactors, which with a mean-field decoupling becomes $\langle n_{i}^{\alpha}n_{j}^{\gamma}\rangle\approx p_{i}^{\alpha}p_{j}^{\gamma}$. Inserting the resulting approximation for the currents into (10), we obtain the mean-field kinetic equations

[TABLE]

where – within the mean-field approximation –

[TABLE]

These kinetic equations are – in spite of the rather different method of derivation – consistent with the mean-field free energy (3) in the sense that they always decrease it over time, $dF/dt\leq 0$. (This consistency is what requires the approximation step we have taken above, of dropping the restrictions on the sums defining $\Delta H_{ij}^{\alpha\gamma}$.60) This is as one would expect in a closed system, where there are no currents crossing boundaries. 53 We defer the derivation of the result $dF/dt$ to Appendix A, which generalizes similar derivations discussed elsewhere in the literature53 to include the case of direct particle-particle swapping.

From the fact that $dF/dt\leq 0$ it follows that the dynamics leads to a state which minimizes the mean-field free energy. This final state may be the ground state (global minimum) or a metastable state (local minimum). The monotonic decrease of the free energy also implies that the mean-field kinetic equations cannot describe nucleation events. Capturing these would require introducing fluctuations, e.g. by adding Langevin noise to our deterministic eqns (16). (This is further discussed in Section 6.)

4 Early-time spinodal dynamics

In this section, we will present a linearized version of our theory: it describes the growth of small fluctuations around an initial homogeneous state (via spinodal decomposition), within the framework of our mean-field kinetic equations. It will be shown that the maximum spinodal growth rates can be expressed in terms of only three moments of the polydisperse distribution. More importantly, we will use the result for the spinodal growth rates to test Warren’s two-stage hypothesis.

We begin by considering a homogeneous system of overall composition described by a list of densities: $\left\{p^{\alpha}\mid\alpha=1,\dots,M\right\}$. The system is perturbed by small fluctuations of the densities:

[TABLE]

where $\delta^{\alpha}_{i}\ll 1$. As shown in Appendix B, linearization of (16) leads to the following equation:

[TABLE]

where we define the homogeneous mobilities as

[TABLE]

We have also introduced the local chemical potentials $\mu_{i}^{\alpha}=\partial F/\partial p_{i}^{\alpha}$. These are given explicitly by

[TABLE]

for $\alpha=1,\dots,M$, while $\mu_{i}^{0}=0$. These expressions are derived from the free energy expression (3) with the explicit substitution $p_{i}^{0}=1-\sum_{\alpha}p_{i}^{\alpha}$. Finally in eqn (19) we use the discrete Laplacian $\Delta_{\rm d}$, which is defined by

[TABLE]

for any site-dependent quantity $g$.

Note that, so far, the theory applies for completely generic attempt rates $w^{\alpha\gamma}$. Furthermore, if these attempt rates are set such that the right-hand side of eqn (19) contains only the particle-vacancy term (i.e. no direct interchanges of particles are allowed), and one considers only the $M=2$ case, then eqn (19) has the same form as eqn (27) in Plapp and Gouyet’s work.53 The only difference is that their expression for the homogeneous mobility is different from (20), as Plapp and Gouyet used Arrhenius jump rates. Their choice reflects the physical assumption for alloys that the elementary particle moves have an energy barrier effectively equivalent to removing a particle from the system. For the colloidal case Glauber rates are rather more plausible.

To solve the linearized mean-field equations, one exploits that a homogeneous system is invariant under translation with respect to the lattice vectors. Solutions are therefore superpositions of time-dependent Fourier modes

[TABLE]

(as will be clear, the $i$ in the exponents refers to the imaginary unit $i\equiv\sqrt{-1}$). Here $\mathbf{k}$ is the fluctuation wave vector and $\mathbf{x}_{j}$ is the position vector in real space of lattice site $j$. Moreover, $\omega$ is the growth rate of the mode and $\delta p^{\alpha}$ indicates the amplitude of the fluctuation associated with species $\alpha$. By inserting eqn (23) into (19) one finds an eigenvalue equation with eigenvalue $\omega$, with the $\delta p^{\alpha}$ being the components of the corresponding $M$-dimensional eigenvector. [See eqn (49) in Appendix B.] Thus we have a stability spectrum with $M$ branches. A branch can be stable ($\omega$ is negative for all wave vectors) or unstable, with $\omega$ being positive in some range of $|\mathbf{k}|$, typically for small $|\mathbf{k}|$. We will be interested in the maximum growth rate $\omega_{\rm max}$ over all branches and wave vectors, which identifies the dominant growing fluctuation mode. Outside of the spinodal region this maximum growth rate becomes zero because the system is stable to all small fluctuations there.

To be able to evaluate the maximum growth rate we need to make specific assumptions about the attempt rates $w^{\alpha\gamma}$. We will set $w^{\alpha 0}=w_{0}$ and $w^{\alpha\beta}=w_{s}$, for any $\alpha\neq 0$ and $\beta\neq 0$, where $w_{0}$ and $w_{s}$ are constant attempt rates associated with particle-vacancy and particle-particle exchanges, respectively. (The ‘[math]’ subscript is for vacancy, and the ‘$s$’ is for swapping.) In principle, one could imagine that $w^{\alpha\beta}$ might not be the same for all pairs $\alpha$ and $\beta$, or that it depends on the temperature. However, our simple choice for the values of $w^{\alpha\gamma}$ is enough to distinguish between particle-vacancy and particle-particle kinetic mechanisms. Moreover, we will see later that the dependence on temperature would be irrelevant for our purposes. Also, Plapp and Gouyet say in ref. 53 that, in their case with particle-vacancy dynamics only (where $M=2$ and Arrhenius rates are used), numerical results indicate that qualitative phase-separation behaviour is unaffected by the attempt rate ratio $w^{A0}/w^{B0}$ as long as it is not too far from unity.

For the above choice of $w^{\alpha\gamma}$, we show in Appendix B, that the (largest branch of) growth rates can be expressed as

[TABLE]

where $A(\mathbf{k})$ is essentially the Fourier transform of the Laplacian. Equation (24) is valid for any spatial dimension given the appropriate expression for $A(\mathbf{k})$, with e.g. for a two-dimensional lattice $A(\mathbf{k})=-4\sin^{2}\left(k_{x}a/2\right)-4\sin^{2}\left(k_{y}a/2\right)$. The moment densities in the above expressions are, in the discrete representation, $\rho_{n}=\sum\limits_{\alpha}\sigma_{\alpha}^{n}p^{\alpha}$.

In the monodisperse limit ($M=1$), expression (24) does not depend on $w_{s}$. This of course should be so as swaps between identical particles do not change the system configuration and hence cannot contribute to the dynamics. To find the maximum growth rate, one needs to maximize $\omega$ with respect to $\mathbf{k}$. As in the regime of interest the maximum occurs for small wavevectors, this can be done by expanding $\omega$ as a function of $\mathbf{k}$ around $|\mathbf{k}|=0$ up to the fourth order in $|\mathbf{k}|$, and then maximising the resulting expression analytically. One can also carry out numerically the maximization of the full $\omega(\mathbf{k})$, with essentially indistinguishable results (see Figures 4a and 7a), but we use the expansion procedure to obtain a closed-form expression for $\omega_{\rm max}$. (This is nonetheless too long to be displayed here, though.)

Because of the moment structure of (24), our linear theory can be applied to a fully polydisperse (i.e. $M\rightarrow\infty$) system, using an experimentally-reasonable distribution like the Schulz-Gamma form61

[TABLE]

In the following we set the mean interaction strength $\bar{\sigma}=1$ as we did in the binary case. The parameter $g$ controls the polydispersity of the distribution, which is given by $1/\sqrt{g+1}$. This means that e.g. the choice $g=15$ produces a standard deviation of $\sigma$ that is $25\%$ of the mean. With these choices, the moments appearing in the growth rates (24) can be expressed in terms of the density as $\rho_{1}=\bar{\sigma}\rho=\rho$, $\rho_{2}=\bar{\sigma}^{2}\rho[1+1/(g+1)]=\rho(g+2)/(g+1)$. Note that because only moments up to second order appear in our mean-field spinodal rates, other distributions with the same mean and variance would give identical results.

Let us now see what our linear theory says about Warren’s scenario as applied to the spinodal dynamics. Fig. 4a shows the maximum growth rate as a function of the overall density, for reasonably dense systems. Here and in the following we fix the overall timescale by setting $w_{0}=1$. The vertical lines indicate the upper limits of the annealed and quenched spinodal regions, respectively. As expected, the maximum growth rate becomes zero beyond the annealed spinodal, where the system is stable to all density fluctuations. More remarkable is that the maximum growth rate increases only very slowly as density is decreased below this point, and only begins to rise substantially at the quenched spinodal. This is exactly in line with what would be expected from Warren’s two-stage hypothesis: inside the quenched spinodal region, the system has fast (stage 1) spinodal dynamics driven by the instability with respect to density fluctuations. (This corresponds to region R1 in Table 1.) Outside the quenched spinodal, on the other hand, there is no spinodal decomposition in stage 1 (corresponding to regions R2/R3) and the spinodal dynamics is produced by the much slower growth of composition fluctuations in stage 2. To the extent that stage 2 dynamics, which involves local fractionation, is slow compared to stage 1, $\omega_{\rm max}$ should therefore be small between quenched and annealed spinodals, as compared to its values inside the quenched spinodal region. This is exactly what we find.

Graphically, the above reasoning means that $\omega_{\rm max}(\rho)$ should have a kink at the quenched spinodal, where it crosses over from small (stage 2) to large (stage 1) values. The situation in Fig. 4a is quite close to such an ideal two-stage scenario. The kink can be seen more clearly by looking at the second derivative of $\omega_{\rm max}$ with respect to $\rho$, which would be large around a kink. Fig. 4b shows that this second derivative does indeed have a maximum, and this is positioned close to, if not exactly at, the quenched spinodal density. Note that this happens even though our calculation of $\omega_{\rm max}$ did at no point involve any restriction to quenched dynamics, i.e. fixed composition. In other words, the full, unrestricted dynamics of the system nonetheless ‘feels’ the presence of the quenched spinodal. This provides strong support for Warren’s two-stage scenario.

Since the physics of Warren’s two-stage scenario requires the system to be dense, we expect it to be break down at lower densities. To check this, we first investigate the behaviour of $\omega_{\rm max}\left(\rho\right)$ around the lower end of the spinodal region. In Fig. 5, it is clear that the second derivative does not have any signatures around the (lower) quenched spinodal density, instead increasing smoothly towards the annealed spinodal. Secondly, returning to the kink in $\omega_{\rm max}$ at the upper quenched spinodal density, we can consider its dependence on temperature. At higher temperatures, the upper spinodal densities become lower, so that the two-stage scenario should be less pertinent. To check this, we show in Fig. 6 the density where the maximum in the second derivative of $\omega_{\rm max}$ for a range of temperatures. We can see that the density at the maximum moves away from the quenched spinodal curve and towards the annealed spinodal as $T$ increases, as expected. More usefully, we can read off from the figure that the two-stage scenario gives a good account of the position of the kink for densities above $\rho\approx 0.9$, so for the given polydispersity this is the threshold where the system becomes sufficiently dense to make fractionation slow.

We now perform a second test that detects directly whether the behaviour we are seeing – namely the near-kink in $\omega_{\rm max}(\rho)$ – is in fact due to fractionation being slow. We do this by turning on direct particle swaps, using a nonzero swap attemp rate $w_{s}$. Fractionation is then possible even in dense systems, without relying on mediation by rare vacancies. The signatures of the two-stage scenario that we have found should therefore disappear as we increase $w_{s}$. Indeed, Fig. 7 shows that for $w_{s}=0.5$, $\omega_{\rm max}$ increases smoothly as density is decreased from the annealed spinodal, rather than remaining small until the quenched spinodal. Likewise the second derivative of $\omega_{\rm max}\left(\rho\right)$ is now featureless around the quenched spinodal and just increases gradually towards its value at the annealed spinodal. Essentially, this is evidence that the two-stage scenario has been destroyed.

So far we have focussed on a system with fixed polydispersity of 25%. As Warren’s argument does not rely on specific features of the $\sigma$-distribution, one would however expect qualitatively similar results also for other polydispersities. To assess this quantitatively one needs to account for the fact that the separation in density between quenched and annealed spinodals grows with polydispersity. This can be done by considering the density difference between the maximum of $\omega_{\rm max}^{\prime\prime}$ and the quenched spinodal, normalized by the difference between the annealed and quenched spinodal densities. When this ratio is $\ll 1$, the kink of $\omega_{\rm max}$ is close to the quenched spinodal as the two-stage scenario predicts. Upon varying temperature, density and polydispsersity we find (data not shown) that the ratio is indeed dependent mostly on density and largely independent of polydispersity, becoming small at high densities as it should.

Summarizing the discussion in this section thus far, our mean-field theory for the spinodal dynamics of polydisperse colloids provides strong support for Warren’s two-stage hypothesis, in the appropriate regime of high densities. It is worth noting that previous support for the two-stage scenario, both in Warren’s original paper48 and in the study by Pagonabarraga and Cates,51 was obtained in the context of polymers, whereas here we have polydispersity in the interaction strength in a context that is more easily connected with the physics of colloids. This indicates that Warren’s scenario may be a general feature of the non-equilibrium dynamics of dense fluid mixtures.

We next briefly discuss the temperature dependence of the maximum growth rates, returning to the physical value $w_{s}=0$. The maximum growth rates are then directly proportional to $w_{0}$. So far we have taken constant $w_{0}=1$, though the results of Fig. 6 would remain the same for any temperature-dependent $w_{0}$: this is because we always consider the second derivative $\omega_{\rm max}^{\prime\prime}$ with respect to density, at fixed temperature. With constant $w_{0}$, we find that the typical maximum growth rates (say, in the middle of the spinodal region) increase towards low temperatures as $\sim 1/T$. This temperature dependence comes from the use of Glauber rates: linearizing the factors of $\exp(\Delta H/T)$, with energy changes from particle moves that are of temperature-independent magnitude, gives a factor of $1/T$ in all growth rates. For quantitative modelling of spinodal growth rates one would therefore want to choose a $w_{0}$ that goes to zero as $T\to 0$, e.g. $w_{0}\sim T$, which would be consistent with the picture of an underlying diffusive dynamics causing attempted particle moves. This would then give constant maximum growth rates for low $T$.

We conclude this section with a discussion of the amplitudes $\delta p^{\alpha}$ of the mode that grows fastest in the spinodal dynamics, with rate $\omega_{\rm max}$. These amplitudes define the relative strength of density fluctuations for each species. As a vector in density space, they identify the early-time (non-equilibrium) fractionation direction, or simply spinodal direction. As discussed above, combining eqns (19) and (23) leads to an eigenvalue equation with eigenvalue $\omega$, and the $\delta p^{\alpha}$ are the $M$ components of the eigenvector corresponding to $\omega_{\rm max}$. This eigenvector can be calculated by solving the eigenvalue equation numerically, or one can derive a closed-form expression in terms of $\omega$, $\rho$, and the attempt rates. We will denote the angle between the spinodal direction and the dilution line by $\theta_{\rm max}$. This quantity is of interest because in an ideal two-stage scenario, it should be zero inside the quenched spinodal region, where the initial dynamics corresponds to stage 1, i.e. pure density fluctuations; it should then rise as one moves from the quenched spinodal density to the annealed one, where the spinodal dynamics corresponds to stage 2. Our numerical data for $\theta_{\rm max}(\rho)$ (not shown) follow this scenario fairly closely in dense systems, providing further support for Warren’s hypothesis. As our system is not ideal in the sense that fractionation is not infinitely slow, we find inside the quenched spinodal region a $\theta_{\rm max}$ that is constant but not quite zero, indicating that even in stage 1 the dominant spinodal mode contains a fractionating component in addition to pure density fluctuations. Likewise the transition to larger values of $\theta_{\rm max}$ is not a sharp kink but a crossover, though importantly this remains located around the quenched spinodal density. As in the case of $\omega_{\rm max}$, we have tested that when direct particle swaps ($w_{s}>0$) are introduced, this crossover disappears and is replaced by a featureless increase of $\theta_{\rm max}$ across the entire annealed spinodal region. This confirms that the crossover in the physical system ($w_{s}=0$) arises from fractionation being a slow process.

5 Beyond the spinodal regime

In order to see what happens after the end of the spinodal regime, and hence investigate the full phase-separation dynamics, we integrate our mean-field kinetic eqns (16) numerically using a forward Euler method, with periodic boundary conditions. The method returns the evolution of the local densities for all species. For ease of visualization and reduction of computation time, spatial dimension $D=2$ was chosen. Firstly, we examine what happens for binary mixtures ($M=2$), and then extend the analysis to $M>2$. We set the jump attempt rates as before, i.e. $w^{\alpha 0}=w_{0}$ and $w^{\alpha\beta}=w_{s}$, for any $\alpha\neq 0$ and $\beta\neq 0$. More specifically we concentrate on the physical setting $w_{s}=0$, unless otherwise stated. To find the evolution of a system whose overall composition is given by a list of densities $\left\{p^{\alpha}\mid\alpha=1,\dots,M\right\}$, we firstly created an initial homogeneous state defined by $p^{\alpha}_{i}=p^{\alpha}$ for all sites $i$. To trigger the phase separation, we added small fluctuations to the initial state of each species by generating $L^{D}$ random numbers, normally distributed, with mean zero and standard deviation $1\%$. We then subtracted the average of these random numbers, to ensure that the overall density of every species remains unchanged. The time step used was limited by the numerical stability of the algorithm, and typically equal to $0.1$; we checked that this value is small enough to give us effectively the solution of the continuous-time equations by running the numerics for exactly the same initial configuration with a time step smaller by a factor of 5 and verifying that the results were virtually the same. When we used $50\times 50$ lattice sites and $M=2$, for instance, our program (written in C) performed $8,000$ time steps in approximately $100$ seconds (running at 2.6 GHz processor speed), which in many instances was enough time to grow relatively large domains.

5.1 Binary fluids

Fig. 8 shows snapshots of the phase-separation dynamics for $M=2$, with $\sigma_{A}=1+d$ and $\sigma_{B}=1-d$. For these binary case numerics we chose $d=0.25$. Thus the first and second moment densities (normalized by $\rho$), i.e. $\rho_{1}/\rho$ and $\rho_{2}/\rho$, are the same as in the continuous distribution analysis presented in Section 4, where the Schulz parameter was $g=15$; also, the same temperature $T=0.3$ was used. Therefore the spinodal growth rates are also the same as in the fully polydisperse case and we expect, at a given density $\rho$, to see a very similar initial dynamics. We chose $\rho=0.82$, which places the system within the quenched spinodal region. To determine the colour of a site $i$, we used a colour scheme in which the colours red, green, and blue are blended together. The intensity of each of these colours at a given site varies from [math] to $1$. In our scheme, red, green, and blue intensities are given exactly by $1-p_{i}^{A},p_{i}^{0},1-p_{i}^{B}$, respectively. (Remember the notation for the local concentration of vacancies, i.e. $p_{i}^{0}=1-p_{i}^{A}-p_{i}^{B}$.) This leads to the colour key shown in the top-left part of Fig. 8. It is plotted in triangular colour space in $(p^{A},p^{B},1-p^{A}-p^{B})$, dropping the site index $i$. For example, if the concentration of particles of species $A$ at one site is high (low), and the concentration of particles of species $B$ at the same site is low (high), then the site colour will tend towards blue (red); if the concentrations of all species are all low, then the site colour will tend to white.

The snapshots in Fig. 8 show the growth of lighter regions of the system that represent gas bubbles. These bubbles are surrounded by a $B$-rich (red) interface separating them from an $A$-rich (blue) continuous liquid phase. Walking from the centre of a gas bubble along an arbitrary direction, one therefore initially sees low concentrations of both particle types, then a high concentration of the particle species with the smallest $\sigma$, i.e. $B$-particles, and eventually one reaches the bulk liquid that contains predominantly $A$-particles, for which $\sigma$ is the largest. (Video showing the full-time evolution of the phase separation process is provided as an animation in the ESI.†)

Of course as we decrease $\rho$ our numerical results show larger vapour fractions. Changing $\rho$ can also lead to more complicated morphologies such as in Fig. 1, which was generated using $T=0.4$, $p^{A}=p^{B}=0.33$, $d=0.3$, $w_{s}=0$, $L=128$, with the last snapshot taken at $t=6000$. Because the value of $\rho$ here is closer to the critical density one observes bicontinuous domains of gas and liquid. The intuition is the same as in a one-species Ising lattice-gas system, where because of the particle-hole (vacancy) symmetry, in a system at the critical density neither gas or liquid can ‘win’ to form bubbles or droplets. Instead, finger-like bicontinuous structures are formed. Moving away from the critical point, one expects these to survive for a certain time until the system ‘notices’ which phase is going to be the majority phase, and forms bubbles or droplets of the minority phase.

In Fig. 9 we introduce a different representation of the time evolution that will prove to be quite revealing. We show two-dimensional density histograms, in species density space, i.e. in the $p^{A}_{i}$–$p^{B}_{i}$ plane that contains all the possible density combinations an arbitrary site $i$ can have; the physically accessible region in this plane is a triangle bounded by $p^{A}_{i}\geq 0$, $p^{B}_{i}\geq 0$ and $p^{A}_{i}+p^{B}_{i}\leq 1$). What the histogram counts is the number of lattice sites $i$ that have species densities $(p^{A}_{i},p^{B}_{i})$ inside each two-dimensional bin. The results were normalized by the total number of lattice sites $L^{D}$. Such a density histogram can then be viewed as a dynamic analogue of an equilibrium phase diagram as sketched in Fig. 2 above. In a density histogram, the parent phase lies on the dilution line (shown dashed in Fig. 9) as before. The low- and high-density daughter phases calculated from the equilibrium phase diagram lie off this dilution line, with the parent on the connecting tieline; the latter defines the equilibrium fractionation direction. At the temperature we are considering, this equilibrium fractionation direction deviates only slightly from the dilution line. From the histograms, we can clearly see different dynamical regimes: initially, the histogram spreads linearly from the parent along the spinodal direction as expected for spinodal dynamics. As nonlinear effects kick in, a curved path of compositions connecting a gas and liquid phase is then formed. This clearly delineated ‘arc’ contains the compositions of the different parts of the system: as one moves in space from a gas bubble into the bulk liquid, one passes through a series of compositions within the interface between these two phases.

Beyond this generic structure, there are several interesting observations we can make from Fig. 9. The density histograms reveal that gas-liquid interfaces are strongly fractionated, with the arc being well away from the dilution line, in the $B$-rich part of the density plane. Physically, the reason is that $B$-particles have smaller $\sigma_{B}$ and hence interact more weakly; they therefore pay a smaller energy penalty for sitting at an interface, where they have fewer neighbouring particles. Interfaces also have a well-defined sequence of density combinations as can be seen from the fact that the gas-liquid arc is quite narrow.

For further analysis it is useful to switch to the two-dimensional density histogram representation in Fig. 10. This shows the same data as in Fig. 9c but now seen from the top, with different heights corresponding to different colours. The peak in the high-density region of the histogram is also marked; this gives the majority composition of the bulk liquid at that time instant. It is interesting to observe that this peak, having started out at the parent composition and ‘walked’ along the spinodal direction initially, does not subsequently move straight away towards the liquid equilibrium composition. Instead, the liquid phase composition stays away from its equilibrium optimum for a long time. In fact, for the dynamics shown in Fig. 10 the liquid peak moves away from its equilibrium point for a long transient period beyond the initial spinodal decomposition dynamics. This arises because the gas-liquid interfaces are strongly enriched in $B$-particles, leaving an unusually $A$-rich bulk liquid. Of course at very long times the equilibrium prediction and the dynamics must eventually agree, and we have verified that the density histogram peaks then indeed centre on the calculated equilibrium compositions while the arc with the interface compositions contains only a small ($\sim W/L$, where $W$ is the interfacial width) fraction of probability. (See also the evolution of two-dimensional density histograms in the animation in the ESI.†)

A final, intriguing feature of the dynamics we observe is the inhomogeneity of the bulk liquid: in the last two snapshots in Fig. 8 one can clearly see well-defined liquid regions that are unusually enriched in $B$-particles. In the density histograms of Figs. 9c and 10 and also Fig. 11a, which shows results at a higher temperature, these regions manifest themselves as an ‘arm’ at high density that is quite distinct from the arc arising from gas-liquid interfaces.

Looking at Fig. 8 carefully, one notices that the origin of the $B$-rich liquid regions lies in the evaporation of gas bubbles. As these shrink, so do their interfaces, eventually forming dense patches. Because the interfaces are strongly fractionated, these dense patches are strongly enriched in $B$-particles. While the density of the patches can rapidly equilibrate to the bulk liquid—as shown by the fact that the arm almost coincides with a line of constant total number density $p^{A}+p^{B}$, as seen from Figs. 9c, 10, and 11a—it requires inter-diffusion of particles to equilibrate their composition. Hence the composition heterogeneities formed by these patchese are unusually long-lived. We thus have here, in the long-time dynamics, another striking manifestation of Warren’s hypothesis, in its general form which says that equilibration of composition is slow in dense systems. In order to check this interpretation we turned on $w_{s}$; as expected, the arm then disappears, and in the real-space images the liquid is clearly homogeneous (Fig. 11b).

5.2 Multicomponent fluids

Finally, we show now that a similar two-dimensional density histogram analysis can be performed even when one considers arbitrary $M$. A full density histogram in $(p_{i}^{1},\ldots,p_{i}^{M})$-density space would be $M$-dimensional, so one needs to project down to a manageable low-dimensional representation. The obvious choice for the new histogram axes are (the local versions of) low-order moment densities, which can be used for arbitrary $M$. We choose the local analogues of $\bar{\sigma}\rho_{0}-\rho_{1}$ and $\rho_{0}$. The reason for this choice is that in the bidisperse case ($M=2$),

[TABLE]

while $\rho_{0}=p^{A}+p^{B}$, so that the new histogram axes are just a rotated version (by 45∘) of the ones we have used so far.

One can now use these quantities to plot two-dimensional density histograms for an arbitrary polydisperse system. It turns out that all the features that had been previously observed in the $(p^{A},p^{B})$-plane density histogram (i.e.  the clearly delineated curved arc of interfaces compositions, the ‘arm’, etc.) remain qualitatively identical in this new coordinate system. Fig. 12 shows our results for systems with two, three, and four species. For $M=3$, we used $\sigma\in\{1-d,1,1+d\}$, with relative densities (composition) given by $\{\frac{1}{4},\frac{1}{2},\frac{1}{4}\}$, whereas for $M=4$ we chose $\sigma\in\{1-d,1-d/2,1+d/2,1+d\}$ with composition $\{\frac{1}{6},\frac{1}{3},\frac{1}{3},\frac{1}{6}\}$. We then identified values of $d$ for each $M$ that give the same set of moment densities ($\rho_{0}$, $\rho_{1}$, and $\rho_{2}$) as we had in Section 5.1. (Thus the spinodal growth rates are also all the same and hence results at the same $t$ are comparable—but similar results were found using different parameter sets chosen by the same method.) Even though these histograms are now projections of $M$-dimensional histograms to two dimensions, we still get thin arcs between daughter phases, showing that there is still a well-defined sequence of compositions in the interfaces. We have performed checks over a larger parameter range, where we find that long-lived composition heterogeneities also appear in situations where one expects a pronounced slowing down of fractionation, exactly as we saw for $M=2$; correspondingly, they disappear (data not shown) when direct particle swaps are turned on.

The fact that the density histograms for different $M$ in Fig. 12 are so similar is quite remarkable. This similarity is surprising because the mean-field dynamical equations (16) do not in general reduce to closed dynamical equations for local moments like $\rho_{0}$ and $\rho_{1}$, due to nonlinear dependences on particle size $\sigma$ in the Glauber rates. Nonetheless our numerical results suggest that degrees of freedom not captured by these moments only influence the dynamics weakly, so that they could provide a useful way of thinking about the dynamics even for $M\to\infty$.

6 Conclusions

In this work we investigated the dynamics of how colloidal polydisperse systems phase-separate, by introducing new kinetic equations based on the Polydisperse Lattice-Gas model.35 As a baseline for our mean-field approach to the dynamics we calculated the mean-field equilibrium phase diagram of the model, including cloud and shadow curves as well as spinodals. To test Warren’s two-stage scenario,48 we obtained both the annealed and quenched versions of the phase diagram, the latter referring to a system that can only change its density but has fixed composition. We analysed the linearized dynamical equations to understand the dynamics of spinodal decomposition, and found clear evidence in support of Warren’s proposal. For the late-stage dynamics, we introduced a two-dimensional density histogram method that allows fractionation effects in the phase-separation dynamics to be clearly visualized. This revealed strongly fractionated interfaces between gas and liquid. It also helped us to detect the existence of long-lived composition heterogeneities in the bulk liquid, which are a further manifestation of the fact that fractionation dynamics is slow as suggested by Warren. This prediction may be amenable to experimental verification in dense colloidal systems. The whole analysis was performed for an arbitrary number of particle species, although much of it was presented in binary mixture context for the sake of simplicity and ease of visualization.

Our main assumptions were that the dynamics can be described by a kinetic lattice model, and that a mean-field approximation is at least qualitatively accurate. One could ask whether non-mean-field effects in the equilibrium phase diagram of the PLG might cause significant changes in the phase-separation dynamics. This could be tested by deploying higher order approximations beyond our dynamical mean-field theory (DMFT), such as the Path Probability Method (PPM).62 One might expect that the PPM would not necessarily lead to new qualitative outcomes for the analysis presented here, since no differences were observed in a previous comparison between PPM and DMFT, though in the somewhat different context of relaxation dynamics in porous materials.63 Direct Kinetic Monte Carlo (KMC) simulations could be used to directly probe the dynamics of the PLG model and so assess the effects of our mean-field approximation. Such simulations were performed in ref. 64 for similar systems consisting of two species plus vacancies, where $p^{A}=p^{B}$ and the Hamiltonian is given by

[TABLE]

with $\alpha=A,B$. This setting includes the PLG Hamiltonian (for $M=2$) given by eqn (1) but allows more generic interactions $\epsilon_{\alpha\beta}$. Also using a mean-field phase diagram as a reference, the authors of ref. 64 investigated segregation kinetics with particle-particle swaps forbidden (corresponding to $w_{s}=0$ in our notation) or highly energetically suppressed. They focussed on the existence of different domain growth morphologies in various parameter regimes. However, their choices for the interaction strengths $\epsilon_{\alpha\beta}$ were rather different from ours. Our separable assignment $\epsilon_{\alpha\beta}=\sigma_{\alpha}\sigma_{\beta}$ makes liquid-gas phase separation the dominant physical process, for which we then study the effects of polydispersity. On the other hand, the authors of ref. 64 look at much less attractive (and even repulsive) $AB$ interactions, always with $\epsilon_{AA}=\epsilon_{BB}$. This leads to distinct physical effects, including the condensation of vacancies at interfaces between $A$- and $B$-rich phases.

KMC simulations could also be used to fit our model parameters to experimental systems. Ref. 65 explains that one way of obtaining a physically reasonable value for the particle-vacancy jump attempt rate $w_{0}$ is by comparing between estimated ($\hat{D}$) and experimental ($D_{s}$) self-diffusion coefficients. If $\hat{D}$ is obtained from KMC simulations in dimensionless units, the desired correspondence would be $D_{s}=\hat{D}a^{2}w_{0}$, where $a$ is the lattice spacing. Therefore, fixing the value of $a$, an estimate of $w_{0}$ can be obtained. (In principle, one could try to develop a similar scheme to obtain a value for $w_{s}$.)

Coming to the limitations of the PLG model itself, one obvious shortcoming is the fact that a lattice model cannot capture gradual increases in density that are possible in an off-lattice setting, where even in an already fairly dense system collective motion can reduce the typical distance between particles. In a lattice model, on the other hand, density increases have to come from localized filling in of vacancies, or equivalently vacancies diffusing away. As we explain below the presence of vacancies also enables fractionation. Because both relaxation of density and of composition are then tied to the presence of vacancies, it is clear that the latter cannot become arbitrarily slow compared to the former. In a continuum model, however, one would expect that large pockets of free volume that are required for the interchange of particles of different species become rare quickly at high density while collective motion to relax density still remains possible. Therefore, if anything, we expect the slowing down of fractionation compared to density relaxation at high densities to be more pronounced than in our lattice model. Hence the effects we observe should be stronger in more realistic, off-lattice models, getting closer to an ideal two-stage scenario for the phase-separation kinetics.

We illustrate how vacancies enable fractionation in Fig. 13, which shows a sketch of a vacancy-mediated interchange between two particles of different species. A sequence of 4 particle moves within the small lattice portion shown ($2\times 2$ lattice sites) is sufficient to interchange the blue and red particle in the lattice row above the initial vacancy. (Note that in this sequence the two blue particles have also interchanged their positions, but as they are indistinguishable this is immaterial.) Therefore, as fractionation requires interdiffusion of particles of different species, it will be able to proceed locally as long as vacancies are present. As a particle swap can be accomplished with a moderate number of elementary particle moves, fractionation cannot become arbitrarily slow compared to density relaxation–which locally requires a single particle move–as claimed above. This is the likely reason why we do not see an ‘ideal’ two-stage scenario, including e.g. the fact that at early times the spinodal dynamics within the quenched spinodal region already has a component of fractionation rather than consisting of the growth of pure density fluctuations.

Of course off-lattice models at the particle level are difficult to deal with theoretically, so before studying phase-separation dynamics one would aim to derive from them approximate representations in terms of a continuum field theory such as the so-called Model B or rather a polydisperse variant thereof. 66, 67 For a monodisperse system and with stochastic fluctuations neglected, this would give an equation of the (generalized) Cahn–Hilliard type68

[TABLE]

where $F$ is the free energy functional, $\rho\left(\mathbf{r},t\right)$ is the continuous density field, and $\mathscr{M}\left(\rho\left(\mathbf{r},t\right)\right)$ is a density-dependent mobility; within a simple approximation scheme69 one finds $\mathscr{M}\left(\rho\left(\mathbf{r},t\right)\right)=\rho\left(\mathbf{r},t\right)$.

In the polydisperse case, the analogue of eqn (28) involves a mobility matrix with $M\times M$ entries. (In fact, our kinetic eqn (16) can be cast in the form of a discrete Cahn–Hilliard equation, generalized to polydisperse fluids and inhomogeneous mobilities; in the linear version, eqn (19), the generalized Cahn–Hilliard form is clear.) Determining such an entire mobility matrix from an off-lattice model is a considerable challenge, also because the results might differ for size versus interaction polydispersity (similarly to the differences found in the mobility coefficients between systems with length and chemical polydispersity, in ref. 51). One route to deriving such Cahn-Hilliard-like equations would be dynamical density functional theory (DDFT); it should be possible to adapt this to a fully polydisperse scenario and then use it to test Warren’s two-stage scenario.70, 71

A further theoretical challenge is the incorporation of stochastic effects, in order to be able to describe nucleation and growth dynamics. The simplest way of achieving this would be to add Langevin noise to either our lattice dynamics (16) or a continuum description like eqn (28). This however requires a consistent, quantitative way of adding multiplicative Langevin noise to the deterministic equations. The need for a quantitative accurate stochastic description lies in the fact that one wants to compare equilibration timescales between systems placed in different regions of the phase diagram, which means one has to assess the competition between nucleation and growth dynamics on the one hand and spinodal decomposition on the other. For instance, for systems placed within R2 (see Table 1) one would like to compare the timescale for nucleation and growth in stage 1 (which should have a fast intrinsic time scale because it does not require fractionation, but could be slowed down by large nucleation barriers) to spinodal dynamics in stage 2 (intrinsically slow because it requires fractionation, but not slowed down by activation barriers). A model being able to describe quantitatively both types of phase ordering dynamics would clearly be valuable here. This could potentially be obtained from the lattice-based theory presented here via systematic coarse-graining. 72 Alternatively, fluctuating hydrodynamics73, 74, 75, 76 is a possible avenue for deriving models incorporating stochasticity, though whether this can be implemented for polydisperse dynamics and would give a reasonable quantitative account of the physics remains an open question.

In addition to stochastic effects we have also neglected hydrodynamic interactions due to the solvent. This is commonly done throughout the literature on polydisperse dynamics, but in the case of dense systems, which are especially pertinent in the context of the two-stage phase-ordering scenario, they may play an important role.70

Our work in this paper could be extended to investigate the relaxation dynamics of polydisperse fluids in porous materials, generalising previous work done by Peter Monson and others. 77, 78, 79 In fact, although the dynamical mean-field theory calculations developed in those papers are restricted to $M=1$ and $M=2$, they are very similar to the ones presented here. For application to porous materials, we would mainly need to adapt our approach to cover given pore geometries. In this scenario, a natural question would be ‘what is the impact of slow fractionation on the relaxation dynamics in porous materials?’. This is potentially one of the simplest future research topics to pursue if one uses our current theory as a starting point.

Finally, we point out that one could use the framework developed here to investigate yet other problems, e.g. the dynamics of polydisperse wetting or the phase-separation behaviour of polydisperse systems of active particles. It would also be interesting to develop perturbative or scaling approaches for the dynamics of systems with small polydispersity, and to consider alternative ways of obtaining mean-field approximations, for example by explicitly taking the limit of high spatial dimension or considering a Kac-like setup with long-range interactions.80, 81, 36, 82

7 Acknowledgments

PdC acknowledges financial support from CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico – Brazil. PS acknowledges the stimulating research environment provided by the EPSRC Centre for Doctoral Training in Cross-Disciplinary Approaches to Non-Equilibrium Systems (CANES, EP/L015854/1).

Appendix A Monotonic decrease of the free energy

Here we show that in the kinetic PLG model with both particle-vacancy and direct particle swaps (with Glauber-like rates, arbitrary $M$, arbitrary overall composition, and arbitrary attempt rates) the free energy obeys $dF/dt\leq 0$. As we will see, the derivation is not as straightforward as in the case of a binary mixture without direct particle swapping, which was discussed in ref. 53. The beginning is identical, though. As the free energy $F$ is a function of all dynamical variables $p^{\alpha}_{i}$ we can write

[TABLE]

which, combined with the definitions for the chemical potential and for the current, can be rewritten as

[TABLE]

Using the fact that $J_{ij}^{\alpha}=-J_{ji}^{\alpha}$ we can write

[TABLE]

Defining a two-species current $J^{\alpha\gamma}_{ij}$ in the obvious way using

[TABLE]

transforms this further into

[TABLE]

The sum can be extended to $\alpha\geq 0$ because $\mu^{0}_{i}=0$, $\forall i$. Also note that $J^{\alpha\gamma}_{ij}$ is antisymmetric under interchange of $\alpha$ and $\gamma$. After adding the same expression with $\alpha$ and $\gamma$ swapped and using this antisymmetry, we have

[TABLE]

However, the term in square brackets can be written as $\ln(C^{\alpha}_{i}C^{\gamma}_{j})-\ln(C^{\alpha}_{j}C^{\gamma}_{i})$, where we have defined the local chemical activity $C^{\alpha}_{i}$ as $C^{\alpha}_{i}=\exp(\mu^{\alpha}_{i}/T)$. This has the same sign as $J^{\alpha\gamma}_{ij}$ because the two-species current can be written as a positive quantity multiplied by $C^{\alpha}_{i}C^{\gamma}_{j}-C^{\alpha}_{j}C^{\gamma}_{i}$. (This fact can be seen from the mean-field expression for the energy change, eqn (17), which enters the rates via $\exp(\Delta H_{ij}^{\alpha\gamma}/T)=p_{i}^{\alpha}p_{j}^{\gamma}C_{i}^{\gamma}C_{j}^{\alpha}/(p_{i}^{\gamma}p_{j}^{\alpha}C_{i}^{\alpha}C_{j}^{\gamma})$.) Therefore, each link $\langle i,j\rangle$ and each pair of $(\alpha,\gamma)$ makes a negative contribution to $dF/dt$. One can also phrase the last step in the other direction, i.e. by writing

[TABLE]

and checking that the site-dependent mobility $\mathscr{M}^{\alpha\gamma}_{ij}$ is positive. This makes sense because a particle swap is driven by how the difference in $\alpha$-chemical potential between sites $i$ and $j$ compares to the corresponding difference in $\gamma$-chemical potential. Note that while the factor $1/2$ in eqn (34) may looks somewhat unexpected, antisymmetry of the current $J^{\alpha\gamma}_{ij}$ again shows that $(\alpha,\gamma)$ and $(\gamma,\alpha)$ make the same contribution, so one could eliminate the $1/2$ again and restrict the summation to ordered pairs of species labels, say $\alpha<\gamma$.

Appendix B Derivation of growth rates

In this appendix we derive the linearized mean-field equations of motion and derive from them the expressions for the spinodal growth rates. Consider a homogeneous system of overall composition defined by a list of species densities $\left\{p^{\alpha}\mid\alpha=1,\dots,M\right\}$. The system is perturbed by small fluctuations of the densities:

[TABLE]

where $\delta^{\alpha}_{i}\ll 1$. In a linear expansion we can write eqn (17) as

[TABLE]

Now remember our kinetic eqns (16)

[TABLE]

Using a simplified notation in which we drop the superscripts $\alpha$ and $\gamma$, the expression between the curly brackets can be written as

[TABLE]

where $A_{ij}=w^{\alpha\gamma}p_{i}^{\alpha}p_{j}^{\gamma}$ and \Delta E_{ij}=\bigl{\langle}\Delta H_{ij}^{\alpha\gamma}\bigr{\rangle}. As a consequence of linearization, the arguments in the exponentials are also small; hence, each sigmoid linearizes as $1/\left(1+\exp\left(x\right)\right)=\left(1/2\right)\left(1-x/2+\dots\right)$, and then we can write (38) in linearized version as

[TABLE]

where we have used the fact that $\Delta E_{ij}=-\Delta E_{ji}$. Now eqn (16) becomes

[TABLE]

As \bigl{\langle}\Delta H_{ij}^{\alpha\gamma}\bigr{\rangle} is already linear in $\delta$’s, the last factor can be replaced by $2p^{\alpha}p^{\gamma}$. Linearizing also the first term then gives for the time evolution of the fluctuations

[TABLE]

The first part can be rewritten as

[TABLE]

where $\Delta_{\rm d}$ is the discrete Laplacian, defined for any site-dependent quantity $g$ as $\Delta_{\rm d}g_{i}=\sum\limits_{j\in\partial i}\left(g_{j}-g_{i}\right)$. For the second part, one can obtain

[TABLE]

where $z$ is the lattice coordination number and we used the fact that

[TABLE]

which can be easily derived in a few steps. Bearing in mind that summations over $\beta$ do not include $\beta=0$ and defining a homogeneous mobility $\mathscr{M}^{\alpha\gamma}=\frac{w^{\alpha\gamma}}{2T}p^{\alpha}p^{\gamma}$, we have altogether

[TABLE]

Note that the second term in the curly braces here comes from the identity $\delta_{i}^{0}=-\sum_{\beta}\delta_{i}^{\beta}$, which is a consequence of the local hard core constraint (2).

Alternatively, using the local chemical potential (21) linearized about the homogeneous state,

[TABLE]

We can now proceed to solve the linearized dynamical equations of motion. As a homogeneous system is invariant under translation with respect to the lattice vectors, solutions of the linearized equations are superpositions of solutions of this form

[TABLE]

where it will be clear from the context that the letter $i$ here refers to the imaginary unit $i\equiv\sqrt{-1}$. Also, $\mathbf{k}$ is the fluctuation wave vector and $\mathbf{x}_{j}$ is the position in real space of a lattice site $j$. Moreover, $\omega$ is the spinodal decomposition growth rate and $\delta p^{\alpha}$ is the amplitude of the spinodal fluctuation associated with species $\alpha$. Inserting this form of the solution into our eqn (46) will give us an equation for $\omega$. In order to have a useful relation, however, the application of the discrete Laplacian operator $\Delta_{\rm d}$ to a $\delta_{j}^{\alpha}$ needs to be translated into some commutative operation so that we can deal with the $\exp\left[i\mathbf{k}\cdot\mathbf{x}_{j}+\omega t\right]$ properly, as follows:

[TABLE]

where (in the case of a $2$-dimensional system) we have $A=A(\mathbf{k})=-4\sin^{2}\left(k_{x}a/2\right)-4\sin^{2}\left(k_{y}a/2\right)$. Then this $A(\mathbf{k})$, or just $A$, replaces $\Delta_{\rm d}$, and eqn (46) can be written as an equation that no longer involves site-dependent quantities

[TABLE]

As we did in the main text, let us now set $w^{\alpha 0}=w_{0}$ and $w^{\alpha\beta}=w_{s}$, for any $\alpha\neq 0$ and $\beta\neq 0$, where $w_{0}$ and $w_{s}$ are constant attempt rates associated with particle-vacancy and particle-particle exchanges, respectively. We next derive two equations involving $\omega$. The first is obtained by summing up eqn (49) over all species $\alpha$. This gives

[TABLE]

where $\delta\rho\equiv\delta\rho_{0}$ and

[TABLE]

[TABLE]

are the moment densities (and their fluctuations), and we have used $\epsilon_{\alpha\beta}=\sigma_{\alpha}\sigma_{\beta}$. For the case with a continuous polydisperse attribute, the summations are integrals. The second equation is obtained by multiplying eqn (49) by $\sigma_{\alpha}$ and then summing it over all species. The result is

[TABLE]

By using that $\delta\rho$ and $\delta\rho_{1}$ should not both be nonzero in the spinodal dynamics and combining the two equations, one can now find $\omega(\mathbf{k})$. One obtains two non-trivial branches $\omega(\mathbf{k})$; we picked the largest one, which is eqn (24). There are a further $M-2$ eigenvalues (and eigenvectors) with $\delta\rho=\delta\rho_{1}=0$, leading to $M-2$ trivial branches such that all have identical negative eigenvalues. Explicitly, these are given by $\omega=\frac{A}{2}\left[w_{0}(1-\rho)+w_{s}\rho\right]$ but as we are looking for the fastest growing modes, they can be ignored.