# Optimal packing of polydisperse hard-sphere fluids

## I Introduction

Synthetic colloids are never perfectly monodisperse. Often, this polydispersity is a drawback—for instance, polydispersity is a problem in the preparation of high-quality colloidal crystals, that are needed in photonic bandgap materials. However, occasionally, polydispersity is desirable, because it allows us to achieve material properties that cannot be realized with monodisperse colloids. For instance, monodisperse colloidal systems can fill at most 74.05 % of space in the crystalline phase (regular close packing) and some 63% in the liquid/glassy state (random close packing). In contrast, colloids with a properly chosen particle-size distribution can be made essentially space filling, both in the crystalline solid (Appolonian packing) and in the liquid. In practice, perfect space filling structures are never achieved because this requires an infinite number of (predominantly small) particles per unit volume. Here, we consider a somewhat simpler problem, namely the filling of a given volume $`V`$ by a fixed number of particles $`N`$, that occupy a prescribed total volume fraction $`\eta `$. We assume that the particles are free to exchange volume. As we have fixed both the number and the total volume of the particles, the average volume per particle is fixed—it defines the natural length-scale in the model. Clearly, the Helmholtz free energy of the system will depend on the non-fixed particle-size distribution. The distribution, however, is restricted by the two constraints of fixed number of particles and fixed total volume. We define the optimal size distribution to be the one that minimizes the Helmholtz free energy under both constraints. In the Monte Carlo simulations (performed in the isothermal-isobaric ensemble) that we report in this paper, we study the density dependence of the particle-size distribution. We compare the simulation results for the size distribution with an analytical estimate that is obtained by solving the Percus-Yevick (PY) equation for an $`N`$-component hard-sphere mixture. Not surprisingly, the PY equation works very well at low densities. However, the theory breaks down at a surprisingly low density ($`\eta `$0.26). Of course, the fact that an approximate theory fails at a given density, does not imply that there is anything special going on in the system at that density. Yet, our simulations indicate that there is—the size distribution that was initially uni-modal, becomes bi-modal. We present a theoretical argument supporting this scenario.

The remainder of this paper is organized as follows: in section II we describe the constant-pressure Monte Carlo simulations. The Percus-Yevick expression for the free energy of a system of polydisperse hard spheres with variable size distribution is discussed in section III and IV. Section V is devoted to the derivation of analytical results turning useful when interpreting the simulation data. The mechanism behind the transition from uni-modal to bi-modal size distribution is also discussed.

## II Monte Carlo Simulations

Monte Carlo simulations were performed in the isothermal-isobaric (constant $`NPT`$) ensemble . This means that the number of particles $`N`$, the pressure $`P`$ and the temperature $`T`$ are fixed. We attempt three distinct types of trial moves. We change the positions of the particles and allow the volume of simulation box to fluctuate, in order to equilibrate with respect to the applied pressure. Since we do not expect any crystalline order at low pressures, a cubic box shape is maintained. The third type of move is the one related to sampling the polydispersity of the system. To this end, we select two particles at random, between which we exchange an amount of volume drawn uniformly from the interval $`[\mathrm{\Delta }V_{max},\mathrm{\Delta }V_{max}]`$ (Fig. 1). The maximum volume change $`\mathrm{\Delta }V_{max}`$ was chosen such that the acceptance of a volume exchange move is between 35 and 50 %. The relative frequency of the three moves is given by $`N:1:N/2`$. The initial configurations are made by $`N`$ monodisperse spheres on a simple cubic lattice.

Simulations were performed for system sizes of $`N`$=512 or 1000, and several different reduced pressures $`P^{}=(k_BT)^1P\sigma ^3`$, where $`\sigma `$ is the diameter of particles, and $``$ denotes an average over particle-size distribution. As $`\sigma ^3`$ is fixed, we choose $`\sigma ^3^{1/3}`$ to define the unit of length that we will use in the remainder of this paper. We use $`k_BT`$ as our unit of energy. All other units that we need, follow from these definitions. The equation of state \[$`P^{}`$ as a function of $`\rho ^{}(\rho \sigma ^3)`$\] and the particle size distribution function were determined in the simulations. The results for the equation of state are shown in Fig. 2. The simulation data have been collected in Table I. At low pressures the particle-size distribution function is a single-peaked function with its maximum at $`v=0`$ (Figs. 3, 4 and 5). At higher pressures (typically, $`P^{}>2.0`$), the particle-size distribution develops a second peak. Actually, this second peak is quite small (i.e. only a small fraction of all particles becomes “large”). However, these particles contribute appreciably to the total volume fraction (Fig. 6). Depending on the pressure this contribution can get as large as 75%.

The formation of big particles in these MC simulations is a rather slow “dynamical” process. In order to speed up calculations, we did additional simulations, in which we started with a bi-disperse distribution, with one or several big particles containing 99% of the total volume occupied by the spheres, surrounded by a sea of small particles containing the remaining volume. In the 512 particle system, only one or two big particles remain for the lower pressures ($`P^{}=2.5`$ and $`P^{}=3.0`$). For higher pressures the number of big particles can stabilize at higher values as well. For the 1000 particle system the maximum number of big particles observed at the lower pressures is three. In addition, the size of these big particles is not the same. It is not clear whether this suggests a further possible fractionation or that it is a consequence of the slow equilibration and that one or more of the big particles are still shrinking.

Below, we discuss these simulation results in the context of the relevant theoretical predictions, but first let us stop to make certain considerations on the ideal entropy associated with this system.

## III Ideal entropy of a polydisperse system

Strictly speaking, the ideal entropy of a polydisperse system is infinite . In a multicomponent system such an entropy is exactly given by

$$Nk_B\underset{i}{}w_i\mathrm{ln}(\mathrm{\Lambda }_i^3\rho w_i),$$
(1)   
where $`w_i`$ is the molar fraction of species $`i`$ and $`\mathrm{\Lambda }_i`$ its thermal wavelength. The usual path towards the entropy of a polydisperse system  (or for that matter, towards the entropy of continuous signals in Information Theory ) is to classify the species into “boxes” according to a certain property which distinguishes them (diameter, volume, molecular weight, etc). If $`x`$ denotes such a property species $`i`$ will denote the box having particles with $`x`$ between $`i\mathrm{\Delta }x`$ and $`(i+1)\mathrm{\Delta }x`$, for a given $`\mathrm{\Delta }x`$ which defines the boxes. If $`W(x)`$ denotes the probability density of a particle having the value $`x`$ for that property, then $`w_i=W(x_i)\mathrm{\Delta }x`$ where $`x_i`$ is a typical value of the $`i`$th box. Then Eq. (1) adopts the form

$$Nk_B\underset{i}{}W(x_i)\mathrm{\Delta }x\mathrm{ln}(\mathrm{\Lambda }_i^3\rho W(x_i)/\mathrm{\Delta }x).$$
(2)   
The entropy of the polydisperse system is obtained from Eq. (2) taking the limit $`\mathrm{\Delta }0`$, and we can clearly see that, besides obtaining the usual expression 

$$S_{\mathrm{poly}}=NK_B𝑑xW(x)\mathrm{ln}[\mathrm{\Lambda }(x)^3W(x)],$$
(3)   
there is also a divergent $`\mathrm{ln}\mathrm{\Delta }x`$, which is simply taken as a “constant” ignorable term.

But expression (3) is not well defined. Suppose we simply change coordinates to label the species from $`x`$ to $`y`$ (say, from the diameter to the volume). Then the probability density in the new variable will be $`\stackrel{~}{W}(y)=W(x)\left|\frac{dx}{dy}\right|`$. It is straightforward to show that in the new labeling the entropy becomes

$$S_{\mathrm{poly}}=NK_B𝑑y\stackrel{~}{W}(y)\mathrm{ln}[\mathrm{\Lambda }(x(y))^3\stackrel{~}{W}(y)]+Nk_B𝑑y\stackrel{~}{W}(y)\mathrm{ln}\left|\frac{dy}{dx}\right|,$$
(4)   
which is different of what we would have obtained had we began with the labeling $`y`$.

This is a well known fact in Information Theory . In the study of fluid equilibria of quenched polydisperse systems this fact turns out to be irrelevant because the new term simply adds the same constant to both sides of the equilibrium equations . However when studying annealed polydispersity this results tells us that the labeling is crucial and has to be dictated by the physical process underlying the polydispersity. In our case the Monte Carlo movements described in Section II are a large scale description of a hypothetical microscopic system of tiny particles of exactly the same size distributed among $`N`$ aggregates of a variable number of particles. The constant volume constraint would correspond to the conservation of the number of tiny particles, and the natural labeling of the aggregates would be the number of tiny particles which form it. As this number is proportional to the volume of the aggregate, in the continuum description is the volume, instead of the diameter, what turns out to be the natural labeling.

Notice that we could have described another model in which the tiny particles aggregated into spherical surfaces. In that case it would be the total surface what would be preserved and the natural labeling of aggregates would be their respective surfaces. As we will discuss in the conclusions, the physics of this model would be completely different.

## IV Percus-Yevick Theory

The Percus-Yevick equation for an $`n`$-component hard-sphere mixture can be solved analytically, to yield the following equation of state :

$$\frac{\pi }{6}P^{}=\frac{\xi _0}{1\xi _3}+\frac{3\xi _1\xi _2}{(1\xi _3)^2}+\frac{3\xi _2^3}{(1\xi _3)^3}.$$
(5)   
The $`j`$-th moment of the particle-size distribution $`\xi _j`$ is defined as

$$\xi _j=\frac{\pi }{6}\underset{i}{}\rho _i\left(\frac{6v_i}{\pi }\right)^{j/3},$$
(6)   
where $`\rho _i=N_i/V`$, the index $`i`$ is used to denote the different particle species, and $`v_i`$ is the volume of the $`i`$-th species.

Equation (5) is also valid for a continuous size distribution, in which case the sum in Eq. (6) is replaced by an integral. The corresponding expression for the chemical potential of a species with radius $`R`$ is

$$\mu ^{}=\mathrm{ln}\left[\rho \mathrm{\Lambda }^3W(v)\right]\mathrm{ln}(1\xi _3)+\frac{6\xi _2R}{(1\xi _3)}+\frac{12\xi _1R^2}{(1\xi _3)}+\frac{18\xi _2^2R^2}{(1\xi _3)^2}+\frac{4\pi }{3}P^{}R^3$$
(7)   
where $`\mathrm{\Lambda }`$ is the de-Broglie thermal wavelength, $`\sqrt{h^2/(2\pi mkT)}`$, and $`W(v)`$ is the probability density to find a particle with a volume around $`v=(4\pi /3)R^3`$. The pressure $`P^{}`$ is given by Eq. (5).

In an (NPT) description, the Gibbs free energy of the system fulfilling the constraints, must be at a minimum. The conservation of the number of particles and of the solid volume fraction, imply that $`W(v)`$ must be of the form:

$$W(v)=\mathrm{exp}\left\{\underset{i=0}{\overset{3}{}}\alpha _iR^i\right\},$$
(8)   
where

$`\alpha _1`$ $`=`$ $`{\displaystyle \frac{6\xi _2}{1\xi _3}},`$ (9)   
$`\alpha _2`$ $`=`$ $`12({\displaystyle \frac{\xi _1}{1\xi _3}}+{\displaystyle \frac{3\xi _2^2}{2(1\xi _3)^2}}).`$ (10)   
The coefficients $`\alpha _0`$ and $`\alpha _3`$ are determined by the constraints that the number of particles and the solid volume fraction are fixed. Note that all $`\xi _i`$ $`(i=1,2,3)`$ are positive. Moreover, $`\xi _3`$ is equal to the volume fraction $`\eta `$, and is therefore necessarily less than one. Hence, $`\alpha _1`$ and $`\alpha _2`$ are always negative. The last coefficient, $`\alpha _3`$, should be negative or zero, because otherwise the particle-size distribution cannot be normalized. Since $`\alpha _1`$, $`\alpha _2`$ and $`\alpha _3`$ are always negative the Percus-Yevick equation predicts that $`W(v)`$ is a monotonically decreasing function of $`v`$. This implies that the size-distribution given by Eq. (8) can never be bi-modal.

Note that these conclusions also hold for the more accurate equation of state of Mansoori et al. . This equation adds an extra term to the pressure given in Eq. (5) depending on $`\xi _3`$. Thus Eq. (7) is the same with the new expression for $`P`$—which turns out to be irrelevant because the $`R^3`$ term is controlled by the Lagrange multiplier associated with the constraint on the total solid particle volume. In other respect, the analysis for the Mansoori equation-of-state is identical to that for PY.

In practice, we solve Eq. (8) numerically. To this end, we represent $`W(v)`$ as a histogram. Initially, the value of $`W(v)`$ in the different bins is assigned an arbitrary non-negative value, compatible with the constraint that $`W(v)`$ is normalized and that $`v`$ is fixed. We fix the density at the desired value. We determine the optimal $`W(v)`$ using the following scheme: we select a bin (say $`i`$) at random and change the value of $`W(v_i)`$ by a random amount $`\mathrm{\Delta }W`$, distributed uniformly in the interval $`[\mathrm{\Delta }W_{max},\mathrm{\Delta }W_{max}]`$. We first check if the new value $`W(v_i)`$ is non-negative. If it is, we satisfy the constraints by scaling the width of all bins and the height of the function by two appropriately chosen factors. We then compute all moments $`\xi _j`$, the pressure and the free energy, and we check if the Helmholtz free energy is smaller than the previous one. If it is, we accept the new value for $`W(v_i)`$, otherwise we reject it. We repeat the procedure until the free energy no longer decreases. We have verified that $`W(v)`$ is indeed of the form given by Eqs. (8) through (10). Figure 4 shows a comparison of the PY estimate for $`W(v)`$, determined in this way, with the results of the full Monte Carlo simulations. We find that $`\alpha _3`$ is a monotonically increasing function of density. A comparison between simulation and PY theory for the equation of state is shown in Fig. 2. Note that, in this figure, the PY solution terminates at pressure $`P_c1.34`$ (the cross in Fig. 2). This is the point where $`\alpha _3`$ becomes zero. Beyond this point we can no longer find a solution for $`W(v)`$ that is of the form given by Eq. (8). In Appendix A, we consider the breakdown of the PY theory in more detail and obtain the packing fraction beyond which the PY approximation breaks down: $`\eta _c=0.260`$, the corresponding pressure being $`P_c^{}=1.343`$. This breakdown of the PY equation at a relatively low density is surprising, as the PY equation works well up to quite high densities for fixed particle size distributions . That it breaks down regardless the accuracy of the equation of state can be inferred from the fact that Mansoori et al.’s equation of state undergoes exactly the same breakdown, though for slightly different values of $`\eta `$ and $`P`$. Besides, from the analysis that we have carried out, it can be seen that a similar breakdown will appear for any other theory yielding an equation of state depending only on $`\xi _i`$, $`i=0,1,2,3`$.

## V Analytical results

In this section, we derive a theoretical bound for the pressure of the polydisperse system, providing the equation of state at the packing fraction where the size distribution becomes bi-modal. To this end, we work in the ($`NVT`$) ensemble and take advantage of the extremality of the Helmholtz free energy under the constraints of constant $`N`$ and $`\eta `$: the “grand potential”

$$=\{W\}_0W(v)𝑑v_1vW(v)𝑑v,$$
(11)   
where $`_0`$ and $`_1`$ are Lagrange multipliers, has to be minimum for the optimal size distribution. In the above relation, the free energy functional $``$ can be cast into the usual ideal and excess contributions

$$\{W\}=Nk_BT𝑑vW(v)\left[\mathrm{ln}\left(\mathrm{\Lambda }^3\rho W(v)\right)1\right]+_{\text{excess}}\{W\}.$$
(12)   
We attempt the following change in the system: the volume of a given particle $`v_0`$ is increased by an amount $`\delta v_0`$, before a rescaling of all volumes by a factor $`\lambda `$ (e.g. $`v\lambda v`$) such that the overall volume change vanishes. This imposes:

$$\lambda =1\frac{\delta v_0}{Nv}+𝒪\left([\delta v_0]^2\right).$$
(13)   
The effect of the expansion of particle $`v_0`$ on the size-distribution can be written:

$$\delta W(v)=\frac{1}{N}\left[\delta \left(vv_0\delta v_0\right)\delta \left(vv_0\right)\right],$$
(14)   
where $`\delta (\mathrm{})`$ denotes the Dirac distribution. The scaling procedure affects $`W`$ according to

$`\delta W(v)`$ $`=`$ $`{\displaystyle \frac{1}{\lambda }}W\left({\displaystyle \frac{v}{\lambda }}\right)W(v)`$ (15)   
$`=`$ $`{\displaystyle \frac{\delta v_0}{Nv}}{\displaystyle \frac{d[vW]}{dv}}+𝒪\left([\delta v_0]^2\right).`$ (16)   
The corresponding variation of the ideal contribution to $``$ reads

$$\delta _{\text{id}}=𝑑v\frac{\delta _{\text{id}}\{W\}}{\delta W(v)}\delta W(v),$$
(17)   
with the functional derivative

$$\frac{\delta _{\text{id}}\{W\}}{\delta W(v)}=Nk_BT\mathrm{ln}[\mathrm{\Lambda }^3W(v)].$$
(18)   
We then get the entropic term

$`\delta _{\text{id}}`$ $`=`$ $`k_BT{\displaystyle 𝑑v\mathrm{ln}[\mathrm{\Lambda }^3W]\left[\delta \left(vv_0\delta v_0\right)\delta \left(vv_0\right)+\frac{\delta v_0}{v}\frac{d[vW]}{dv}\right]}`$ (19)   
$`=`$ $`k_BT\left\{{\displaystyle \frac{W^{}(v_0)}{W(v_0)}}+{\displaystyle \frac{1}{v}}\right\}\delta v_0,`$ (20)   
where $`W^{}`$ is the derivative of $`W`$.

The variation of the excess free energy reduces to the reversible work needed to perform the transformation, and is derived in appendix B:

$$\frac{\delta W_{\text{rev}}}{\delta v_0}=\rho k_BT𝑑vW(v)g\left(\frac{\sigma _0+\sigma }{2}\right)\left(1+\frac{\sigma }{\sigma _0}\right)^2\frac{P_{\text{excess}}}{\eta },$$
(21)   
where $`g(\sigma _0/2+\sigma /2)`$ denotes the radial distribution function evaluated at contact between species of diameters $`\sigma _0`$ (having volume $`v_0`$) and $`\sigma `$ (having volume $`v`$). When $`\sigma _0\sigma `$, we can replace the density at the surface of particle $`v_0`$ by that at a planar wall, and Eq. (21) becomes

$$\delta W_{\text{rev}}=\left[P\frac{P_{\text{excess}}}{\eta }\right]\delta v_0.$$
(22)   
In this limit $`v_0v`$

$$\frac{\delta }{\delta v_0}=\frac{\delta }{\delta v_0}=k_BT\left\{\frac{W^{}(v_0)}{W(v_0)}+\frac{1}{v}\right\}+P\frac{P_{\text{excess}}}{\eta }.$$
(23)   
For the optimal size distribution, $`\delta F`$ vanishes so that

$$k_BT\frac{W^{}(v_0)}{W(v_0)}+\frac{1\eta }{\eta }\left[\frac{2\rho k_BT}{1\eta }P\right]=0.$$
(24)   
Assuming $`W(v)`$ to be a normalizable distribution, $`W^{}(v)`$ has to be negative for large arguments, which sets the upper bound:

$$P<\mathrm{\hspace{0.17em}2}\frac{\rho k_BT}{1\eta }\text{or}P^{}<\frac{12\eta }{\pi (1\eta )}$$
(25)   
for the rescaled pressure. For low packing fractions ($`\eta <\eta _c`$) where the PY solution is available, the above inequality is fulfilled (Fig. 2). At the threshold $`\eta =\eta _c`$ where the second polydispersity peak appears, $`W^{}`$ changes sign which means

$$P_c^{}=\frac{12\eta _c}{\pi (1\eta _c)}.$$
(26)   
The above relation is remarkably well obeyed within the PY approximation (see the data of section IV or figure 2: the PY expression crosses the line given by Eq. (25) exactly at $`\eta _c`$). For $`\eta >\eta _c`$, the upper bound is violated by the simulation results reported in Table I and figure 2. However, the data suggest non negligible finite-size effects: increasing $`N`$ shifts the pressure closer to the theoretical bound. Besides, starting from bi-disperse initial conditions (cf the procedure described in section II), supposed to be closer to the expected polydispersity, has the same effect. According to expressions (23) and (24), the violation of Eq. (25) means that $`\delta =\delta <0`$ for $`\delta v_0>0`$, so that the biggest particle tends to expand. Its growth is however necessarily limited by the length $`L`$ of the simulation box. This is supported by the observation that, even for the largest system investigated ($`N=10^3`$ particles), the size of the biggest particle obtained is determined by $`L`$ ($`\sigma _{\text{biggest}}>L/3`$, irrespective of the packing fraction, see for example figure 6). This suggests that system sizes that would presumably allow the system to reach thermodynamic equilibrium (and fulfill inequality (25)) are numerically out of reach. Consequently, the question of the extensivity/intensivity of the number of large particles cannot be addressed by simulations; a theoretical investigation seems to require the detailed knowledge of the interfacial free energy between “large” and “small” species. At this stage, we cannot tell whether a true phase transition is associated with the occurrence of the second peak in the particle-size distribution. The available data certainly do not rule out this possibility.

Finally, the integration of Eq. (24) yields the tail of the optimal size-distribution:

$$\mathrm{ln}W(v)v/v\text{for}vv.$$
(27)   
Equation (8) obeys this relation, which cannot be tested against simulation results because of the lack of statistics for very large particles (not more than 5 in a typical run).

## VI Conclusion

At this stage, we can only speculate what will happen at larger $`N`$ and/or larger densities. Conceivably, once the volume-fraction of the large particles exceeds a certain threshold, proliferation of still larger particles can occur, and so on, until eventually an ”Appolonian” packing of the liquid is achieved. The theoretical analysis of this scenario is non-trivial, as the small particles now induce attractive depletion forces between the large particles. Unfortunately, the systems that we can conveniently study by simulation are too small to allow us to investigate this regime.

We stress that the specific model we have chosen to study is somewhat arbitrary. For instance, rather than fixing the number of particles, one might have chosen to fix the total surface area of the particles. The latter constraint would be logical if one aims to model the size distribution of droplets covered with a fixed amount of surfactant. In addition, we assume that the surface free energy of the spheres is negligible. Again, this constraint can be removed. We hope that the rather surprising results of the present study will stimulate research into these related models.

## Acknowledgments

We thank Bela Mulder for critical discussions. The work of the FOM Institute is part of the research program of FOM and is made possible by financial support from the Netherlands Organization for Scientific Research (NWO). Junfang Zhang acknowledges financial support from Chinese government. José Cuesta participated in this work during a stay at AMOLF financed by the project no. PR95-558 from the Dirección General de Enseñanza Superior (DGES). His work is also part of the DGES project no. PB96–0119. He also wants to thank the hospitality found at AMOLF, specially from Daan Frenkel’s and Bela Mulder’s groups.

## A

Let us consider the range of densities where a solution of the Percus-Yevick equation is possible. As stated in section IV, it is essential that $`\alpha _3`$ be non-positive. Hence, the pressure at which $`\alpha _3=0`$ defines the end point of the theory. To locate this point, consider the form of the solution at the point where $`\alpha _3=0`$. Then the distribution function reduces to

$`W(v)=\mathrm{exp}\left\{{\displaystyle \underset{i=0}{\overset{2}{}}}\alpha _iR^i\right\},`$where $`R`$ hereafter denotes a reduced radius measured in units of $`\sigma ^3^{1/3}`$. Using the two constraints for normalization and for the average volume of the particles, we can express the coefficients $`\alpha _0`$ and $`\alpha _2`$ in terms of $`\alpha _1`$. If we combine Eqs. (9) and (10) to eliminate $`\xi _3`$, we obtain

$`\alpha _2`$ $`=`$ $`2\alpha _1{\displaystyle \frac{\xi _1}{\xi _2}}{\displaystyle \frac{1}{2}}\alpha _1^2`$ (A1)   
$`=`$ $`\alpha _1{\displaystyle \frac{R}{R^2}}{\displaystyle \frac{1}{2}}\alpha _1^2,`$ (A2)   
where the ratio of the moments $`\xi _1`$/$`\xi _2`$ does not depend explicitly on the density $`\rho `$ or on $`\alpha _0`$, but it only depends on $`\alpha _1`$ itself. In the second line we have used the definition

$$R^i=_0^{\mathrm{}}R^iW(v)4\pi R^2𝑑R.$$
(A3)   
But we have another relation between the moments of the particle-size distribution: partial integration of $`R^nW(v)𝑑v`$ yields

$$\left[R^{n+1}W(v)\right]_0^{\mathrm{}}=_0^{\mathrm{}}\left(2\alpha _2R^{n+2}+\alpha _1R^{n+1}+(n+1)R^n\right)W(v)𝑑R=0$$
(A4)   
for $`n0`$. This leads to the identities

$`R^3`$ $`=`$ $`{\displaystyle \frac{\alpha _1}{2\alpha _2}}R^2{\displaystyle \frac{4}{2\alpha _2}}R,`$ (A5)   
$`R^2`$ $`=`$ $`{\displaystyle \frac{\alpha _1}{2\alpha _2}}R{\displaystyle \frac{3}{2\alpha _2}}R^0.`$ (A6)   
This allows us to write $`R`$ and $`R^2`$ as a function of $`\alpha _1`$, $`\alpha _2`$, $`R^0`$ and $`R^3`$; i.e.

$`R`$ $`=`$ $`{\displaystyle \frac{3\alpha _1R^0+4\alpha _2^2R^3}{\alpha _1^28\alpha _2}},`$ (A7)   
$`R^2`$ $`=`$ $`{\displaystyle \frac{12R^02\alpha _1\alpha _2R^3}{\alpha _1^28\alpha _2}}.`$ (A8)   
But $`R^0=1`$ and $`R^3=1/8`$; hence

$$\frac{R}{R^2}=\frac{3\alpha 1+\frac{1}{2}\alpha _2^2}{12\frac{1}{4}\alpha _1\alpha _2}.$$
(A9)   
Substituting this expression in Eq. (A1), we can eliminate $`R/R^2`$ to obtain

$$\alpha _2+\frac{1}{2}\alpha _1^2=\alpha _1\frac{3\alpha 1+\frac{1}{2}\alpha _2^2}{12\frac{1}{4}\alpha _1\alpha _2}.$$
(A10)   
For what follows, it is convenient to introduce two new variables $`f`$ and $`c`$

$`\alpha _1`$ $`=`$ $`fc,`$ (A11)   
$`\alpha _2`$ $`=`$ $`f^2;`$ (A12)   
then

$$R^i=\frac{1}{f^{i+3}}R^i_{f=1}.$$
(A13)   
We can use this relation to express $`f`$ as function of $`c`$. We use the fact that the average volume per particle is fixed, to rewrite

$`{\displaystyle \frac{R^3}{R^0}}`$ $`=`$ $`{\displaystyle \frac{1}{8}}={\displaystyle \frac{1}{f^3}}{\displaystyle \frac{R^3_{f=1}(c)}{R^0_{f=1}(c)}}`$ (A14)   
$``$ $`{\displaystyle \frac{1}{f^3}}Y(c),`$ (A15)   
where the second line defines the function $`Y(c).`$ Hence

$$f=2Y^{1/3}(c)$$
(A16)   
As $`Y(c)`$ can be expressed explicitly in terms of error functions, we now know $`f`$ as a function of $`c`$. Equation (A11) allows us to express both $`\alpha _1`$ and $`\alpha _2`$ as explicit functions of $`c`$. We can then use Eq. (A10) to determine $`c`$ numerically. We find that this equation has a unique solution. Once the value of $`c`$ has been determined, we know $`\alpha _1,\alpha _2`$ and $`\alpha _0`$ (from the normalization condition). Eq. (9) finally yields the packing fraction beyond which the PY approximation breaks down: $`\eta _c=0.260198`$. The corresponding pressure is $`P_c^{}=1.343442`$.

## B

We first note that the reversible work done by an operator rescaling both particle volumes ($`v_i\lambda v_i,i`$) and container volume ($`VV^{}=\lambda V`$), is

$$\delta W_1=P_{\text{ideal}}\delta V=\rho k_BT\delta V,\text{where}\delta V=(\lambda 1)V.$$
(B1)   
In the transformation, the total volume $`V_p`$ of the particles changes according to

$$\frac{\delta V_p}{V_p}=\frac{\delta V}{V}\delta V_p=\eta \delta V.$$
(B2)   
Keeping the particle volumes fixed and going back to the original container volume ($`V^{}V^{}/\lambda `$) requires the reversible work

$$\delta W_2=P(VV^{})=P\delta V.$$
(B3)   
It is then straightforward to obtain the work associated with a rescaling of particle volumes at constant accessible volume $`V`$:

$$\delta W_{v\lambda v}=\delta W_1+\delta W_2=\frac{P_{\text{excess}}}{\eta }\delta V_p,$$
(B4)   
valid for all polydispersities $`W(v)`$.

In the remainder, the shall derive the work needed to grow a particle $`v_0`$ by an amount $`\delta v_0`$ ($`\delta v_0=\pi \sigma _0^2\delta \sigma _0/2`$). We assume the normalization $`W𝑑v=1`$ to hold. Consider species having volumes between $`v`$ and $`v+\delta v`$ (diameters between $`\sigma `$ and $`\sigma +\delta \sigma `$). They exert a pressure $`\rho k_BTW(v)dvg(\sigma _0/2+\sigma /2)`$ on particle $`v_0`$, involving the radial distribution function at contact between species $`\sigma _0`$ and $`\sigma `$. For the above pair $`(\sigma _0,\sigma )`$, the excluded volume sphere has diameter $`\sigma _0+\sigma `$, and sweeps a volume

$$\delta V_{\text{sweep}}=\pi (\sigma _0+\sigma )^2\frac{\delta \sigma _0}{2}=\left(1+\frac{\sigma }{\sigma _0}\right)^2\delta v_0$$
(B5)   
during the growth of particle $`v_0`$. Summing over all species $`v`$, the work performed by the operator takes the form

$$\delta W_{\text{growth }v_0}=\rho k_BT\delta v_0𝑑vW(v)g\left(\frac{\sigma _0+\sigma }{2}\right)\left(1+\frac{\sigma }{\sigma _0}\right)^2.$$
(B6)   
For the size modification considered in section V, the global volume change of the particles vanishes, such that $`\delta V_p+\delta v_0=0`$. Summing the contributions arising from Eqs. (B4) and (B6), the results of Eq. (21) is recovered:

$$\frac{\delta W_{\text{rev}}}{\delta v_0}=\rho k_BT𝑑vW(v)g\left(\frac{\sigma _0+\sigma }{2}\right)\left(1+\frac{\sigma }{\sigma _0}\right)^2\frac{P_{\text{excess}}}{\eta }.$$
(B7)   
Note that a similar argument can be invoked to compute the reversible work needed to rescale all particle diameters, at constant $`V`$. After integration of Eq. (B6) over all $`W(v_0)dv_0`$, we obtain:

$$\delta W_{v\lambda v}=\rho k_BT𝑑v𝑑v^{}W(v)W(v^{})g\left(\frac{\sigma +\sigma ^{}}{2}\right)(\sigma +\sigma ^{})^2\sigma \frac{\delta V_p}{\sigma ^3}.$$
(B8)   
Inserting this result into Eq. (B4) provides the equation of state for a polydisperse fluid of hard spheres

$$\frac{P}{\rho k_BT}=1+\eta 𝑑v𝑑v^{}W(v)W(v^{})g\left(\frac{\sigma +\sigma ^{}}{2}\right)(\sigma +\sigma ^{})^2\frac{\sigma }{\sigma ^3}.$$
(B9)   
For a monodisperse fluid, $`W(v)=\delta (vv_0)`$ and we recover the well known relation

$$\frac{P}{\rho k_BT}=1+4\eta g(\sigma _0).$$
(B10)   
FIGURE CAPTIONS

1. Schematic drawing of the Monte Carlo trial used to sample the polydispersity. Of two randomly chosen particles the volume of one is increased, while the other is decreased in volume by the same amount.
2. The equation of state of the polydisperse system. The solid line is the PY prediction, which can not be extended beyond the cross. The circles correspond to a system initially prepared monodisperse, while the squares are final values in which we initially started with one big particle. The solid squares and the open symbols are from a 512 and 1000 particle system respectively. The dashed line is the upper bound (26) for the pressure.
3. The Monte Carlo results for the distribution of particle volumes, $`W(v)`$, as a function of $`v`$, for several reduced pressures for a 1000 particle system. For $`P^{}>2.0`$ the distribution develops a second peak (with statistical noise) at much larger volumes (note the change in scale on the right-hand side of the figure). Although there are only one or several of these big particles, they can contribute over 30% of the total volume of all particles. The diameters can get larger than a third of the length of the simulation box.
4. Comparison of the numerical results (dashed line) for the particle-size distribution $`W(v)`$, with the corresponding prediction of the PY theory (solid line). These results were obtained at a relatively low pressure ($`P^{}=0.5`$), which corresponds to a volume fraction $`\eta =0.151`$. Note that at this density, simulation and theory are in quite good agreement.
5. Snapshot of a typical configuration at reduced pressure $`P^{}=1.0`$ and volume fraction $`\eta =0.244`$.
6. Snapshot of a typical configuration at reduced pressure $`P^{}=3.0`$ and volume fraction $`\eta =0.425`$. In this case there are two big particles.