It was recently shown that the real part of the frequency-dependent fluidity for several glass-forming liquids of different chemistry conforms to the prediction of the random barrier model (RBM) devised for ac electrical conduction in disordered solids [Bierwirth et al., Phys. Rev. Lett. 119, 248001 (2017)]. Inspired by these results, we introduce a crystallization-resistant modification of the Kob–Andersen binary Lennard-Jones mixture for which the results of extensive graphics-processing-unit-based molecular-dynamics simulations are presented. We find that the low-temperature mean-square displacement is fitted well by the RBM prediction, which involves no shape parameters. This finding highlights the challenge of explaining why a simple model based on hopping of non-interacting particles in a fixed random energy landscape with identical minima can reproduce the complex and highly cooperative dynamics of 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.1 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, ⟨Δr2(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 ⟨Δr2(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 kBT 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.4 The only requirement for this to apply is that the energy barrier probability distribution is continuous at the percolation threshold.4 Physically, the response is universal because the dynamics at extreme disorder is dominated by percolation in the 3d random energy landscape.4
In a simple analytical approximation, the frequency-dependent diffusion coefficient of the RBM defined by is predicted to be the solution of , in which is a scaled frequency and .10 This is derived by combining the Alexander–Orbach conjecture that the percolation cluster has harmonic dimension 4/3 (independent of dimension)11 with the effective-medium approximation applied to diffusion on the percolation cluster.10 The quoted equation provides an excellent fit to computer simulations of the RBM,10 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: .10 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 expression12
According to mode-coupling theory, the exponent b is non-universal.13 Tokuyama has discussed common features of the MSD of different models,14 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) simulations17 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 , 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.12 The trick is to have a strong AB non-ideal (non-Lorentz–Berthelot) attraction impeding phase separation. The parameters of the KA model are12 σ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.20 The mode-coupling temperature (the temperature at which idealized mode-coupling theory based on higher-temperature data predicts a diverging relaxation time16) was estimated to be Tc = 0.435.12 As computers became faster, it eventually became possible to investigate the model below Tc (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 < Tm = 1.03.24
Although the strong AB attraction impedes phase separation, the supercooled KA system eventually does crystallize by phase separating into a pure A phase.24 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.
(a) Pair potentials in units of εAA plotted as a function of pair distance in units of σAA. Dashed lines: The standard Kob–Andersen (KA) pair potential with shifted potential cutoffs at r = 2.5 σαβ. Full lines: The modified binary Lennard-Jones (mBLJ) pair potentials, which introduce shifted-force cutoffs at r = 1.5 σAA for the AA and BB interactions and at r = 2.5 σAB for the AB interaction. The significant reduction of the AA attraction obtained in this way suppresses the tendency to phase separate. (b) The red circles show the average crystallization times of the KA model23 and the red dashed curve is a parabolic fit to these data. This figure also shows the simulation times for the mBLJ model at four temperatures: T = 0.37, 0.38, 0.39, and 0.40 (black crosses). Simulation times are scaled to be comparable to the crystallization times of Ref. 23, which used N = 10 000. At T = 0.40, several independent simulations were performed, none of which crystallized. From this fact one deduces a higher estimate of the minimum crystallization time. At each temperature, the black rectangles indicate the estimated range of possible crystallization times for the mBLJ model.
(a) Pair potentials in units of εAA plotted as a function of pair distance in units of σAA. Dashed lines: The standard Kob–Andersen (KA) pair potential with shifted potential cutoffs at r = 2.5 σαβ. Full lines: The modified binary Lennard-Jones (mBLJ) pair potentials, which introduce shifted-force cutoffs at r = 1.5 σAA for the AA and BB interactions and at r = 2.5 σAB for the AB interaction. The significant reduction of the AA attraction obtained in this way suppresses the tendency to phase separate. (b) The red circles show the average crystallization times of the KA model23 and the red dashed curve is a parabolic fit to these data. This figure also shows the simulation times for the mBLJ model at four temperatures: T = 0.37, 0.38, 0.39, and 0.40 (black crosses). Simulation times are scaled to be comparable to the crystallization times of Ref. 23, which used N = 10 000. At T = 0.40, several independent simulations were performed, none of which crystallized. From this fact one deduces a higher estimate of the minimum crystallization time. At each temperature, the black rectangles indicate the estimated range of possible crystallization times for the mBLJ 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, mA = mB = 1, and kB = 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)23 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 108 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.
All-particle MSD as a function of time for the mBLJ system at the four temperatures T = 0.37, 0.38, 0.39, and 0.40. Full lines are best fit to the von Schweidler expression [Eq. (1)].12 At short times, the data follow the line of slope 2 expected from ballistic behavior; at long times, the slope is unity reflecting diffusive behavior. The liquid dynamics of the mBLJ model is close to that of the KA model, which is illustrated by the blue dashed line giving the MSD of the KA model at T = 0.40.
All-particle MSD as a function of time for the mBLJ system at the four temperatures T = 0.37, 0.38, 0.39, and 0.40. Full lines are best fit to the von Schweidler expression [Eq. (1)].12 At short times, the data follow the line of slope 2 expected from ballistic behavior; at long times, the slope is unity reflecting diffusive behavior. The liquid dynamics of the mBLJ model is close to that of the KA model, which is illustrated by the blue dashed line giving the MSD of the KA model at T = 0.40.
At very short times, one finds the well-known free-particle ballistic MSD ∝ t2, 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.1 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.
(a) All-particle MSD at the four temperatures plotted as a function of time scaled by the diffusion coefficient D, showing data corresponding to t > 20 in MD units. The data superpose not just at long times but also in the transition region. The short-time “plateau” regions, however, change with temperature. The inset shows results for longer times, demonstrating that the diffusive behavior is given by the red dashed curve. (b) The same data as in (a) but scaled further on both axes by a squared empirical length c, showing a near-perfect collapse, i.e., demonstrating time–temperature superposition. The black line is the von Schweidler fit to the T = 0.40 data.
(a) All-particle MSD at the four temperatures plotted as a function of time scaled by the diffusion coefficient D, showing data corresponding to t > 20 in MD units. The data superpose not just at long times but also in the transition region. The short-time “plateau” regions, however, change with temperature. The inset shows results for longer times, demonstrating that the diffusive behavior is given by the red dashed curve. (b) The same data as in (a) but scaled further on both axes by a squared empirical length c, showing a near-perfect collapse, i.e., demonstrating time–temperature superposition. The black line is the von Schweidler fit to the T = 0.40 data.
All-particle MSD of true and inherent dynamics (full lines and plus symbols) for four temperatures compared to the RBM prediction [Eqs. (2) and (3)], respectively. The fit parameters were determined as described in the text. The agreement with the inherent dynamics demonstrates that the RBM, despite having no shape parameters, provides a good representation of the liquid dynamics. Adding a constant to represent the contribution from cage-rattling gives a very good agreement with the true MSD.
All-particle MSD of true and inherent dynamics (full lines and plus symbols) for four temperatures compared to the RBM prediction [Eqs. (2) and (3)], respectively. The fit parameters were determined as described in the text. The agreement with the inherent dynamics demonstrates that the RBM, despite having no shape parameters, provides a good representation of the liquid dynamics. Adding a constant to represent the contribution from cage-rattling gives a very good agreement with the true MSD.
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 FRBM(t), the RBM prediction thus becomes
with three parameters c, α, and β. The long-time limit ⟨Δr2(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,27 which were later termed “inherent states” by Stillinger and Weber.28 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 RI(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.29 We run the same data analysis program on both the true configurations R(t) and the quenched inherent configurations RI(t). Note that RI(t) in the course of time jumps discontinuously from one constant vector to another. Below Tc, the dynamics separates into oscillations around inherent states and transitions between these,29 as predicted by Goldstein.27 The point is that the effect of oscillations is removed by considering the inherent dynamics. Thus, the inherent MSD, , 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.1
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.1 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.34 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.
Testing the RBM separately for each particle species with (a) showing the A particle MSD vs the RBM prediction and (b) showing the same for the B particles.
Testing the RBM separately for each particle species with (a) showing the A particle MSD vs the RBM prediction and (b) showing the same for the B particles.
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.8 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).
APPENDIX: NUMERICAL APPROXIMATION TO THE RBM MEAN-SQUARE DISPLACEMENT AS A FUNCTION OF TIME
The Random Barrier Model (RBM) was solved in the real-Laplace-frequency domain on a cubic lattice with a box distribution of energy barriers.10 The frequency-dependent diffusion coefficient D(s) for β = 320 is shown in Fig. 6 as symbols,10 where β is the maximum barrier over temperature for the box distribution of barriers. At low temperatures, the shape of D(s) for the RBM becomes universal, i.e., independent of temperature and energy barrier distribution. In the frequency range shown, the data in Fig. 6 are a good representation of the universal RBM prediction.
Symbols: Numerical solution of the frequency-dependent diffusion coefficient D(s) for the random barrier model at β = 320, where s is the real Laplace frequency. The box distribution was used for the energy barriers, and β is the inverse temperature times the maximum barrier. The green dashed line is the fit according to Eq. (A1). The lower panel shows the deviation between fit and numerical data.
Symbols: Numerical solution of the frequency-dependent diffusion coefficient D(s) for the random barrier model at β = 320, where s is the real Laplace frequency. The box distribution was used for the energy barriers, and β is the inverse temperature times the maximum barrier. The green dashed line is the fit according to Eq. (A1). The lower panel shows the deviation between fit and numerical data.
To facilitate transformation from the real-Laplace-frequency domain to the time domain over the many decades involved, we fit to the following function:
in which .
The fitting was performed after taking the logarithm of both the numerical data and the fitting function, resulting in the following parameters (a1, …, a10) = (−1.1914, 11.2368, −34.2903, 26.6019, 47.0002, −96.2905, 77.4671, −27.0535, 7.725 35, −0.178 844). The resulting fit is shown as the green dashed line in Fig. 6, and the corresponding error is shown in the lower panel.
Using (see the main text), Eq. (A1) implies by inverse Laplace transformation