Glass Transitions Structure of Glasses Theory and Modeling of Glass Transitions

The most conspicuous aspect of the phenomenon referred to as “the glass transition” is the immense increase in relaxation time of a supercooled liquid upon a modest reduction of its temperature. It is not uncommon to observe experimentally an increase of $10-14$ orders of magnitude in the measured viscosity. The phenomenon of slowing down is correlated with the dynamic heterogeneity in super-cooled liquids which is seen in both experiments and simulations [1]. Dynamic heterogeneity indicates the existence of a dynamical lengthscale which is revealed by studying multi-point correlation functions [2, 3, 4, 5, 6, 7, 8, 9]; more particles move in a correlated way when the temperature of the supercooled liquid is reduced. Although this is an interesting feature of supercooled liquids, it is still not clear to what extent dynamic correlations are the consequence or the primary origin of slow dynamics [8].

On the other hand a major theoretical question that is still not fully answered is whether the tremendous increase in relaxation time is accompanied by a corresponding increase in a typical (static) lengthscale [10, 11, 12, 13] similarly to what is observed in critical phenomena. The identification of a candidate lengthscale characterizing the glass transition attracted considerable amount of attention, with two candidates showing the most promise. The first method is based on the “point-to-set” (PTS) length [14, 15]. This length allows one to probe the spatial extent of positional amorphous order; it was measured in several numerical simulations  [16, 17, 18, 19, 20] and shown to grow mildly in the (rather high) temperature regime investigated.

Another, superficially unrelated way to define a static scale (hereafter denoted as $\xi$) was announced in Ref. [21] and employed further in Ref. [22]. The relation between these two differently defined lengthscales was studied carefully in Ref. [23] with the conclusion that they are in fact in full agreement, except that the point-to-set lengthscale is accessible at higher temperatures whereas the length $\xi$ is more accurately evaluated at lower temperatures.

The starting point for the definition of the latter scale is the fact that the low frequency tail of the density of states (DOS) of amorphous solids reflects the excess of plastic modes which do not exist in the density of states of purely elastic solids [24, 25]. This excess of modes is sometimes referred to as the ‘Boson Peak’ [26]. Here and below a ‘mode’ refers to an eigenfunction of the underlying Hessian matrix calculated at the local minimum of the potential energy function $U$ (the so called ‘inherent state’). The Hessian matrix is defined as ${\cal H}_{ij}^{\alpha\beta}=\frac{\partial^{2}U}{\partial x_{i}^{\alpha}\partial x_{j}^{\beta}}$ where $x_{i}^{\alpha}$ denotes the $\alpha^{th}$ component of coordinate of particle $i$. Recently [27] it was proposed that the eigenvalues $\{\lambda_{i}\}_{i=1}^{dN}$ (with $d$ being the space dimension and $N$ the number of particles) appear in two distinct families in generic amorphous solids, one corresponding to eigenvalues of the Hessian matrix that are only weakly sensitive to external strains; the other group consists of eigenvalues that go to zero at certain values of the external strain, thus leading to a plastic failure. The first group of modes is decently described by the Debye model of an elastic body, but this is not the case for the second group corresponding to the density of plastic modes. A simple model for the full density of states was suggested [27, 28] as a sum of the Debye and the plastic modes in the form

where the pre-factor $B(T)$ depends on the temperature, $\lambda_{D}\simeq\mu\rho^{2/d-1}$ is the Debye cutoff frequency and $\mu$ is the shear modulus. Particular models for the function $f_{\rm pl}\left(\frac{\lambda}{\lambda_{D}}\right)$ were presented in [28]. For our purposes here it is only important to understand that this function is a partial characterization of the degree of disorder which decreases upon approaching the glass transition.

The idea to determine the static typical scale is that the minimal eigenvalue $\lambda_{\rm min}$ observed in a system of $N$ particles will be determined by either the first or the second term in Eq. (1). For a system large enough, local disorder will be irrelevant in determining $\lambda_{\rm min}$, and it will be determined by the Debye contribution. For small systems the opposite is true. Thus there exists a system size where a cross-over occurs. This cross-over is interpreted in terms of a typical lengthscale $\xi$ separating correlated disorder from asymptotic elasticity.

The calculation of the length $\xi$ at various T requires equilibrating the supercooled liquid. For this reason, until now one could only compute $\xi$ over a range that corresponded to a change in the $\alpha$ relaxation time $\tau$ of up to $4-5$ orders of magnitude. In this Letter we explore a “swap Monte Carlo” [29] technique to reach a measurement of $\xi$ over a range that corresponds to a change in $\tau$ of between 12 and 16 orders of magnitude (depending of the form of the extrapolation or the functional dependence of $\tau$ on $\xi$ considered). The system under study is a ternary mixture of $N$ point particles made of $15\%$ particles A, $30\%$ particles B and $55\%$ particles C, which interact via a pure repulsive soft sphere potential similar to that in Ref. [29]:

Here, $\sigma_{ij}=\sigma_{i}+\sigma_{j}$, where $\sigma_{i}$ is determined by the particle type, and $r_{ij}\equiv|{\bm{r}}_{i}-{\bm{r}}_{j}|$. The energy scale $\epsilon$ as well as the particle mass $m$ are taken as unity. Time units are given by $\sqrt{m\,l_{0}^{2}/\epsilon}$, where $l_{0}$ gives the unit of length. The coefficients $c_{2l}$ are chosen in such a way that the potential and its first and second derivatives vanish at the cutoff $x_{c}=1.25$. We choose $\sigma_{A}$, $\sigma_{B}$ and $\sigma_{C}$ such that $\sigma_{A}/\sigma_{B}=\sigma_{B}/\sigma_{C}=1.25$, with $\sigma_{C}=0.424741977632123\,l_{0}$. This guarantees that the average value of the diameter of the dominating $\left(\frac{\sigma_{ij}}{r_{ij}}\right)^{12}$ soft sphere part of the potential is unity. The simulations are performed in the $NVT$ ensemble at a density $\rho=1.1$ with periodic boundaries in 3D. The use of three components and the choice of the compositions and scales are designed to avoid the crystallization observed in the original binary system which was reported elsewhere  [30, 31]. This system does not show a tendency to crystallize at any of the temperatures examined for the time spans considered, as evidenced in the profile of the radial distribution functions and by the absence of abrupt transitions in the potential energy time series.

The Swap Monte Carlo method combines standard Monte Carlo moves with particle swaps. Both types of moves are accepted or rejected according to the standard Metropolis algorithm (detailed balance holds for both types of moves under this criterion). While particle swaps are accepted at lower rates as the temperature is decreased, they nonetheless greatly accelerate equilibration by providing a way for the particles to break away from the cages in which they would be otherwise stuck, at least as long as the acceptance rate is of the order of $10^{-4}$ or even $10^{-5}$. An illustrating comparison between the regular and swap Monte Carlo methods is provided in Fig. 1, where it is shown that Swap Monte Carlo is able to relax the system at such low temperatures for which regular Monte Carlo shows no perceptible departure from the initial configuration after millions of sweeps. These results are in qualitative agreement with those observed for a binary mixture of soft particles in Ref. [29].

In order to study the typical lengthscale characterizing the glass transition, we need to define the range of temperatures at which we are able to equilibrate our system. For this task we consider molecular dynamics and Swap Monte Carlo. Regular Monte Carlo is not considered in the following for its relative inefficiency with respect to Swap Monte Carlo, as illustrated above. The upper panel of Fig. 2 shows the low T limit for which equilibrium can be attained via molecular dynamics, as evidenced by the deviation of the low T potential energy from the extrapolated supercooled liquid curve. The initial configurations are prepared by equilibrating a fluid with molecular dynamics at $T=0.5$ for $1000$ time units, and then cooling it down to the target temperature with a rate of $3.33\times 10^{-4}$, a procedure that is also followed to generate the initial conditions for the Swap Monte Carlo simulations. The potential energy was measured across a time window of $50000$ time units after an equilibration time of the order of ${10^{7}}$ time units. Despite the enormous computational efforts this required (it took several weeks of CPU time), it becomes impossible to equilibrate the system for $T\leq 0.26$, as the conspicuous departure from the supercooled liquid curve indicates. Additionally, during the same time window in which the potential energy was measured, we computed the self-intermediate scattering function,

where ${\bm{r}}_{i}(t)$ is the position of the $i^{th}$ particle at time $t$ and ${\bm{q}}$ is a wave vector. We choose the length of ${\bm{q}}$ to correspond to the first peak of the structure factor, and we average across 50 orientations (as orientation should not be of relevance in an isotropic system). Every time $F_{s}(q,t)$ decays to $1/e$, a ‘relaxation instance’ is said to occur, and $F_{s}(q,t)$ is restarted to $1$ (the current configuration is taken as $t=0$ reference configuration for the next decay). The relaxation time $\tau$ reported below in Figs. 4 and 5 is computed as the average duration of a relaxation instance.

As for the Swap Monte Carlo results, the arrival to equilibrium was assessed by testing the attainment of a canonical energy distribution as described below (cf. Fig. 2 lower panel). This was done for temperatures as low as ${T=0.22}$, the equilibration time for even lower temperatures being unreasonably long. Additionally we include ${T=0.20}$ in the equilibration analysis to illustrate the failure of equilibration. The resulting potential energies are shown in Fig. 2 upper panel and can be compared to those of molecular dynamics.

For the study of the static lengthscale below it is crucial to ascertain that our system is indeed equilibrated for $T\geq 0.22$. For this purpose, we follow Ref. [32] and collapse the probability distribution of the potential energy $P(U,T)$ for different temperatures on the probability distribution of a reference temperature ${T_{0}=0.30}$. The energy distribution is rescaled according to

which for equilibrated systems rescales the distribution to a canonical one at a reference temperature $T_{0}$. The rescaled $P_{T_{0}}(U,T)$ for $T=\left\{0.20,0.22,0.24,\ldots,0.40\right\}$ with $T_{0}=0.30$ are shown in the lower panel of Fig. 2. The excellent collapse confirms that the Swap Monte Carlo technique allows us to equilibrate the system at $T\geq 0.22$. In order to strengthen the evidence for attaining equilibrium for all $T\geq 0.22$ we have checked that our system reached detail balance for all these temperatures. We determined the Gaussian pdf for energy fluctuations at each temperature $T$, and computed the variance of this pdf. Detailed balance is guaranteed when the variance is $\sigma^{2}=kT^{2}C_{V}$ where $C_{V}$ is the specific heat. Computing $C_{V}$ from $\partial\langle U\rangle/\partial T$ we ascertained agreement to high precision, cf. inset in Fig. 2 upper panel. For lower temperatures even the present technique fails to equilibrate the system within reasonable time (see inset of Fig. 2, lower panel). While a decrease in temperature from $T=0.28$ to $T=0.22$ may not seem impressive to the uninitiated eye, we will show now that it has a major effect on the typical length.

The Swap Monte Carlo method is used to obtain a representative sampling of configuration space for a given temperature $T$ in the range for which equilibration is attainable. The decay of the intermediate scattering function allows us to focus on sufficiently uncorrelated configurations (i.e., configurations separated by one relaxation instance). These configurations are then quenched to zero temperature using conjugate gradient minimization [33], landing them on an inherent state. We accumulated $1500$ inherent states for each temperature with systems of sizes ranging from $N=954$ to $N=20000$. The minimal eigenvalue of the Hessian matrix is calculated at each of these inherent states using the Lanczos algorithm [34].

Note that in systems of finite size both $\lambda_{\rm min}$ and $\lambda_{D}$ are fluctuating from one amorphous realization to the other, and we therefore consider their mean over realizations, denoted as $\langle\lambda_{\rm min}\rangle$ and $\langle\lambda_{D}\rangle$.

To proceed we integrate the density of states (Eq. (1)) from zero to the average of the first eigenvalue. This integral picks up (on the average) one out of $dN$ eigenvalues. Because we integrate up to the average $\langle\lambda_{\rm min}\rangle/\langle\lambda_{\rm D}\rangle$ we do not get exactly 1 but ${\cal O}(1)$ in a finite system. (In an infinite system the exact $1$ is recaptured). Explicitly, $Nd\int_{0}^{\frac{\langle\lambda_{\rm min}\rangle}{\langle\lambda_{\rm D}\rangle}}P\left(x\right)dx={\cal O}(1)$. Introducing Eq. 1 into the integral and doing some simple algebra one can show

where ${\cal G}\left(\frac{\langle\lambda_{\rm min}\rangle}{\langle\lambda_{\rm D}\rangle}\right)\equiv\int_{0}^{\frac{\langle\lambda_{\rm min}\rangle}{\langle\lambda_{\rm D}\rangle}}f_{pl}(x)dx$ and $\tilde{A}$ and $\tilde{B}$ are the same as $A$ and $B$ up to absorbed $T$-independent constants. Changing slightly this equation, writing $\xi^{d}(T)\equiv\frac{1}{\rho\tilde{B}(T)}$, one obtains the following scaling function

where ${\mathcal{F}}\equiv{\mathcal{G}}^{-1}$. Since the function ${\mathcal{G}}$ is monotonically increasing, its inverse is well defined. The typical scale $\xi(T)$ is calculated by demanding that all the data obtained for different system sizes and temperatures should collapse into a master curve just by appropriately choosing the $\xi(T)$.

The raw results for the present method, before rescaling, are shown in the upper panel of Fig. 3. The same results after rescaling by $\xi(T)$ are shown in the lower panel of the same figure, and the resulting typical scale as a function of $T$ is shown in in the inset. It is evident that the addition of the three (lowest temperature) points to the typical scale $\xi(T)$, results in an increase of $\xi(T)$ by a substantially larger factor than that reached by direct molecular dynamics simulations [22].

The Hessian-based analysis becomes less reliable at higher temperatures due to anharmonic effects. Thus for higher temperatures we employ yet another approach to determine the lengthscale $\xi$, extending its range of validity even further. In Ref. [23], it was argued that the point-to-set method offers a computation of the lengthscale at high temperatures. It was also found that finite size effects of the relaxation times $\tau$ at different temperatures are perfectly correlated to the point-to-set lengthscale. This was tested using the lengthscale obtained in Ref.  [23] for the Kob-Andersen model system. Since a finite size scaling analysis of the relaxation time as in Ref. [22] is significantly easier than calculating the point-to-set lengthscale, we employ the former to extract the static lengthscale at higher temperatures as explained next.

In the upper panel of Fig. 4, we show the system size dependence of the relaxation time measured by molecular dynamics, rescaled by the value for the largest $N$ considered. The scattered data themselves are a first indication that even at higher temperatures the effect of the static lengthscale is substantial. In the lower panel we have collapsed the data by appropriately choosing the lengthscale.

In both the minimal eigenvalue and the finite size $\tau$ scaling methods the lengthscales are defined up to a constant factor. For the latter method we choose this factor so that the lengthscale is unity at the highest temperature considered. The factor for the minimal eigenvalue method is chosen appropriately, considering that the results obtained from both methods should match at some intermediate temperature range. In Fig. 5 upper panel, we show the combined results of both methods. One can clearly see that the static lengthscale varies by more than 500%, which is quite remarkable.

Presently we relate the results of the measured lengthscale $\xi(T)$ to the $\alpha$ relaxation time $\tau$ as measured from the intermediate scattering function for temperatures $T\geq 0.28$ (and to extrapolations for lower $T$). Our values of $\tau(T)$ were measured by averaging over at least fifty relaxation times. Traditionally the relaxation times are fitted to the Vogel-Fulcher-Tamman (VFT) formula: $\tau=\tau_{0}\exp\left(\frac{A}{T-T_{0}}\right)\ .$ We fitted the three parameters in this formula to our measured relaxation times in the range of temperatures $T\in[0.28,0.50]$, as shown in Fig. 5 lower panel, continuous black line. While this formula fits well the measured relaxation time, it predicts a divergence of the relaxation time at $T_{0}=0.23$. Since we succeed to equilibrate the system at $T=0.22$ we propose that this formula is untenable for the present system, ruling out the VFT fit.

Another proposed fit formula is the so called modified Bässler proposition  [35], $\tau=\tau_{0}\exp\left[B\left(\frac{T_{1}}{T}-1\right)^{2}\right]\ ,$ which is again a three parameter fit. Fitting again to our available relaxation times we find the dashed line in Fig. 5 lower panel. The small deviation of this fit at temperatures above $T_{1}=0.48$ is expected, and has been attributed to the facilitated particle motion at those high temperatures  [35]. However since the relaxation time of $4.6\times 10^{5}$ expected by this fit at $T=0.26$ is within the reach of our simulation (spanning $10^{7}$ time units), and we were not able to equilibrate the system by molecular dynamics at that temperature (see upper panel of Fig. 2), we can rule out also this extrapolation. Indeed, the authors in Ref. [35] predict the existence of a lower bound for the range of validity of this fit providing a reason for the inadequacy of this extrapolation.

Finally we refer to the fit formula proposed in Ref. [22]:

which is a 2-parameter fit since $\xi(T)$ is determined. Note that in previous publications it was unclear whether the length scale $\xi$ in the exponent in Eq. 7 should be raised to power $d$ or not. In Fig. 5 we also show the predictions of the same formula with $\xi$ raised to power 1/2 and 2. The former is very close to the Bässler prediction, and the latter predicts relaxation times that are two long. The fit with 1/2 can again be dismissed since it predicts a relaxation time that could be reached by MD simulations where we know that equilibration was not achieved. The fit with 2 overestimates the relaxation time even at $T=0.28$ where MD data exists. Assuming that the Eq. 7 is applicable, and determining the parameters from the measured relaxation times, we obtain the curve shown in a continuous grey line in Fig. 5 lower panel. While this is not rigorously justified at this time, the reader should note that at the lowest temperature one predicts $\tau=1.64\times 10^{15}$, which if correct extends the analysis presented here to a range of $15$ orders of magnitude in relaxation times.

In conclusion, we have presented a calculation of the typical static lengthscale $\xi$ in a range of temperatures that spans the interval $T\in[0.22,0.50]$, with an increase in $\xi$ by a factor of about 5, being unprecedented at this point in time. Whichever model relating this lengthscale to the relaxation time $\tau$ one takes, this increase in $\xi$ is associated with an increase in $\tau$ by a factor between $10^{12}$ and $10^{16}$, meaning that we are able to offer an analysis that is competitive with the best experiments.