Sedimentation path theory for mass-polydisperse colloidal systems

Both polydispersity and the presence of a gravitational field are inherent to essentially any colloidal experiment. While several theoretical works have focused on the effect of polydispersity on the bulk phase behavior of a colloidal system, little is known about the effect of a gravitational field on a polydisperse colloidal suspension. We extend here sedimentation path theory to study sedimentation-diffusion-equilibrium of a mass-polydisperse colloidal system: the particles possess different buoyant masses but they are otherwise identical. The model helps to understand the interplay between gravity and polydispersity on sedimentation experiments. Since the theory can be applied to any parent distribution of buoyant masses, it can be also used to study sedimentation of monodisperse colloidal systems. We find that mass-polydispersity has a strong influence in colloidal systems near density matching for which the bare density of the colloidal particles equals the solvent density. To illustrate the theory, we study crystallization in sedimentation-diffusion-equilibrium of a suspension of mass-polydisperse hard spheres.


I. INTRODUCTION
A certain degree of polydispersity in e.g. the size and the shape of the particles is inherent to all natural colloids.Even tough modern synthesis techniques allow the preparation of almost monodisperse colloidal particles [1][2][3][4], a small degree of polydispersity is unavoidable.Understanding bulk phase equilibria in polydisperse systems is a significant challenge [5].Polydispersity alters the relative stability between bulk phases [6][7][8][9][10].Phases that are metastable in the corresponding monodisperse system can become stable due to polydispersity.Examples are the occurrence of hexatic columnar [11] and smectic phases [12] in polydisperse discotic liquid crystals, as well as macrophase separation in diblock copolymer melts [13].The opposite phenomenon can also occur.For example, crystallization in a suspension of hard-spheres is suppressed above a terminal polydispersity [14][15][16].Also, fractionation into several phases appears if the degree of polydispersity is high enough [17][18][19].A smectic phase of colloidal rods is no longer stable above a terminal polydispersity in the length of the particles [20].Dynamical processes such as shear-induced crystallization [21] are also affected by polydispersity.During drying, a strong stratification occurs in polydisperse colloidal suspensions [22,23], and the dynamics of large and small particles is different if the colloidal concentration is large enough [24,25].
Sedimentation-diffusion-equilibrium experiments are a primary tool to investigate bulk phenomena in colloidal suspensions.However, the effect of the gravitational field on the suspension is far from trivial [26][27][28][29][30] and it needs to be understood in order to draw correct conclusions about the bulk [31].Gravity adds another level of complexity to the already intricate bulk phenomena of a polydisperse suspension.To understand the interplay between sedimentation and polydispersity, we introduce here a mass-polydisperse colloidal suspension: a collection of colloidal particles with the same size and shape (and also identical interparticle interactions) but with buoyant masses that follow a continuous distribution.Since the interparticle interactions are identical, mass-polydispersity does not have any effect in the bulk phase behaviour.Hence, our model isolates the effects of a gravitational field on a polydisperse colloidal system from the effects that shape-and size-polydispersity generate in bulk.
We formulate a theory for mass-polydisperse colloidal systems in sedimentation-diffusion-equilibrium.The theory is based on sedimentation path theory [32,33] which incorporates the effect of gravity on top of the bulk description of the system.Sedimentation path theory uses a local equilibrium approximation to describe how the chemical potential of a sample under gravity changes with the altitude.So far, sedimentation path theory has been used to study sedimentation in colloidal binary mixtures [28,[31][32][33][34][35][36][37][38][39].In this work, we extend sedimentation path theory to mass-polydisperse systems.Using statistical mechanics, we obtain the exact expression for the sedimentation path of the mass-polydisperse suspension combining the individual paths of all particles in the distribution.We use a model bulk system to illustrate and highlight the key concepts of the theory, such as the construction of the sedimentation path and that of the stacking diagram (which is the analogue of the bulk phase diagram in sedimentation).The theory is general and can be applied to any colloidal system in sedimentationdiffusion-equilibrium.Moreover, the theory contains the description of a monodisperse system as a special limit (delta distribution of the buoyant masses).As a proof of concept, we study sedimentation of a suspension of masspolydisperse hard-spheres with different buoyant mass distributions.We find that mass polydispersity plays a major role in systems near density matching.For example, near density matching the packing fraction and the height of the sample at which crystallization is observed in sedimentation-diffusion-equilibrium are strongly influenced by the details of the mass distribution.

A. Bulk
We use classical statistical mechanics to describe the thermodynamic bulk equilibrium of our masspolydisperse colloidal system.The term bulk refers here to an infinitely large system in which boundary effects can be neglected and that is not subject to any external field.The particles differ only in their buoyant masses.Since the buoyant mass does not play any role in bulk, the bulk phenomenology of our model is identical to that of a monocomponent system in which only one buoyant mass is present.Only when gravity is incorporated into both systems the buoyant mass becomes a relevant parameter and the behaviour of the mass-polydisperse and the monodisperse colloidal systems will differ from each other.
The total Helmholtz free energy F is the sum of the ideal and the excess contributions, i.e.F = F id + F exc .In a mass-polydisperse system, the free energy is a functional of ρ m , the density distribution of species with buoyant mass m.For simplicity, we work with a scaled, dimensionless, buoyant mass m = m b /m 0 , where m b is the actual buoyant mass of a particle and m 0 is a reference buoyant mass.Sensible choices relate m 0 to e.g. the average buoyant mass of the distribution or its standard deviation.The concrete definition of m 0 is given below in each considered system.
The ideal contribution to the free energy is a functional of ρ m and is given exactly by where k B is the Boltzmann's constant and T is the absolute temperature.Without loss of generality we measure ρ m relative to the thermal de Broglie wavelengths with reduced Planck's constant .Note that the value of Λ m does not play any role here since altering Λ m simply adds a term to the free energy that is proportional to the total number of particles with buoyant mass m.Such term can be reinterpreted as a change of the origin of the chemical potential of the species with buoyant mass m.
The integration over m in Eq. ( 1) reflects the fact that due to the mass-polydispersity, the buoyant mass is a continuous variable.For the shake of simplicity, we omit the positional argument r in the density distribution as well as its corresponding space integral that appear in bulk-phases with positional order such as crystalline phases.
The ideal free energy, Eq. ( 1), accounts for the entropy of mixing of our mass-polydisperse system.The overall density across all species ρ follows directly from the density distribution of buoyant masses ρ = dm ρ m . ( Since the interparticle interaction is independent of the buoyant masses of the particles, only the density across all species ρ enters into the excess (over ideal) free energy.Hence, the excess free energy functional must satisfy The grand potential is also a functional of ρ m given by where µ m is the chemical potential of the species with buoyant mass m.In equilibrium Ω[ρ m ] is minimal w.r.t. the mass-density distribution, i.e.
The Euler-Lagrange equation associated to Eq. ( 5), see derivation in Appendix A, reads where µ is the chemical potential of a monodisperse system with overall density ρ, see Eq. ( 2).Hence, it follows from Eq. ( 6) that the density of particles with buoyant mass m can be written as Integrating Eq. ( 7) over m on both sides, and using Eq. ( 2) on the left hand side, leads to Since ρ = 0, we obtain which constitutes an exact analytic expression for the chemical potential of the monodisperse bulk system in terms of the chemical potentials of the individual species µ m in the mass-polydisperse system.In a monodisperse system there exists only a single species and Eq. ( 10) holds trivially.

B. Particle Model
To proceed we need the bulk equation of state (EOS) of the monodisperse colloidal system, ρ EOS (µ).Given an interparticle interaction potential, several methods can be used to obtain the corresponding bulk EOS.These include e.g.density functional theory [40], liquid state integral equation theory [41][42][43], computer simulations [44][45][46] and empirical expressions [47][48][49].Here, and with the only purpose of illustrating our theory we use a model (fabricated) EOS that contains two phase transitions, see Fig. 1(a).Our model EOS satisfies both the ideal gas limit and also the close packing limit characteristic of systems with hard core interactions here η EOS is the packing fraction (percentage of volume occupied by the particles) according to the EOS and η cp is the close packing fraction.Such EOS could represent e.g. a lyotropic colloidal system with two first-order bulk phase transition, say isotropic-nematic and nematicsmectic.
Apart from the model EOS, we also illustrate and validate the theory by studying sedimentation of a suspension of hard-spheres.We use the analytical EOS proposed by Hall [50], which describes the liquid (L) and solid crystalline (S) phases of a hard sphere system.The Hall EOS was originally formulated using the compressibility factor as a function of the density.Following Ref. [51], we numerically integrate the analytical Hall EOS to obtain the chemical potential as a function of the density, see Fig. 1(b) for a graphical representation.It is sufficient to fix ρ EOS (µ) up to an arbitrary additive constant in µ.Hence, for convenience, we choose µ = 0 as the chemical potential at the liquid-solid first order phase transition.

C. Sedimentation
To incorporate gravity into our theory, we extend sedimentation path theory [32,52] as formulated for finite height samples [31,33] to include mass-polydispersity.As often done in colloidal sedimentation, we assume that all horizontal slices of a sample in sedimentation-diffusionequilibrium can be described as a bulk equilibrium state, and also that they are independent of each other.This local-equilibrium approximation is justified if the correlation lengths are small compared to the gravitational lengths ξ m = k B T /(m b g), which is the case in many colloidal systems.Here g is the acceleration of gravity.
We work in units of the thermal energy k B T , the gravitational constant g, and the reference mass m 0 for ease of comparability between different systems.Using m 0 we define a reference gravitational length ξ = k B T /(m 0 g), which acts as our fundamental length scale.
We treat the slices for each elevation z as a bulk system with local chemical potentials for each species µ m given by Here µ 0 m is the chemical potential of the species with buoyant mass m at elevation z = 0.The set of constant offsets µ 0 m in µ m (z) is a priori unknown and must be determined via an iterative numerical procedure to match the prescribed mass-resolved density distribution ρ m .Returning to the discussion about the thermal wavelengths, altering the value of Λ m would only introduce a constant term ln(Λ m ) in Eq. ( 6) that can be reabsorbed in Eq. ( 13) as a shift of the chemical potential µ m via the offset µ 0 m .The offsets µ 0 m depend therefore on the choice of Λ m .However, the sedimentation profiles ρ m (z) remain unchanged, since µ 0 m are determined to match the prescribed density distribution.
Equation ( 13) is the sedimentation path [31][32][33]52] of the species with buoyant mass m.It hence describes how the chemical potential of each species varies linearly with z in the range 0 ≤ z ≤ h, with h the sample height.The local chemical potential for each species either decreases (m b > 0) or increases (m b < 0) with the elevation z, depending on the sign of the buoyant mass.
The sedimentation path of each species µ m (z) is just a straight line, see Fig. 2(a), as in the case of monodisperse systems.Next, we combine all paths at each elevation z to obtain an effective chemical potential µ eff (z).Inserting µ m (z) in Eq. ( 13) into Eq.( 10) yields the sedimentation path of a mass-polydisperse system Equation ( 14), which has the form of a LogSumExp function, describes how the effective chemical potential of the mass-polydisperse system varies vertically along the sample in sedimentation-diffusion-equilibrium.We give an example of µ eff (z) in Fig. 2(b).The sedimentation path is obtained from the set of µ m (z) in Fig. 2(a) via Eq.( 14).The sedimentation path is no longer a straight line even though the individual paths for each species are lines.Since (i) the logarithm is a concave function, (ii) the scalars exp(βµ 0 m ) are positive, and (iii) the exponential is a convex function, it follows that µ eff (z), as given by Eq. ( 14), is a convex function of the elevation z.This is a strong constraint on the possible shapes of µ eff (z).It means that (i) µ eff (z) can have only one minimum and also that (ii) the local maxima of µ eff (z) in the interval 0 ≤ z ≤ h are either z = 0, z = h, or both of them.As we discuss below, the extrema of the path µ eff (z) are important because they determine the layers of different bulk phases that form in the sample.
Via the equation of state ρ EOS (µ) for the bulk density we then obtain the density profile across all species at elevation z from Eq. ( 14).
The density of species with buoyant mass m at elevation z follows then by inserting Eqs. ( 13) to (15) into Eq.( 7) The average density of particles with buoyant mass m in a sample with height h is then given by The value of ρm is also the density of particles with buoyant mass m in the initial distribution, i.e. before the particles sedimented and equilibrated.
The average packing fraction is where v 0 is the particle volume.
The parent distribution, which gives the overall probability of finding a particle with buoyant mass m anywhere in the sample can be obtained as with N = hA dm ρm = A dm h 0 dz ρ m (z) the total number of particles, and A being the area of a cross section of the sample.Both η and f P (m) are directly comparable with experimental results, since η is the concentration of particles in the stock solution (before sedimentation) and f P (m) describes the mass-polydispersity of the particles, normalized by the total concentration.
From the definitions ( 18) and ( 19) we can get back the average density of specie m via To obtain the sedimentation-diffusion-equilibrium of a mass-polydisperse colloidal system, we start prescribing the sample height h, the average packing fraction of the sample η, and the parent distribution f P (m).An illustrative parent distribution that contains particles with both positive and negative buoyant masses is shown in Fig. 2(c).These initial conditions are sufficient to find the as yet undetermined offsets on the chemical potential for each species µ 0 m , see Eq. ( 13).We discretize f P (m) and then numerically determine µ 0 m via a least square algorithm which iteratively solves for the prescribed ηf P (m) in a sample of height h.With the offsets µ 0 m we calculate the corresponding µ eff (z) via Eq.( 14).Next we obtain ρ(z) and ρ m (z) via Eqs.( 15) and ( 16), respectively.The profiles ρ(z) and ρ m (z) determine both η and f P (m) via Eqs.( 18) and ( 19) respectively.The least square algorithm finds then the offsets that minimize the difference to the prescribed (target) values of η and f P (m).
For example, we show in Fig. 2(d) the offsets µ 0 m corresponding to the distribution prescribed in Fig. 2(c).We discretize in m, and hence the number of input variables ρm and unknown variables µ 0 m is the same.The selfconsistency problem of finding µ 0 m is therefore well defined.The set of sedimentation paths µ m (z) in Fig. 2(a) are obtained with the offsets calculated in Fig. 2(d).The effective sedimentation path µ eff (z), see Fig. 2(b), of the mass-polydisperse system follows then from the set of paths for each species µ m (z).
The sedimentation path of the mass-polydisperse system determines the stacking sequence, i.e. the set of layers of bulk phases that are observed in the sample under gravity.Every time the path crosses the coexistence chemical potential of a bulk transition, an interface between the coexisting phases appears in the cuvette.By looking at the crossings between the sedimentation path and the bulk binodals we determine the stacking sequence and the position of the interfaces between stacks.For example, the sequence corresponding to the path in Fig. 2(b) is BABC (from top to bottom of the sample).
Extended Gibbs phase rule.Given the convexity properties of the sedimentation path, recall our discussion below Eq. ( 14), we conclude that the maximum number of layers that can appear in a sedimented sample of a mass-polydisperse system is 2n b − 1, with n b the number of different stable phases in bulk.This corresponds to the stacking sequence of a mass-polydisperse suspension with positive and negative buoyant masses in which all phases occur repeatedly except the middle layer, which corresponds to the bulk phase stable at low chemical potential.In our model EOS, the stacking sequence with the maximum number of layers is CBABC, for which the sedimentation path is similar to the one in Fig. 2(b) but extended such that it reenters the C region at high elevations.
If the parent distribution contains only buoyant masses of the same sign, then the maximum number of layers in a stacking sequence is simply n b , the number of stable bulk phases.

D. Stacking diagram
Different sedimentation paths can give rise to distinct stacking sequences.The set of all possible stacking sequences can be represented in a stacking diagram.In binary mixtures, the sedimentation paths of both species vary linearly with z.In mass polydisperse systems, we average the linear local chemical potentials µ m (z), Eq. ( 13), of all species together, according to Eq. ( 14), and obtain a non-linear effective chemical potential µ eff (z).Even though the sedimentation paths are no longer straight lines, the same ideas as in the case of binary mixtures [31,32] apply for the construction of the stacking diagram.In short, we must find all the sedimentation paths that constitute a boundary between two or more stacking sequences in the stacking diagram.Examples of such paths are shown in Fig. 3(a).The boundary paths are the sedimentation paths µ eff (z) that either end [paths 1 and 4 in Fig. 3(a)], start (paths 2 and 5), or are tangent (paths 3 and 6) to a bulk binodal.These paths are a boundary between two or more stacking sequences since an infinitesimal change of the path in general alters the stacking sequence.Without gravity (i.e. in bulk) the mass-polydisperse system behaves like a monocomponent system, since the interparticle interaction potential is independent of the buoyant mass.Thus, in bulk, there is only a single relevant chemical potential.In the chemical potential vs. height plane, the bulk transitions are simply horizontal lines independent of z, see Fig. 3(a).Hence, given that the sedimentation path is convex, a path tangent to a bulk binodal is also a path for which the minimum coincides with the chemical potential of the bulk transition, like e.g.paths 3 and 6 in Fig. 3(a).For other types of bulk phase coexistence such as critical and triple points, the procedure to find the boundary paths is the same as the one just described for a bulk binodal.
Next we find the total density profile ρ(z) and the average packing fraction η corresponding to each of the boundary sedimentation paths via Eqs.( 15) and ( 18), respectively.To obtain the full stacking diagram, we repeat the procedure for every sample height h ranging from zero to the desired maximal sample height.This provide us with the stacking diagram in the (experimentally relevant) plane of average packing fraction η and sample height h, see Fig. 3(b).Each point in the stacking diagram represents one sedimentation path and it hence represents one specific sample in sedimentation-diffusionequilibrium.
For each bulk phase transition there can be at most three boundary lines in the stacking diagram, so-called sedimentation binodals [31][32][33].The sedimentation binodals corresponding to the paths that either start or end at the binodal are always present independently of the parent distribution and the sample height.On the other hand, the sedimentation binodal corresponding to paths tangent to the bulk transition appears if and only if the sedimentation path presents a minimum at intermediate values of z.It follows from Eq. ( 14) that a minimum in µ eff (z) not located at the bottom (z = 0) or the top (z = h) of the sample can appear only if the parent distribution contains both positive and negative buoyant masses.Even in that case, there might be sample heights for which the path does not have a minimum at intermediate elevations.
In our illustrative example, there are two bulk phase transition (A − B and B − C), see Fig. 1(a), and the parent distribution is made of particles with positive and negative buoyant masses, see Fig. 2(c).The stacking diagram contains six sedimentation binodals, see Fig. 3(b).For sample heights h/ξ 0.4 only two types of sedimentation binodals can be observed.In this low height regime, we cannot find sedimentation paths tangent to the binodal since µ eff (z) does not have a minimum at intermediate elevations 0 < z < h.
Within our local equilibrium approximation, in the limit h → 0 the sedimentation path reduces to a point and hence the stacking diagram reduces to the bulk phase diagram.In a real system, confinement and surface effects such as wetting and layering will become relevant in the limit of short sample heights.2(c).The coexistence chemical potentials for the A-B (µcoex = 0) and for the B-C (βµcoex = 0.25) bulk transitions are indicated by solid-black horizontal lines.There are six sedimentation paths labeled from 1 to 6.The average packing fraction η of each sample is such that the corresponding path either ends at (solid paths 1 and 4), starts at (dashed paths 2 and 5) or is tangent to (dotted paths 3 and 6) a bulk binodal (horizontal lines).The points where the paths touch the coexistence bulk chemical potential are marked by solid circles.(b) Stacking diagram in the plane of average packing fraction η/ηcp and sample height h/ξ for the model EOS in Fig. 1(a) and parent distribution as in Fig. 2(c).The position of the six boundary sedimentation paths in (a) is marked in (b) using the corresponding labels 1 to 6.The sedimentation binodals of paths that end, start, and are tangent to the bulk binodals are indicated with solid-, dashed-and dottedlines, respectively.The stacking sequences are labeled from the top of the sample to the bottom.Each point in the stacking diagram is a sample in sedimentation.The sketch shows the stacking sequence BABC and relative layer thicknesses of the sample with η/ηcp = 0.45 and h/ξ = 3.5 (indicated by a black circle).Mass-monodisperse system.Our method to construct the stacking diagram for mass-polydisperse systems contains as a limiting case the monodisperse system.In a monodisperse system all particles possess the same buoyant mass.Hence, Eq. ( 14) reduces to Thus, as expected, the sedimentation path of a monodisperse system is the segment of a line, linear in z.In the stacking diagram, only the sedimentation binodals of paths that start, i.e. µ eff (h) = µ coex , or end, i.e. µ eff (0) = µ coex , at the bulk binodal (given by µ coex ) appear.The sedimentation path of a monodisperse system can never have a minimum at intermediate elevations.

III. RESULTS
We next apply our theory to the arguably best studied colloidal system to date: hard spheres.We study sedimentation of a mass-polydisperse hard sphere system using the Hall equation of state [50], represented in the plane of µ and η in Fig. 1(b), to describe the bulk of the system.

A. Species-resolved probability distributions in mass-polydisperse systems
The imposed parent distribution of the masspolydisperse system, f P (m), describes the probability of finding a particle with a certain buoyant mass m anywhere in the system.Experimentally, this corresponds to the stock solution.After letting the dispersion settle under gravity to reach sedimentation-diffusion-equilibrium, a height-dependent density profile develops.The overall probability distribution integrated over the whole sample is still f P (m) since particles are conserved.However, at each horizontal slice the mass composition is generally different from f P (m).One expects e.g.heavier particles to concentrate next to the bottom of the sample as compared to lighter particles.Sedimentation path theory allows to carry out a detailed study of the mass distribution along the sample.
We study first a mass-polydisperse dispersion of hard spheres with only positive buoyant mass.The parent distribution is a Gaussian centered around m = 1 and cut at m = 0 and m = 2, i.e.only buoyant masses in the range 0 ≤ m ≤ 2 are allowed.The mean packing fraction is η/η cp = 0.6.Under gravity, the sample develops the stacking sequence: top liquid and bottom solid (LS).We show the probability f (m, z) of finding a particle with buoyant mass m at elevation z in Fig. 4(a).The probability distribution f m (z) for a fixed buoyant mass m and resolved in z, as well as the probability distribution f z (m) for a fixed z resolved in m are shown in Figs.4(b) and 4(c), respectively.The distributions f m (z) and f z (m) correspond to vertical and horizontal slices of the full distribution f (m, z), respectively.The distributions f z (m) are shifted and skewed, Fig. 4(c), as compared to the parent distribution f P (m) (black-dashed line) which is symmetric w.r.t.m = 1.As expected, heavier particles are more frequently found at the bottom of the sample.This becomes more apparent when we look at f m (z).Fig. 4(b).There is a clear depletion of lighter particles from the bottom of the sample.Interestingly, the probability distribution along z of particles with m 1.01 is not monotonically increasing towards the bottom of the sample, but has a maximum up to 0.5h above the bottom.Lighter particles are displaced by heavier particles from the bottom as a result of a balance between only two contributions: the gravitational energy and the entropy of mixing.The excess free energy does not play a role in determining the relative position of the particles according to their buoyant masses.Interchanging heavier for lighter particles and vice versa does not alter the overall density, and thus the excess free energy F exc [ρ], which is a functional of only the overall density ρ, is not affected.
We also show in Figs.4(d), 4(e), and 4(f), the massand height-resolved probability distributions of a sample with a parent distribution containing both positive and negative buoyant masses.The parent distribution is a Gaussian centered around m = 0.03 and cut at m = ±1.9.The initial packing fraction is η/η cp = 0.7 and the stacking sequence is SLS.The liquid-solid interfaces occur at elevations z/h = 0.25 and 0.8 and are visible as discontinuities of the distribution functions.On the top (bottom) of the sample particles with negative (positive) buoyant masses are more frequently found.This is visible in Fig. 4(f) as a shift toward negative or positive buoyant masses of the distributions belonging to the solid crystalline layers.

B. Mass-polydispersity close to density matching
In density matching colloidal experiments, the mass density of the colloidal particles is very close to the mass density of the solvent.If the density match between particle and solvent is perfect, the buoyant mass of the colloids vanishes and therefore gravity has no effect on the sample.This, in principle, would allow to carry out a direct comparison between bulk phenomena and sedimentation experiments.In practice, however, preparing experimentally a perfect density matching solution is challenging.Density matching is typically achieved by combining solvents with different mass densities in the correct proportions to match the mass density of the particles [53][54][55][56].To sterically stabilize the colloidal particles, they are frequently coated with a polymer layer of a different density than that of the particle core [57][58][59][60][61]. Due to the polydisperse nature of most colloidal systems, the effective particle density (including both the core and the coating layer) can vary between the particles.As a result, not all the particles in the solution can have neu- tral buoyancy.The buoyant mass of the particles falls within a range roughly centered around neutral buoyancy.In general there will be particles that have either slightly positive or slightly negative buoyant masses.We will see here that small deviations from density matching can have a strong effect on sedimentation-diffusionequilibrium experiments.
We model a system close to density matching by a parent Gaussian distribution f P (m) roughly centered around a buoyant mass m = 0, as shown in Fig. 5(a).We study four different cases, with the mean of the Gaussian m sightly shifted in the range of ±0.02, which is approximately 10% of their standard deviations.The distributions are cut at m = ±0.95around their respective mean.
The stacking diagram for the case m = 0 is shown in Fig. 5(b).Near density matching, the sedimentation paths are rather horizontal and sensitive to the precise form of the parent distribution.Hence, the small deviations between the (imposed) target and the (actual) numerical parent distributions that arise in the iterative procedure due to numerical inaccuracies, can have a noticeable effect.This is the reason behind the scattered data points (symbols) in the sedimentation binodals of Fig. 5(b).With a symmetrical parent distribution around m = 0 (i.e.m = 0) neither particles with positive nor with negative buoyant mass are favored.Thus, only symmetric stacking sequences (with respect to the midpoint of the sample z = h/2) occur, namely L, S and SLS.Asymmetric sequences like LS or SL do not appear.
The situation is different for m = −0.02,where particles with negative buoyant mass that cream up are predominant, see Fig. 5(c).Consequently, we also observe the stacking sequence SL, with the denser, solid phase, The position of the sedimentation binodals for the cases m = −0.02and 0.02 are identical, but the associated stacking sequences are inverted.This was expected, since changing from m = −0.02 to 0.02 is equivalent inverting the direction of gravity and thus interchanging the meaning of top and bottom of the sample.This is also the reason why we observe the stacking sequence LS in Fig. 5(e) in the region occupied by SL in Fig. 5(c).The case m = 0.01 shows the same characteristics as m = 0.02, but with the position of the sedimentation binodals roughly rescaled in the h/ξ axis by a factor of 1/2, which is the ratio between the mean of the corresponding parent distributions.Most notable, there is a qualitative difference between the case m = 0 and any other parent distribution considered, namely the lack of asymmetric stacking sequences such as LS and SL.Mass-polydispersity plays therefore an important role in colloidal suspensions close to density matching and even small deviation from density matching can have drastic effects on the stacking diagram.

C. Mass-polydispersity away from density matching
Not all types of parent distributions are as sensitive to mass-polydispersity as those representing a system near density matching.In many cases the stacking diagram is robust against perturbations of the parent distribution.To show this, we construct here four classes of parent distributions and calculate the corresponding sedimentation paths.The sedimentation paths are quite similar within each class.We hence can conclude that the corresponding stacking diagrams are also alike.Recall that the stacking diagram is constructed from the set of special paths, µ eff (z), that either start at, end at, or are tangent to the bulk binodal (see Fig. 3).Thus, if two systems share similar paths for a range of packing fractions and sample heights, then the stacking diagrams will also be similar.
The four classes of parent distributions and the corresponding sedimentation paths are shown in Fig. 6.We construct several distributions within each class by varying a control parameter.In Fig. 6(a) we increase the mass-polydispersity by interpolating between unimodal and bimodal Gaussian distributions.In Fig. 6(b) we increase the variance of a Gaussian distribution.In Fig. 6(c) we vary the skewness of the distribution while keeping the first and the second moment unaltered.In all cases the distributions contain only positive masses and varying the control parameter has little effect on the sedimentation paths, even when we e.g.drastically increase the degree of mass-polydispersity (second moment of the distribution).The corresponding sedimentation paths, shown in Fig. 6 panels (e) to (g), deviate only slightly from a straight line with a slope given by the mean buoyant mass of the distribution.Hence, in sedimentation-diffusion-equilibrium, In all cases we use the Hall EOS for hard spheres, the packing fraction is η/ηcp = 0.5, and the sample height is h/ξ = 80.The dashed lines in panels (e) to (g) are the linear trends of the corresponding sedimentation paths (displaced vertically for a better visualization).The slopes (indicated next to each dashed line) coincide in all three cases with the mean mass of the corresponding family of distributions, which are 1 (a), 1 (b) and 0.6 (c).mass-polydisperse systems in which only positive or negative buoyant masses are present are similar to a reference monodisperse system.(Recall that µ(z) = µ 0 − m b gz for a monodisperse system.)The monodisperse reference system has the same particle mass as the mean of the mass distribution of the mass-polydisperse We also consider a class of parent distributions with both and negative buoyant masses, where we increase the variance of a bimodal distribution, see Fig. 6(d).to the presence of buoyant masses with different sign, the suspension does not behave like a monodisperse system under gravity, and hence µ eff (z) is not close to a straight line, see Fig. 6(h).Still the increase in the degree of mass-polydispersity does not affect the behaviour of the system strongly since the paths do not deviate much from each other.

IV. SUMMARY AND CONCLUSIONS
Sedimentation path theory [32,33] was initially developed to study sedimentation-diffusion-equilibrium of binary mixtures.The theory describes systems that are in equilibrium under the presence of a gravitational field and therefore cannot be used to describe non-equilibrium phenomena such as drying [62,63] or systems that get arrested due to e.g. the formation of glasses [58] and nonequilibrium gels [64].Depending, among other factors, on the the buoyant mass of the colloids, the experimental equilibration times can vary from a few hours to several months [28].We have extended here sedimentation path theory to deal with mass-polydisperse colloidal systems, i.e. the particles are identical except for the value of their buoyant masses.We derived an exact equation for the sedimentation path of the mass-polydisperse system, Eq. ( 14), that combines all the sedimentation paths of the individual species.The resulting equation has the structure of the LogSumExp function, often used in machine learning algorithms for its smooth approximation to the maximum function [65].Adding mass polydispersity to a binary mixture is in principle a straight forward extension of the present work.
In bulk, mass-polydispersity has no effect on the phase behaviour.Hence, our mass-polydisperse model allows us to highlight the interplay between polydispersity and gravity, eliminating by construction the complex effects that shape-and size-polydispersity generates in bulk [5, 7-9, 12, 14, 16, 19, 20, 66].Beyond its fundamental interest, a mass-polydisperse system can be specifically realized experimentally by e.g.synthesizing core-shell nanoparticles [67][68][69][70] with the same overall size but different relative size between the core and the shell.In addition, if the degree of size-polydispersity is small, then mass-polydispersity is likely the dominant effect in sedimentation-diffusion-equilibrium.This is particularly relevant in colloidal suspensions near density matching [53][54][55][56], in which mass-polydispersity has a big effect on the stacking diagram under gravity: two mass distributions that are only slightly different can give rise to topologically different stacking diagrams containing different stacking sequences.
Granular media is a related system in which the particles can be polydisperse [71][72][73][74].It would be interesting to analyse the effects of mass-polydispersity in granular systems.For example, phase separation induced by mass-polydispersity might occur in vibrated monolayers [75,76] of granular systems.
Despite the relevance of the hard-sphere model in soft matter, there are only a few experiments on the sedimentation-diffusion-equilibrium of (quasi) hardspheres.Moreover, colloids with a relatively large buoyant mass are often used [57,59] and the sample height is not used as control parameter.A systematic experimental study of the stacking diagram of hard spheres would be valuable.
In bulk, it is sometimes possible to approximate the free energy of a polydisperse system using only a finite number of moments of the parent distribution [77,78].In a similar way, using the first moment of the parent distribution it is possible to obtain a reasonable approximation for the effective sedimentation path of the masspolydisperse system, and hence an approximated stacking diagram.
Polydispersity in the size of the particles affects the bulk behaviour of the suspension and therefore also the sedimentation-diffusion-equilibrium.For example, van der Kooij et al. studied the sedimentation of polydisperse colloidal platelets [26].Their particle distribution contained platelets of different sizes but only positive buoyant masses.By changing the overall packing fraction, they found a striking inversion of the stacking sequence from the expected top isotropic and bottom nematic, IN , to top nematic and bottom isotropic, N I. Due to the geometric properties of the sedimentation path of a mass-polydisperse system, such inversion the sequence cannot occur in a mass-polydisperse system that contains particles with only positive (or only negative) buoyant masses.correctly pointed out in Ref. [26], the inversion must therefore be a consequence of interplay between gravity and size-polydispersity.Sedimentation path theory could be applied on top of a bulk theory for size-polydisperse systems in order to describe such effects.

Figure 1 .
Figure 1.Packing fraction ηEOS relative to close packing ηcp as a function of the scaled chemical potential βµ for (a) our model equation of state, and (b) the Hall equation of state [50] for hard spheres.Our model EOS (a) contains three different bulk phases named A, B and C which could correspond to e.g. the isotropic, the nematic, and the smectic phases of a lyotropic liquid crystal.The Hall EOS (b) describes the liquid (L) and the solid crystalline (S) phases of a hard-sphere system.The vertical dotted lines indicate the chemical potentials of the different bulk phase transitions.Without loss of generality, we have translated the origin of chemical potential such that it coincides with the chemical potential of (a) the A-B, and (b) the L-S transitions.

Figure 2 .
Figure 2. (a) Local chemical potentials βµm(z) as a function of the elevation z scaled with the sample height h.Each sedimentation path varies linearly with z and it is colored according to the buoyant mass of the species m (top-left color bar).(b) Effective chemical potential βµ eff (z) as a function of the scaled elevation z/h.The non-linear sedimentation path of the mass-polydisperse system (b) is the result of combining the sedimentation paths of each species (a) according to the LogSumExp structure in Eq. (14).The sedimentation paths are defined in the interval 0 ≤ z ≤ h, and the sample height is h = 2ξ.(c) Imposed parent distribution fP(m), and plots of the parent distribution scaled with several values of the average packing fraction η (see color bar).The imposed fP(m) and the value of η in (c) fix the offsets βµ 0 m for all buoyant masses m, which are shown in panel (d).The pink arrows illustrate the iterative procedure to find the effective sedimentation path: for a fixed distribution fP(m) and packing fraction η (c), we give an initial guess for the offsets µ 0 m (d), calculate the individual paths µm(z) (a), and combine them to get the effective path µ eff (z) (b).Using the effective path we obtain the density profile and then the resulting distribution of particles and the average packing fraction.With this information, we readjust the offsets µ 0 m

Figure 3 .
Figure 3. (a) Sedimentation paths in the plane of effective chemical potential βµ eff as a function of the elevation z/h for samples of height h = 1.7ξ and a parent distribution like in Fig.2(c).The coexistence chemical potentials for the A-B (µcoex = 0) and for the B-C (βµcoex = 0.25) bulk transitions are indicated by solid-black horizontal lines.There are six sedimentation paths labeled from 1 to 6.The average packing fraction η of each sample is such that the corresponding path either ends at (solid paths 1 and 4), starts at (dashed paths 2 and 5) or is tangent to (dotted paths 3 and 6) a bulk binodal (horizontal lines).The points where the paths touch the coexistence bulk chemical potential are marked by solid circles.(b) Stacking diagram in the plane of average packing fraction η/ηcp and sample height h/ξ for the model EOS in Fig.1(a) and parent distribution as in Fig.2(c).The position of the six boundary sedimentation paths in (a) is marked in (b) using the corresponding labels 1 to 6.The sedimentation binodals of paths that end, start, and are tangent to the bulk binodals are indicated with solid-, dashed-and dottedlines, respectively.The stacking sequences are labeled from the top of the sample to the bottom.Each point in the stacking diagram is a sample in sedimentation.The sketch shows the stacking sequence BABC and relative layer thicknesses of the sample with η/ηcp = 0.45 and h/ξ = 3.5 (indicated by a black circle).

Figure 4 .
Figure 4. (a) Probability of finding a particle with buoyant mass m at elevation z relative to the sample height h in a hard sphere system modeled using the Hall EOS.The light contour lines in (a) indicate points (m, z) at which the chemical potential µm(z) [see Eq. (13)] is equal to the bulk liquid-solid coexistence chemical potential.The white arrow indicates the position of the liquid (L) -solid crystalline (S) interface.(b) Vertical slices of panel (a) for fixed buoyant mass m given by the colorbar.The inset is a sketch of the sample.(c) Horizontal slices of panel (a) for fixed elevation z given by the colorbar.The imposed parent distribution fP(m) is a Gaussian with standard deviation 0.4 and centered around m = 1, i.e. the reference buoyant mass m0 is the mean of the parent distribution (see black dashed-line).The sample has a height h = 80ξ with gravitational length ξ and a packing fraction η = 0.6ηcp, relative to the close packing fraction ηcp.Panels (d-f) are the same as panels (a-c), but for a Gaussian with standard deviation 0.6 centered around 0.03 as the parent distribution fP(m) (black dashed line), slightly favoring particles with positive buoyant mass, sample height h = 120ξ and packing fraction η = 0.7ηcp.The crosses in (b) and (e) indicate the position of the local maxima in the probability distribution fm(z) along elevation z for fixed buoyant mass m.

Figure 5 .
Figure 5. (a) Four parent fP Gaussian distributions with slightly shifted mean m in the range 0 ± 0.02, as indicated.The standard deviation is 0.2 in all cases, and the distributions are cut at m = ±0.95around their respective mean.The corresponding stacking diagrams for a mass-polydisperse system of hard spheres in the plane of average packing fraction η (relative to close packing ηcp) and sample height h (relative to the gravitational length ξ) are shown in (b) for the parent distributions with m = 0, in (c) for m = −0.02, in (d) for m = 0.01, and in (e) for m = 0.02.The sedimentation binodals of paths that end, start, and are tangent to the bulk binodals are indicated with solid-, dashed-and dotted-lines, respectively.The symbols are the data points.The sedimentation binodals of the case m = 0 are shown for reference in all the stacking stacking diagrams.Note that for the case m = 0 the sedimentation binodals of paths that either start or end at the bulk transition coincide since the parent distribution is symmetrical around m = 0.The bulk system exhibits liquid (L) and solid crystalline (S) phases.The stacking sequences are labeled from the top to the bottom of the sample.

Figure 6 .
Figure 6.Families of parent distributions fP(m) in the form of (a) the sum of two Gaussians which move apart symmetrically around the buoyant mass m = 1, (b) Gaussian with increasing standard deviation from 0.1 to 0.5, (c) χ 2 -distribution with the degree of freedom increasing from 3 to 13, and (d) sum of two Gaussians with mean values at m = −0.5 and m = 1.5 and increasing standard deviation from 0.1 to 0.5.The distributions in (c) have mean 0.6 and standard deviation 1.The effective sedimentation paths βµ eff as a function of the scaled elevation z/h corresponding to the families of distributions in panels (a) to (d) are shown in panels (e) to (h), respectively.In all cases we use the Hall EOS for hard spheres, the packing fraction is η/ηcp = 0.5, and the sample height is h/ξ = 80.The dashed lines in panels (e) to (g) are the linear trends of the corresponding sedimentation paths (displaced vertically for a better visualization).The slopes (indicated next to each dashed line) coincide in all three cases with the mean mass of the corresponding family of distributions, which are 1 (a), 1 (b) and 0.6 (c).