The Early Crystal Nucleation Process in Hard Spheres shows Synchronised Ordering and Densification

We investigate the early part of the crystal nucleation process in the hard sphere fluid using data produced by computer simulation. We find that hexagonal order manifests continuously in the overcompressed liquid, beginning approximately one diffusion time before the appearance of the first `solid-like' particle of the nucleating cluster, and that a collective influx of particles towards the nucleation site occurs simultaneously to the ordering process: the density increases leading to nucleation are generated by the same individual particle displacements as the increases in order. We rule out the presence of qualitative differences in the early nucleation process between medium and low overcompressions, and also provide evidence against any separation of translational and orientational order on the relevant lengthscales.


I. INTRODUCTION
Homogenous nucleation in suspensions of hard spherical particles is a topic of particular interest within the wider study of nucleation phenomena. This system is often chosen for study because it is arguably the simplest possible off-lattice nucleation process to treat theoretically, and is also fairly convenient experimentally. For the theoretician this system offers the benefit that temperature does not enter into the phase diagram, which is determined purely by the configurational entropy, and also that interactions are extremely short ranged thus providing a system in which a large number of particles can be simulated inexpensively. For the experimentalist, a suspension of colloidal hard spheres has the advantage that individual particles of approximately 500nm can be resolved using confocal microscopy, and that the diffusive timescales (in a suitably viscous solvent) are in the manageable range of milliseconds to minutes. Hard spheres have only two equilibrium phases: a fluid and an FCC crystal. The HCP crystal is metastable with respect to FCC by a tiny amount due to the minutely larger phonon entropy of the FCC lattice 1 . There is no distinction between liquid and gas in the fluid phase due to the lack of any attractive interaction. Length and time measurements from hard sphere systems are typically presented using dimensionless units such that r * = r/σ and t * = t/ σ m/k B T where σ is the particle diameter, m is the mass, T is the temperature and k B is Boltzmann's constant. Considering only the appearance of structural order, the nucleation process has been shown to begin with the formation of diffuse and ramified nascent nuclei which have a mixture of liquidlike and hexagonal internal order (with the hexagonal order provided by a mixture of FCC and RHCP) without measurable BCC or icosahedral structure 2,3 .
The question of the mechanism of homogenous nucleation in hard sphere suspensions remains controversial despite a considerable amount of work from both the-ory and experiment (see 4 for a recent review). Previous simulation studies of crystallization in hard spheres have shown that dense but disordered structures form before the appearance of dense and ordered crystallites 5 , and strong theoretical arguments have been made that in general a long-wavelength density fluctuation should precede the formation of an ordered cluster 6,7 . Experimental studies of near-monodisperse hard sphere systems using time-resolved light scattering have also shown evidence for a two-stage nucleation process 8 , although this was initially interpreted as due to fractionation of similar but non-identical particles delaying the formation of the fully ordered crystallite.
A position diametrically opposite to the idea that density fluctuations are prior to structural ordering has also been taken 9 , arguing from constant pressure Brownian Monte-Carlo simulations that the pre-critical nucleus is in fact considerably less dense than the equilibrium solid, and that structural order reaches a higher proportion of its final value more quickly than density does. The method of comparing changes in the density and in the structural order by renormalising them to their equilibrium values in the crystal does not address the question of which quantity has the greatest determining role in selecting the fate of a pre-critical nucleus, as a small early fluctuation in one quantity may nonetheless be a cause of nucleation while large late ones in another quantity may be purely a consequence. To resolve this ambiguity we here collect and compare simulation trajectories not only of nucleation events, but also of order parameter fluctuations not resulting in nucleation of the solid phase.

II. SIMULATION METHOD
The sphere dynamics were evolved using event driven molecular dynamics (EDMD) 10-13 at fixed particle number, volume and temperature. The particle number arXiv:1503.07732v1 [cond-mat.soft] 26 Mar 2015 was set to N = 1000 spheres, and the volume fraction was chosen as η = 0.525. At this packing fraction the chemical potential difference per particle between the metastable fluid and the stable crystalline state is The FRESHS rare event sampling system 15 was used to generate and collect multiple EDMD runs to make Forward Flux Sampling (FFS) calculations. FFS is an established simulation method for calculating reaction pathways for non-equilibrium rare events 16 , which operates by 'splitting' trajectories of a stochastic system as the state evolves with respect to some progress coordinate, enabling an importance sampling to be carried out on trajectories which move towards the rare event of interest. To enforce stochasticity on the otherwise deterministic trajectories generated by EDMD, thus enabling otherwise identical configurations to split into different trajectory paths, velocities were randomised at time intervals of 0.1 σ m/k B T giving Brownian dynamics over reduced times above 0.1. The progress coordinate X was defined following ten Wolde et al. 17 as the size of the largest cluster of connected particles having solid-like bond order 18 , with interfaces placed at integer multiples of 2 particles. Statistics were obtained by averaging over 120 independent FFS calculations. The overcompressed fluid configurations did not show pre-existing crystallites that might have been created during the preparation process.

III. RESULTS
Averaging over all FFS runs, a nucleation rate of 3.8(±2.5) × 10 −12 was found, lying within the spread of results from previous calculations at this overcompression 4 . Trajectories from the FFS runs were then classified according to those which returned to the quiescent state at some point after forming an initial cluster of two or more ordered particles and those which progressed to an irreversibly large cluster size. The timezero for each trajectory was set as the first formation of a cluster of two particles. Averages were prepared by taking a histogram over time of the appropriate quantity. For successful trajectories, averages were weighted by the nucleation rate resulting from that individual FFS calculation. For unsuccessful trajectories the weight was set by the partial rate to the highest interface on X to be reached.
Density fluctuation magnitudesρ k at integer wavenumber k were defined as the root of the sum of squared magnitudes for the Fourier components f k of the particle distribution at wavelength L/k: wavelengths L, L/2, L/3 and L/4 are shown where L = 9.991σ is the length of the cubic simulation box.
Examining the time series of the different fluctuations ( fig. 1), the differences in density fluctuations at long wavelengths between systems which are destined to nucleate and those which are not can be observed long before the first structural fluctuation. The time course of these slow density fluctuations follows three stages, with only the third stage corresponding to the appearance and increase of structural order.
The first stage of density fluctuation in successfully nucleating trajectories begins earlier than the earliest times for which trajectory information was collected (t * = −20). At the long wavelengths k = 1, 2, 3 an increased density fluctuation is visible indicating an imbalance of mass towards one half of the simulation box and thus a local enhancement of the chemical potential driving nucleation. At the shortest wavelength examined (L/4 = 2.498σ, comparable to the diameter between opposite first peaks in the fluid radial distribution function 19 ) greater smoothness is instead observed during this stage, indicating a net decrease of coherent liquid-like order coincident to the bulk density imbalance. In the second stage (t * = −8 to t * = 0) the long-range density fluctuations at k = 1, 2, 3 disappear, being replaced by a large-amplitude, short-wavelength density fluctuation of enhanced liquid-like order within which the hexagonally ordered nucleus will form. In the third stage, from t * = 0 onwards, crystalline order appears and the density fluctuations follow the growth and compaction of the crystal phase.

IV. DISCUSSION
We have shown here that long wavelength density fluctuations prefigure the development of structural order in the hard sphere system at a low overcompression, by a large time interval of over 20 σ m/k B T . This phenomenon has previously been repudiated by other researchers, however we can offer as an explanation for this that the onset of structural ordering is accompanied by a relaxation of these density fluctuations, allowing them to be overlooked if only the time immediately preceding the formation of a critical cluster is examined.
In the dilute limit the expected lifetime of diffusively generated density fluctuations is given by Fick's laws of diffusion as being proportional to 1 D L k 2 , where D is the particle self-diffusion coefficient and L k is the wavelength. Away from the dilute limit the correct value of D becomes wavelength and concentration dependent due to collective effects, however the expectation remains that density fluctuations should be increasingly slow as wavelength increases. This idea is supported by the long duration of the density fluctuations observed to drive nucleation events.
Although not all splitting-type sampling methods for the study of rare event processes necessarily require a Figure 1. Cluster size (top) and density fluctuation amplitudes (lower four) for trajectories which nucleated (red) and those which reached the committor then returned to a cluster size of 2 or smaller (black). A variation in the size of density fluctuations between those trajectories which will nucleate and those which will not is manifest at very early times. The zero of the time axis is set at the first appearance of a cluster of size 2.
Markovian dynamics on the progress coordinate 20 , all such approaches are more effective with the benefit of a clear understanding of the reaction pathway, which we suggest a consideration of the density fluctuation dynamics might provide. Given that the progress of the density fluctuation dynamics prior to nucleation is nonmonotonic and wavelength dependent, as well as potentially varying with the overcompression, a comprehensive picture of this aspect of nucleation phenomena will require considerable further study.