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 the 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 also be used to study the 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 the size and the shape of the particles, for example, is inherent to all natural colloids. Even though modern synthesis techniques allow the preparation of almost monodisperse colloidal particles,^{1–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–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–16} Also, fractionation into several phases appears if the degree of polydispersity is high enough.^{17–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–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 behavior. 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} that 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–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 analog of the bulk phase diagram in sedimentation). The theory is general and can be applied to any colloidal system in sedimentation–diffusion–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 mass-polydisperse 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.

## II. THEORY

### A. Bulk

We use classical statistical mechanics to describe the thermodynamic bulk equilibrium of our mass-polydisperse 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 behavior 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 $\Lambda m=2\pi \u210f2/(mbkBT)$ 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 sake 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

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 with respect to the mass-density distribution, i.e.,

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

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–43} computer simulations,^{44–46} and empirical expressions.^{47–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

where $\eta 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 nematic–smectic.

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–diffusion–equilibrium 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, $\mu m0$ is the chemical potential of the species with buoyant mass *m* at elevation *z* = 0. The set of constant offsets $\mu m0$ 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 $\mu m0$. The offsets $\mu m0$ depend therefore on the choice of Λ_{m}. However, the sedimentation profiles *ρ*_{m}(*z*) remain unchanged, since $\mu m0$ are determined to match the prescribed density distribution.

Equation (13) is the sedimentation path^{31–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(\beta \mu m0)$ 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)–(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 $\rho \u0304m$ 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\u222bdm\rho \u0304m=A\u222bdm\u222b0hdz\rho m(z)$ the total number of particles, and *A* being the area of a cross section of the sample. Both $\eta \u0304$ and *f*_{P}(*m*) are directly comparable with experimental results, since $\eta \u0304$ 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.

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 $\eta \u0304$, 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 $\mu m0$ [see Eq. (13)]. We discretize *f*_{P}(*m*) and then numerically determine $\mu m0$ via a least square algorithm that iteratively solves for the prescribed $\eta \u0304fP(m)$ in a sample of height *h*. With the offsets $\mu m0$, 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 $\eta \u0304$ 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 $\eta \u0304$ and *f*_{P}(*m*).

For example, we show in Fig. 2(d) the offsets $\mu m0$ corresponding to the distribution prescribed in Fig. 2(c). We discretize in *m*, and hence the number of input variables $\rho \u0304m$ and unknown variables $\mu m0$ is the same. The self-consistency problem of finding $\mu m0$ 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 following Eq. (14), we conclude that the maximum number of layers that can appear in a sedimented sample of a mass-polydisperse system is 2*n*_{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, 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, 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 $\eta \u0304$ 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 $\eta \u0304$ 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–diffusion–equilibrium.

For each bulk phase transition, there can be at most three boundary lines in the stacking diagram, so-called sedimentation binodals.^{31–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.

**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

### A. Species-resolved probability distributions in mass-polydisperse systems

The imposed parent distribution of the mass-polydisperse 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 us 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 $\eta \u0304/\eta 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) that is symmetric with respect to *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 toward the bottom of the sample, but has a maximum up to 0.5*h* 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(f), the mass- and 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 $\eta \u0304/\eta 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–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–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 neutral 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–diffusion–equilibrium 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\u0304$ slightly shifted in the range of ±0.02, which is ∼10% of their standard deviations. The distributions are cut at *m* = ±0.95 around their respective mean.

The stacking diagram for the case $m\u0304=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\u0304=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, such as *LS* or *SL*, do not appear.

The situation is different for $m\u0304=\u22120.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, on top of the sample.

In Figs. 5(d) and 5(e), we show the stacking diagram for the remaining cases $m\u0304=0.01$ and 0.02, respectively. For comparison, we show always the sedimentation binodals of the buoyant neutral suspension with $m\u0304=0$. The position of the sedimentation binodals for the cases $m\u0304=\u22120.02$ and 0.02 are identical, but the associated stacking sequences are inverted. This was expected, since changing from $m\u0304=\u22120.02$ to 0.02 is equivalent to 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\u0304=0.01$ shows the same characteristics as $m\u0304=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\u0304=0$ and any other parent distribution considered, namely the lack of asymmetric stacking sequences, such as *LS* and *SL*. Mass-polydispersity therefore plays 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, 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(e)–(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, 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 system.

We also consider a class of parent distributions with both positive and negative buoyant masses, where we increase the variance of a bimodal distribution [see Fig. 6(d)]. Due 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 behavior 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 non-equilibrium gels.^{64} Depending, among other factors, on 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 straightforward extension of the present work.

In bulk, mass-polydispersity has no effect on the phase behavior. Hence, our mass-polydisperse model allows us to highlight the interplay between polydispersity and gravity, eliminating by construction the complex effects that shape-polydispersity 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–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, mass-polydispersity is likely the dominant effect in sedimentation–diffusion–equilibrium. This is particularly relevant in colloidal suspensions near density matching^{53–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–74} It would be interesting to analyze 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) hard-spheres. Moreover, colloids with a relatively large buoyant mass are often used^{57,59} and the sample height is not used as a 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 mass-polydisperse system and hence an approximated stacking diagram.

Polydispersity in the size of the particles affects the bulk behavior 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, *NI*. Due to the geometric properties of the sedimentation path of a mass-polydisperse system, such inversion of the sequence cannot occur in a mass-polydisperse system that contains particles with only positive (or only negative) buoyant masses. As correctly pointed out in Ref. 26, the inversion must therefore be a consequence of the 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.

## ACKNOWLEDGMENTS

We acknowledge useful discussions with Peter Sollich and Stefan Egelhaaf. This work was supported by the German Research Foundation (DFG) via Project No. 436306241.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Tobias Eckert**: Conceptualization (equal); Investigation (lead); Methodology (lead); Software (lead); Writing – original draft (equal); Writing – review & editing (equal). **Matthias Schmidt**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Daniel de las Heras**: Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX: EULER–LAGRANGE EQUATION

Carrying out the functional derivative with respect to *ρ*_{m} in Eq. (5) for the ideal free energy contribution to Ω yields

For the excess contribution, we find

where we have used the functional chain-rule and also the definition of the overall density [Eq. (2)] to calculate the functional derivative

Hence, introducing the excess chemical potential *μ*_{exc} = *δF*^{exc}[*ρ*]/*δρ* in Eq. (A2), it follows

Here, *μ* = *μ*_{exc} + *k*_{B}*T* ln(*ρ*) is the total chemical potential, including the ideal contribution *k*_{B}*T* ln(*ρ*), of the corresponding monodisperse system with the same overall density *ρ* as the mass-polydisperse system.

For the last contribution to Ω[*ρ*_{m}] in Eq. (4), we get