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

TL;DR

This paper demonstrates that the mean-square displacement in a modified Lennard-Jones glass-forming liquid aligns with the random barrier model, revealing simple underlying dynamics in complex glassy systems.

Contribution

Introduces a crystallization-resistant Lennard-Jones mixture and shows its low-temperature dynamics fit the RBM without shape parameters.

Findings

Mean-square displacement matches RBM predictions

Simple hopping model explains complex glass dynamics

GPU simulations validate the model's applicability

Abstract

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 [S. P. Bierwirth \textit{et al.}, Phys. Rev. Lett. {\bf 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 (GPU)-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 can reproduce the complex and highly cooperative dynamics of glass-forming liquids.

Equations

$$⟨\Delta r^2(t)⟩=r_0^2+a(6Dt{)}^b+6Dt.$$

$$⟨\Delta r^2(t)⟩=c{F}_{\mathrm{RBM}}(\alpha t)+\beta$$

$$⟨\Delta r_I^2(t)⟩=c{F}_{\mathrm{RBM}}(\alpha t).$$

$$\tilde{F}(\tilde{\omega })={\tilde{F}}_{\mathrm{RBM}}(\tilde{\omega })+i\tilde{\omega }\beta$$

$$\frac{D(\tilde{s})}{D_0}=1+\sum _{j=1}^{10}{a}_j{\tilde{s}}^{j/10}$$

$$⟨\Delta r^2(\tilde{t})⟩=6D_0\tilde{t}\left(1+\sum _{j=1}^{10}\frac{a_j{\tilde{t}}^{1-j/10}}{\Gamma (2-j/10)}\right)$$

Full text

From experimental data for nine glass-forming liquids, Gainaru and coworkers recently demonstrated a striking universality of the real part of the frequency-dependent fluidity Bierwirth et al. (2017). The fluidity is defined as $1/\eta(\omega)$ in which $\omega$ is the angular frequency and $\eta(\omega)$ 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 Dyre (1985, 1988); Schrøder and Dyre (2000), which 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 see whether the RBM does describe the particle dynamics.

The relaxation time increases dramatically when a liquid is supercooled and approaches the glass transition Brawer (1985); Angell (1995); Debenedetti and Stillinger (2001); Dyre (2006); Berthier and Biroli (2011). A standard probe of the dynamics is the single-particle mean-square displacement (MSD) as a function of time, $\langle\Delta r^{2}(t)\rangle$, in which $\Delta 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 $\langle\Delta r^{2}(t)\rangle=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 like oxide glasses, polymers, amorphous semiconductors, etc Dyre (1985, 1988). In the extreme-disorder limit, i.e., when $k_{B}T$ 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 Schrøder and Dyre (2000). The only requirement for this to apply is that the energy barrier probability distribution is continuous at the percolation threshold Schrøder and Dyre (2000). Physically, the response is universal because the dynamics at extreme disorder is dominated by percolation in the 3d random energy landscape Schrøder and Dyre (2000).

In a simple analytical approximation the frequency-dependent diffusion coefficient of the RBM defined by $D(\omega)\equiv(-\omega^{2}/6)\int_{0}^{\infty}\langle\Delta r^{2}(t)\rangle\exp(i\omega t)dt$ is predicted to be the solution of $\ln\tilde{D}=(i\tilde{\omega}/\tilde{D})^{2/3}$ in which $\tilde{\omega}$ is a scaled frequency and $\tilde{D}\equiv D(\omega)/D(0)$ Schrøder and Dyre (2008). This is derived by combining the Alexander-Orbach conjecture that the percolation cluster has harmonic dimension $4/3$ (independent of dimension) Alexander and Orbach (1982) with the effective-medium approximation applied to diffusion on the percolation cluster Schrøder and Dyre (2008). The quoted equation provides an excellent fit to computer simulations of the RBM Schrøder and Dyre (2008), except at the lowest frequencies where the transition to a frequency-independent diffusion constant is better described by the solution of the following equation $\ln\tilde{D}=(i\tilde{\omega}/\tilde{D})(1+(8/3)i\tilde{\omega}/\tilde{D})^{-1/3}$ Schrøder and Dyre (2008).

MSD simulation data are conveniently fitted to the von Schweidler empirical expression Kob and Andersen (1995a)

[TABLE]

According to mode-coupling theory, the exponent $b$ is non-universal Götze (2008). Tokuyama has discussed common features of the MSD of different models Tokuyama (2011), but to the best or 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. Bierwirth et al., 2017, 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 Ninarello et al. (2017); Ozawa et al. (2019). 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-art graphics-processing unit (GPU) simulations Bailey et al. (2017) to access equilibrium dynamics at very low temperatures. The duration of the longest simulation was four months, which with traditional 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 $\varepsilon$ is a characteristic energy and $\sigma$ a characteristic length: $v(r)=4\varepsilon\left[\left({r}/{\sigma}\right)^{-12}-\left({r}/{\sigma}\right)^{-6}\right]$ Lennard-Jones (1924); Hansen and McDonald (2013). 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 Kob and Andersen (1995a). The trick is to have a strong AB non-ideal (non-Lorentz-Berthelot) attraction impeding phase separation. The parameters of the KA model are Kob and Andersen (1995a) $\sigma_{\rm BB}/\sigma_{\rm AA}=0.88$, $\sigma_{\rm AB}/\sigma_{\rm AA}=0.8$, $\varepsilon_{\rm AB}/\varepsilon_{\rm AA}=1.5$, and $\varepsilon_{\rm BB}/\varepsilon_{\rm AA}=0.5$. The KA model quickly became the standard model for simulations of viscous liquid dynamics Ashwin and Sastry (2003). The mode-coupling temperature (the temperature at which idealized mode-coupling theory based on higher-temperature data predicts a diverging relaxation time Ozawa et al. (2019)) was estimated to be $T_{c}=0.435$ Kob and Andersen (1995a). As computers became faster, it eventually became possible to investigate the model below $T_{c}$, see e.g. Refs. Ashwin and Sastry, 2003 and Coslovich et al., 2018. The KA model crystallizes in very lengthy simulations Toxvaerd et al. (2009); Ingebrigtsen et al. (2019); in fact, at the standard density 1.2 the KA liquid is supercooled whenever $T<T_{m}=1.03$ Pedersen et al. (2018). Although the strong AB attraction impedes phase separation, the supercooled KA system eventually does crystallize by phase separating into a pure A phase Pedersen et al. (2018).

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\,\sigma_{\rm AA}$ for the AA and BB interactions Allen and Tildesley (1987); Toxvaerd and Dyre (2011). 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 in 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 Allen and Tildesley (1987); Toxvaerd and Dyre (2011), 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 $\rho\equiv N/V=1.2$; the temperature $T$ was controlled by a Nose-Hoover thermostat. Unless otherwise noted, results are reported in standard MD-units defined by $\sigma_{AA}=1$, $\epsilon_{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 performed for this study. Figure 1(b) shows as a function of temperature the average crystallization time of the original KA model (red circles) Ingebrigtsen et al. (2019), 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^{8}$ MD time units, corresponding to 0.2 milliseconds in argon units. The data are for $N=8000$ particles; size independence was checked by simulating also $N=27000$ particles at $T=0.40$, which gave indistinguishable results. The dashed blue line gives 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 liquid dynamics.

At very short times one finds the well-known free-particle ballistic MSD $\propto t^{2}$ 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 MSD $\propto t$. 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 Dyre (2006); Berthier and Biroli (2011). 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 of the axes, as found in many experiments Bierwirth et al. (2017). 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 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 include 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_{\rm RBM}(t)$, the RBM prediction thus becomes:

[TABLE]

with three parameters $c,\,\alpha,\,\beta$. For dimensional reasons, there must be at least two parameters, a length and a time. The long-time limit $\langle\Delta r^{2}(t)\rangle=6Dt$ results in $c\alpha=D$. 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 model 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 Goldstein (1969), which were later termed “inherent states” by Stillinger and Weber Stillinger and Weber (1983). If the $N$ particle coordinates are collected into a single $3N$-dimensional vector denoted by $\mathbb{R}$, one can distinguish between the “true” Newtonian dynamics $\mathbb{R}(t)$ and its corresponding quenched “inherent” dynamics $\mathbb{R}_{\rm 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 Schrøder et al. (2000). We run the same data analysis program on both the true configurations $\mathbb{R}(t)$ and the quenched inherent configurations $\mathbb{R}_{\rm I}(t)$. Note that $\mathbb{R}_{\rm 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 Schrøder et al. (2000), as predicted by Goldstein Goldstein (1969). The point is that the effect of oscillations is removed by considering the inherent dynamics. Thus the inherent MSD, $\langle\Delta r_{I}^{2}(t)\rangle$, can be compared directly to the RBM prediction without additional constants:

[TABLE]

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) (full green line) agrees nicely with the data.

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 (Eq. (3)). Subsequently, the plateau parameter $\beta$ 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 (black dashed line). This is comparable to the increase in fluidity observed in experimental data Bierwirth et al. (2017).

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 coworkers Bierwirth et al. (2017). Both assumptions are highly nontrivial. The first assumption is supported by our simulation results presented above. The second assumption, the Stokes-Einstein assumption, is definitely challenged in glass-forming liquids Hodgdon and Stillinger (1993); Tarjus and Kivelson (1995); Douglas and Leporini (1998); Debenedetti and Stillinger (2001); Sengupta and Karmakar (2014). 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 $\beta$ not necessary in the analysis of experimental data in Ref. Bierwirth et al., 2017? Letting tilde denote a suitable scaling, we get from Eq. (2):

[TABLE]

For the real part of the fluidity, which was the quantity investigated in Ref. Bierwirth et al., 2017, there is no contribution from the plateau parameter. Referring to the simulation data in Fig. 4, one can easily imagine the plateau parameter $\beta$ to be non-universal, leading to a universal inherent MSD but a non-universal full MSD. Via 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, $\eta(\omega)=1/F(\omega)$, 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 complexity may be replaced by randomness in the high-dimensional configuration space Wolynes (1992). 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 addresses also the challenge that a given particle does not experience a frozen landscape. Figure 4 showed that the all-particle MSD is fitted well by the RBM prediction. Figure 5(a) shows results for the same analysis restricted to the A particles, while (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 they 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 landscape. Because the best fit to the RBM is found for the all-particle MSD, we speculate that the “correct” explanation of why the RBM works should refer to the high-dimensional configuration space, not to 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 Dyre (2006). The latter views the metastable equilibrium supercooled liquid as behaving more like “a solid that flows” than like an ordinary liquid Dyre (2006); the fact that particle motion in a highly viscous liquid in the present work has been shown to be much like particle motion in an 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. Bierwirth et al., 2017 that are also well described by the RBM. Our results leaves 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 and perhaps most importantly: How general are our findings?

I Appendix

The Random Barrier Model (RBM) was solved in the real-Laplace-frequency domain on a cubic lattice with a box distribution of energy barriers Schrøder and Dyre (2008). The frequency-dependent diffusion coefficient $D(s)$ for $\beta=320$ is shown in Fig. 6 as symbols Schrøder and Dyre (2008), where $\beta$ is the maximum barrier over the temperature. 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 is a good representation of the universal RBM prediction.

To facilitate transformation from the real-Laplace-frequency domain to the time domain over the many decades involved, we fit to the following empirical fitting function in which $\tilde{s}\equiv s/D(s)$

[TABLE]

The fitting was performed after taking the logarithm of both the numerical data and the fitting function, resulting in the following parameters $(a_{1},...,a_{10})$ = (-1.1914, 11.2368, -34.2903, 26.6019, 47.0002, -96.2905, 77.4671, -27.0535, 7.72535, -0.178844). The resulting fit is shown as the green dashed line in Fig. 6 and the corresponding error is shown in the lower panel.

Using $D(s)\equiv(s^{2}/6)\int_{0}^{\infty}\langle\Delta r^{2}(t)\rangle\exp(-st)dt$ (see the main text), Eq. (5) implies

[TABLE]

This is the equation used to represent the RBM in Fig. 4 and Fig. 5, in which a further empirical rescaling of $\tilde{t}$ is introduced.