# Freezing line of polydisperse hard spheres via direct-coexistence simulations

Antoine Castagnède, Laura Filion, Frank Smallenburg

► To cite this version:

Antoine Castagnède, Laura Filion, Frank Smallenburg. Freezing line of polydisperse hard spheres via direct-coexistence simulations. The Journal of Chemical Physics, 2025, 163 (13), pp.134121. 10.1063/5.0281621. hal-05380366

HAL Id: hal-05380366

https://hal.science/hal-05380366v1

Submitted on 24 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Freezing line of polydisperse hard spheres via direct-coexistence simulations

Antoine Castagnède,<sup>1, a)</sup> Laura Filion,<sup>2</sup> and Frank Smallenburg<sup>1</sup>

<sup>1)</sup>Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France

<sup>2)</sup>Soft Condensed Matter and Biophysics, Debye Institute for Nanomaterials Science, Utrecht University, Utrecht, Netherlands

In experimental systems, colloidal particles are virtually always at least somewhat polydisperse, which can have profound effects on their ability to crystallize. Unfortunately, accurately predicting the effects of polydispersity on phase behavior using computer simulations remains a challenging task. As a result, our understanding of the equilibrium phase behavior of even the simplest colloidal model system, hard spheres, remains limited. Here, we present a new approach to map out the freezing line of polydisperse systems that draws on direct-coexistence simulations in the semi-grand canonical ensemble. We use this new method to map out the conditions where a hard-sphere fluid with a Gaussian size distribution becomes metastable with respect to partial crystallization into a face-centered-cubic crystal. Consistent with past predictions, we find that as the polydispersity of the fluid increases, the coexisting crystal becomes increasingly size-selective, exhibiting a lower polydispersity and larger mean particle size than the fluid phase. Finally, we exploit our direct-coexistence simulations to examine the characteristics of the fluid-crystal interface, including the surface stress and interfacial absorption.

# I. INTRODUCTION

When studying colloidal self-assembly, we often assume that we are working with systems that are either monodisperse, or consist of a limited number of distinct species, each with a well-defined size, shape, and other characteristics. In practice, however, colloidal particles inevitably exhibit a finite degree of polydispersity: variations of the particle characteristics around their mean value. This is particularly crucial in crystallization studies, where polydispersity can hinder the formation of well-ordered structures, or even induce formation of complex crystal phases¹⁻⁴.

Unfortunately, we still have a limited understanding of the effect of polydispersity on the phase behavior of even the simplest of colloidal model systems: hard spheres. Colloidal hard spheres are arguably one of the most fundamental model systems in soft-matter physics, as their simplicity makes them extremely amenable for study in theory, simulations, and experiments⁵. As a result, significant efforts have been made to explore the effects of size polydispersity on e.g. the dynamics⁶⁻⁸, phase behavior³,⁹⁻¹⁵ and crystallization³,⁴,¹⁶⁻²³ of hard spheres, as well as the behavior of polydisperse hard-sphere mixtures near a wall²⁴,²⁵. Despite these efforts, our knowledge of the equilibrium crystallization behavior is still incomplete. Specifically, simulation studies that focus on drawing equilibrium phase boundaries for polydisperse hard-sphere mixtures⁹⁻¹⁵ typically consider only crystallization into face-centered cubic (FCC) crystals – the crystal phase that is known to be stable for monodisperse hard spheres. However, a number of recent simulation studies have revealed that at relatively high polydispersities (≥ 10%) more complex crystal structures can spontaneously form, including Laves phases²⁻⁴, whose stability has not yet been explored in detail.

The reason for our lack of clear results for this fundamental model system is the fact that determining phase boundaries for polydisperse systems is notoriously difficult¹⁵. In polydisperse systems coexisting phases can have significantly different distributions of the relevant particle properties, a phenomenon known as fractionation. As a result, determining phase boundaries requires considering a vast space of possible mixtures of particles with different size distributions, significantly complicating any computational methodology.

Past studies have addressed this problem in a variety of ways, often making use of semi-grand canonical ensemble simulations, where the particles can fluctuate in size but the total number of particles stays fixed¹¹,¹³,¹⁴,²⁶. Using this ensemble, methodologies such as Gibbs-Duhem integration¹¹,¹²,¹⁴ or phase-switch Monte Carlo¹³,²⁷,²⁸ provide a route towards locating coexisting state points. Alternatively, theoretical methods can be used to estimate phase boundaries based on free-energy calculations that build on approximate equations of state for the fluid and crystal phases⁹,¹⁰. These methods have given us a solid picture of the phase behavior of hard spheres in the low-polydispersity regime, but often struggle at higher polydispersities.

Here, we propose a new approach for predicting fluid-solid phase coexistences based on direct-coexistence simulations in the semi-grand canonical ensemble. In particular, we build on our recent work on direct-coexistence simulations for monodisperse systems²⁹, and combine it with efficient event-driven simulations in which particles can fluctuate in size based on an externally imposed relative chemical potential. Using this approach, we map out the fluid-FCC freezing line (i.e. the fluid cloud curve) as a function of polydispersity, and examine the behavior of the particle size distribution across the interface. Our methodology can be readily extended to other model systems and more complex crystal phases, opening the door to new experiments.

<footnote>a) Author to whom correspondence should be addressed: antoine.castagnede@universite-paris-saclay.fr</footnote>
to a more complete exploration of polydisperse phase behavior in the future.

## II. METHODS

The main aim of this work is to determine the thermodynamic conditions for the fluid-solid phase coexistence of a polydisperse hard-sphere system, where the particle size distribution in the fluid phase is characterized by a Gaussian distribution

$$ \mathcal{P}(\sigma) = \frac{1}{p\bar{\sigma}\sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\sigma - \bar{\sigma}}{p\bar{\sigma}}\right)^2\right). \quad (1) $$

Here, $p$ is the polydispersity of the parent fluid phase, defined as the ratio of the standard deviation of the particle diameter $\sigma$ to its mean value $\bar{\sigma}$. Note that fixing the size distribution of the fluid phase does not fix the overall size distribution in the coexisting system. Since the crystal phase is generally expected to have a different size distribution (and polydispersity) from the fluid, the overall size distribution of any coexisting state also depends on the fraction of volume taken up by each phase. In this work, we specifically focus on the identification of the freezing line of a fluid with fixed size distribution, also known as the fluid cloud curve⁹. The corresponding set of coexisting crystal phases is known as the shadow curve.

When two polydisperse phases coexist, thermodynamic equilibrium requires that their temperature $T$, pressure $P$, and absolute chemical potential $\mu(\sigma)$ are equal in the two phases. The core idea of the strategy we propose here is to look at a series of “trial” fluids at different densities, and to check whether there is a polydisperse crystal phase which satisfies these criteria. To narrow down the search for this crystal phase, we note that a weaker (and insufficient) requirement for such a phase is that it has the same temperature, pressure, and relative chemical potential $\delta\mu(\sigma)$. Here, $\delta\mu(\sigma) = \mu(\sigma) - \mu(\sigma_{\text{ref}})$ represents the difference in chemical potential between particles of size $\sigma$ and particles of an arbitrary reference size $\sigma_{\text{ref}}$. Since the combination of $T$, $P$ and $\delta\mu(\sigma)$ uniquely define a crystal state point, there is only one trial crystal state that can coexist with a given trial fluid. The final requirement for an equilibrium coexistence is then the constraint that the absolute chemical potentials of the two phases are also equal. This can be checked in a direct-coexistence simulation by confirming that the trial crystal can coexist with the initial trial fluid without deforming. For this, we use the trick of Ref. 29 of checking for crystal deformation by monitoring the pressure component $P_{zz}$ along the direction perpendicular to the interface between the two phases.

Finding the equilibrium freezing point for a polydisperse system with a given polydispersity $p$ and temperature $T$ therefore comes down to the following steps:

* Choose a series of “trial” fluid densities $\rho_f^{\text{trial}}$, and measure the corresponding pressures $P_f^{\text{trial}}$ and relative chemical potentials $\delta\mu_f^{\text{trial}}(\sigma)$.

* For each $\rho_f^{\text{trial}}$, determine the density (and corresponding lattice spacing) of the corresponding trial crystal phase that could coexist with this fluid. In other words, find the equilibrium bulk crystal that is stable at $P_f^{\text{trial}}$ and $\delta\mu_f^{\text{trial}}(\sigma)$.

* For each $\rho_f^{\text{trial}}$, perform a direct-coexistence simulation in the semi-grand canonical ensemble between the trial fluid and trial crystal phase. If the coexistence remains stable, check for deformation of the crystal phase by measuring $P_{zz}$.

After these steps, the identification of the equilibrium freezing line is done by determining the value of $\rho_f^{\text{trial}}$ where $P_{zz}$ in the direct-coexistence simulation is equal to $P_f^{\text{trial}}$. Note that the above steps ensure that the fluid phase in the direct coexistence simulation has the same pressure and relative chemical potential as the trial fluid, which fixes its size distribution to be equal to that of the trial fluid as well.

In the following subsections, we outline each of these steps in more detail.

### A. Chemical potential of the fluid phase

We start by measuring the bulk properties of our trial fluids. To this end, we perform event-driven molecular dynamics (EDMD) simulations (based on Ref. 30) of polydisperse hard-sphere systems in the canonical ensemble. We use system sizes of both 16000 and 32000 particles, and simulate for a total simulation time of $1.5 \times 10^5 \tau$ (1/3 of it being dedicated to equilibration). Here, our time unit is $\tau = \sqrt{\beta m \sigma^2}$, where $\beta = 1/k_B T$ and $m$ is the particle mass, which is taken to be the same for all particles. Hence, $\tau$ corresponds to the time it takes a typical particle in free flight to move over a distance roughly equivalent to its own diameter. For the sake of reproducibility and to avoid significant deviations from the expected size distribution, we chose to generate particle sizes deterministically using the inverse of the cumulative probability distribution function. Further details on this can be found in the supplementary material (SM).

From each of these simulations, our aim is to determine both the pressure $P$ and the relative chemical potential $\delta\mu(\sigma)$. We obtain the pressure via the standard virial expression for hard-sphere systems³¹:

$$ \beta P_{kl}/\rho = 1 + \frac{\sum_{ij} m v_{ij}^k r_{ij}^l}{N(t_b - t_a)}, \quad (2) $$

where $P_{kl}$ indicates the kl-component of the pressure tensor (with $k, l \in \{x, y, z\}$), $\rho = N/V$ is the number density, and the sum is taken over all collisions in a time interval
$[t_a; t_b]$. Additionally, $r_{ij}^k$ and $v_{ij}^k$ are the k-components of the relative distance and velocity vectors between the colliding particles i and j, respectively.

For polydisperse mixtures, the total chemical potential can be expressed as the sum of two contributions, namely the mixing part $\mu_{\text{mix}}(\sigma)$ and the configurational part $\mu_{\text{conf}}(\sigma)$, which are both functions of the particle size. The former is an intrinsic property of the particle size distribution, given by $\beta\mu_{\text{mix}}(\sigma) = -\log(P(\sigma)\lambda)$, with $\lambda$ an arbitrary length scale analogous to the thermal wavelength. The derivative of the latter, $\mu'_{\text{conf}}(\sigma)$, can be linked to the internal energy of the system following (see SM for details)

$$ \mu'_{\text{conf}}(\sigma_i) = \left\langle \frac{\partial U}{\partial \sigma_i} \right\rangle, \qquad (3) $$

where $\langle \dots \rangle$ denotes an ensemble average, $U$ is the potential energy of the system, and $\sigma_i$ refers to the size of particle i.

As we show in the SM, Eq. 3 can then be interpreted in terms of an average force experienced by the particle radius $R_i$:

$$ \mu'_{\text{conf}}(\sigma_i) = -\langle f_{R_i} \rangle / 2, \qquad (4) $$

where a positive $f_{R_i}$ corresponds to an outward force on the particle surface. In practice, for hard spheres, the forces $f_i$ on any particle's surface are due to collisions, and always directed radially inwards, such that $f_{R_i}$ is always negative (and hence $\mu'_{\text{conf}} > 0$). We can write the ensemble average:

$$ \langle f_{R_i} \rangle = \frac{-1}{t_b - t_a} \sum_k m |(\delta\mathbf{v}_i)_k|, \qquad (5) $$

where $(\delta\mathbf{v}_i)_k$ is the change in velocity of particle i during collision k, and the sum is taken over all collisions in the time interval $[t_a, t_b]$.

Combining $\mu'_{\text{conf}}(\sigma)$ and $\mu'_{\text{mix}}(\sigma)$, we can calculate the derivative of the total chemical potential, $\mu'(\sigma)$, which can be fitted and integrated to obtain $\delta\mu(\sigma)$ – the function needed to perform semi-grand canonical simulations.

The typical behavior of the surface forces is shown in Fig. 1a for systems with polydispersity $p = 0.06$ (see particle size distribution in Fig. 1b). We observe that the forces become monotonically stronger as the density increases – which makes sense when considering that higher densities will lead to a higher collision rate. Similarly, the forces are stronger for larger particle sizes, which can be understood from the fact that larger particles will experience more collisions on their larger surface. For all of our data, a quadratic fit in $\sigma$ is in good agreement with the measured forces. To obtain a consistent fit across all investigated densities, we fit the forces with a polynomial function of the form

$$ \beta f(\sigma, \rho\bar{\sigma}^3)\bar{\sigma} = \sum_{i=0}^{2} \sum_{j=0}^{2} c_{i,j} (\rho\bar{\sigma}^3)^i \left(\frac{\sigma}{\bar{\sigma}}\right)^j. \qquad (6) $$

The outcome of this fit $f(\sigma, \rho\bar{\sigma}^3)$ is shown in dashed lines in Fig. 1a, and is in good agreement with the computed forces. This procedure allows us to explore systems with any density in this range without having to recalculate the forces. We have confirmed that surface forces measured in larger systems (up to 32000 particles) are consistent with our predictions in order to rule out significant finite-size effects.

We show the different contributions to the chemical potential difference $\delta\mu(\sigma, \rho\bar{\sigma}^3)$ in Fig. 1c for a system with polydispersity $p = 0.06$ and density $\rho\bar{\sigma}^3 = 0.96$. Fig. 1d shows the total chemical potential difference for a range of densities for this same system. Note that in Fig. 1d we plot $-\delta\mu(\sigma)$ as this corresponds to the external potential field $V$ that we impose in the semi-grand-canonical simulations, described in the next section.

## B. Mechanical equilibrium with the crystal phase

Now that we have access to the behavior of the chemical potential difference as a function of system size, $\delta\mu(\sigma)$, we turn to simulations of hard spheres in the semi-grand canonical ensemble. In a similar approach to the continuous time swap approach developed in Ref. 32, we consider a semi-grand ensemble simulation in which particles have three positional degrees of freedom supplemented by a fourth one: their radius $R_i = \sigma_i/2$. The Hamiltonian for this system can then be written as:

$$ H = \sum_{i} \left[ \frac{\mathbf{p}_{i}^{2}}{2m} + \frac{\varphi_{i}^{2}}{2M} + V(R_{i}) \right] + \sum_{i<j} U_{ij}(\mathbf{r}_{i}, \mathbf{r}_{j}, R_{i}, R_{j}), \qquad (7) $$

where $\mathbf{p}_i$ is the translational momentum of particle i, $m$ is the (translational) particle mass, $\varphi_i$ is the momentum associated with its radius $R_i$, $M$ is the corresponding mass, $V(R_i) = -\delta\mu(\sigma_i = 2R_i)$ is the external field controlling the particle size, and $U_{ij}$ represents the pair interaction (which is either 0 or $+\infty$ for hard spheres).

The Hamiltonian in Eq. 7 naturally gives rise to equations of motion for the positions and radii of the particles. Note that since our interaction potential is discontinuous, forces between the particles only occur at contact, and the collision rules can be derived by regarding the hard-sphere interaction as the limiting case of a sharply repulsive continuous interaction (see SM). These collision rules can be straightforwardly implemented in an EDMD simulation. However, the motion of the radial degree of freedom is additionally affected by the continuous field $V(R_i)$, which is not directly compatible with the event-driven nature of EDMD simulations. Fortunately, this field can be efficiently handled with the event-driven Monte Carlo (EDMC) approach introduced by Peters and De With$^{33,34}$. In this event-driven approach, continuous interactions (or in this case, fields) are handled stochastically. Specifically, whenever we predict future collisions for a particle, we also consider collisions with the field $V(R_i)$, which are determined based on a maxi-
<figure><img src="image_2.png" /><figcaption>FIG. 1. a) Surface forces $f_i = -f_{\sigma_i}$ measured for a 16000 particle system with Gaussian particle size distribution and polydispersity $p = 0.06$, for a range of system densities $\rho\bar{\sigma}^3$ (top to bottom: 0.94, 0.95, 0.96, 0.97, and 0.98). Shading corresponds to one standard error and dashed black lines represent the fits $f(\sigma, \rho\bar{\sigma}^3)$ we employed in subsequent simulations. b) Particle size distribution for a system with polydispersity $p = 0.06$. c) Configurational contribution (dashed line), mixing contribution (dotted line), and total chemical potential difference for the above system at density 0.96. d) Total chemical potential difference for the systems shown in the left panel.</figcaption></figure>

momentum increase $\Delta V$ drawn from a Boltzmann distribution (see SM for details). Collisions with the field $V(R_i)$ simply reverse the associated momentum $\wp_i$ for the particle under consideration.

We now use the resulting semi-grand canonical simulations to identify the trial crystal phase that could coexist with a given trial fluid. Note that for a given $\rho_f^{\text{trial}}$, we already know the pressure $P^{\text{trial}}$ and the associated chemical potential function $\delta\mu(\sigma)$. The coexisting crystal phase must match the fluid in both of these quantities. To find this coexisting crystal, we perform simulations of FCC crystals at different densities $\rho_\chi$ but at the same relative chemical potential $\delta\mu^{\text{trial}}(\sigma)$, and measure the pressure in each simulation. We then locate the density $\rho_\chi^{\text{trial}}$ such that $P_\chi(\rho_\chi^{\text{trial}}) = P^{\text{trial}}$. Figure 2 shows this procedure for a system with polydispersity $p = 0.06$ and trial fluid density $\rho_f^{\text{trial}}\bar{\sigma}^3 = 0.9618$. We used crystal bulk systems with 16384 particles set up on a FCC lattice, which were equilibrated for $25 \times 10^3\tau$ and sampled over an additional $25 \times 10^3\tau$. We find that the pressure is a monotonically increasing function of $\rho_\chi$ (as expected), and simply locate the intersection of $P_\chi(\rho_\chi^{\text{trial}})$ and $P^{\text{trial}}$ by fitting the former with a polynomial and numerically locating the crossing point based on the fit.

### C. Semi-grand direct-coexistence simulations

At this point, for each given trial fluid, we have identified a crystal phase that could coexist with it: the two phases have the same temperature, pressure, and relative chemical potential. However, we have not yet confirmed that this coexistence is indeed a stable equilibrium, as there is no guarantee that the trial fluid and the associated crystal have the same absolute chemical potential. To address this, we now perform semi-grand direct-coexistence simulations following the methodology developed in Ref. 29. We initialize a crystal slab with density $\rho_\chi^{\text{trial}}$ in an elongated box in the $z$ direction with aspect ratio 1 : 3. The global density inside the box is set to $(\rho_f^{\text{trial}} + \rho_\chi^{\text{trial}})/2$, and approximately half of it is initially filled with crystal. The rest of the box is filled with smaller fluid particles to match the desired global density. These particles are then rapidly grown to their target size distribution and the system is equilibrated for $10^5\tau$. A typical snapshot of a coexisting system is shown in Fig. 3a.

We note that by using this choice for the global density, we place our coexisting system essentially at the midpoint of the tie-line that connects the cloud and shadow point corresponding to the two coexisting phases. This is similar in spirit to the approach used in Ref. 36 for determin-
<figure><img src="image_3.png" /><figcaption>FIG. 2. Crystal pressure as a function of crystal density obtained from semi-grand simulations of a bulk crystal of 16384 particles. The imposed chemical potential difference corresponds to that of a trial fluid system with polydispersity $p = 0.06$ at density $\rho_{f}^{trial}\bar{\sigma}^{3} = 0.9618$. Pressure of the trial fluid is reported on the horizontal dotted line. Standard errors are shown for all points but are typically found to be smaller than the point size. The dashed line corresponds to the quadratic fit employed to determine the crystal density at mechanical equilibrium.</figcaption></figure>

ing the cloud point in a polydisperse lattice gas. Rather than direct coexistence simulations, Ref. 36 used grand-canonical ensemble simulations and iterative histogram reweighting approaches to locate coexistence where one of the two fluid phases has the desired size distribution. While this approach would be impractical for our system (e.g. because particle insertions are not effective in dense or crystalline phases), this has the similar effect of setting up a coexistence situation where each phase occurs with equal weight in order to locate the cloud point.

To determine whether the obtained coexistence is thermodynamically stable, we check whether the fluid and crystal phase remain at their initial state point by monitoring the pressure tensor. If the chemical potential is not equal in the two phases, their densities will adjust to reach thermodynamic equilibrium. This results in a deformation of the crystal along the long axis of the box (as the lattice spacings in the other directions are fixed by the box) and a change in overall density in the fluid. This deformation will inevitably change the pressure component $P_{zz}$ (along the long axis of the box), which is homogeneous throughout the system²⁹. In contrast, at equilibrium coexistence, both phases keep their initial density, and $P_{zz}$ will be equal to $P^{trial}$, resulting in an unstrained crystal coexisting with the trial fluid, at equal temperature, pressure, and chemical potential. Moreover, since both the pressure and relative chemical potential of the fluid phase in this coexistence simulation match those of our bulk trial fluid, its size distribution necessarily also matches that of the trial fluid. In contrast, the crystal phase will generally have a significantly different size distribution. Hence, the phase coexistence we found corresponds exactly to the fluid cloud point and

the associated FCC crystal shadow phase.

Fig. 3b shows the result of this procedure for the $p = 0.06$ system, repeated for a range of trial fluids. Specifically, we plot the measured value of $P_{zz}$ (green points) as a function of the trial pressure $P^{trial}$, and numerically locate the crossing point with the line $P_{zz} = P^{trial}$. This yields the freezing point for this polydisperse mixture. Note that while significant noise is visible in Fig. 3b, the effect of this noise on the resulting coexistence pressure is weak: we estimate the coexistence pressure to be $\beta P^{coex}\bar{\sigma}^{3} = 13.37(1)$.

Additionally, we plot the components of the pressure tangential to the interfaces ($P_{xx}$ and $P_{yy}$), which are related to the stress on the simulation box caused by the interfaces²⁹. As expected from the choice of crystal lattice orientation inside the box (square plane facing the fluid), we recover the equality between the $x$ and $y$ components of the pressure tensor.

### III. RESULTS AND DISCUSSION

To study the full extent of the fluid-FCC freezing lines in hard-sphere mixtures with Gaussian polydispersity, we repeat the procedure outlined in the previous sections on systems with different polydispersities ranging from 2% to 14%. The main results for the phase behavior are shown in Fig. 4. Specifically, in Fig. 4a we report the fluid-FCC coexistence pressure as a function of polydispersity. As expected, increasing polydispersity destabilizes the crystal, resulting in higher coexistence pressures. The same trend is seen in the packing fractions $\eta$ of the two coexisting phases (Fig. 4b), which both monotonically increase with increasing $p$. Intriguingly, plotting the number density $\rho = N/V$ of the two coexisting phases (Fig. 4c) shows a different perspective. While the number density of the coexisting fluid increases with polydispersity, the corresponding crystal density decreases. This is indicative of the fact that with increasing polydispersity, the mean size of the particles also increases, as shown in Fig. 5a. This results in a larger lattice spacing and hence a lower number density. At a polydispersity close to $p \simeq 9\%$, the number densities of the fluid and crystal phase cross. This renders the precise measurement of the coexistence points more difficult. Specifically, when the density difference between the two phases is small, fluctuations in the amount of crystal phase in the direct-coexistence simulations grow larger, making it more likely that the crystal fully melts or starts percolating the box. Hence, in this regime we cannot keep the coexistence stable over long enough simulation times to accurately measure the pressure tensor. Due to this, we were unable to accurately determine the coexistence conditions at $p = 10\%$ and omit this point from Fig. 4. At even higher polydispersities ($p \ge 12\%$) we can once again obtain a stable coexistence. In this regime, the number density in the FCC phase is lower than in the fluid, consistent with earlier predictions by
<figure><img src="image_4.png" /><figcaption>FIG. 3. a) Snapshot of a fluid with polydispersity $p = 0.06$ in coexistence with a crystal. For visualization purposes, particles found to be in a local crystalline environment using ten Wolde's local bond order parameter$^{35}$ $q_6$ are displayed in red, while particles found to be fluid-like are displayed in blue. b) Pressure tensor components $P_{xx}$, $P_{yy}$, and $P_{zz}$ as a function of trial pressure $P^{trial}$ obtained from semi-grand direct-coexistence simulations for a 16384 particle system in which the fluid has polydispersity $p = 0.06$. Error bars show one standard error. These systems were constructed such that the fluid-crystal interface corresponds to the square-lattice plane of an FCC crystal.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIG. 4. a) Coexistence pressure as a function of the polydispersity of the fluid phase. b) Freezing line (blue) and corresponding shadow curve (red) in the polydispersity-packing fraction plane. c) Freezing line and shadow curve in the polydispersity-density plane. Data presented here corresponds to a fluid in coexistence with a FCC crystal, the fluid-crystal interface corresponds to the square lattice plane of the FCC crystal. Note that error bars are shown on all points but are smaller than the point size (typically on the order of 0.1%). All lines correspond to the functional forms reported in the main text. Data for the monodisperse case was obtained from Ref. 29.</figcaption></figure>

Wilding and Sollich$^{13}$.

well-fitted by the following functions:

$$ \beta P^{\text{coex}} \bar{\sigma}^3 = 11.5646 + 433.63950p^2 + 2308.37288p^3 - 19711.1098p^4 + 19077.0257p^5, \quad (8) $$

$$ \eta_f = 0.4918 + 4.90278p^2 + 23.48901p^3 - 375.76734p^4 + 1021.79611p^5, \quad (9) $$

$$ \eta_{\chi} = 0.5432 + 0.31861p^2 + 60.26721p^3 - 361.66200p^4 + 545.04336p^5, \quad (10) $$

$$ \rho_f \bar{\sigma}^3 = 0.93918 + 5.63121p^2 + 93.78531p^3 - 1645.24628p^4 + 8878.27295p^5 - 18195.4842p^6, \quad (11) $$

$$ \rho_{\chi} \bar{\sigma}^3 = 1.0375 - 14.54067p^2 + 184.78167p^3 - 1074.12989p^4 + 1577.06258p^5. \quad (12) $$

Empirically, we find that the coexistence conditions are
<figure><img src="image_6.png" /><figcaption>FIG. 5. Mean particle size (a) and polydispersity (b) of the crystal phase as a function of the fluid polydispersity $p$. The dashed lines indicate the corresponding mean size and polydispersity for the fluid phase for comparison.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Fluid cloud curve and the associated crystal shadow curve in the packing fraction - polydispersity plane. Results from this work are shown as blue points and lines. We compare these to data points obtained from Refs. 9, 10, and 13, which correspond respectively to a symmetric triangular, Schultz, and top-hat particle size distribution of the parent phase.</figcaption></figure>

Here, $\eta_f$ and $\eta_\chi$ are the packing fractions of the fluid and FCC phase at coexistence, and $\rho_f$ and $\rho_\chi$ are the corresponding number densities.

Taking a closer look at the behavior of the composition of the crystal line along the shadow curve, we observe that while the mean particle size rises significantly with increasing fluid polydispersity, the polydispersity of the crystal saturates, approaching values around 5.6% (Fig. 5b), consistent with past predictions using various parent size distributions$^{9–12}$.

To compare to past predictions in more detail, we plot in Fig. 6 the fluid cloud curve and the corresponding FCC shadow curve in the packing fraction - polydispersity plane. We compare these curves to the results from Refs. 9, 10, and 28, which correspond respectively to hard-sphere fluids with a symmetric triangular, Schultz, and top-hat particle size distribution. Note that in this representation, we use the polydispersity of the crystal phase for the crystal branch. Despite the differences in size distribution, we observe good agreement between all sets of data. Nonetheless, our new results deviate more from the data of Refs. 9, 10, and 13 than the latter three sets differ between each other. In fact, the deviations already start in the limit of zero polydispersity, where the distribution shape should not matter. We attribute this deviation to inaccuracies in the equations of state used in Refs. 9, 10, and 13. Specifically, for the fluid phase, these works use the Boublik, Mansoori, Carnahan, Starling and Leland (BMCSL) expression$^{37,38}$, which reduces to the standard Carnahan-Starling expression$^{39}$ in the monodisperse case. For monodisperse hard spheres, Carnahan-Starling is known to underestimate the pressure at high densities$^{40}$. From our own pressure measurements (not shown), the BMCSL similarly underestimates the pressure of polydisperse mixtures. Such deviations would cause a shift in the freezing transition to higher packing fractions. This is consistent with the data of Refs. 9, 10, and 13, who report the monodisperse freezing and melting packing fractions at approximately 0.494 and 0.545, respectively, which are high compared to modern estimates of approximately 0.492 and 0.543$^{5,29}$.

We expect the phase behavior reported in Fig. 4 to be sensitive to finite-size effects which will typically serve to help stabilize the crystal phase. The data we report in Fig. 4 is obtained for systems with $N = 16000$ particles and the FCC lattice orientation in the box was chosen so that the square crystal plane is facing the fluid. To examine the strength of the finite-size effects, we repeat our calculations for different system sizes (up to $N = 64000$) for both the same orientation of the crystal phase (denoted FCC$_a$) and one where the hexagonal planes are aligned along one of the long walls of the box (FCC$_b$, see Fig. 7a). For small system sizes, we find variations in our predictions of the coexistence pressure for $p = 0.06$ as we report in Fig. 7b. Additional data on these different system sizes is provided in the SM. Nevertheless, for the largest systems considered, coexistence pressure predictions are indistinguishable for the two crystal orientations we consider. This is consistent with the fact that the crystal orientation should not affect coexistence conditions in the thermodynamic limit.

We now turn to the characterization of the fluid-crystal interface at the freezing points determined above.

First, we measure the surface stresses at the fluid-crystal interface. We can write the xx-component of the surface stress as$^{41}$

$$f_{xx} = \frac{1}{2} L_z (P_{zz} - P_{xx}), \qquad (13)$$
<figure><img src="image_8.png" /><figcaption>FIG. 7. a) Schematic representation of the choice of crystal orientation in the simulation box. Note that the FCC_b orientation used in this work differs from the one used in Ref. 29 by 90 degree rotation in the xz plane. b) Coexistence pressure estimates as a function of system size for a system with polydispersity $p = 0.06$, and for the two different FCC crystal lattice orientations shown in the top panel. Lines are guides to the eye. Error estimates for $P^{coex}$ are determined from the noise in the individual pressure measurements (e.g. those shown in Fig. 3).</figcaption></figure>

where $L_z$ denotes the length of the simulation box in the elongated $z$ direction. The expression for the $yy$ component is analogous. The stress tensor is measured directly from semi-grand simulations at the coexistence points determined above, for systems of 16000 particles equilibrated for $10^5 \tau$. These quantities are reported in Fig. 8, along with $\frac{1}{2}(f_{xx} + f_{yy})$, which for the specific case of the square symmetry of the FCC_a interface corresponds to the scalar value of the surface stress.

For the FCC_a interface, increasing polydispersity rapidly causes the surface stress to depart from the known value for the monodisperse case$^{41,42}$. This is not the case for the FCC_b interface, for which the surface stresses remain relatively constant in the low polydispersity regime, and close to the monodisperse values we evaluated from additional direct-coexistence simulations. Upon increasing polydispersity, we notice an upward trend for both the FCC_a and FCC_b interfaces. For the two cases, the surface stress becomes positive at polydispersities around 0.06 and 0.1, respectively. These measurements are particularly sensitive to the precise determination of the coexistence point, which translates into the error bars reported in Fig. 8. While the sign change

<figure><img src="image_9.png" /><figcaption>FIG. 8. Surface stresses measured at coexistence for polydisperse systems in the case of a FCC_a (a) and FCC_b (b) fluid-crystal interface plane. Error bars correspond to twice the standard error obtained from 10 to 20 independent runs at the coexistence we determine for 16384 particles systems. Note that the actual error is likely to be larger, due to systematic errors introduced by inaccuracies in our determination of the coexistence point. Points for the monodisperse case were obtained from Ref. 41 for FCC_a and from additional monodisperse direct-coexistence simulations for FCC_b. All lines are guides to the eye.</figcaption></figure>

of the surface stress $f$ may at first glance appear surprising, one can recall the Shuttleworth equation$^{43,44}$ that links it to the surface free energy $\gamma$:

$$f_{ij} = \delta_{ij}\gamma + \frac{\partial\gamma}{\partial u_{ij}}, \qquad (14)$$

with $u_{ij}$ the surface strain tensor and $\delta_{ij}$ the Kronecker delta. For monodisperse hard spheres, the surface stress of the investigated interfaces is negative. As $\gamma$ must be positive, this negative sign indicates that the second term in Eq. 14 is negative: stretching the interface results in a lower $\gamma$. The increase in $f$ with increasing polydispersity can then be attributed to two possible effects. First, the rise in coexistence pressure with increasing polydispersity likely gives rise to higher surface free energies $\gamma$. Second, the ability of particle sizes to adapt to the structure of the interface may make $\gamma$ less dependent on the strain applied to the interface, which would lower the magnitude of the (negative) second term. A more careful analysis of the surface free energy would be required to disentangle these two effects.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Scan of particle size distributions across the FCCa fluid-crystal interface for a system of 16000 particles where the fluid has polydispersity $p = 0.06$ (a) and $p = 0.12$ (b). The dashed black lines correspond to the respective size distributions of the reference fluid and crystal phases. Note that the range of sizes particles are allowed to explore in simulation far exceeds the range shown in this figure.</figcaption></figure>

Next, we examine how the distribution of particles varies in space as we move from one phase to the other through the interface. From the semi-grand simulations we performed at the coexistence point, we find the interface to be highly diffusive and very mobile. Hence, to be able to characterize the interfacial structure, we shift the position of all particles in each snapshot to correct for the motion of the interfaces. To this end, we monitor the mean displacement of the crystal in the elongated $z$-direction and shift the particles accordingly. However, this trick does not prevent the crystal slab from growing or shrinking. As a second measure to correct for shifts in the positions of the two individual interfaces, we measure the ten Wolde $q_6$ bond orientational order parameters³⁵ as a function of $z$ and locate one interface at half the height of a fitted sigmoid function.

We measure the size distribution as a function of the distance to the interface along the long axis of the simulation box. Fits to these size distributions are shown in Fig. 9 for systems with polydispersity $p = 0.06$ and $p = 0.12$ in the FCCa orientation. These measurements corroborate what we previously reported in Fig. 4: as we progress from the fluid phase to the FCC crystal phase, the particle size distributions become increasingly more peaked and their means shift to larger sizes. In other words, we see a crystal structure preferentially made of larger particles, and with less size diversity compared to what can be found in the fluid, consistent with what was

predicted in Refs.³,⁹⁻¹⁵. As we show the same range for both scans in $z$ in Fig. 9, we also note that the composition of the system can undergo drastic changes over a small distance. Closer inspection of the size distributions in the crystal bulk and at the interface reveal that they deviate from the expected Gaussian behavior by being skewed towards the smaller sizes.

Finally, our direct-coexistence simulations give us access to the interfacial excess absorption of particles of a given size $\sigma_i$. This quantity reflects the excess amount (which can be either positive or negative) of a component – here a particle of size $\sigma_i$ – present in the system with respect to a reference system. Following Gibbs formalism, we interpret our system as consisting of two bulk phases separated by a dividing surface which has no volume but can have particles absorbed to it. Hence, the total volume of the system is given by

$$V = V_f + V_{\chi}, \qquad (15)$$

were the subscripts $f$ and $\chi$ refer respectively to the fluid and crystal phase. The total particle size distribution can be written as

$$N(\sigma) = N_f(\sigma) + N_{\chi}(\sigma) + N_I(\sigma), \qquad (16)$$

where $N(\sigma)$ denotes the total number of particles of a given size $\sigma$, and the subscript $I$ refers to the interface. As our dividing surface, we choose the equimolar surface,
<figure><img src="image_11.png" /><figcaption>FIG. 10. Surface excess amount of particles with size $\sigma_i$ for an equimolar Gibbs dividing surface in a system with 16384 particles at coexistence. Here the crystal lattice orientation is such that the square plane of the FCC crystal faces the fluid (FCC$_a$). Polydispersities are reported on each line, all lines are vertically shifted for readability.</figcaption></figure>

which satisfies

$$ \int d\sigma N_I(\sigma) = 0. \qquad (17) $$

We set up our coexistence simulations so that the global density of the box is precisely halfway in the coexistence region. For the equimolar surface, this implies that $V_F = V_X = V/2$. The excess absorption at the interface $N_I(\sigma)$ is then given by

$$ N_I(\sigma) = N(\sigma) - \frac{V}{2}(\rho_f P_f(\sigma) - \rho_\chi P_\chi(\sigma)), \quad (18) $$

where $\rho_f$ and $\rho_\chi$ are the coexistence densities of the fluid and crystal phase, respectively, and $P_f$ and $P_\chi$ are the corresponding particle size distributions. This quantity, normalized by the total interface area A, is reported in Fig. 10 for the range of polydispersities we investigated and for the FCC$_a$ interface. We observe a significant change in the qualitative behavior of the surface excess which is strongly correlated with the crossing of the coexistence density lines we show in Fig. 4. For the lower polydispersity regime ($p \le 0.08$, leaving out $p = 0.02$), we observe a deficit of small particles and an excess of

large particles at the interface, both of which are growing in amplitude with $p$. For the higher polydispersity regime ($p \ge 0.12$), i.e. above the crossing of the density lines, this trend shifts entirely. We now see a net excess of both very small and very large particles, accompanied by a large deficit of particles with intermediate-to-large sizes. Interestingly, at low polydispersities the picture is different for the FCC$_b$ interface (see SM), where we observe a deficit of the smaller particles and an excess of the larger ones from $p = 0.02$ and up. We note that these measurements are likely to be sensitive to small inaccuracies in our determination of the coexistence conditions, as well as to long-time fluctuations in the amount of crystal present in the simulation box. For reference, we show in the SM the composition of all investigated systems at coexistence as well as the composition of the reference fluid and crystal phases.

## IV. CONCLUSION

In conclusion, we have introduced a novel method for locating phase boundaries in polydisperse systems, based on direct-coexistence simulations. Using this approach, we have mapped out the freezing line of polydisperse hard spheres with a Gaussian size distribution, finding good agreement with past predictions. Additionally, we have used the direct-coexistence simulations to probe the properties of the fluid-crystal interface. We found that the surface stress increases as a function of polydispersity, so much as to become positive for $p \ge 0.06$, suggesting that particles at the interface adapt their sizes to reduce interfacial stress. This idea is corroborated by the observation of size-selectivity of excess particles associated with the interface.

While the discussion presented here is limited to coexistence between a Gaussian polydisperse fluid and a FCC crystal, a natural next step would be to use the same method to examine the stability of the more exotic phases that were seen to form spontaneously in highly polydisperse mixtures$^{2–4}$. Additionally, it would be interesting to examine how the choice of particle size distribution in the fluid phase impacts the emerging crystal phase.

We emphasize that the methodology presented here is not limited to hard spheres. Simulations in the canonical and semi-grand ensemble can be performed on systems interacting via any potential, and the only quantities we measure are particle forces (in the canonical ensemble) and pressure tensors (in the semi-grand ensemble), which are standard quantities of interest. The main caveat is that – as seen in our case for polydispersity $p = 0.10$, the direct-coexistence method struggles in systems where a stable interface is hard to form, which occurs when the density gap between the phases is very small. Additionally, in systems where dynamics are slow, equilibration of the coexisting system can take a long time, and hence be computationally expensive. Finally, we emphasize
size that the method proposed here is aimed at locating coexistences where the size distribution in one of the two phases (in this case the fluid) is known. This makes the method ideal for locating cloud and shadow curves, but unsuitable for finding the coexistence that would result from phase separation in the canonical ensemble, where the total size distribution (of both phases and the interface combined) would be fixed. For addressing the latter question, more complex simulation approaches or theoretical approximations are likely required, such as those developed in Refs. 15, 28, 45, and 46. Despite these limitations, we believe that our approach provides an elegant and straightforward method for finding cloud and shadow curves in polydisperse systems.

## V. SUPPLEMENTARY MATERIAL

Additional details for the generation of the deterministic particle size distribution as well as derivations for the effect of particle size on the free energy can be found in the supplementary material available online. We also provide the derivation for the equations of motions in the semi-grand canonical ensemble, as well as implementation details for the EDMC approach. An overview of the effects of system size on the determination of coexistence conditions is also provided, along with additional details on the phase composition of the system at coexistence, and the excess surface absorption for the $FCC_b$ interface.

## VI. ACKNOWLEDGMENTS

We thank Marjolein de Jager for helpful discussions. AC and FS acknowledge the use of the CERES high performance computer cluster at the Laboratoire de Physique des Solides of Université Paris-Saclay (Orsay, France). AC and FS acknowledge funding from the Agence Nationale de la Recherche (ANR), grant ANR-21-CE30-0051. LF acknowledges funding from the Vidi research program with project number VI.VIDI.192.102 which is financed by the Dutch Research Council (NWO).

## VII. AUTHOR DECLARATIONS

### A. Conflict of Interest

The authors have no conflicts to expose.

## VIII. DATA AVAILABILITY

A version of the semi-grand canonical EDMC simulation code for polydisperse hard spheres is available online at https://github.com/FSmallenburg/. The data that support the findings presented in this study are available

in a data package on Zenodo⁴⁷. The latter also includes additional surface forces data measured for a wide range of systems.

## IX. REFERENCES

¹B. Cabane, J. Li, F. Artzner, R. Botet, C. Labbez, G. Bareigts, M. Sztucki, and L. Goehring, “Hiding in plain view: colloidal self-assembly from polydisperse populations,” Phys. Rev. Lett. 116, 208001 (2016).

²P. K. Bommineni, N. R. Varela-Rosales, M. Klement, and M. Engel, “Complex crystals from size-disperse spheres,” Phys. Rev. Lett. 122, 128005 (2019).

³D. Coslovich, M. Ozawa, and L. Berthier, “Local order and crystallization of dense polydisperse hard spheres,” J. Phys. Condens. Matter 30, 144004 (2018).

⁴B. A. Lindquist, R. B. Jadrich, and T. M. Truskett, “Communication: From close-packed to topologically close-packed: Formation of laves phases in moderately polydisperse hard-sphere mixtures,” J. Chem. Phys. 148 (2018).

⁵C. P. Royall, P. Charbonneau, M. Dijkstra, J. Russo, F. Smallenburg, T. Speck, and C. Valeriani, “Colloidal hard spheres: Triumphs, challenges, and mysteries,” Rev. Mod. Phys. 96, 045003 (2024).

⁶R. P. Sear, “Molecular dynamics of a dense fluid of polydisperse hard spheres,” J. Chem. Phys. 113, 4732–4739 (2000).

⁷D. Orsi, A. Fluerasu, A. Moussaïd, F. Zontone, L. Cristofolini, and A. Madsen, “Dynamics in dense hard-sphere colloidal suspensions,” Phys. Rev. E 85, 011402 (2012).

⁸R. C. Roberts, N. Marioni, J. C. Palmer, and J. C. Conrad, “Dynamics of polydisperse hard-spheres under strong confinement,” Mol. Phys. 118, e1728407 (2020).

⁹M. Fasolo and P. Sollich, “Equilibrium phase behavior of polydisperse hard spheres,” Phys. Rev. Lett. 91, 068301 (2003).

¹⁰M. Fasolo and P. Sollich, “Fractionation effects in phase equilibria of polydisperse hard-sphere colloids,” Phys. Rev. E 70, 041410 (2004).

¹¹P. G. Bolhuis and D. A. Kofke, “Monte Carlo study of freezing of polydisperse hard spheres,” Phys. Rev. E 54, 634 (1996).

¹²D. A. Kofke and P. G. Bolhuis, “Freezing of polydisperse hard spheres,” Phys. Rev. E 59, 618 (1999).

¹³N. B. Wilding and P. Sollich, “Phase behavior of polydisperse spheres: Simulation strategies and an application to the freezing transition,” J. Chem. Phys. 133 (2010).

¹⁴P. G. Bolhuis and D. A. Kofke, “Numerical study of freezing in polydisperse colloidal suspensions,” J. Phys. Condens. Matter 8, 9627 (1996).

¹⁵P. Sollich and N. B. Wilding, “Crystalline phases of polydisperse spheres,” Phys. Rev. Lett. 104, 118302 (2010).

¹⁶S. Auer and D. Frenkel, “Quantitative prediction of crystal-nucleation rates for spherical colloids: A computational approach,” Annu. Rev. Phys. Chem. 55, 333–361 (2004).

¹⁷H. J. Schöpe, G. Bryant, and W. Van Megen, “Small changes in particle-size distribution dramatically delay and enhance nucleation in hard sphere colloidal suspensions,” Phys. Rev. E 74, 060401 (2006).

¹⁸E. Zaccarelli, C. Valeriani, E. Sanz, W. C. K. Poon, M. E. Cates, and P. N. Pusey, “Crystallization of hard-sphere glasses,” Phys. Rev. Lett. 103, 135704 (2009).

¹⁹P. K. Bommineni, M. Klement, and M. Engel, “Spontaneous crystallization in systems of binary hard sphere colloids,” Phys. Rev. Lett. 124, 218003 (2020).

²⁰S. Martin, G. Bryant, and W. Van Megen, “Crystallization kinetics of polydisperse colloidal hard spheres: Experimental evidence for local fractionation,” Phys. Rev. E 67, 061405 (2003).

²¹R. A. LaCour, T. C. Moore, and S. C. Glotzer, “Tuning stoichiometry to promote formation of binary colloidal superlattices,” Phys. Rev. Lett. 128, 188001 (2022).
²²R. Evans and C. Holmes, “Diffusive growth of polydisperse hard-sphere crystals,” Phys. Rev. E **64**, 011404 (2001).

²³R. Kurita, D. B. Ruffner, and E. R. Weeks, “Measuring the size of individual particles from three-dimensional imaging experiments,” Nat. Commun. **3**, 1127 (2012).

²⁴M. Buzzacchi, I. Pagonabarraga, and N. B. Wilding, “Polydisperse hard spheres at a hard wall,” J. Chem. Phys. **121**, 11362–11373 (2004).

²⁵I. Pagonabarraga, M. Cates, and G. Ackland, “Local size segregation in polydisperse hard sphere fluids,” Phys. Rev. Lett. **84**, 911 (2000).

²⁶D. A. Kofke and E. D. Glandt, “Nearly monodisperse fluids. i. Monte Carlo simulations of lennard-jones particles in a semigrand ensemble,” J. Chem. Phys. **87**, 4881–4890 (1987).

²⁷N. Wilding and A. Bruce, “Freezing by Monte Carlo phase switch,” Phys. Rev. Lett. **85**, 5138 (2000).

²⁸N. B. Wilding, “Solid-liquid coexistence of polydisperse fluids via simulation,” J. Chem. Phys. **130** (2009).

²⁹F. Smallenburg, G. Del Monte, M. de Jager, and L. Filion, “A simple and accurate method to determine fluid–crystal phase boundaries from direct coexistence simulations,” J. Chem. Phys. **160**, 224109 (2024).

³⁰F. Smallenburg, “Efficient event-driven simulations of hard spheres,” Euro. Phys. J. E **45**, 22 (2022).

³¹B. J. Alder and T. E. Wainwright, “Studies in molecular dynamics. ii. behavior of a small number of elastic spheres,” J. Chem. Phys. **33**, 1439 (1960).

³²L. Berthier, E. Flenner, C. J. Fullerton, C. Scalliet, and M. Singh, “Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids,” J. Stat. Mech.: Theory Exp. **2019**, 064004 (2019).

³³E. A. Peters and G. de With, “Rejection-free Monte Carlo sampling for general potentials,” Phys. Rev. E **85**, 026703 (2012).

³⁴A. Castagnède, L. Filion, and F. Smallenburg, “Fast event-driven simulations for soft spheres: from dynamics to laves phase nucle-

ation,” J. Chem. Phys. **161**, 024116 (2024).

³⁵P.-R. Wolde *et al.*, “Simulation of homogeneous crystal nucleation close to coexistence,” Faraday Discuss. **104**, 93–110 (1996).

³⁶M. Buzzacchi, P. Sollich, N. B. Wilding, and M. Müller, “Simulation estimates of cloud points of polydisperse fluids,” Physical Review E—Statistical, Nonlinear, and Soft Matter Physics **73**, 046110 (2006).

³⁷T. Boublík, “Hard-sphere equation of state,” J. Chem. Phys. **53**, 471–472 (1970).

³⁸G. Mansoori, N. F. Carnahan, K. Starling, and T. Leland Jr, “Equilibrium thermodynamic properties of the mixture of hard spheres,” J. Chem. Phys. **54**, 1523–1525 (1971).

³⁹N. F. Carnahan and K. E. Starling, “Equation of state for nonattracting rigid spheres,” J. Chem. Phys. **51**, 635–636 (1969).

⁴⁰S. Pieprzyk, M. N. Bannerman, A. C. Brańka, M. Chudak, and D. M. Heyes, “Thermodynamic and dynamical properties of the hard sphere system revisited by molecular dynamics simulation,” Phys. Chem. Chem. Phys. **21**, 6886–6899 (2019).

⁴¹M. de Jager, C. Vega, P. Montero de Hijes, F. Smallenburg, and L. Filion, “Statistical mechanics of crystal nuclei of hard spheres,” J. Chem. Phys. **161**, 184501 (2024).

⁴²R. L. Davidchack and B. B. Laird, “Simulation of the hard-sphere crystal–melt interface,” J. Chem. Phys. **108**, 9452–9462 (1998).

⁴³R. Shuttleworth, “The surface tension of solids,” Proc. Phys. Soc. **63**, 444 (1950).

⁴⁴N. Di Pasquale and R. L. Davidchack, “Shuttleworth equation: A molecular simulations perspective,” J. Chem. Phys. **153**, 154705 (2020).

⁴⁵P. Sollich and M. E. Cates, “Projected free energies for polydisperse phase equilibria,” Phys. Rev. Lett. **80**, 1365 (1998).

⁴⁶P. Sollich and N. B. Wilding, “Polydispersity induced solid–solid transitions in model colloids,” Soft Matter **7**, 4472–4484 (2011).

⁴⁷A. Castagnède, L. Filion, and F. Smallenburg, “Data package for ‘freezing line of polydisperse hard spheres via direct-coexistence simulation’,” (2025), https://doi.org/10.5281/zenodo.15470568.
# Supplementary Material: Freezing line of polydisperse hard spheres via direct coexistence simulations

Antoine Castagnède,<fnref n="1" />,ª) Laura Filion,<fnref n="2" /> and Frank Smallenburg<fnref n="1" />

## I. GENERATING A DETERMINISTIC GAUSSIAN SIZE DISTRIBUTION

We consider the freezing line of polydisperse mixtures in which the particle sizes follow a Gaussian distribution:

$$ \mathcal{P}(\sigma) = \frac{1}{p\bar{\sigma}\sqrt{2\pi}} \exp\left[-\frac{1}{2}\left(\frac{\sigma - \bar{\sigma}}{p\bar{\sigma}}\right)^2\right]. \quad (1) $$

Here, $\bar{\sigma}$ is the mean particle diameter, and $p$ is the polydispersity, defined as the ratio between the standard deviation and the mean of the distribution.

To avoid any effects of randomness in the size distribution of the particles on simulations in the canonical ensemble, it is convenient to generate particle sizes deterministically. If we are generating $N$ particles according to a size distribution $P(\sigma)$, we can do this by generating sizes using the inverse cumulative distribution function $\Pi^{-1}(x)$, defined as the inverse of the cumulative distribution function

$$ \Pi(\sigma) = \int_{-\infty}^{\sigma} dx \mathcal{P}(x). \quad (2) $$

We can then choose

$$ \sigma_i = \Pi^{-1}\left(\frac{i-1/2}{N}\right), \quad (3) $$

as our set of particle sizes, where $i = 1, ..., N$. For a Gaussian distribution with polydispersity $p$, this results in

$$ \Pi(\sigma) = \frac{1}{2} \left( 1 + \text{Erf} \left( \frac{\sigma - \bar{\sigma}}{p\bar{\sigma}\sqrt{2}} \right) \right) \quad (4) $$

$$ \Pi^{-1}(x) = \bar{\sigma} \left( 1 + p\sqrt{2} \text{Erf}^{-1}(2x-1) \right). \quad (5) $$

Note that due to the finite number of particles, this results in a distribution with a slightly lower polydispersity than $p$ for small systems. However, for a system of $N = 1000$ particles the relative deviation in $p$ is already on the order of ~ 0.1%, and it further decreases (approximately $\propto 1/N$) for larger system sizes. Hence, we consider this deviation negligible.

## II. POLYDISPERSE MIXTURES IN THE CANONICAL ENSEMBLE

In the canonical ensemble, a polydisperse system of particles is characterized by a fixed number of particles $N$, volume $V$, temperature $T$, and a size distribution $\{\sigma\}$. We can write the associated Helmholtz free energy as:

$$ F(N, V, T, \{\sigma\}) = F_{\text{mix}} + F_{\text{conf}} \quad (6) $$

$$ F_{\text{conf}}(N, V, T, \{\sigma\}) = -k_B T \log Z_{\text{conf}}(N, V, T, \{\sigma\}) \quad (7) $$

$$ Z_{\text{conf}}(N, V, T, \{\sigma\}) = \frac{1}{\Lambda^3 N!} \int d\mathbf{r}^N \exp(-\beta U(\mathbf{r}^N, \{\sigma\})). \quad (8) $$

Here, $F_{\text{mix}}$ and $F_{\text{conf}}$ are the mixing and configurational parts of the free energy, respectively, $k_B$ is Boltzmann's constant, and $\beta = 1/k_B T$. Additionally, $U(\mathbf{r}^N, \{\sigma\})$ is the energy of a given configuration, and $\Lambda$ is the thermal wavelength.

For a macroscopic system ($N \to \infty$) with a size distribution $P(\sigma)$, we can define the mixing free energy by discretizing the size distribution<fnref n="1" />:

$$ \beta F_{\text{mix}} = N \sum_{i} x_i \log x_i, \quad (9) $$

where

$$ x_i = \mathcal{P}(i\lambda)\lambda \quad (10) $$

is the fraction of particles within some small interval $\lambda$ around a discretized size $\sigma_i^{\text{disc}} = i\lambda$. In the limit of small intervals, this can also be written as:

$$ \beta F_{\text{mix}} = N \int d\sigma \mathcal{P}(\sigma) \log(\mathcal{P}(\sigma)\lambda). \quad (11) $$

Importantly, this mixing entropy is independent of the choice of $\lambda$ up to a constant shift in the free energy per particle, which will not influence phase behavior<fnref n="1" />.

The chemical potential of a particle of size $\sigma_i$ is then defined as

$$ \mu(\sigma_i) = \left( \frac{\partial F}{\partial N_i} \right)_{V,T,N_i \neq i} \quad (12) $$

$$ = \mu_{\text{mix}}(\sigma_i) + \mu_{\text{conf}}(\sigma_i). \quad (13) $$

Here, the mixing term $\mu_{\text{mix}}(\sigma_i)$ is given by

$$ \mu_{\text{mix}}(\sigma_i) = k_B T \log x_i \quad (14) $$

$$ = k_B T \log(\mathcal{P}(\sigma_i)\lambda). \quad (15) $$

<footnote n="1">¹) Université Paris-Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France</footnote>

<footnote n="2">²) Soft Condensed Matter and Biophysics, Debye Institute for Nanomaterials Science, Utrecht University, Utrecht, Netherlands</footnote>

<footnote>ª) Author to whom correspondence should be addressed: antoine.castagnede@universite-paris-saclay.fr</footnote>
Note that similar to the free energy, a change in $\lambda$ only results in a trivial shift in the chemical potential, which again does not influence phase behavior. For the purpose of this paper, we will simply set $\lambda$ equal to the average particle size in the fluid phase $\bar{\sigma}$.

### III. MEASURING RELATIVE CHEMICAL POTENTIALS

We would like to measure the relative chemical potential function $\delta\mu(\sigma)$ in a simulation in the canonical ensemble. Here, $\delta\mu(\sigma) = \mu(\sigma) - \mu(\sigma_{ref})$ is the chemical potential of a particle of size $\sigma$ relative to that of a particle of an arbitrary reference size $\sigma_{ref}$. To this end, we first consider how changing the size of a single particle changes the canonical free energy of the system.

#### A. Effect of changing one particle's size on the free energy

We can now consider the effect of a change to the size of particle $i$, such that its size changes from $\sigma_i$ to $\sigma_i + \epsilon$. This affects both the mixing free energy and the configurational part. The total change in free energy can be written as:

$$ \Delta F = \Delta F_{mix} + \Delta F_{conf}. \quad (16) $$

Note that $\Delta F$ can also be regarded as the effect of taking out a particle of size $\sigma_i$ and putting in a particle of size $\sigma_i + \epsilon$. In other words,

$$ \Delta F = \mu(\sigma_i + \epsilon) - \mu(\sigma_i). \quad (17) $$

For the mixing free energy, we can simply consider removing a particle of size $\sigma_i$ and inserting a particle of size $\sigma_i + \epsilon$. Using the mixing part of the chemical potential, this yields:

$$ \beta\Delta F_{mix} = \beta\mu_{mix}(\sigma_i + \epsilon) - \beta\mu_{mix}(\sigma_i). \quad (18) $$

For the configurational part of the free energy:

$$ \beta\Delta F_{conf} = -\log \frac{\int d\mathbf{r}^N \exp(-\beta U(\mathbf{r}^N, \{\tilde{\sigma}\}))}{\int d\mathbf{r}^N \exp(-\beta U(\mathbf{r}^N, \{\sigma\}))}, \quad (19) $$

where $\{\tilde{\sigma}\}$ is the set of particle sizes after the change. We can rewrite the numerator:

$$ \begin{aligned} \beta\Delta F_{conf} &= -\log \frac{\int d\mathbf{r}^N \exp(-\beta\Delta U) \exp(-\beta U(\mathbf{r}^N, \{\sigma\}))}{\int d\mathbf{r}^N \exp(-\beta U(\mathbf{r}^N, \{\sigma\}))}, \\ &= -\log \langle \exp(-\beta\Delta U) \rangle_{\{\sigma\}}, \end{aligned} \quad (20) $$

where $\Delta U$ is the change in potential energy due to the change in the size of particle $i$, and the subscript on the angular brackets indicates that the ensemble average is taken at the original distribution of sizes $\{\sigma\}$. For hard spheres, $\Delta U$ is either infinite or zero, depending on

whether the particle size change creates an overlap, and hence $\exp(-\beta\Delta U)$ evaluates to either 0 or 1.

We now consider an infinitesimal change to the size of particle $i$, by taking the limit $\epsilon \to 0$. From Eq. 17, we obtain:

$$ \frac{\partial F}{\partial \sigma_i} = \mu'(\sigma_i) = \delta\mu'(\sigma_i). \quad (21) $$

Additionally, from Eqs. 18 and 20, we obtain:

$$ \frac{\partial F}{\partial \sigma_i} = \mu'_{mix}(\sigma_i) + \left\langle \frac{\partial U}{\partial \sigma_i} \right\rangle. \quad (22) $$

Since $\mu'_{mix}$ is known from the size distribution, we can determine the derivative of the chemical potential $\delta\mu'(\sigma_i)$ for any particle size, by measuring the ensemble average $\langle \frac{\partial U}{\partial \sigma_i} \rangle$ for different particle sizes $\sigma_i$. Subsequently, $\delta\mu(\sigma_i)$ can be trivially obtained by integrating $\delta\mu'(\sigma_i)$. For continuous potentials, $\langle \frac{\partial U}{\partial \sigma_i} \rangle$ can be measured straightforwardly. However, since our system of interest consists of hard spheres, the energy of the system is not a continuously differentiable function of the particle sizes, and slightly more care needs to be taken to measure $\delta\mu$.

#### B. Measuring $\delta\mu(\sigma)$ for hard spheres

The hard-sphere pair potential depends on the particle sizes as

$$ U_{ij}(\mathbf{r}_{ij}, \sigma_i, \sigma_j) = \tilde{U}^{\text{HS}} \left( r_{ij} - \frac{\sigma_i}{2} + \frac{\sigma_j}{2} \right). \quad (23) $$

Hence, it is more convenient to work in terms of the radius $R_i = \sigma_i/2$ of the particle than the diameters:

$$ \frac{\partial U}{\partial \sigma_i} = \frac{1}{2} \frac{\partial U}{\partial R_i} = -\frac{1}{2} \langle f_{R_i} \rangle, \quad (24) $$

where $f_{R_i}$ is the force experienced by the radius of particle $i$. We can write this average as a time-average:

$$ f_{R_i} = \frac{-1}{t_b - t_a} \int_{t_a}^{t_b} dt \sum_j \frac{\partial U_{ij}}{\partial R_i} \quad (25) $$

$$ = \frac{-1}{t_b - t_a} \int_{t_a}^{t_b} dt \sum_j |\mathbf{f}_{ij}| \quad (26) $$

$$ = \frac{-1}{t_b - t_a} \sum_c m |(\delta \mathbf{v}_i)_c|, \quad (27) $$

Here, $\mathbf{f}_{ij}(t) = -\frac{\partial U_{ij}}{\partial \mathbf{r}_{ij}}$ is the force on particle $i$ due to particle $j$. We have used the fact that the hard-sphere potential is repulsive (i.e. $\frac{\partial U_{ij}}{\partial \sigma_i} \ge 0$) and acts along the radius of each particle. In the last step, $(\delta \mathbf{v}_i)_c$ is the change in velocity of particle $i$ during collision $c$, $m$ is mass of a particle (taken equal for all particles), and the sum is taken over all collisions in the time interval $[t_a, t_b]$.
#### IV. POLYDISPERSE MIXTURES WITH DYNAMIC PARTICLE SIZES

We are interested in the equilibrium behavior of polydisperse mixtures in three dimensions. For our direct coexistence simulations, it is helpful to consider these systems in the semi-grand-canonical ensemble, i.e. at constant number of particles $N$, volume $V$, temperature $T$, and relative size-dependent chemical potential $\delta\mu(\sigma)$. Analogous to the “continuous time swap” approach of Ref. 2, we can consider the semigrand ensemble as a description of a system where each particle has three positional degrees of freedom and one for its radius. The Hamiltonian for this system can then be written as:

$$H = \sum_i \left[ \frac{p_i^2}{2m_i} + \frac{\wp_i^2}{2M} + V(R_i) \right] + \sum_{i<j} U_{ij}(\mathbf{r}_{ij}, R_i, R_j), \quad (28)$$

where $\mathbf{p}_i$ is the translational momentum of particle $i$ with mass $m$, $\wp_i$ is the momentum associated with its radius $R_i$, $M$ is the corresponding mass, $V(R_i) = -\delta\mu(\sigma_i = 2R_i)$ is the external field controlling the particle sizes, $U_{ij}$ represents the pair interaction, and $\mathbf{r}_{ij}$ is the vector connecting particles $i$ and $j$.

The canonical partition function for this system can then be written as

$$Z(N, V, T, \delta\mu(\sigma)) = \frac{1}{h^{4N}N!} \int d\mathbf{r}^N d\mathbf{p}^N d\sigma^N d\wp^N \exp(-\beta H), \quad (29)$$

where $h$ is a discretization constant used to ensure $Z$ is dimensionless, and $\beta = 1/k_B T$ with $k_B$ Boltzmann's constant. Integrating out the (translational and size) momenta yields:

$$Z = \frac{1}{\Lambda^{3N}\lambda^N N!} \int d\mathbf{r}^N d\mathbf{R}^N \exp\left(-\beta \sum_{i<j} U_{ij}\right) \times \exp\left(\beta \sum_i \delta\mu(R_i)\right), \quad (30)$$

where $\Lambda = h/\sqrt{2\pi mk_B T}$ is the usual thermal wavelength, $\lambda = h/\sqrt{2\pi M k_B T}$ is its analogue for the size degree of freedom. Note that apart from the irrelevant prefactor $1/\lambda^N$, this corresponds to the semigrand-canonical partition function that one would obtain for a system which can swap particles (at fixed total $N$) with an external particle reservoir where the chemical potential of a particle of size $\sigma$ is given by $\delta\mu(\sigma) + c$, where $c$ is an arbitrary constant. Hence, a simulation in the canonical ensemble governed by the Hamiltonian in Eq. 28 will correctly sample the semigrand ensemble governed only by the particle interactions and an imposed relative chemical potential $\delta\mu(\sigma)$.

For systems with continuous interaction potentials, the Hamiltonian in Eq. 28 gives rise to equations of motion that can be directly implemented in standard (time-driven) molecular dynamics schemes. However, for hard

spheres, the interaction potential is discontinuous, and hence we will use event-driven simulations instead.

This poses two problems. First, we will need to derive a collision rule for our variable-size particles, which will update not only the translational velocities of the colliding particles, but also their size velocities $\dot{\sigma}_i$. Second, we need to predict and handle “collisions” with the external potential $V(\sigma_i)$, which will affect only the radial velocities of the particles. We will discuss these two topics in the following two subsections.

#### A. Pair collision rules

For a pair potential of the form $U_{ij} = U(r_{ij} - R_i - R_j)$, the Hamiltonian in Eq. 28 gives rise to the following equations of motion:

$$\dot{\mathbf{r}}_i = \frac{\partial H}{\partial \mathbf{p}_i} = \mathbf{p}_i/m_i \quad (31)$$

$$\dot{R}_i = \frac{\partial H}{\partial \wp_i} = \wp_i/M \quad (32)$$

$$\dot{\mathbf{p}}_i = -\frac{\partial H}{\partial \mathbf{r}_i} = \sum_{i \neq j} \mathbf{f}_{ij} \quad (33)$$

$$\dot{\wp}_i = -\frac{\partial H}{\partial R_i} = -V'(R_i) + \sum_{i \neq j} \mathbf{f}_{ij} \cdot \hat{\mathbf{r}}_{ij}. \quad (34)$$

During a pair collision, the changes in $\dot{\mathbf{p}}_i$ and $\dot{\wp}_i$ are fully dominated by the pairwise forces, and hence the effects of the external field $V$ can be ignored. To determine the collision rule, we can regularize the interaction as

$$U_{ij}(\mathbf{r}_i, \mathbf{r}_j, R_i, R_j) = -\frac{\xi(\mathbf{r}_{ij}, R_i, R_j)}{\epsilon} \Theta(-\xi(\mathbf{r}_{ij}, R_i, R_j)), \quad (35)$$

where

$$\xi(\mathbf{r}_{ij}, R_i, R_j) = |\mathbf{r}_{ij}| - R_i - R_j \quad (36)$$

is the surface-to-surface distance between the spheres, $\Theta$ is the Heaviside step function, and $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$. This continuous pair potential represents a short-range repulsion for overlapping particles ($\xi < 0$), which converges to the hard-sphere potential in the limit $\epsilon \to 0$. The equations of motion for the velocities during the collision then become:

$$\dot{\mathbf{p}}_1 = \Theta(-\xi) \frac{\hat{\mathbf{r}}_{12}}{\epsilon} \quad (37)$$

$$\dot{\mathbf{p}}_2 = -\Theta(-\xi) \frac{\hat{\mathbf{r}}_{12}}{\epsilon}. \quad (38)$$

$$\dot{\wp}_i = -\frac{\Theta(-\xi)}{\epsilon} \quad (39)$$

During the collision, i.e. during the period where $\xi < 0$, we obtain:

$$\dot{\xi} = \mathbf{v}_{12} \cdot \hat{\mathbf{r}}_{12} - \frac{\wp_1 + \wp_2}{M} \quad (40)$$

$$\ddot{\xi} = \frac{1}{\epsilon} \left( \frac{1}{m} + \frac{1}{M} \right) + \mathbf{v}_{12} \cdot \frac{\partial \hat{\mathbf{r}}_{12}}{\partial t} \quad (41)$$
In the limit of small $\epsilon$, we can neglect the last term in $\ddot{\xi}$. The acceleration is then constant during the collision, and we can solve for the total duration of the collision:

$$ \Delta t = \frac{\dot{\xi}(t = t_{col})}{\ddot{\xi}(t = t_{col})} = \epsilon \frac{b(t_{col})}{\frac{1}{m} + \frac{1}{M}}, \quad (42) $$

where

$$ b(t) = \mathbf{v}_{12}(t) \cdot \hat{\mathbf{r}}_{12}(t) - \frac{\varphi_1(t) + \varphi_2(t)}{M} \quad (43) $$

is the surface-to-surface velocity of the two particles along the $\hat{\mathbf{r}}_{12}$ direction, which we evaluate at the time $t_{col}$ where the collision starts. The change in momenta is then given by

$$ \Delta \mathbf{p}_1 = \frac{b(t_{col}) \hat{\mathbf{r}}_{12}}{\frac{1}{m} + \frac{1}{M}} \quad (44) $$

$$ \Delta \mathbf{p}_2 = - \frac{b(t_{col}) \hat{\mathbf{r}}_{12}}{\frac{1}{m} + \frac{1}{M}} \quad (45) $$

$$ \Delta \varphi_i = \frac{b(t_{col})}{\frac{1}{m} + \frac{1}{M}} \quad (46) $$

Note that this result is independent of $\epsilon$ and hence also holds in the limit of $\epsilon \to 0$, where the collision becomes a hard-core collision. This gives us a set of collision rules for polydisperse hard spheres of variable size. These rules conserve the translational momentum and the total kinetic energy (translational + radial) of the system.

The introduction of the radial “mass” $M$ introduces a free parameter into the simulation, given by $\alpha = m/M$. Low $\alpha$ (high $M$) will result in relatively slower motion of the particle sizes. The choice of $\alpha$ should not affect any static equilibrium properties (including phase behavior) of the system, but will influence dynamics. As the growing and shrinking of particles is not physical anyway, we arbitrarily set $\alpha = 1$.

Finally, we note that predicting the collision time for hard spheres with a constant radial velocity is largely analogous to predicting collisions of hard spheres of constant size. Specifically, it involves solving a quadratic equation

$$ |\mathbf{r}_{ij}(t)|^2 = \sigma_{ij}(t)^2, \quad (47) $$

with $\sigma_{ij}(t) = R_i(t) + R_j(t)$ the time-dependent contact distance between the two spheres. Since both sides of the equation can be written as second-order polynomials of $t$, this can be analytically solved for $t$, yielding:

$$ t_{\text{col}} - t_{\text{current}} = \frac{-b - \sqrt{b^2 - a(r_{ij}^2 - \sigma_{ij}^2)}}{a}, \quad (48) $$

$$ a = (v_{ij}^2 + \dot{\sigma}_{ij}^2), \quad (49) $$

where all time-dependent quantities are evaluated at the current simulation time $t_{current}$. Note that future collisions only occur when $b < 0$ or $a < 0$, and additionally the discriminant (inside the square root) is larger than zero. Otherwise, no future collisions are expected.

## B. Collision rules with the field $V(R)$

Traditionally, event-driven simulations are commonly used for potentials and energy fields that are piecewise constant functions, such as hard spheres or other shapes, square-well interactions, and systems of these particles interacting with with hard walls or external stepwise potential fields. For these types of systems, the momentum of a particle only changes at discrete events in time. However, in the canonical ensemble, continuous potentials can also be addressed using an event-driven approach. Specifically, Peters and De With³ have introduced a rejection free event-driven Monte Carlo (EDMC) algorithm for simulating particles with continuous interaction potentials. In this approach, the interaction between pairs of particles gives rise to stochastic collisions, with the collision distance determined based on a random number drawn from a Boltzmann distribution at the imposed temperature $T$. Between these collisions, the particles experience free flight as normal in an event-driven molecular dynamics (EDMD) simulation. As a result, it has been shown that the EDMC approach can be an efficient alternative to conventional MD simulations in the canonical ensemble, assuming the interaction potential is short-ranged and easily invertible⁴.

Here, we adapt this approach to handle the effects of the field $V(R)$ that controls the size distribution of the particles. Specifically, whenever we predict the next collision for a given particle, we also consider the possibility of colliding with the field $V(R)$. To this end, we draw a random energy $\Delta V$ from a Boltzmann distribution, and we find the radius $R_{col}$ such that the total increase in energy encountered between the current particle radius $R(t)$ and $R_{col}$ is equal to $\Delta V$. As proven in Ref. 3, this approach will correctly sample the canonical ensemble.

In practice, our function $V(R)$ is taken to be a cubic polynomial (which we find fits our data for $\delta\mu(\sigma)$ excellently, see main text). Additionally, we only allow $R$ to vary within a region $[R_{min}, R_{max}]$ of permitted values, chosen large enough to ensure that the probability of hitting either edge is negligible. Moreover, we find that $V(R)$ is always either a monotonically decreasing function of $R$ or has a single minimum at $R_{well}$ (and no maximum) within the allowed region. We therefore have to consider three cases (see Fig. 1):

1. The particle’s radius is going “uphill” (i.e. $V(R)$ is increasing as a function of time). An example is shown in Fig. 1a for a shrinking particle. In this case, the particle is allowed to continue going uphill until an amount of energy equal to $\Delta V$ has been expended. In other words, the radius at which the collision with the field occurs is such that $V(R_{col}) = V_{col} = V(R) + \Delta V$, with $R$ the current particle radius. If this radius lies outside the interval $[R_{min}, R_{max}]$, the collision occurs at the edge of this interval, instead. For a growing particle whose radius starts off larger than $R_{well}$, the procedure is
<figure><img src="image_12.png" /><figcaption>FIG. 1. Illustration of the prediction of particle-field collisions in the semigrand-canonical ensemble simulations. The solid red lines depict the field $V(R)$ felt by the particle's radius. The black dot and arrow denote the particle's current radius and whether the particle is currently growing or shrinking. The blue dot denotes the predicted value of $R$ at the moment of collision (which may depend on the stochastically chosen value of $\Delta V$). The three depicted cases correspond to a) uphill motion, b) downhill motion without a minimum, and c) downhill motion with a minimum.</figcaption></figure>

exactly the same.

2. The particle's radius is going “downhill” (i.e. $V(R)$ is decreasing as a function of time), and will reach an edge of $[R_{\min}, R_{\max}]$ before encountering a local minimum. For our systems, this only applies to growing particles (Fig. 1b). The field collision is then scheduled at $R_{\text{col}} = R_{\text{max}}$.

3. The particle's radius is going downhill towards the minimum in $V(R)$ at $R_{\text{well}}$. In this case, we ignore the downhill part of the trajectory, and we find the collision radius corresponding to $V(R_{\text{col}}) = V_{\text{col}} = V(R_{\text{well}}) + \Delta V$. In other words, $\Delta V$ now measures the permitted increase in energy starting from the minimum³ (Fig. 1c). Again, if the collision radius falls outside $[R_{\min}, R_{\max}]$, the collision occurs at the edge of this interval instead.

Finally, for all three cases described above, once $R_{\text{col}}$ is determined, the associated collision time can be trivially determined as $(R_{\text{col}} - R)/\dot{R}$, making use of the fact that $\dot{R}$ is constant between collisions.

Note that since our expression for $V(R)$ is a cubic polynomial with real coefficients, solving $V(R_{\text{col}}) = V_{\text{col}}$ can be done analytically. This is done by transposing $V(R)$ into a depressed cubic of the form $t^3 + pt + q = 0$ (with $p$ and $q$ constants), and expressing the real roots using hyperbolic or trigonometric functions depending on their multiplicity. For more complex forms of $V(r)$, numerical root-finding methods or lookup tables would likely be necessary.

V. SYSTEM SIZE EFFECTS

We show in Fig. 2 the effect of system size on the determination of the coexistence pressure for hard spheres with a polydispersity of $p = 0.06$. As one might expect, larger systems lead to less noise. However, the obtained coexistence pressures are not strongly affected by system size.

<figure><img src="image_13.png" /><figcaption>FIG. 2. Pressure tensor components $P_{xx}$, $P_{yy}$, and $P_{zz}$ as a function of trial pressure $P^{\text{trial}}$ obtained from semi-grand direct coexistence simulations for system with 16000 particles, polydispersity $p = 0.06$, and varying system size. Error bars show one standard error. These systems were set so the fluid-crystal interface corresponds to the square lattice plane of an FCC crystal.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>FIG. 3. Surface excess concentration of particles per unit area with size $\sigma_i$ for an equimolar Gibbs dividing surface in a system with 16384 particles at coexistence. Here the crystal orientation is such that the hexagonal FCC planes are parallel to the elongated dimension of the simulation box (FCC$_b$). Polydispersities are reported on each line, all lines are vertically shifted for readability.</figcaption></figure>

## VI. PHASE COMPOSITION AND EXCESS SURFACE ABSORPTION AT COEXISTENCE

Fig. 3 shows the excess surface absorption for the FCC$_b$ interface. Additionally, we show in Fig. 4 the particle size distribution at coexistence in the fluid phase, FCC phase, and coexisting system (FCC$_a$ interface) for each of the investigated polydispersities. We note that the behavior of the coexisting system with FCC$_b$ interface is qualitatively very similar.

## VII. REFERENCES

¹Daan Frenkel. Why colloidal systems can be described by statistical mechanics: some not very original comments on the gibbs paradox. *Mol. Phys.*, 112(17):2325-2329, 2014.

²Ludovic Berthier, Elijah Flenner, Christopher J Fullerton, Camille Scalliet, and Murari Singh. Efficient swap algorithms for molecular dynamics simulations of equilibrium supercooled liquids. *J. Stat. Mech.: Theory Exp.*, 2019(6):064004, 2019.

³Elias AJF Peters and G de With. Rejection-free Monte Carlo sampling for general potentials. *Phys. Rev. E*, 85(2):026703, 2012.

⁴Antoine Castagnède, Laura Filion, and Frank Smallenburg. Fast event-driven simulations for soft spheres: from dynamics to laves phase nucleation. *J. Chem. Phys.*, 161(2):024116, 2024.
<figure><img src="image_15.png" /><figcaption>FIG. 4. Particle size distribution of direct coexistence simulations (FCC<sub>a</sub> interface, solid green lines) for the range of polydispersities explored in this work. In each direct coexistence simulation, the global density was chosen to be exactly halfway between the fluid and crystal coexistence densities. Also shown are the compositions of the bulk coexisting fluid (dashed orange lines) and crystal (dashed blue lines) phases.</figcaption></figure>