Solid-like mean-square displacement in glass-forming liquids

From experimental data for nine glass-forming liquids, Gainaru and co-workers recently demonstrated a striking universality of the real part of the frequency-dependent fluidity.¹ The fluidity is defined as 1/η(ω), in which ω is the angular frequency and η(ω) is the complex frequency-dependent linear shear viscosity. The data involved van der Waals, ionic, and hydrogen-bonding liquids, i.e., chemically quite diverse systems. The universal fluidity data were shown to fit well to a prediction based on the random barrier model (RBM). This is surprising in view of the fact that this model has no shape parameters and was devised for describing the completely different hopping conduction in disordered solids. The RBM considers non-interacting particles jumping stochastically on a simple cubic lattice with identical site energies and random energy barriers for nearest-neighbor jumps.^2–4 This is not at all how one thinks about a liquid. To illuminate this puzzling situation, we have carried out extensive computer simulations of a highly viscous model liquid in order to investigate whether the RBM does describe the particle dynamics.

The relaxation time increases dramatically when a liquid is supercooled and approaches the glass transition.^5–9 A standard probe of the dynamics is the single-particle mean-square displacement (MSD) as a function of time, ⟨Δr²(t)⟩, in which Δr(t) is the distance traveled by a given atom or molecule in time t and the angular brackets denote an ensemble average. At long times, the MSD is proportional to time and determines the (self) diffusion coefficient D via ⟨Δr²(t)⟩ = 6Dt. The transition to a linear-time MSD takes place roughly at the time at which the particles on average have moved an interatomic distance. Since all liquids become diffusive at long times, it is the subdiffusive regime that reveals details about the liquid dynamics.

The RBM was devised as an idealized model of ac ionic or electronic hopping conduction in disordered solids such as oxide glasses, polymers, and amorphous semiconductors.^2,3 In the extreme-disorder limit, i.e., when k_BT is much smaller than the relevant energy barriers, the model predicts a universal MSD such that the MSD as a function of time is the same for all barrier distributions except for a scaling of time and space.⁴ The only requirement for this to apply is that the energy barrier probability distribution is continuous at the percolation threshold.⁴ Physically, the response is universal because the dynamics at extreme disorder is dominated by percolation in the 3d random energy landscape.⁴ 

In a simple analytical approximation, the frequency-dependent diffusion coefficient of the RBM defined by $D(ω)≡(−ω2/6)∫0∞⟨Δr2(t)⟩exp(iωt)dt$ is predicted to be the solution of $lnD̃=(iω̃/D̃)2/3$⁠, in which $ω̃$ is a scaled frequency and $D̃≡D(ω)/D(0)$⁠.¹⁰ This is derived by combining the Alexander–Orbach conjecture that the percolation cluster has harmonic dimension 4/3 (independent of dimension)¹¹ with the effective-medium approximation applied to diffusion on the percolation cluster.¹⁰ The quoted equation provides an excellent fit to computer simulations of the RBM,¹⁰ except at the lowest frequencies where the transition to a frequency-independent diffusion constant is better described by the solution of the following generalized equation: $lnD̃=(iω̃/D̃)(1+(8/3)iω̃/D̃)−1/3$⁠.¹⁰ The  Appendix provides a numerical approximation to the RBM mean-square displacement as a function of time.

MSD simulation data are conveniently fitted to the von Schweidler empirical expression¹² 

According to mode-coupling theory, the exponent b is non-universal.¹³ Tokuyama has discussed common features of the MSD of different models,¹⁴ but to the best of our knowledge, the possibility of a universal viscous-liquid MSD has not been considered in the literature. This means that after the publication of Ref. 1, the glass community finds itself in the unusual situation that experiments suggest a more universal behavior than reported in simulations. An important difference between experiments and simulations, however, is that the latter cannot yet cover the long time scales of experiments on highly viscous liquids. Is this why the MSDs reported in simulations, though similar, are not universal? To address this question, one needs a viscous model liquid that is fast and easy to simulate and which does not crystallize, even in extremely long simulations.

Recent exciting developments with swap dynamics have made it possible to generate equilibrium states of liquids with astronomically long relaxation times.^15,16 Unfortunately, probing the alpha relaxation dynamics on these time scales remains out of reach, so for studying the equilibrium dynamics, brute-force molecular dynamics (MD) is still the only available option. We utilize state-of-the-art graphics-processing unit (GPU) simulations¹⁷ to access equilibrium dynamics at very low temperatures. The duration of the longest simulation was four months, which with traditional central-processing unit (CPU) computing would have taken several years.

For almost a century, the standard model in liquid-state theory has been the Lennard-Jones (LJ) system based on the following pair potential $v(r)=4εr/σ−12−r/σ−6$⁠, in which ε is a characteristic energy and σ is a characteristic length.^18,19 The LJ liquid cannot be studied in the supercooled phase because it crystallizes. In 1995, Kob and Andersen proposed a binary LJ system that is easily supercooled. The Kob–Andersen (KA) model is a mixture of 80% large A particles and 20% small B particles.¹² The trick is to have a strong AB non-ideal (non-Lorentz–Berthelot) attraction impeding phase separation. The parameters of the KA model are¹² σ_BB/σ_AA = 0.88, σ_AB/σ_AA = 0.8, ε_AB/ε_AA = 1.5, and ε_BB/ε_AA = 0.5. The KA model quickly became the standard model for simulations of viscous liquid dynamics.²⁰ The mode-coupling temperature (the temperature at which idealized mode-coupling theory based on higher-temperature data predicts a diverging relaxation time¹⁶) was estimated to be T_c = 0.435.¹² As computers became faster, it eventually became possible to investigate the model below T_c (see, e.g., Refs. 20 and 21). The KA model crystallizes in very lengthy simulations;^22,23 in fact, at the standard density 1.2, the KA liquid is supercooled whenever T < T_m = 1.03.²⁴ 

Although the strong AB attraction impedes phase separation, the supercooled KA system eventually does crystallize by phase separating into a pure A phase.²⁴ Is it possible to modify the KA model to make it even less prone to crystallization? We do this by introducing a shifted-force cutoff at r = 1.5 σ_AA for the AA and BB interactions.^25,26 Figure 1(a) shows the original KA pair potentials (dashed lines) compared to the modified binary Lennard-Jones (mBLJ) pair potentials (full lines). The AA attraction of the later is visibly weaker, and the BB attraction has also been weakened. The motivation for using a shifted-force cutoff is that this is known to leave the liquid dynamics almost unchanged,^25,26 thus facilitating a comparison between the original and the modified model.

We performed molecular-dynamics simulations in the NVT ensemble with N = 8000 particles (unless otherwise noted) at the “standard” number density ρ ≡ N/V = 1.2; the temperature T was controlled by using a Nose–Hoover thermostat. Unless otherwise noted, the results are reported in standard MD units defined by σ_AA = 1, ϵ_AA = 1, m_A = m_B = 1, and k_B = 1. The time step was 0.005 in this unit system.

The mBLJ liquid did not crystallize during the months of GPU simulations carried out for this study. Figure 1(b) shows, as a function of temperature, the average crystallization time for the original KA model (red circles)²³ and the total simulation times for the mBLJ system (black crosses). Since the mBLJ system did not crystallize, at each temperature, the total simulation time gives a lower bound on the crystallization time. The black rectangles indicate where the unknown crystallization times are to be found. At T = 0.40, several independent simulations were performed. This includes simulations that were first equilibrated at the lower temperatures (0.37, 0.38, and 0.39, respectively), a procedure known to increase the tendency to crystallize. Based on the data presented, we estimate that the mBLJ liquid has at least a 100 times longer crystallization time than the original KA liquid.

Having modified the original KA model such that crystallization is, in practice, avoided, we turn to studying the supercooled liquid dynamics. Figure 2 shows the mBLJ liquid’s all-particle MSD as a function of time at four temperatures. The figure presents data going to times larger than 10⁸ MD time units, corresponding to 0.2 ms in argon units. The data are for N = 8000 particles; size independence was checked by simulating also N = 27 000 particles at T = 0.40, which gave indistinguishable results. The blue dashed line shows data for the original KA model at T = 0.40, which are close to those of the mBLJ model (blue crosses). This confirms that the two models have very similar dynamics.

At very short times, one finds the well-known free-particle ballistic MSD ∝ t², after which there is a plateau where the MSD is almost constant. This derives from “cage rattling” of the particles in local potential-energy minima, reflecting the fact that a viscous liquid over short time scales is virtually indistinguishable from an amorphous solid. At longer times, the MSD increases, of course, and eventually one finds the standard diffusive behavior with the MSD proportional to time. Note the dramatic slowing down upon cooling: a temperature decrease of 7.5% results in more than one decade’s slowing down. What causes this is, in a nutshell, the mystery of the glass transition.^8,9 The full curves in Fig. 2 are fits to the von Schweidler expression [Eq. (1)].

Figure 3 investigates time–temperature superposition (TTS), i.e., whether data for different temperatures can be made to collapse by scaling the axes, as found in many experiments.¹ In their original paper, Kob and Andersen reported TTS obtained by scaling time with the diffusion coefficient. In Fig. 3(a), we apply the same scaling to our low-temperature MSD data, finding differences in the plateau regime. In Fig. 3(b), we perform a further scaling on both axes by a parameter c that has the dimension of a squared length, thus making both axes dimensionless. The determination of c is described below in connection with Fig. 4. From Fig. 3(b), one concludes that TTS applies.

Next, we compare our low-temperature MSD data to the RBM prediction (Fig. 4). Full lines are the MSD data scaled as in Fig. 3, including now also the ballistic regime. The RBM relates to spatially discrete particle jumps, so the predicted MSD does not reproduce the short-time “cage-rattling” plateau present in any realistic viscous liquid model. Consequently, in order to compare the RBM to the MSD data, we add a constant reflecting the cage-rattling contribution to the MSD. If the universal RBM prediction for the MSD corresponding to unit diffusion coefficient is denoted by F_RBM(t), the RBM prediction thus becomes

with three parameters c, α, and β. The long-time limit ⟨Δr²(t)⟩ = 6Dt results in cα = D. Note that, for dimensional reasons, there must be at least two parameters, a length and a time. Equation (2) is plotted as the green dashed line in Fig. 4. We conclude that despite having just a single shape parameter—compared to the two shape parameters of the von Schweidler expression [Eq. (1)]—Eq. (2) fits the MSD data very well.

In the following, we show that one can go one step further and compare the RBM to the dynamics of the mBLJ model using just the two scaling parameters, D and c, i.e., without any shape parameters. Already in 1969, Goldstein recognized the significance of potential-energy minima,²⁷ which were later termed “inherent states” by Stillinger and Weber.²⁸ If the N particle coordinates are collected into a single 3N-dimensional vector denoted by R, one can distinguish between the “true” Newtonian dynamics R(t) and its corresponding quenched “inherent” dynamics R_I(t). As illustrated in the inset of Fig. 4, the latter is arrived at by quenching configurations from an equilibrium MD simulation to their inherent states.²⁹ We run the same data analysis program on both the true configurations R(t) and the quenched inherent configurations R_I(t). Note that R_I(t) in the course of time jumps discontinuously from one constant vector to another. Below T_c, the dynamics separates into oscillations around inherent states and transitions between these,²⁹ as predicted by Goldstein.²⁷ The point is that the effect of oscillations is removed by considering the inherent dynamics. Thus, the inherent MSD, $⟨ΔrI2(t)⟩$⁠, can be compared directly to the RBM prediction without additional constants

In Fig. 4, the inherent MSD is plotted as crosses for all four temperatures. It obeys time–temperature superposition, and despite having no shape parameters, the RBM prediction [Eq. (3)] (the full green line) describes the data well.

In the above analysis, the scaling parameter c was determined by minimizing the root-mean-square difference between the inherent MSD and the RBM prediction. Subsequently, the plateau parameter β was fitted for the true MSD [Eq. (2)]. The true MSD can be fitted directly to Eq. (2), but using the resulting c parameters for the inherent MSD results in a considerably worse fit. This reflects the fact that the inherent contribution to the MSD is small at short times. Note, however, that the inherent MSD at short times is more than a factor 100 larger than the extrapolation of the diffusive regime (the black dashed line). This is comparable to the increase in fluidity observed in experimental data.¹ 

We turn now to the relation between our results and the experimental findings. If the RBM describes the liquid MSD and if one assumes that the macroscopic shear viscosity controls the microscopic frictional forces via the Stokes–Einstein relation between diffusion coefficient and viscosity, the frequency-dependent fluidity is proportional to the RBM universal prediction as found for nine liquids by Gainaru and co-workers.¹ Both assumptions are highly nontrivial, though. The first assumption is supported by our simulation results. The second assumption, the Stokes–Einstein assumption, is definitely challenged in glass-forming liquids.^7,30–33 We shall not discuss this further here, but note that for the arguments presented, it is enough that the frequency-dependent viscosity and diffusivity are inversely proportional—the proportionality coefficient may well depend on temperature.

Why was the plateau parameter β not necessary in the analysis of experimental data in Ref. 1? Letting tilde denote a suitable scaling, we get from Eq. (2) the following equation:

For the real part of the fluidity, which was the quantity investigated in Ref. 1, there is no contribution from the plateau parameter. Referring to the simulation data in Fig. 4, one can easily imagine the plateau parameter β to be non-universal, leading to a universal inherent MSD but a non-universal full MSD. By Eq. (4), this would lead to a universal real part of the fluidity, but a non-universal imaginary part. This may explain why fluidity universality was not noted before: The more commonly studied frequency-dependent viscosity, η(ω) = 1/F(ω), is not universal.

Why do particles in a viscous liquid move like in a disordered solid? The liquid is disordered, of course, but the more or less random potential-energy landscape experienced by any given particle changes with time. This argument refers to three dimensions. Taking a more abstract approach, it has been argued that the complexity may be replaced by randomness in the high-dimensional configuration space.³⁴ This is similar to the philosophy of statistical mechanics, which ignores the extreme complexity of a many-body system and models it by a probability distribution. This way of thinking about the problem also addresses the challenge that a given particle does not experience a frozen landscape. Figure 4 shows that the all-particle MSD is fitted well by the RBM prediction. Figure 5(a) shows the results of the same analysis restricted to the A particles, while Fig. 5(b) shows it for the B particles. The fit to the RBM is not as good for the A particles as for all particles (Fig. 4), indicating that the dynamics of the A particles by themselves is not accurately modeled by the RBM. Interestingly, the B particles are fitted better by the RBM. The B particles are smaller than the A particles and move faster, with a characteristic time for transition to diffusive dynamics that is about 1/10 of that of the A particles. This means that the B particles to a higher extent than the A particles move in a frozen three-dimensional landscape. Alternatively, because the best fit to the RBM is found for the all-particle MSD, one may speculate that the “correct” explanation of why the RBM works should refer to the high-dimensional configuration space, not 3d space.

Our findings focus on the shape of the MSD and do not have direct consequences for understanding the relaxation time’s temperature dependence as modeled in, e.g., the Adam-Gibbs, random first order transition, or shoving models.⁸ The latter views the metastable equilibrium supercooled liquid as behaving more like “a solid that flows” than like an ordinary liquid. The fact that particle motion in a highly viscous liquid in the present work has been shown to be much like particle motion in a disordered solid is in line with this view.

In summary, we have shown that the MSD of the KA model modified to avoid crystallization follows the zero-shape-parameter RBM prediction. This is consistent with the experimental findings of Ref. 1 that are well described by the RBM. Our results leave important questions for future works: Why do the mBLJ particles move as if they were hopping in a random solid with identical energy minima? How can one justify using a Stokes–Einstein argument for converting MSD to the frequency-dependent fluidity? Finally, how general are our findings?

This work was supported by the VILLUM Foundation (Grant Nos. 00016515 and 00023189).