{'fraction_non_alphanumeric': 0.029944396265877672, 'fraction_numerical': 0.010916696424016731, 'mean_word_length': 4.762492651381541, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hard sphere crystal nucleation rates: Reconciliation of simulation and experiment

Over the past two decades, a large number of studies addressed the topic of crystal nucleation in suspensions of hard spheres. The shared result of all these efforts is that, at low super-saturations, experimentally observed nucleation rates and numerically computed ones differ by more than ten orders of magnitude. We present precise simulation results of crystal nucleation rate densities in the meta-stable hard sphere liquid. To compare these rate densities to experimentally measured ones, we propose an interpretation of the experimental data as a combination of nucleation and crystal growth processes (rather than purely the nucleation process). This interpretation may resolve the long standing dispute about the differing rates.

Precise predictions of crystallization and solidification kinetics are highly sought after in scientific fields like atmospheric physics, metallurgy and chemical engineering. Arguably, the simplest model to study the liquid-to-crystal transition is the hard sphere fluid. Following the seminal computer simulation study by Alder and Wainwright in 1955 [1] and the resourceful experiments by Pusey and van Megen in 1987 [2], the crystal nucleation process in suspensions of colloidal hard spheres has been the topic of extensive research in the past decades. Strikingly, despite the simplicity of the system, there are pronounced differences between theoretically predicted and experimentally measured nucleation rate densities, as remarked by Auer and Frenkel in 2001 [3]. This observation has been confirmed repeatedly, and many authors emphasized their astonishment about the unexplained discrepancy of 10 orders of magnitude between theory and experiment [4][5][6][7][8].
As causes for the discrepancy on the experimental side have been suggested the size polydispersity of the colloidal particles [3], their electrostatic charge [9], hydrodynamic coupling between the colloids due to the solvent [10][11][12], sedimentation in the samples [6,8], and the fact that it is difficult to control the volume fraction [13][14][15]. However, although these effects are relevant to a certain degree, neither resolved the disagreement. On the theoretical side, the choice of sampling technique has been suggested [4], the choice of reaction coordinate [16,17] and the validity of the assumption, that the nucleation process were Markovian [18,19]. Again, neither of these could explain a discrepancy of such a staggering magnitude.
The central "bone of contention" is shown in fig. 1, which presents the crystal nucleation rate density as a function of the volume fraction of spheres. As unit of length we used the particle diameter σ and as unit of time the long-time self-diffusion time τ L = (6D S L ) −1 , where D S L is the long-time self-diffusion coefficient. This choice of units is common in the literature, because it allows for a consistent comparison of the rather varied experimental realizations of the system. The striking point about Figure 1 is the difference in slope between the experimental data and the numerical data. While the effects listed above could explain uncertainties in the location of some of the curves, none explained the significant difference in slope.
To solve this issue, we investigated nucleation and crystal growth by means of an efficient event driven molecular dynamics simulation code (EDMD), which we developed based on the algorithm proposed by Bannerman in 2014 [26]. In this letter, we first present a precise computation of the homogeneous nucleation rate density.  [3,17,[20][21][22][23][24][25], as well as this work (black symbols, see also Table I and Table II).
which hard sphere nucleation experiments are evaluated and use the Johnson-Mehl-Avrami- Kolmogorow model (JMAK model) to compare our data to the experiments. In particular, we take into account the fact that polycrystalline clusters form by a successive twinning process. We conclude that the experimental data presented in fig. 1 is not the crystal nucleation rate but the rate at which mono-crystalline domains form. Thus we argue that the discrepancy between experiment and theory is due to a misinterpretation of the experimental data.
The simulated system consists of up to a few million mono-or polydisperse hard particles enclosed by cubic periodic boundary conditions. The particles' (average) diameter σ, mass m, and the thermal energy k B T are chosen as units resulting in a unit of time of δt = m/(k B T ) σ. For nucleation measurements, the system is rapidly compressed from the stable liquid state to a meta-stable volume fraction.
We monitor the crystallinity of the system in terms of the q6q6-bond-order parameters [27,28], details are given in [19]. As the main observable we measure the number of particles within the largest crystalline cluster. The time between the volume quench and the first nucleation event (the "nucleation time") is determined by following the trajectory of the largest cluster within a crystallized system back to the point in time when its size last crossed the size of the average largest cluster of the meta-stable fluid. As the nucleation time increases drastically for low supersaturation, it is not feasible to run a set of simulations in which the complete ensemble ends in the crystalline state. Therefore, a definition of the nucleation rate, that does not depend on a fully crystallized simulation ensemble, is required.
Here, we define the nucleation rate via the assumption of exponentially distributed nucleation times. For each measured distribution, we use the maximum likelihood estimator (M.L.) of a censored exponential distribution [29].
An example of a nucleation time distribution with the corresponding maximum likelihood estimation is shown in Figure 2. The example is measured in a mono-disperse system with a volume fraction of η = 0.534. We observe, that the simple exponential model describes the distribution well, indicating that the assumption of a time-independent nucleation rate is justified.  To determine the nucleation rates, we ran 500 trajectories of mono-disperse samples  we carried out 100 simulations, each consisting of 16 384 particles with 5% Gaussian polydispersity. The analysis of the poly-disperse samples relies on the absence of competing phase transitions, but, as remarked by Bolhuis and Kofke [30] fractionation effects only occur above volume fractions of 5 %. As we have not seen any delay in the nucleation times after the quench, our simulations confirm that fractionation is not significant. For poly-disperse systems, it does not make sense to plot the rates versus the plain volume fraction, as the coexistence region is shifted with respect to the one of mono-disperse hard spheres. In principle, it would be desirable to plot all data as a function of the difference in chemical potential between the meta-stable liquid and the stable crystal. However, as it is time-consuming to compute this difference, we use instead which is close to the mapping given by the chemical potential, much simpler and frequently used in experimental analyses. The poly-disperse transition points η freeze and η melt are taken from ref. [30].
Our results, Tab. I and Tab. II, agree well with older, less accurate simulation results. We conclude that the mismatch between experiments and simulations is not due to the choice of simulation method. Thus we revisit the analysis of the experimental data.
In fig. 3 we show a slice through a crystallite obtained in our simulations. The colours encode crystalline structures, particles in the liquid are not shown. Clearly the cluster is poly-crystalline and twinned. Other snapshots (see SI) show icosahedral domains in addition.
We observed only very few mono-crystals. Poly-crystalline clusters produced by a successive twinning process have also been observed by O'Malley and Snook [32]. We will argue in the following that this observation contradicts one of the basic assumptions made in the analysis of the experimental data, and that this is the cause of the discrepancy in nucleation rates. In ref. [23,[33][34][35] the Bragg scattering data was analysed as follows: The structure factor of a crystallizing sample was measured as a function of time. Then the structure factor of the liquid phase was subtracted in order to obtain only the signal that stemmed from the crystalline components. The degree of crystallinity χ(t) was then determined based on the area under the main Bragg peak of the crystalline structure factor, and the average linear size L(t) of the crystallites was determined based on the width of the peak at half maximum (according to the Scherrer equation). Finally, the number of crystallites was computed as To our knowledge, this type analysis was first presented for hard spheres by Harland and van Megen in the 1990s, and all analyses carried out in more recent work follow their protocol. Harland and van Megen assumed that L(t) was the size of individual crystals rather than domains, although domains produce a Bragg peak of the same shape. We do not know why they excluded the possibility of poly-crystalline growth.
To check whether our simulation data is consistent with the published Bragg scattering signals, we have produced "scattering data" from simulated configurations, i.e. we computed the radial distributions functions of the crystallizing systems and Fourier transformed them.
Then we subtracted the liquid background. The resulting structure factor is shown in fig. 4 (solid line). For comparison we show the experimental structure factors from Schöpe [34]  (dashed line) and Iacopini [35] (dash-dotted line). As we do not know the crystallinity exactly, we scaled the structure factors to match the heights of the first Bragg peaks, thus only the relative intensities between the peaks are meaningful. Furthermore, we multiplied the q-data of ref. [34] by a factor of 1.01 and the one of ref. [35] by 0.97, which is compatible with the uncertainty of measured sphere radii. The resulting signals are nearly identical.
Hence we conclude that our simulations, which show poly-crystalline growth, are consistent with the Bragg scattering signals from experiments in the literature. Regarding the small angle scattering experiments (SAS), we cannot show a direct comparison, because our simulations do not allow us to access sufficiently large length scales.
Thus we base our arguments on ref. [22]: The SAS signal stems from spatial variations in the refractive index, which are caused by spatial variations in the colloid concentration.
Ackerson and Schätzel assumed that the signal was produced by the drop in concentration at the crystal-liquid interface [22]. However, they observed that the signal continued to evolve in time after they expected crystal growth to have stopped. They attributed this phenomenon to coarsening, i.e. to the dynamics of the grain boundaries, at which there is also a drop in concentration. According to Babinet's theorem the crystal-liquid interface and the amorphous slit between two crystalline domains produce the same scattering pattern.
Hence it makes sense that crystallites embedded in a liquid and grain boundaries in the crystal produce similar SAS signals. We conclude that the SAS signals are also consistent with the assumption of poly-crystalline growth.
To interpret the nucleation measurements, we use the JMAK-model [38][39][40][41], which describes a crystallization process with stochastically distributed, independent nucleation events and spherically growing clusters. Here we consider the case of a time-independent homogeneous nucleation rate density, κ hom , and a time-independent radial growth rate of the clusters, c (for a justification of the latter assumption see SI). According to the JMAK-model, the fraction of the system's volume which is filled by crystals increases as (See SI for a brief derivation of this equation.) Thus the characteristic filling time t * p at which a fraction p of the whole system has undergone crystallization is given by For the following analysis we need to define a time that can be compared to the experimental induction time. We define a characteristic time t * at which crystallization is in full progress As experiments, simulations and the JMAK-model alike show a steep increase of the crystallinity after an initially slow growth phase, the specific choice of p does not influence t * by much. Therefore, the results obtained by setting it to the computationally convenient value p = 1 − exp − π 3 do not differ much from those obtained with other values of p. Prior to events in which clusters merge, the total number of mono-crystalline domains is given by the sum of the homogeneously nucleated clusters and of the mis-oriented domains grown on the cluster surface. In Equation 5 we write this relation in terms of rates, with the total mono-crystalline domain rate density, κ mono , and the rate density for the growth of mis-oriented domains on the interface, κ domain . (We average over the rates for events on differently orientated crystalline planes and for different twinning processes.) where A crys is the interfacial area between the crystallites and the liquid and V is the system

Inserting Equation 6 into Equation 5 we obtain
To compare κ mono to the experimental data, we evaluate Equation 7 at the characteristic time t * from Equation 4.
We determined the value of κ domain by visual inspection of our simulation snapshots to be κ domain = 5 · 10 −4 (σ 2 τ L ) −1 (see SI for a detailed description of the procedure).
We suppose that the rate that has been determined in Bragg-scattering experiments is κ mono rather than the bare homogeneous nucleation rate κ hom . In Figure 5 we show the same experimental data as Figure 1, but all simulation data is now given in terms of κ mono .
The ten orders of magnitude gap is closed and there is remarkable agreement between the simulations and the experiments.
The agreement is particularly clear, if the poly-dispersity of the experimental samples is taken into account. The larger the poly-dispersity, the lower the rate, i.e. the simulation data on the top of the bundle of lines, which stems from mono-disperse systems and the experimental data on the bottom, which stems from samples with more than 5 % polydispersity, are consistent with our interpretation.
In summary, we analysed and characterized the nucleation and crystallization process of the meta-stable hard sphere system in detail. We calculated nucleation rate densities The gap between experiment and simulation is closed. seems to result mostly from interpreting the shape of scattering signals as originating from full clusters rather than mono-crystalline domains.  The crystal in Fig. 6 is clearly twinned. In Fig. 7 and Fig. 8, there is a mixture of various twinned structures. Each cluster we observed had grown from one initial nucleus.
We hardly ever observed the incorporation of an independently nucleated cluster into the growing cluster. Some domains formed by secondary nucleation on the cluster surface rather than by twinning, but this process was rare, too. Fig. 9 shows a mono-crystalline cluster.   As shown in the example snapshots, the domains can easily be distinguished by eye (if only the order of magnitude of their rate of formation is of interest). While the fcc/hcp type clusters mostly contain fewer than ten domains, the icosahedral clusters typically incorporate more than ten domains for clusters with a few hundred thousand particles. We counted Here, A cluster (t) is the surface area of the cluster and c is the radial growth rate. We assumed that κ domain is independent of time, which neglects the fact, that there might be an impact from the cluster's varying curvature on the twinning process. As we are not studying a large range of curvatures, this assumption is reasonable. Also, we assumed a constant spherical radial growth rate c, which is justified by our analysis presented below in the section regarding the crystal growth rate. The number of particles in the cluster N cluster is linked to the cluster growth time T and the radius of the cluster R(t) by: Inserting the relation given in eqn. B2 into eqn. B1 we obtain for κ domain Regarding the dependence of κ domain on the supersaturation of the liquid, we analysed the mono-disperse systems, ranging from η = 0.531 to η = 0.54. We could not find a definitive trend for κ domain . While this is partly due to our approximative method of determining κ domain , it is also reasonable to assume that κ domain does not depend sensitively on the volume fraction of the system. Most of the domains are formed by successive twinning processes. The energy cost associated with the formation of a twinned domain stems from the cost of formation of a grain boundary, and the kinetics of twinning is driven by the elastic properties of the crystal. Both, the elastic constants as well as the cost of a grain boundary, depend only very weakly on the total volume fraction of the system. The dependence is due to the pressure applied on the crystal by the liquid, thus it is of a much weaker nature than e.g. the dependence of the nucleation rate on the supersaturation. We therefore carried out our analysis with a rate which is independent of the density of the system: Appendix C: Crystal growth rate We calculated shape descriptors based on the tensor of gyration of the crystallites and found the crystallites to be mostly spherical for particle numbers above a few thousand.
Furthermore, we calculated power law fits to the particle number growth law of the largest cluster within the system. A quantitative evaluation of the power law coefficient, n, for the 5% Gaussian poly-disperse system at η = 55.4%, resulted in n = 2.97 ± 0.03, consistent with the theoretical value for constant radial growth being n = 3.
Qualitative deviations from this behaviour were not observed in any system, except for fluctuations at cluster sizes below 1000 particles. Therefore, we concluded that constant radial growth is the dominant growth mode in the system.
The cluster growth rates were measured by linear regressions to the third root of the particle number in the largest cluster. In Fig. 13 we show the data of a poly-disperse seeded system as an example.
While we are confident that our measurement is accurate, we also note, that the radial growth rate enters the calculation of the crystallization rate only as a fourth root. Thus, to shift the rate densities presented in the letter by one order of magnitude, the value of c would have to vary by a factor of 10,000.  Cluster growth rate, measured in a 5% Gaussian poly-disperse system of 500 000 particles, with an initial seed of about 7000 particles.
volume V is dN cluster dt = V κ hom .
At time t, a cluster which was formed at a time t has grown to a volume V cluster (t, t ) = 4π 3 c 3 (t − t ) 3 .
Thus in the "ideal case", in which the clusters do not overlap, the total volume filled by crystals at time t, V cryst, ideal (t), would be given by i.e. per infinitesimal time interval, the volume filled by crystals grows as However, in a real system, the volume available for nucleation and growth shrinks while the crystalline clusters grow. Therefore all processes need to be weighted with the probability of still finding non-transformed space. If the actually crystallized volume is V cryst (t), then the remaining volume fraction is given by 1 − Vcryst(t) V . Expressed in the differential form: