# Polydisperse Hard Spheres

###### tags: `PhD projects`

---

### Owners

Marjolein, Antoine, Frank, Laura

---

## Coexistence points (Antoine)

Coexistence conditions obtained from direct coexistence semi-grand canonical simulations of systems of 16000 particles, with a crystalline slab configuration and a $(100)$ crystal plane facing the fluid.

| $p$ | $\beta P_\textrm{coex} \bar{\sigma}^3$ | $\rho_F \bar{\sigma}^3$ |
| ------ | ----- | ----- |
| 0.02 | 11.742 | 0.9418 |
| 0.04 | 12.367 | 0.9510 |
| 0.06 | 13.379 | 0.9643 |
| 0.08 | 14.778 | 0.9802 |

Composition of the crystal phase at coexistence for 16000 particles systems:

| $p$ | $\langle R \rangle_\chi / \bar{\sigma}$ | $SD(R)_\chi / \bar{\sigma}$ |
| ------ | ----- | ----- |
| 0.02 | 0.50069 | 0.00980 |
| 0.04 | 0.50271 | 0.01836 |
| 0.06 | 0.50594 | 0.02488 |
| 0.08 | 0.51061 | 0.02876 |
| 0.12 | 0.52543 | 0.02883 |
| 0.14 | 0.53609 | 0.02660 |

---

## Equilibrium nuclei (Antoine)

### Turning a monodisperse HS nucleus into a polydisperse one

Every attempt to re-grow a monodisperse system into a 4% polydisperse one lead to it fully melting. I tried different approaches, such as preferentially putting larger particles inside the nucleus, but that requires for the system to grow from too small particle sizes, which makes the nucleus melt immediately.
I am discarding this route to make a polydisperse nucleus and preffering the second one:

### Partial melting of a chunck of polydisperse FCC crystal

This is the route used by Marjolein to produce monodisperse nuclei.
I am now looking at 4% polydisperse systems for packing fractions $\in [ 0.545 0.550 ]$ and trying to melt them in EDMD with size-bias swaps.

Nothing conclusive so far -- box too tight?

### Umbrella Sampling - Particle size distrib. of critical nuclei

I have obtained free energy barriers for 4% and 6% polydisperse systems and I am still running some biased simulations at the top of the barriers to improve statistics.

Here are the average trends for $\langle\sigma(r)\rangle$ and $\text{SD}(\sigma)_\text{nuc}$ for a 4% polydisperse system. It is not clear from these early results that we recover the trends Marjolein observed for dynamical nuclei. I should also emphasize that we are looking at relatively small nuclei, with $n*$ reported below.


> 4% polydisperse

| $\eta$ | 0.544 | 0.545 | 0.546 | 0.547 | 0.548 | 0.549 | 0.550 |
| ----- | ----- | ----- | ----- | ----- | ----- | ----- | ----- |
| $n^*$ | 97 | 84 | 86 | 77 | 76 | 64 | 54 |

---

## Brute force nucleation (Marjolein)

I will look at nucleation of 4% and 6% polydisperse systems. The 8% polydisperse system is "too glassy" and will probably take too long to nucleate (if it even nucleates...).

### Nuclation rate
The simulations for 4% polydispersity are all finished. I considered the packing fraction range [0.544, 0.550], and performed 100 brute force simualtions of max $10^6\tau$ for each packing fraction. I stop the simulations once the nucleus reaches a size larger than 800 particles. For packing fractions [0.545,0.550], all 100 simulations managed to nucleate. For 0.544, 92/100 simulations nucleated.

The simulations for 6% polydispersity are still running. I considered the packing fraction range [0.558, 0.562], and again performed 100 brute force simulations of max $10^6\tau$ for each packing fraction. For packing fractions 0.562 and 0.561, all 100 simulations managed to nucleate. For 0.560, 0.559, and 0.558, we have that, respectively, 97/100, 89/100, and 51/100 simulations nucleated.

Underneath are the (preliminary) results on the nucleation rate. To compute the rates, I used
$$
J = \frac{1}{\langle t_\text{nuc}\rangle V},
$$
where $V$ is the volume of the box and $\langle t_\text{nuc}\rangle$ is the average time it takes the system to nucleate. This includes the time of the unsuccessful events and for the successful events it takes the last time for which the nucleus last has a size of 100 particles. Note that particles are classified as crystal when they have 7 or more crystal-like ten Wolde bonds ($d_6>0.7$).
I scale the rates with the long time diffusion time and compare them to the nucleation rates of monodisperse HS (from Frank) and to the rates from **[Wohler & Schilling, PRL 128 (2022)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.128.238001)** (from Tab. I and II, see further down below).


Underneath, I rescaled the packing fraction with the freezing packing fractions from Antoine. Note that the trends of the polydisperse systems do not collapse onto the trend of monodisperse HS.
(Note that for this I used that $\eta_F^\text{mono}=0.491753$, $\eta_F^\text{4%}=0.500316$, and $\eta_F^\text{5%}=0.506$. This last value I obtained from Fig. 3 in **[Bolhuis & Kofke, PR E 54 (1996)](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.54.634)**.)


Lastly, note that additionally replacing $\langle\sigma\rangle$ with $\sigma^\text{eff}$, which is given by $\sigma^\text{eff}=\langle\sigma\rangle(\eta_F^\text{HS}/\eta_F)^{1/3}$, will not improve the mapping. Since, the freezing packing fraction increases with increasing polydispersity, $\sigma^\text{eff}/\langle\sigma\rangle<1$. This will shift the rates in the wrong direction (down instead of up). And it will also be quite a small shift (3-6%).

> (Antoine)
> The following rates were obtained for 4% polydisperse systems from umbrella sampling simulations, for system sizes such that the critical nucleus does not exceed 10% of the total system size:
>
> | $\eta$ | $J \bar{\sigma}^5 / (6 D_L)$ |
> | ----- | -----|
> | 0.5275 | 1.85e-23 |
> | 0.5280 | 8.14e-23 |
> | 0.5300 | 5.19e-20 |
> | 0.5380 | 2.00e-12 |
> | 0.5450 | 2.65e-08 |
>
> *19/02/25:* Will update shortly with other nuc. rates predictions for 4% and 6% systems.

The rates from **[Wohler & Schilling, PRL 128 (2022)](https://doi.org/10.1103/PhysRevLett.128.238001)** come from the tables:

#### Other scaling options


For the above we used the freezing and melting packing fractions **[Bolhuis & Kofke, PR E 54 (1996)](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.54.634)** for 5% (as Schilling did). However, the melting packing fraction in particular is much higher (0.550) than what I would expect from Antoine's values for 4% and 6% (between 0.5385 and 0.5357). If we take a melting packing fraction of 0.537 for 5%, then we get the figure below. Note that Antoine's values show a coexistence region that shrinks with increasing polydispersity, which does not "help" in this effective scaling of the packing fracting.

#### Literature overview of nucleation rates

The results of Auer & Frenkel on monodisperse HS and 5% polydisperse (Gaussian) HS are reported in quite a few papers and with some varying values... Underneath is an overview of these different values.

Let's focus on the monodisperse HS first.
| $\eta$ | Auer & Frenkel (2001) </br>$\log_{10}(I\sigma^5/D_0)$ | Auer & Frenkel (2004) </br>$\log_{10}(I\sigma^5/D_0)$ |
| -------- | -------- | -------- |
| 0.5207 | -18.45 | -19.3 |
| 0.5277 | -12.96 | -13.5 |
| 0.5343 | -8.61 | -9.14 |

The 2nd column is from **[Auer & Frenkel, Nature 409 (2001)](https://doi.org/10.1038/35059035)** (read off Fig. 2).
The 3rd column is from **[Auer & Frenkel, JCP 120 (2004)](https://doi.org/10.1063/1.1638740)** (read from Tab. I). Note that this paper contains the same results as reported in **[Auer & Frenkel, Annu. Rev. Phys. Chem. 55 (2004)](https://doi.org/10.1146/annurev.physchem.55.091602.094402)**.

We can correct the above nucleation rates (i.e. use $D_L$ instead of $D_0$ and include a factor 1/6), such that we can compare them to other papers that report Auer & Frenkel's values. For this, I use $D_L(\eta)=D_0(1-\eta/0.58)^{2.6}$.

| $\eta$ | Auer & Frenkel (2001) </br>$\log_{10}(I\sigma^5/6D_L)$ | Auer & Frenkel (2004) </br>$\log_{10}(I\sigma^5/6D_L)$ | Filion et al (2004) </br>$\log_{10}(I\sigma^5/6D_L)$ | Wohler & Schilling (2022) </br>$\log_{10}(I\sigma^5/6D_L)$ |
| ------ | ------ | ------ | ------ | ------ |
| 0.5207 | -16.65 | -17.5 | -17.1 | -17.9 |
| 0.5277 | -11.02 | -11.6 | -11.1 | -12.3 |
| 0.5343 | -6.52 | -7.05 | -6.56 | -7.89 |

The 4th column is from **[Filion et al, JCP 133 (2010)](https://doi.org/10.1063/1.3506838)** (read off Fig. 8).
The 5th column is from **[Wohler & Schilling, PRL 128 (2022)](https://doi.org/10.1103/PhysRevLett.128.238001)** (read off Fig. 1).

Let us now turn to the 5% polydisperse HS.
| $\eta$ | Auer & Frenkel (2001) </br>$\log_{10}(I\sigma^5/D_0)$ | Auer & Frenkel (2004) </br>$\log_{10}(I\sigma^5/D_0)$ |
| -------- | -------- | -------- |
| 0.535 | -22.0 | -24.4 |
| 0.537 | -17.9 | -18.5 |
| 0.541 | -14.1 | -15.1 |
| 0.545 | -11.8 | -12.7 |
| 0.553 | -8.5 | -8.9 |

The 2nd column is from **[Auer & Frenkel, Nature 409 (2001)](https://doi.org/10.1038/35059035)** (read off Fig. 2), and the 3rd column is from **[Auer & Frenkel, JCP 120 (2004)](https://doi.org/10.1063/1.1638740)** (read off Fig. 11).

Again, if we correct these nucleation rates (i.e. use $D_L$ instead of $D_0$ and include a factor 1/6), then we can compare them to other papers that report Auer & Frenkel's values. For this, I again use $D_L(\eta)=D_0(1-\eta/0.58)^{2.6}$. This is gives similar results to Laura's version of Auer & Frenkel's values; yet, notice that they agree less well than the results of monodisperse HS.
| $\eta$ | Auer & Frenkel (2001) </br>$\log_{10}(I\sigma^5/6D_L)$ | Auer & Frenkel (2004) </br>$\log_{10}(I\sigma^5/6D_L)$ | Filion et al (2004) </br>$\log_{10}(I\sigma^5/6D_L)$ | Wohler & Schilling (2022) </br>$\log_{10}(I\sigma^5/6D_L)$ |
| ------ | ------ | ------ | ------ | ------ |
| 0.535 | -19.9 | -22.3 | -19.2 | -21.4 |
| 0.537 | -15.8 | -16.3 | -15.1 | -17.3 |
| 0.541 | -11.8 | -12.8 | -11.2 | -13.4 |
| 0.545 | -9.4 | -10.3 | -8.9 | -11.1 |
| 0.553 | -5.8 | -6.2 | -5.2 | -7.5 |

Underneath is an overview of the various figures in papers containing the nucleation rates.
| **[Auer & Frenkel, Nature 409 (2001)](https://doi.org/10.1038/35059035)** (Fig. 2) </br>Plots: $\log_{10}(I\sigma^5/D_0)$ | **[Auer & Frenkel, JCP 120 (2004)](https://doi.org/10.1063/1.1638740)** (Fig. 11) </br>Plots: $\log_{10}(I\sigma^5/D_0)$ |
| -------- | -------- |
|  |  |
| **[Auer & Frenkel, Annu. Rev. Phys. Chem. 55 (2004)](https://doi.org/10.1146/annurev.physchem.55.091602.094402)** (Fig. 2) </br>Plots: $\log_{10}(I\sigma^5/D_0)$ | **[Filion et al, JCP 133 (2010)](https://doi.org/10.1063/1.3506838)** (Fig. 8) </br>Plots: $\log_{10}(I\sigma^5/6D_L)$ |
|  |  |
| **[Wohler & Schilling, PRL 128 (2022)](https://doi.org/10.1103/PhysRevLett.128.238001)** (Fig. 1) </br>Plots: $\log_{10}(I\sigma^5/6D_L)$ | **[Zaccarelli et al, PRL 103 (2009)](https://doi.org/10.1103/PhysRevLett.103.135704)** (Fig. 5c) </br>Plots: $\tau_\text{nuc}/\tau_d$ </br>(0%=circles, 5%=downward triangles) |
|  |  |

Note that the range of supersaturations studied by **[Zaccarelli et al, PRL 103 (2009)](https://doi.org/10.1103/PhysRevLett.103.135704)** is higher than what we studied, i.e. $\eta\in[0.540,0.580]$ for 0% and $\eta\in[0.560,0.580]$ for 5%.

**[Gispen & Dijkstra, JCP 159 (2023)](https://doi.org/10.1063/5.0165159)**

#### Comparison with Auer & Frenkel's rates

Note that corrected nucleation rates for monodisperse HS of Auer & Frenkel (2001) now agree with Laura's version of their rates (red and blue open diamonds). For the 5% polydisperse HS, they almost agree (red and blue open triangles).

#### Better freezing packing fraction for 5% polydispersity

It bothered me that I used the freezing packing fraction from **[Bolhuis & Kofke, PR E 54 (1996)](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.54.634)** for 5% polydispersity, as our direct coexistence simulations are probably more accurate. That's why I tried fitting the trend in Antoine's freezing packing fractions. Underneath, you see the freezing packing fraction as a function of the polydispersity. The red dots are Antoine's data points and the dashed line is just the interpolation between the points obtained by Mathematica.


I think this interpolation looks great. So, if we use it for 5%, then we find $\eta_F=0.50496$. Although this only differs 0.001 from the value (i.e. 0.506) of **[Bolhuis & Kofke, PR E 54 (1996)](https://journals.aps.org/pre/pdf/10.1103/PhysRevE.54.634)**, it has a signficant effect on the mapping. As you can see below, with the old value of 0.506, the trend of 5% collapses approximately on top of the trend of 4%. However, with the new value of 0.50496, the trend of 5% lies neatly between the trends of 4% and 6%.

**OLD: Mapping using $\eta_F^\text{5%}=0.506$**

**NEW: Mapping using $\eta_F^\text{5%}=0.50496$**

#### Fitting the nucleation rates
In order to figure out a good fitting function for the nucleation rates, I first focussed on monodisperse HS, as we know everything about them. Since we actually know the height of the nucletion barrier, $\Delta G$, as a function of the fluid density from our paper on HS nuclei, we can just see what the dependense of $\Delta G$ on the fluid density is.
Below, I show the height of the nucleation barrier obtained from integrating $\Delta N$ (blue dots). Note that it is a log-plot. The dashed lines are fits of the form
$$
y = \frac{a}{(x-b)^c},
$$
where $a$, $b$, and $c$ are the fit parameters. The black dashed line takes all the data points into account; therefore, putting too much emphasis on correctly fitting the points close to the freezing density. The green dashed line takes the points with $\rho\sigma^3>0.965$ into account, and the red dashed line the points with $\rho\sigma^3>0.99$. Although the red dashed line is slightly better at capturing the tail of the trend, the green dashed line is better at capturing the entire trend. Moreover, the green dashed line makes a better estimation of the freezing density (fit parameter $b$).
Thus, for fitting the nuclation rates, I will take $c=2$.
Note that this makes sense if you use that $\Delta G=\frac{16\pi\gamma^3}{3(\Delta P)^2}$ and assume that, close to the coexistence density, $|\Delta P|$ is linear in $(\rho-\rho_F)$ and $\gamma$ is approximately constant.


Thus, to fit our nucleation rates, I use the fitting function
$$
\log\left(J\bar{\sigma}^3\tau_D\right) = \log(A) -\frac{a}{(\eta-\eta_F)^2},
$$
where $A$ and $a$ are fit parameters, and $\eta_F$ is the freezing packing fraction of the system under consideration. Note that $A$ represents the kinectic prefactor of the nucleation rate. The resulting fits are shown below.


Just out of curiousity, if we use a power 4 instead of 2, the fits look like this


So based on the point from umbrella sampling (black square), I would say that indeed the fit with $(\eta-\eta_F)^2$ is better than with $(\eta-\eta_F)^4$.

### Onset of nucleation

To study the onset of nucleation, I focus on the composition of the (successful) nuclei. For every run, I determine what the successful nucleus is and calculate the average and standard deviation of the particle diameter of the particles inside the nucleus, $\langle\sigma\rangle_\text{nuc}$ and $\text{SD}(\sigma)_\text{nuc}$.

By plotting these as a function of time for each nucleation event, I quickly saw that the nuclei generally contain more larger particles. Even starting from its very first occurance. Furthermore, the standard deviation of the particle diameter is generally smaller than in the fluid. Underneath are some examples for 4% polydispersity at $\eta=0.550$.

|  |  |
| -------- | -------- |
|  |  |

Initially, I wanted to try and overlay all nucleation events to get an average trend over time. However, the time it takes for the nucleus to grow from a few particles to 800 particles differs significantly over the different nucleation events. I will still think about a good way of overlaying the various nucleation events. But for the moment, I just made some density plots of the nucleus size vs $\langle\sigma\rangle_\text{nuc}$. Note that these only include the data from the succesful nuclei.
The black dashed line indicates $\langle\sigma\rangle_\text{nuc}/\langle\sigma\rangle=1.0$, and the color gradient goes from dark blue to dark yellow to light yellow. White means no data points.

| Mean diameter | Standard deviation |
| -------- | -------- |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |
|  |  |

We clearly see that the nuclei generally contain more larger particles and with a smaller spread in particle size. This is even more so for 6% polydispersity.

| Mean diameter | Standard deviation |
| -------- | -------- |
|  |  |

Next, I determined the average trend as a function of the nucleus size. Underneath are the resulting trends for the different packing fractions. The top row is for 4% and the bottom row for 6%. Note that, to reduce the noise of the trends a bit, I used a bin size of 25 particles for the nucleus size.
For $\langle\sigma\rangle_\text{nuc}$, it is difficult too say whether to degree of supersaturation has an effect. However, for $\text{SD}(\sigma)_\text{nuc}$ in the 4% system, increasing the supersaturation seems to shift the trend to slightly lower $\text{SD}(\sigma)_\text{nuc}$. In any case, the supersaturation does not have a large effect and we clearly see that the initial nuclei prefer larger particles.

| Average trend $\langle\sigma\rangle_\text{nuc}$ | Average trend $\text{SD}(\sigma)_\text{nuc}$ |
| -------- | -------- |
|  | |
|  |  |

Instead of studying the size distribution inside the nuclei, we can also study the onset of nucleation by looking at the size distribution of the particles inside the region where the nucleation event is about to happen. For each event, I pick this region using the center of mass of the critical nuclei and take a sphere of radius $2.2\langle\sigma\rangle$ around it. In practice this region contains around 45 particles. I then select the final snapshot before nucleation kicks off. I then determine the size distribution using the particles inside this region of the final snapshot before nucleation kicks off of each nucleation event.

Underneath are the resulting distributions for each supersaturation studied. The thick black line indicates the size distribution of the entire system and the thin blue lines indicate the distributions for the different packing fractions. The blue vertical dashed lines indicate the mean of theses distributions (these lines fall on top of each other as they approximately all have the same mean). For the 6% polydisperse system, we clearly see that the distributions are shifted towards slightly larger particles. For the 4% system, the shift is slighly smaller, but there is still a significant shift.
For 6%, the mean of the distributions is $1.0072(4)\sigma$. For 4%, the mean of the distributions is $1.0032(3)\sigma$.

| 4% | 6% |
| -------- | -------- |
|  |  |

I also looked at the radial composition of the nuclei. Even though I do not have that much data (because it is not in equilibrium), I still get surprisingly okay results. Underneath you can see the radial profiles of the average particle diameter (top) and density (bottom) for 4% polydispersity at a global packing fraction of 0.549. The gray area indicates the rough size of a nucleus of 800 particles and the dashed line in the bottom plot indicates the global density.
To obtain these trends, I first compute the average radial profiles for each nucleation event. I do this by considering the snapshots of that event for which the nucleus has a size $600<N_\text{nuc}<820$. Next, I take the average over all nucleation events.


I did this for all supersaturations of both 4% and 6% polydispersity. Underneath are the results. We clearly see that the core of the nucleus contain on average more larger particles. However, notice that this does not extend to the surface of the nucleus. Although the profiles suggest that the fluid--crystal interface is located at around $r/\langle\sigma\rangle\approx5$, the trend of $\langle\sigma(r)\rangle$ decays back to the global $\langle\sigma\rangle$ before that. This suggests that the initial (precritical?) nuclei strongly prefer larger particles, whereas the nuclei in the growth stage just "eat up" any surrounding fluid particles.
Note that a sphere of radius $2.5\langle\sigma\rangle$ to $3.0\langle\sigma\rangle$ contains approximately 65 to 115 particles. So the trend in $\langle\sigma(r)\rangle$ indeed suggest that the preference for larger sized particles significantly decreases once the nucleus has reached a critical size (around 60-100 according to Antoine for the 4% state points studied by us).

| 4% | 6% |
| -------- | -------- |
|  |  |

Additionally, I did the same analysis but then for nuclei in the range $100<N_\text{nuc}<200$. Underneath are the results. The gray area indicates the rough size of a nucleus of 200 particles. For these nuclei, we clearly see that the density and $\langle\sigma(r)\rangle$ follow similar radial decays.

| 4% | 6% |
| -------- | -------- |
|  |  |