Arrhenius temperature dependence of the crystallization time of deeply supercooled liquids

Usually, supercooled liquids and glasses are thermodynamically unstable against crystallization. Classical nucleation theory (CNT) has been used to describe the crystallization dynamics of supercooled liquids. However, recent studies on overcompressed hard spheres show that their crystallization dynamics are intermittent and mediated by avalanche-like rearrangements of particles, which largely differ from the CNT. These observations suggest that the crystallization times of deeply supercooled liquids or glasses cannot be described by the CNT, but this point has not yet been studied in detail. In this paper, we use molecular dynamics simulations to study the crystallization dynamics of soft spheres just after an instantaneous quench. We show that although the equilibrium relaxation time increases in a super-Arrhenius manner with decreasing temperature, the crystallization time shows an Arrhenius temperature dependence at very low temperatures. This is contrary to the conventional formula based on the CNT. Furthermore, the estimated energy barrier for the crystallization is surprisingly small compared to that for the equilibrium dynamics. By comparing the crystallization and aging dynamics quantitatively, we show that a coupling between aging and crystallization is the key for understanding the rapid crystallization of deeply supercooled liquids or glasses.


I. INTRODUCTION
Classical nucleation theory (CNT) effectively describes first-order phase transition dynamics.This theory predicts that nucleation first takes place, where sufficiently large droplets of the ordered phase are created in the sea of the disordered phase and then the droplets grow steadily [1][2][3].The nucleation is driven by thermal fluctuations and the nucleation time τ N , the inverse of the nucleation rate, is given by where k B is Boltzmann's constant and T is the temperature.The constant τ t is the characteristic time scale of the transport process and ∆G is the free energy barrier to form the critical droplet, which is given by the maximum value of the droplet formation energy as a function of the droplet size.In its simplest version, the droplet formation energy is given by the combination of the free energy gain due to the low energy ordered phase and the surface tension between two phases.The CNT has been successfully applied to a wide variety of systems undergoing first-order phase transitions.In the Ising model, for example, the direct numerical estimate of τ N was shown to accurately be described by Eq. ( 1) [4][5][6].
For the crystallization of slightly supercooled liquids, the formula Eq. ( 1) is known to work well [7,8], although the first-principle theoretical prediction of τ t and ∆G is difficult and still under debate [9,10].As liquids are further supercooled, the glass transition is approached, * yukitakaha@g.ecc.u-tokyo.ac.jp which results in a dramatic slowing of the particle transport process.To effectively take this into account in the CNT, τ t is conventionally approximated by the equilibrium relaxation time τ α [11] [12], leading to Notably, this formula predicts that τ N is a nonmonotonic function of temperature.The nucleation time τ N diverges as the temperature increases since the first-order transition is approached and the free energy gain diminishes so that ∆G diverges.On the other hand, τ N drastically increases as the temperature decreases since the glass transition is approached and τ α drastically increases.As a result, τ N has a minimum value at a certain temperature.Reflecting this, the CNT predicts a similar nonmonotonic temperature dependence for the crystallization time τ cry , which is the time required for crystalline regions to occupy a major part of the sample.Indeed, several experimental [13][14][15] and numerical [16,17] studies confirmed that the nucleation and crystallization time has a minimum value at a certain temperature.Despite these successes of the CNT, it has recently been reported that the crystallization of deeply supercooled liquids and glasses is qualitatively different from that in the CNT.In highly overcompressed hard spheres, crystallization proceeds before the transport of particles becomes diffusive [18].The crystallization is caused by small displacements of particles, which are even smaller than the particles' diameter [19].Moreover, avalanchelike crystallization has been reported for well annealed "mature" glasses in hard spheres [20] and pseudo hard spheres [21,22], in which intermittent stochastic rearrangements of particles cause crystallization.Some similarities between avalanche-like crystallization dynamics and aging dynamics have been reported [23].An experiment of a colloidal system also shows observations consistent with these results [24].These phenomena do not appear in the CNT, which assumes equilibrium conditions and a continuum picture.
These peculiar crystallization dynamics are interesting, but quantitative analysis of the crystallization time is still lacking.In particular, the temperature dependence of the crystallization time has not yet been studied in detail.This is partly because the previous studies focused mainly on hard spheres, where the temperature does not play any role.In this paper, we perform molecular dynamics (MD) simulations of soft spheres.We mainly focus on the dynamics of crystallization after an instantaneous quench from a high temperature to low temperatures.We show that, while the equilibrium relaxation time depends on the temperature in a super-Arrhenius manner, the crystallization time shows an Arrhenius-like temperature dependence at very low temperatures.Moreover, the corresponding energy barrier is much smaller than that for the equilibrium dynamics.This clearly contradicts Eq. ( 2).We show that the crystallization process at low temperatures consists of an early process, in which the system falls into the nearest inherent structure, and a late process, in which transitions between inherent structures occur.The crystallization time is dominated by the late process, which is widely different from the continuum picture in the CNT.Furthermore, we show that although this late process has some similarities with the aging dynamics, the crystallization time cannot be directly explained by the characteristic time scales of the aging dynamics.Instead, a coupling between the aging and crystallization dynamics is the key for understanding the rapid crystallization.

II. METHODS
We perform MD simulations of monodisperse and polydisperse soft spheres.The interaction potential is an inverse-power-law (IPL) potential given by where ε is an energy scale; r ij = |r i − r j |, with r i being the position of particle i; and σ ij ≡ (σ i + σ j ) /2, where σ i represents the effective diameter of particle i.
The coefficients A and B force the continuity of the first and second derivatives at the cutoff length r c,ij ≡ 1.5σ ij .
For polydispersity, we introduce a top-hat distribution of particle size with the average particle size σ and the width ∆, in which a non-trivial particle arrest is known not to occur [25].We consider the systems with ∆ = 0.00, 0.12, 0.24, and 0.36 in this work, where ∆ = 0.00 corresponds to the monodisperse system.To approximate a bulk system, we use the periodic boundary condition in a cubic box (volume V ).Mass, length, time, and temperature are measured in units of m, σ, σ m/ε and ε/k B , respectively, where m is the particle mass and k B is the Boltzmann constant.Throughout this paper, the temperature is controlled by the Nosé-Hoover thermostat, where the thermostat mass is tuned so that unphysical temperature oscillation decays sufficiently rapidly even after the instantaneous quench.The time step of the MD simulation is 0.01.The number of particles N is 16384 unless otherwise noted.The packing fraction φ = i πσ 3 i /6) V is fixed to be π/6 for both monodisperse and polydisperse systems.At this density, the freezing temperature is T Freezing = 0.59 for the monodisperse system [26].
To measure the degree of crystallization of the samples, we use the method introduced in Ref. [18,27].First, each particle is assigned a vectorial bond-order parameter d 6 [28].Pairs of neighbors are then identified by using a cutoff distance 1.4σ ij , and a bond of a pair is deemed "solid-like" if the scalar product of their d 6 vectors exceeds 0.7.A particle is labeled solid-like if it has at least six solid connections.Finally, we calculate the crystallinity X(t) of a sample as the proportion of solid-like particles at a given time t.

A. Equilibrium dynamics
To discuss the crystallization dynamics on firm ground, we first clarify the equilibrium dynamics of the model.We equilibrate the system at a target temperature T and polydispersity ∆ and then perform the production runs starting from the equilibrium configurations.We perform 20 independent runs at each state point.The overlap function is then calculated from the time series of the particle configurations, where the bracket represents the average over 20 samples.This calculation is performed with N = 2048.
Figure 1(a) shows the overlap functions for ∆ = 0.36.Clearly, the overlap function decays rapidly at higher T and the relaxation becomes much slower with decreasing T .It shows a two-step relaxation, which is a hallmark of the dynamics of supercooled liquid.We measure the equilibrium relaxation time τ α defined as F O (t = τ α ) = 0.4 for various state points (T , ∆).For some (T, ∆) conditions, we find that the system undergoes crystallization during the equilibration or production runs.Since we focus solely on the equilibrium dynamics here, we do not calculate τ α if X(t = τ α ) ≥ 0.5.We note that the equilibration time is always 40 times larger than the relaxation time τ α .dence, although the relaxation time cannot be calculated with small ∆ and low T because of crystallization.Importantly, the relaxation times for various ∆ collapse onto almost the same curve, meaning that τ α does not strongly depend on ∆.In previous studies, it has been reported that the fragility depends on the polydispersity; the more polydisperse the system, the stronger the dynamics [29].We can indeed find such a tendency, albeit slightly, when we closely compare ∆ = 0.24 and 0.36.However, this effect is minor, as we focus on the small ∆ region, and we can regard the relaxation time as almost independent from ∆.

B. Crystallization time
We now focus on the crystallization dynamics of the monodisperse system.We first prepare equilibrium liquid configurations at T = 3 and then perform MD simulations at the target temperatures starting from these configurations.This corresponds to an instantaneous quench from T = 3 to the target temperatures.We show the re-sults with N = 16384, which is free from a finite size effect [30]. Figure 2(a) shows the time evolution of the crystallinity X(t) of typical runs.At T = 0.34, the crystallinity X(t) fluctuates around X = 0 for a long time, and then it suddenly starts to increase at very large t.This time evolution is qualitatively consistent with the dynamics predicted by the CNT, which consist of rare nucleation and steady growth of the nucleus.On the other hand, X(t) for T = 0.20, 0.06, 0.02 is very different.The crystallinity X(t) shows an increase up to 10% in the short-time region t 10 2 , which we call the early process of the crystallization.Interestingly, X(t) for various temperatures coincides in this time region, meaning that the early process does not strongly depend on the temperature if the temperature is sufficiently low.At larger time t 10 2 , gradual and slow increases are observed, which we call the late process of the crystallization.Unlike the early process, the late process becomes slower and slower with decreasing temperature.At the lowest temperature, X(t) fluctuates around 10% and increases intermittently.Then, it finally grows at very large t.These results mean that the intermittent crystallization takes place not only in overcompressed hard spheres [20] but also in soft spheres at low temperature.These crystallization dynamics clearly differ from that described by the CNT.
To discuss the statistical law that governs the crystallization, we calculate the survival probability [8].We repeat the simulations described above for 100 independent samples and calculate the crystallinity X(t) for each of them.We label the samples as "crystallized" at a given time t if the crystallinity X(t) becomes larger than 0.5 [31].Then, we calculate the survival probability P (t), which is defined as the proportion of noncrystallized samples at time t.Figures 2(b) and 2(c) show the survival probability P (t) at each temperature.Figure 2(c) shows that P (t) does not decrease from 1 in a certain time, called a lag time, and then decays exponentially at larger t.This means that crystallization of a given sample can be seen as a random event with a constant rate, although there is a minimum latency time.Note that similar behavior of P (t) was previously reported for the liquidgas condensation dynamics in the Lennard-Jones system, where the CNT works well [8].In this case, the lag time was identified as the time required for the growth of the nucleus, and the exponential decay time was identified as the nucleation time.
We now extract the characteristic time scales of the crystallization.The simplest definition of the crystallization time is the time at which the survival probability reaches some threshold value.In practice, we define τ cry as P (t = τ cry ) = 0.4, the time at which 60 % of samples become crystallized.Additionally, we can define the lag time τ lag and the exponential decay time τ exp by fitting P (t) to exp((t − τ lag )/τ exp ) [32].These three time scales are plotted in Fig. 3(a).At approximately the freezing temperature, the exponential decay time τ exp is much larger than the lag time τ lag and then τ cry ≈ τ exp .This means that rare nucleation is the rate-controlling process and growth is much quicker.With decreasing temperature, τ exp drops more sharply than τ lag and then τ lag dominates τ cry .In this temperature range, nucleation becomes fast and comparable to the microscopic time scale, which may be called "spinodal nucleation" [27].Finally, at the lowest temperatures, the exponential decay time τ exp increases rapidly, while the lag time τ lag does not, and as a result, τ cry is again dominated by τ exp .Therefore, our results establish that the exponential decay time τ exp is the unique time scale of crystallization in the lowtemperature regime where the intermittent crystallization takes place, and this time scale can be well captured by τ cry .The irrelevance of the lag time suggests that any deterministic process, such as the crystal growth, is absent or at least is not a rate-controlling process in this temperature regime.
Finally, we quantitatively discuss the temperature dependence of the crystallization time.For comparison, the equilibrium relaxation time of polydisperse system ∆ = 0.36 is also plotted in Fig. 3(b).From the viewpoint of Eq. ( 2), τ cry of the monodisperse system should be compared with τ α of the monodisperse system, but this cannot be done because τ α cannot be measured for the monodisperse system at low temperatures due to rapid crystallization.However as shown in Sec.III, τ α does not strongly depend on the polydispersity in our model, so τ α of the monodisperse system can be well approximated by those of polydisperse systems [18].The crystallization time τ cry diverges not only at the freez-ing point T freezing = 0.59 but also at lower temperatures and thus has a minimum at a finite temperature (see Fig. 3(b)).This is qualitatively consistent with the conventional crystallization theory.However, at sufficiently low temperature, the crystallization time becomes much smaller than the equilibrium relaxation time, τ α τ cry .Importantly, the temperature dependence of the crystallization time is Arrhenius-like, while that of the equilibrium relaxation time increases in a super-Arrhenius manner.This is contrary to Eq. ( 2).By fitting the crystallization time and the equilibrium relaxation time with τ ∼ exp(∆E/T ), we estimate the characteristic energy barrier ∆E for these processes.The energy barrier ∆E for the crystallization is 0.16 while ∆E for the equilibrium relaxation of the high-temperature and lowtemperature liquids are 1.0 and 6.5, respectively (see Fig. 3(b)).The energy barrier ∆E for the crystallization is surprisingly small, which is even smaller than that for the equilibrium dynamics of high-temperature liquids.Therefore, we conclude that crystallization in deeply supercooled liquids proceeds by crossing energy barriers that are much smaller than those experienced by the equilibrium liquids.

C. Early process and late process
The early and late processes, which we have divided the crystallization dynamics into above, can be qualitatively distinguished from the viewpoint of the potential energy landscape.We apply the optimization method FIRE [33] to the equilibrium configurations at T = 3 and obtain their inherent structures (ISs).For 20 independent samples, we measure the crystallinity X of these ISs and define their average as X IS .While the crystallinity X of the original equilibrium liquid configurations is almost zero, X IS is approximately 8%.This value is almost equal to the crystallinity reached by the early process (see Fig. 2(a)).This suggests that the early process corresponds to the dynamics in which the system falls into the nearest IS.This is consistent with the observation that the dynamics in the early process hardly depend on the temperature.
The observation above also suggests that the late process of crystallization consists of transitions between the ISs.To directly access this dynamics, we perform simulations that mimic the late process; namely, we perform MD simulations at the target temperatures starting from the ISs.Here, the initial momenta are generated following the Maxwellian distribution at the target temperatures.We perform 20 independent runs, monitor the time evolution of X(t) and analyze the survival probabilities as described in the previous section.By this, we obtain the crystallization time for this dynamics, which we call τ cry,late .As shown in Fig. 3(b), this crystallization time τ cry,late is equal to the original crystallization time τ cry .Therefore, we conclude that the early process (falling into ISs) is not relevant for the crystallization time and that the late process (transitions between ISs) dominates the crystallization time in deeply supercooled liquids.

D. Aging and crystallization
The inequality τ α τ cry implies that the deeply supercooled liquid crystallizes before equilibration.This naturally explains the reason why Eq. ( 2) does not work because it assumes that the transport of particles is controlled by the equilibrium dynamics.This observation suggests that the crystallization time would be better compared with the relaxation time of the system just after an instantaneous quench, which is the aging dynamics.To discuss this possibility, we perform the following analysis.For monodisperse ∆ = 0.00 and polydisperse ∆ = 0.36 systems, we first prepare the equilibrium configurations at T = 3, apply FIRE to obtain the corresponding ISs, and perform MD simulations at target temperatures starting from these ISs as described in the previous section.We then measure the overlap function F O (t) to monitor the relaxation dynamics.This protocol is similar to the one used in the studies of the aging dynamics for zero waiting time, although crystallization is expected to take place during the simulation runs in the monodisperse case.Note that we perform simulations for 20 independent samples and the results are averaged over them.
Figure 4(a) shows the overlap functions F O (t) of the monodisperse and polydisperse systems with T = 0.08.We first focus on the polydisperse case where crystallization does not occur.The correlation function has two properties that are not observed for the equilibrium dynamics.(1) The short-time relaxation (t 10 2 ).The overlap function F O (t) starts to decrease at approximately t = 10.This short-time relaxation is not due to the vibrations of particles but the rearrangements of particles; we checked that the amplitude of the vibrations is too small to affect the overlap function since the target temperature is very low.To extract the time scale τ ini,poly characterizing these earliest rearrangements, we fit F O (t) with exp (−t/τ ini,poly ) in the short-time region.It is clear that the fitting works well in the region F O (t) 0.8, namely for the rearrangement of approximately the first 20% of particles.(2) The long-time relaxation (t 10 2 ).The overlap function F O (t) becomes widely separated from the exponential function and has a tail in the longtime regime.This means that the relaxation of the majority of particles is qualitatively different from that of the first 20% of particles and the former cannot be seen as the repetition of the latter.Note that a similar longtime tail in the aging dynamics was reported in Ref. [34].We introduce the relaxation time in this time regime as F O (t = τ aging,poly ) = 0.4.
The two relaxation times obtained for the polydisperse system, τ ini,poly and τ aging,poly , are plotted in Fig. 4 laxation time in the long-time regime τ aging,poly shows an Arrhenius-like temperature dependence at low temperatures.However, they are still larger than the crystallization time.An Arrhenius fit of this relaxation time yields the energy barrier ∆E = 0.8, which is far from that for crystallization.On the other hand, the relaxation time in the short-time regime τ ini,poly is much smaller than the crystallization time.We find that its temperature dependence is almost linear in 1/T rather than Arrheniuslike, suggesting that the corresponding energy barrier is vanishingly small.Therefore, we conclude that the relaxation times in the aging dynamics do not give a direct explanation of the crystallization time; namely, one cannot improve Eq. ( 2) by replacing τ α with τ ini,poly or τ aging,poly .Finally, we focus on the monodisperse case.The overlap function is shown in Fig. 4(a).We define the relaxation time in the short-time τ ini,mono and long-time regimes τ aging,mono in the same way as the polydisperse case and plot them in Fig. 4(b) [35].For the short-time decay, the overlap function of the monodisperse system is almost the same as that in the polydisperse system.Also, the relaxation times τ ini,mono and τ ini,poly coincide at all the temperatures studied.Therefore, crystallization does not affect the first rearrangements of particles in glasses.However, this relaxation time is widely separated from the crystallization time.On the other hand, in the longtime regime, the overlap function of the monodisperse system decays much faster than that in the polydisperse system.This means that the crystallization itself accelerates the relaxation dynamics.Moreover, the relaxation time τ aging,mono has a similar temperature dependence as the crystallization time.Indeed, a fitting of τ aging,mono with exp(∆E/T ) yields ∆E = 0.16, which coincides with that for τ cry .Therefore, the crystallization itself accelerates the relaxation dynamics, and the resulting relaxation time gives the crystallization time.This suggests that a coupling between the relaxation and crystallization dynamics is the key for understanding the rapid crystallization of glasses.We also find that the crystallinity X(t) begins to increase on this time scale τ aging,mono , which is consistent with the discussion above.
We additionally mention that the equivalence between τ aging,mono and τ cry is not trivial because the former is measured from the overlap function, which does not monitor the crystallinity or other orders.It monitors only the small displacement (up to 0.3) of particles.This observation is consistent with Ref. [18], where the crystallization due to small displacements was reported.

IV. CONCLUSION
In this paper, we used MD simulations to study the dynamics of monodisperse and polydisperse soft spheres, with a particular focus on the crystallization dynamics just after the instantaneous quenching from a high temperature to low temperatures.The equilibrium relaxation time increases in a super-Arrhenius manner with decreasing temperature.Despite this, the crystallization time shows an Arrhenius temperature dependence at very low temperature, and as a result, the crystallization becomes much faster than the equilibrium relaxation.This is contrary to the conventional expression of the crystallization time Eq. ( 2) based on the CNT.Furthermore, the estimated energy barrier for the crystallization is surprisingly small compared to those for the equilibrium dynamics.
To discuss this result, we first performed an energy topographic analysis.This analysis reveals that the crystallization process can be divided into an early process, in which the system falls into the nearest IS, and a late process, in which transitions between ISs occur.We show that the temperature dependence of the crystallization time is dominated by the late process.Second, we analyze the aging dynamics with and without the crystallization process.The aging dynamics consist of short-time relaxation, in which only the fastest particles undergo rearrangements and the correlation function decays ex-ponentially, and the long-time relaxation, in which the majority of particles undergo rearrangements and the correlation function shows a long-time tail.In the case of aging without crystallization, the characteristic time scales of these dynamics are far from the crystallization time.Instead, in aging with crystallization, the relaxation is accelerated by crystallization, and the resulting relaxation time is quantitatively equivalent to the crystallization time.This means that a coupling between aging and crystallization is the key for understanding the rapid crystallization in deeply supercooled liquids or glasses.
Previously, the crystallization of hard spheres at high density was studied, and various peculiar properties of crystallization dynamics have been reported, including the intermittent and avalanche-like crystallization dynamics and similarities between crystallization and aging dynamics [18][19][20][21][22][23]27].Our work is partly consistent with these reports.Intermittent crystallization takes place not only in hard spheres but also soft spheres, suggesting that these peculiar crystallization dynamics are universal in deeply supercooled liquids or glasses.Moreover, our work reveals several new aspects of this crystallization.Most importantly, we could discuss the role of temperature in the crystallization dynamics and find the Arrhenius temperature dependence of the crystallization time with a very small energy barrier.Our findings would be relevant to understanding the crystallization of atomic and molecular glasses, where the temperature is the key control parameter.Since the model studied in this paper is a glass made by instantaneous quench from high temperature, it can be seen as a simplified model of molecular glasses made by vapor deposition at very low temperature, where the cooling rate is extremely large [36].Our results suggest that the crystallization in such systems proceeds very rapidly with an Arrhenius temperature dependence even if they are fragile glass formers.
In this paper, we focus solely on the instantaneous quench of the monodisperse system.It is interesting to extend our work to different preparation protocols and various types of glass formers.Such studies would enable us to understand the role of the maturity of glasses on the crystallization dynamics and the glass-forming ability at very low temperatures.It is also interesting to analyze the intermittency of the crystallization dynamics.In the case of the yielding of sheared glasses, the intermittency and avalanche formation have been studied in detail [37][38][39], but a similar analysis seems to be lacking for systems driven by small thermal agitation.Such studies would provide more microscopic understandings of the crystallization and aging at very low temperatures.

FIG. 2 .
FIG. 2. (a) Time series of crystallinity X(t).The blue shadow represents the average XIS of crystallinity of the inherent structures of the equilibrium liquid configurations with T = 3.00.(b)(c) Survival probability P (t), which is the probability of uncrystallized samples.

FIG. 3 .
FIG. 3. (a) Temperature dependence of τcry, τexp, and τ lag .(b) Temperature dependence of the crystallization time and equilibrium relaxation time.The inverse temperature β ≡ 1/T is used on the horizontal axis.The time scale τ cry,late represents the crystallization time when the initial configurations were the ISs of the liquids with T = 3.0.The three black lines represent the slope of exp(∆E/T ).

16 FIG. 4 .
FIG. 4. (a) Overlap functions of nonequilibrium dynamics with T = 0.080.The black line shows the fitting curve of FO(t) of a monodisperse system with exp(−t/τini,mono).(b) Relaxation times of nonequilibrium dynamics and crystallization time.The inverse temperature β ≡ 1/T is used on the horizontal axis.The two black lines represent the slope of exp(∆E/T ).