We show that the aging dynamics of a strong glass former displays a strikingly simple scaling behavior, connecting the average dynamics with its fluctuations, namely, the dynamical heterogeneities. We perform molecular dynamics simulations of SiO2 with van Beest-Kramer-van Santen interactions, quenching the system from high to low temperature, and study the evolution of the system as a function of the waiting time tw measured from the instant of the quench. We find that both the aging behavior of the dynamic susceptibility χ4 and the aging behavior of the probability distribution P(fs,r) of the local incoherent intermediate scattering function fs,r can be described by simple scaling forms in terms of the global incoherent intermediate scattering function C. The scaling forms are the same that have been found to describe the aging of several fragile glass formers and that, in the case of P(fs,r), have been also predicted theoretically. A thorough study of the length scales involved highlights the importance of intermediate length scales. We also analyze directly the scaling dependence on particle type and on wavevector q and find that both the average and the fluctuations of the slow aging dynamics are controlled by a unique aging clock, which is not only independent of the wavevector q, but is also the same for O and Si atoms.

If a glass-forming liquid is cooled from a high temperature to a low temperature and crystallization is avoided, the relaxation times of the system increase dramatically. Depending on the experimental (or simulation) time accessible in comparison with this growing relaxation time, either a supercooled liquid (in equilibrium) or a glass (out of equilibrium) is observed.1–3 In the non-equilibrium (aging) case, after a temperature quench, the dynamics at low temperature depends on the waiting time tw — the time elapsed since the temperature quench. To investigate this rich dynamics, a large variety of approaches (experiments, computer simulations, and theoretical techniques) have been used and many different systems have been studied. Previous work on the dynamics of both of supercooled liquids and of glasses ranges from small molecules, to polymers, to network glasses, to colloidal glasses, and to granular systems (in the last two cases density is the control parameter, instead of temperature). For reviews we refer the reader to Refs. 1–5.

A common finding of these studies, and in that sense a universal feature, is that the dynamics is spatially heterogeneous, meaning that there are fast and slow regions in space.6–8 One route for probing the extent of universality is to investigate the possible presence of similar scaling behaviors of the dynamical heterogeneities in diverse systems. We take this route in the work presented here. Specifically, we investigate the scaling of the dynamic susceptibility χ4 and of the distribution P(fs,r) of the local incoherent intermediate scattering function fs,r. In previous work, there have been relatively few studies on P(fs,r).9–14 Also, most previous work on χ4 has focused on the dependence of χ4 on the temperature (or density) in the supercooled liquid regime (see Sec. 3.2.4.3 of Ref. 6 and for SiO2 specifically see Refs. 15–18) There have been fewer studies for the aging dynamics, i.e., the dependence of χ4 on the waiting time tw, which we discuss in this paper.13,19–25

Predictions for the scaling of P(fs,r) with respect to tw follow from a theoretical framework for the aging dynamics that explains dynamical heterogeneities in terms of the presence of Goldstone modes associated with a broken continuous symmetry under time reparametrizations.9–11,26–32 To study the dynamical heterogeneity, i.e., the local fluctuations in the relaxation, we focus on a local two-time correlation, the local incoherent intermediate scattering function fs,r(tw, tw + t), which depends on the position r, the waiting time tw, and the time interval t. The Goldstone modes correspond to space dependent shifts of the time variable tϕr(t) such that

fs,r(tw,tw+t)C(ϕr(tw),ϕr(tw+t)).
(1)

Here C(tw, tw + t) corresponds to the global two-time correlation function.9,31 More generally, a simple Landau-theory approximation for the dynamical action predicts that quantities describing fluctuations in the system depend on the waiting time tw and the time interval t essentially only through the global two-time correlation function C(tw, tw + t).9–11 Thus it is expected that the probability distribution P(fs,r(tw, tw + t)) should collapse for different waiting times tw, for (tw, t) pairs chosen such that C(tw, tw + t) is held fixed. This prediction is consistent with spin glass simulation results.9–11 Despite the theory being initially derived for spin glasses, simulation results for structural glasses12,13,20,21,25,30,31 and experimental results for a polymer glass22 find this predicted scaling to hold. The numerical simulations also show that χ4(tw, tw + t) is a product of two factors: a waiting-time dependent scale that grows with tw and a scaling function that depends on (tw, t) only through the value of C(tw, tw + t).

The large variety of structural glass formers can be divided into two broad groups, called fragile and strong glass formers, due to their different dependence of the viscosity (and the relaxation time) on temperature.1–3 All of the previous tests of predictions of the Goldstone mode approach in structural glasses have been for the case of fragile glass formers. It is an open question whether the behavior of dynamical heterogeneity in strong glass formers is also well described by the same theoretical framework. In this paper, we address precisely that question. We present here molecular dynamics simulation results for the network former SiO2 which is a strong glass former. The van Beest-Kramer-van Santen (BKS) potential33 which we use has not only been shown to be an excellent model for real silica34–37 but also previous work of the last 20 years provides us with detailed insight into many of the properties of this system, including its phase diagram,38–45 energy landscape,46–50 structure34,35,38–40,51–53,82 vibrational spectrum,37,54–57 dynamical heterogeneities,18,58–63 and aging.64–68 

Further motivation for the present work is the unexpected similarity that has recently been found between the dynamics of the strong glass former SiO2 and fragile glass formers.67 Whereas in Ref. 67 the microscopic dynamics is studied via single particle jump analysis, we investigate in this paper whether this surprising similarity of strong and fragile glass dynamics also holds true for the scaling of dynamical heterogeneities. We find not only that indeed most results confirm universal dynamics, but we also gain deeper insight into the involved length and time scales. The scaling of P(fs,r) uncovers the importance of intermediate length scales and the scaling of C, χ4, and P(fs,r) all indicates a common aging clock which is the same for Si and O atoms.

To model amorphous SiO2, we used the BKS potential.33 We carried out molecular dynamics (MD) simulations with NSi = 112 silica atoms and NO = 224 oxygen atoms, at a constant volume V=16.920468Å3 which corresponds to a density ρ = 2.323 g/cm3. For further details on the interaction see Ref. 66.

At 6000 K we generated 200 independent configurations (at least 1.63 ns apart) which then were fully equilibrated at initial temperature Ti = 5000 K for 3.27 ns, followed by an instantaneous quench to lower temperature Tf = 2500 K, i.e., below Tc = 3330 K.69 Unique to our simulations is that we applied the Nosé-Hoover temperature bath at Tf only for the first 0.327 ns (NVT) and then continued with constant energy (NVE) for 98.1 ns. We switch to (NVE), i.e., the integration of Newton’s equations, to avoid any algorithm specific influence of a thermostat. The shortest possible application of the (NVT) thus allows us to disturb the dynamics minimally. We confirmed that Tf stays constant and is similar to Tf(t) as shown in Fig. 2 of Ref. 66. In all the following results, we find no indication of a sudden change in the dynamics due to this switch from (NVT) to (NVE), instead all investigated quantities show a smooth transition at the corresponding time. The MD time step was 1.02 fs and 1.6 fs during the (NVT) and (NVE) runs, respectively. In what follows, we analyzed the combined (NVT) and (NVE) simulation runs at Tf.

The main difference between the present simulation and the ones discussed in Refs. 51, 66, and 67 is that our new dataset has increased statistics (200 independent runs instead of 20) and that each NVE run at Tf has a longer duration (98.1 ns instead of 32.7 ns). As described in Sec. III, this increased statistics and the longer simulation runs allowed us to gain insight into the scaling of the two point correlation function Cα significantly beyond the results of Ref. 66. Furthermore, having 200 independent simulation runs gave good enough statistics to make it possible to determine the dynamic susceptibility and the distribution of the local incoherent intermediate scattering function (see Secs. IV and V).

Let us first look at the global generalized incoherent intermediate scattering function

Cα(tw,tw+t,q)=fsα(tw,tw+t,q)
(2)

with

fsα(tw,tw+t,q)=1Nαj=1Nαcosq(rj(tw+t)rj(tw)).
(3)

Here rj(t) is the position of particle j at time t, tw is the waiting time elapsed since the temperature quench from 5000 K to 2500 K, and Nα is the total number of particles of type α (α ∈ {Si, O, all}). The notation indicates an average over wave vectors q of fixed magnitude q and over the 200 independent simulation runs. In Eq. (3) the sum is over particles in the complete simulation box, whereas in Sec. V the sum is only over particles within a local sub-box. We call Cα the “global” incoherent intermediate scattering function to stress this distinction. Error bars for Cα are given by the statistical error of the average over the 200 independent simulation runs. Even though averaging over 200 independent runs allows us to remove a lot of the noise in the results, we further smooth the results by additionally applying a time average. The time average is computed by using logarithmic time bins and by averaging t-values, Cα-values, and error bars ΔCα within the same time bin.

Unique to the present work is that we investigate directly the influence of particle type on scaling. To do so, we distinguish three different cases, labeled by the symbol α. In the case of α = Si, the sum in Eq. (2) is exclusively over Si atoms and in the case of α = O the sum is exclusively over O atoms. In the third case, α = all, the sum is over all particles, i.e., including both Si and O atoms. We use this notation throughout the whole paper, including in our discussion of the dynamic susceptibility χ4 and the probability distribution P(fs,r).

Fig. 1 shows Cα(tw, tw + t, q) for oxygen atoms and for q = 2.7 Å−1. We find that with increasing waiting time tw the correlation function decays more slowly. To quantify this, we define the relaxation time τqα as the time when Cα has decayed to a certain value Ccutα(q),

Cα(tw,tw+τqα,q)=Ccutα(q).
(4)
FIG. 1.

Global generalized incoherent intermediate scattering function at q = 2.7 Å−1 for oxygen atoms, CO(tw, tw + t, q) as defined in Eq. (2), for 12 waiting times tw between 0 ns and 81.9 ns. To avoid cluttering the graph, in this figure and in most other figures where results for different waiting times are compared, statistical error bars are shown for just one of the waiting times. Here and in the following figures, the shown error bars are of the same order as the not shown error bars for different waiting times.

FIG. 1.

Global generalized incoherent intermediate scattering function at q = 2.7 Å−1 for oxygen atoms, CO(tw, tw + t, q) as defined in Eq. (2), for 12 waiting times tw between 0 ns and 81.9 ns. To avoid cluttering the graph, in this figure and in most other figures where results for different waiting times are compared, statistical error bars are shown for just one of the waiting times. Here and in the following figures, the shown error bars are of the same order as the not shown error bars for different waiting times.

Close modal

Instead of the commonly used choice Ccutα=1/e, we adjust to the varying plateau height of Cα(q) for different q. Hence for Ccutα(q), we choose the values listed in Table I, each of which is given by 1/e times the corresponding plateau value of Cα(q).

TABLE I.

Ccutα(q) values.

qSiOAll
1.7 0.323 0.298 0.306 
2.7 0.265 0.217 0.233 
3.4 0.221 0.162 0.182 
4.6 0.144 0.085 0.105 
qSiOAll
1.7 0.323 0.298 0.306 
2.7 0.265 0.217 0.233 
3.4 0.221 0.162 0.182 
4.6 0.144 0.085 0.105 

The resulting relaxation times for oxygen atoms are shown in Fig. 2 as functions of the waiting time tw. As in Ref. 66, different regimes for tw can be identified from Figs. 1 and 2. In Ref. 66 it was found that for small waiting times, tw ≲ 0.1 ns, Cα(tw, tw + t) does not form a plateau, for intermediate tw a plateau is formed, and time superposition applies, and for sufficiently large waiting times Cα becomes tw-independent, i.e., equilibrium is reached. The increased statistics of the present simulations allow the identification in Fig. 2 of the transition from small to intermediate tw as a change in the tw dependence of the relaxation times τα(tw) from non-power law to power law behavior. Fig. 2 also shows that τqα(tw) does not reach a plateau, i.e., the waiting times are not long enough to reach equilibrium.

FIG. 2.

Relaxation times for oxygen atoms τqO(tw) as defined in Eq. (4). For clarity, the data for q = 4.6 Å−1, 3.4 Å−1, 2.7 Å−1, and 1.7 Å−1 have been shifted by factors of 1, 3, 9, 27, respectively. The lines are power law fits. The fitted exponents μ are shown in the inset, both for the case of O atoms and for the case of all atoms.

FIG. 2.

Relaxation times for oxygen atoms τqO(tw) as defined in Eq. (4). For clarity, the data for q = 4.6 Å−1, 3.4 Å−1, 2.7 Å−1, and 1.7 Å−1 have been shifted by factors of 1, 3, 9, 27, respectively. The lines are power law fits. The fitted exponents μ are shown in the inset, both for the case of O atoms and for the case of all atoms.

Close modal

Let us next investigate further the dynamics for intermediate waiting times tw ≳ 0.1 ns. The inset in Fig. 2 shows the power law fit exponents μ as functions of wave vector q. We find that within the error bars, μ seems to be independent of q for q ≥ 2.7 Å−1. Similar results for τ(tw) have been found experimentally for a metallic glass70 and for a colloidal glass.71 The inset in Fig. 2 also shows that μ is independent of the particle type α. (The particle type is indicated by a square for α = O and by a rhombus for α = all.) This independence of α is rather surprising, since Horbach and Kob had found that the dynamics of silicon and oxygen atoms is very different for temperatures below 3330 K.35 Saksaengwijit and Heuer72 relate this decoupling of silicon and oxygen dynamics to rotational processes.

We interpret the α-independence of μ as an evidence for the existence of a common “single aging clock” in the system, despite the different dynamics of Si and O atoms. This allows us to analyze Si and O atoms together (α = all). To directly test the hypothesis of a single aging clock, we generalize an approach introduced by Kob and Barrat.73 They had investigated the q-dependence of Cα via a parametric plot of Cα(q2) versus Cα(q1) for various tw. Whereas for their system, a binary Lennard-Jones system, they found no data collapse,73 for our system, SiO2, it was found in Ref. 66 that data collapse indeed happens. This indicates that Cα(tw,tw+t,q)=Cα(z̃(tw,t,α),q).66 In other words, for each particle type α, there is a unique q-independent aging clock represented by z̃(tw,t,α). Notably, this seems the only non-universal result we encounter in our comparison of the dynamics of fragile and strong glass formers. However, Kob and Barrat included the full range of waiting times, whereas we included only tw ≥ 0.14 ns. When we include in a parametric plot of Cα(q2) versus Cα(q1) also too short waiting times tw < 0.14 ns, we obtain no data collapse for SiO2 similar to the results of Kob and Barrat. According to Fig. 4 of Ref. 73, their smallest three or four waiting times correspond to too short waiting times. When we exclude these waiting times in their Fig. 8, their results for the binary Lennard-Jones system show data collapse similar to our results for SiO2.

We now address the question of the existence of a common aging clock for Si and O atoms, in other words, whether or not the function z̃(tw,t,α) is independent of α. The question is therefore whether it is true that

Cα(tw,tw+t,q)C(tw,tw+t,q,α)=?C(z(tw,t),q,α).
(5)

To answer this question, we investigate directly the α-dependence via a parametric plot of CO(tw, tw + t, q) versus CSi(tw, tw + t, q) as shown in Fig. 3 for q = 2.7 Å−1 (and in the inset for q = 1.7 Å−1). We find almost perfect data collapse for all investigated q and conclude that Eq. (5) is correct. Hence the (tw, t)-dependence is solely governed by one function z(tw, t), i.e., an “inner aging clock” which is not only q-independent but also the same for different particle types.

FIG. 3.

To study directly the dependence of the global incoherent intermediate scattering function Cα on the particle type α, we show here two parametric plots of CO versus CSi for various waiting times tw: one for q = 2.7 Å−1 (main panel) and another for q = 1.7 Å−1 (inset).

FIG. 3.

To study directly the dependence of the global incoherent intermediate scattering function Cα on the particle type α, we show here two parametric plots of CO versus CSi for various waiting times tw: one for q = 2.7 Å−1 (main panel) and another for q = 1.7 Å−1 (inset).

Close modal

This may appear surprising at first, since, as mentioned before, it is known that in SiO2 the oxygen atoms have a faster dynamics than the silicon atoms.35,72 Our results do not contradict this statement. To illustrate how to reconcile a common clock and yet different dynamics of Si and O atoms, we show in Fig. 4 (tw, tw + t)-pairs for fixed Cα. For example, the blue circles of the bottom curve were obtained by finding for each tw the corresponding tw + t for which CSi(tw, tw + t, q = 2.7 Å−1) = 0.575 with 1% accuracy. For equilibrium dynamics, this curve would be trivial: it would be the set of (tw, tw + t)-pairs for a certain constant value of t. For aging dynamics, the curve is non-trivial: t changes as tw changes. Nevertheless, we obtain the same non-trivial curve for O atoms (blue triangles) for fixed CO = 0.463. Similarly we obtain identical curves for Si atoms (circles) and oxygen atoms (triangles) for different choices of CSi and CO (upper three curves), hence a common clock for Si and O atoms. The fact that the dynamics of O atoms is faster than the dynamics of Si atoms is reflected by the constant value of CO being lower than the corresponding constant value of CSi on the same curve.

FIG. 4.

This figure illustrates the compatibility of Si and O atoms having both different speeds as well as a common clock. For each tw the corresponding (tw + t) was determined to obtain the specified Cα = Cα(tw, tw + t, q = 2.7 Å−1) with 1% accuracy. The solid lines are for the guidance of the eye. The common clock is apparent by identical curves for silicon (circles) and oxygen (triangles). The fact that the dynamics of O atoms is faster than the dynamics of Si atoms is reflected by the constant value of CO being lower than the corresponding constant value of CSi on the same curve.

FIG. 4.

This figure illustrates the compatibility of Si and O atoms having both different speeds as well as a common clock. For each tw the corresponding (tw + t) was determined to obtain the specified Cα = Cα(tw, tw + t, q = 2.7 Å−1) with 1% accuracy. The solid lines are for the guidance of the eye. The common clock is apparent by identical curves for silicon (circles) and oxygen (triangles). The fact that the dynamics of O atoms is faster than the dynamics of Si atoms is reflected by the constant value of CO being lower than the corresponding constant value of CSi on the same curve.

Close modal

In this section, we study the dynamic susceptibility, χ4α, which is a four-point correlation function. χ4α quantifies thermal fluctuations of the incoherent intermediate scattering function. Following the notation of Berthier,17χ4α is defined74 to be

χ4α(tw,tw+t,q)=Nαfsα2fsα2,
(6)

where fsα(tw,tw+t,q), Nα, and are as defined at the beginning of Sec. III. To obtain error bars for χ4α, we divide the 200 independent simulation runs into 20 subsets each of 10 independent simulation runs and compute a value χ4(α,i) for each subset i, with i = 1, …, 20. The error bar of χ4α is the standard deviation of the mean over those 20 independent χ4(α,i) values. As in the case of Cα, we further smooth the results for χ4α by applying a time average using logarithmic time bins. This means that for each logarithmic bin, all unsmoothed data (t, χ4, Δχ4), which occur at a time t within the specified bin time window, are averaged to (t¯,χ4¯,Δχ4¯).

Fig. 5 shows the resulting dynamic susceptibility for oxygen atoms and q = 2.7 Å−1. Since χ4α is a four-point correlation function, it is a measure of dynamic heterogeneities. The dynamic susceptibility is small both for very short times and for very large times t and it has a maximum χmaxα at an intermediate time tmaxα. This maximum can be interpreted as a maximal number of particles in a dynamically correlated region. For a thorough discussion of this maximum and its scaling dependence on temperature and system size in the case of a supercooled liquid and a dense granular system, see Refs. 16, 17, and 75. We investigate here instead the aging dynamics and thus the dependence of χ4α on waiting time tw. To quantify the dependence of this maximum on tw, q, and α, we show the peak height χmaxα and the peak position tmaxα as functions of tw in Figs. 6 and 7, respectively. Similar to previous results in fragile glass formers,13,20,21,23–25χmaxα and tmaxα increase with increasing tw. It is possible that χmaxα reaches a plateau for large tw, but the noise in the results is too large to allow for any definite conclusions to be drawn. As above, the dependence of tmaxα(tw,q) on tw is consistent with the presence of two regimes: tw ≲ 0.1 ns and tw ≳ 0.1 ns. As in the case of τqα(tw), in the longer time regime, the time scale tmaxα(tw,q) has a power law dependence on tw, with exponents μmax which are again independent of particle type, as shown in the inset of Fig. 7.

FIG. 5.

Dynamic susceptibility χ4O(tw,tw+t,q) as defined in Eq. (6) for q = 2.7 Å−1, for the same waiting times tw as in Fig. 1.

FIG. 5.

Dynamic susceptibility χ4O(tw,tw+t,q) as defined in Eq. (6) for q = 2.7 Å−1, for the same waiting times tw as in Fig. 1.

Close modal
FIG. 6.

Peak value, χmaxα, of χ4α(tw,tw+t,q) (see Fig. 5). Here χmaxα as a function of waiting time tw is shown for oxygen atoms, i.e., α = O. We find similar results for α = Si and α = all.

FIG. 6.

Peak value, χmaxα, of χ4α(tw,tw+t,q) (see Fig. 5). Here χmaxα as a function of waiting time tw is shown for oxygen atoms, i.e., α = O. We find similar results for α = Si and α = all.

Close modal
FIG. 7.

Peak position tmaxO of the maximum of χ4O(tw,tw+t,q) (see Fig. 5). For clarity, the data for q = 4.6 Å−1, 3.4 Å−1, 2.7 Å−1, and 1.7 Å−1 have been shifted by factors of 1, 3, 9, 27, respectively. The lines are power law fits tmaxOtwμmax with exponents as shown in the inset.

FIG. 7.

Peak position tmaxO of the maximum of χ4O(tw,tw+t,q) (see Fig. 5). For clarity, the data for q = 4.6 Å−1, 3.4 Å−1, 2.7 Å−1, and 1.7 Å−1 have been shifted by factors of 1, 3, 9, 27, respectively. The lines are power law fits tmaxOtwμmax with exponents as shown in the inset.

Close modal

We next address the question of how the dynamic susceptibility scales with respect to waiting time. As described in the Introduction, numerical simulations for fragile glasses find13,20,21,25 a scaling behavior of the (tw, t) dependence of the dynamic susceptibility given by

χ4α(tw,tw+t,q)=χ40(tw,q,α)ϕ(Cα(tw,tw+t,q),q,α).
(7)

A similar more general result follows from Ref. 76. Without loss of generality, we choose χ40(tw,q,α) in Eq. (7) to be the maximum height χmaxα. To test the validity of Eq. (7), we plot in Fig. 8χ4α/χmaxα as a function of 1Cα for the case of oxygen atoms (α = O) and q = 2.7 Å−1 (for details see endnote77). Fig. 8 shows indeed data collapse within the error bars. To quantify how good this data collapse is, we determine the left crossing point of a horizontal line at 0.6 in Fig. 8, i.e., (1 − Cα)cross is defined as the smaller of the two solutions of the equation

χ4α/χmaxα((1Cα)cross)=0.6.
(8)
FIG. 8.

We show here χ4O/χmaxO as a function of (1 − CO) for different waiting times tw. The data collapse confirms Eq. (7), i.e., the normalized dynamic susceptibility depends on tw only via the global incoherent intermediate scattering function Cα. We find similar data collapse for all other investigated q and α.

FIG. 8.

We show here χ4O/χmaxO as a function of (1 − CO) for different waiting times tw. The data collapse confirms Eq. (7), i.e., the normalized dynamic susceptibility depends on tw only via the global incoherent intermediate scattering function Cα. We find similar data collapse for all other investigated q and α.

Close modal

Fig. 9 shows the resulting (1 − Cα)cross as a function of tw for α = O, for q = 1.7 Å−1, 2.7 Å−1, 3.4 Å−1, and 4.6 Å−1. To a first approximation, (1 − Cα)cross is waiting time independent, consistent with Eq. (7). However, a slight systematic increase of (1 − Cα)cross with tw appears to be present, which may be evidence for the existence of small corrections to Eq. (7).

FIG. 9.

(1 − CO)cross as defined in Eq. (8) is plotted as a function of waiting time tw.

FIG. 9.

(1 − CO)cross as defined in Eq. (8) is plotted as a function of waiting time tw.

Close modal

We now investigate the dependence of χ4 on wave vector q and particle type α. We use Eq. (5) to rewrite Eq. (7) as

χ4α(tw,tw+t,q)=χmax(tw,q,α)ϕ(C(z(tw,t),q,α),q,α).
(9)

The dependence on time t enters into Eq. (9) only via C and therein only through z(tw, t), where z is independent of q, as it has been shown in Ref. 66, and independent of α, as shown above in Fig. 3. Therefore we have

χ4α(tw,tw+t,q)=χmax(tw,q,α)χˆ(z(tw,t),q,α),
(10)

where χˆ(z,q,α) is a tw-independent function in the sense that all (t, tw)-dependences enter only via z. As a consequence of z being independent of q, we obtain data collapse for a parametric plot of χ4α/χmaxα for q = q2 as a function of χ4α/χmaxα for q = q1, with q2q1, and similarly as a consequence of z being independent of α we obtain data collapse for a parametric plot of (χ4O/χmaxO) versus (χ4Si/χmaxSi) (see the  Appendix).

The common aging clock for Si and O atoms allows for the analysis of Si and O atoms together, thus leading to the data collapse of χ4all/χmaxall(1Call) for different tw, which is shown in Fig. 10 for q = 2.7 Å−1 and quantified for all investigated q values via (1 − Call)cross in the inset of Fig. 10. This common aging clock might be the reason why a data collapse was also found in previous work on fragile glass formers, in which case different particle types were analyzed together.13,20,21

FIG. 10.

Similar to Fig. 8, we show here χ4α/χmaxα as a function of (1 − Cα) for different waiting times tw but now for α = all, i.e., both Si- and O-atoms were included in the analysis. The inset shows the corresponding (1 − Call)cross as a function of tw for all investigated q values.

FIG. 10.

Similar to Fig. 8, we show here χ4α/χmaxα as a function of (1 − Cα) for different waiting times tw but now for α = all, i.e., both Si- and O-atoms were included in the analysis. The inset shows the corresponding (1 − Call)cross as a function of tw for all investigated q values.

Close modal

In Sec. IV, we found scaling for the dynamic susceptibility, which can be thought of as a measure of the thermal fluctuations of the global incoherent intermediate scattering function fsα(tw,tw+t,q). In this section, we present results on the probability distribution for the local coarse grained intermediate scattering function

fs,rα(tw,tw+t,q)=1Nrαrj(tw)Brcosqrj(tw+t)rj(tw),
(11)

where the sum is over particles of type α which are at time tw within a local sub-box Br.78 By contrast, in Eq. (3) the sum is over all particles in the system. Our definition of fs,r is identical to the definition of Cr in Refs. 9–11, 26, and 28. We choose a different notation here to emphasize the fact that fs,r is not an ensemble-averaged quantity. By definition fs,r = 1 for t = 0. Relaxation in a region corresponds to the decay of the value of fs,r from 1 to 0. Spatial fluctuations of fs,r quantify dynamical heterogeneities: the “slow” regions have values of fs,r that remain non-negligible for a longer time, and “fast” regions correspond to local values of fs,r that decay more rapidly towards 0.

In the following, we determine the probability distribution P(fs,rα(tw,tw+t,q)) of the local correlations fs,rα. As described in Sec. I, a Landau-theory approximation for spin glasses9–11,26–29,32 predicts for this distribution P(fs,rα) that all (tw, t)-dependences are solely governed by Cα(tw, tw + t, q). This is rather surprising, since the prediction is for the full distribution P(fs,rα) of these local fluctuations, yet Cα(tw, tw + t, q) is not only a scalar but also a global quantity which is equal to the average

Cα(tw,tw+t,q)=fsα(tw,tw+t,q).
(12)

The theory therefore predicts that if (tw, t) pairs are chosen such that Cα(tw, tw + t, q) is fixed, the corresponding P(fs,rα) should be tw independent. This data collapse has been confirmed for spin glasses9,11 and for fragile structural glass formers.12,13,21,30 In this section, we investigate the scaling of P(fs,rα(tw,tw+t,q)) for our SiO2 simulation data, i.e., for a strong glass former.

As will be shown below, we find that the goodness of the scaling depends on the involved length scales via q and via the chosen size of the local sub-box Br. Since the specifics of the analysis influence the chosen length scales, we include in the following all necessary details. Our procedure for the choice of sub-box size and the corresponding set of sub-boxes within the complete simulation box of length L = 16.9205 Å is as follows: We first divide the simulation box into very small sub-boxes of length L/M. The length of a sub-box Br is then an integer b times this very small sub-box, i.e., L × b/M. The average number of particles in Br is therefore 〈Nr〉 = Nα(b/M)3. We present in this paper results for 〈Nr〉 ≈ 5 and 〈Nr〉 ≈ 40 as listed in Table II. The distribution P(fs,rα(tw,tw+t,q)) is then obtained for a specific set of (tw, t, q, α) and for a specific simulation run via measurements of fs,rα for all M3 possible Br and for all q-vectors of magnitude q. To obtain the M3 possible sub-boxes Br, the sub-box is shifted in the three directions and periodic boundary conditions were used.

TABLE II.

Table of analyzed (M, b)-values which specify the local sub-box Br.

αMbNrMbNr
Si 14 5.1 14 10 40.8 
5.22 41.8 
All 5.25 42.0 
αMbNrMbNr
Si 14 5.1 14 10 40.8 
5.22 41.8 
All 5.25 42.0 

To test whether the (tw, t) dependence of P(fs,rα) is governed by Cα(tw, tw + t, q), we use the same approach as in previous work.12,13,21,30 We choose a fixed value Cfix of the global Cα, and for each waiting time tw we determine the time tfix such that Cα(tw, tw + tfix, q) = Cfix (see endnote79). For these time pairs (tw, tfix) we then determine for each independent simulation run c = 1, 2, …, 200, a distribution Pc(fs,rα(tw,tw+tfix,q)), using for Eq. (11) the positions ri(tw)c and ri(tw+tfix)c. Please note that the parameters tw, tw + tfix, and q are common to all simulation runs. For each bin of the distribution, we obtain the mean Pc(fs,rα(tw,tw+tfix,q)) over the 200 independent simulation runs. The error bars represent the standard deviation of this mean over runs. The result of this computation is shown in Fig. 11 for the case of α = O and 〈Nr〉 = 5.1.

FIG. 11.

Distribution of the local incoherent intermediate scattering function for oxygen atoms, P(fs,rO), for q = 2.7 Å−1 and 〈Nr〉 = 5.1.

FIG. 11.

Distribution of the local incoherent intermediate scattering function for oxygen atoms, P(fs,rO), for q = 2.7 Å−1 and 〈Nr〉 = 5.1.

Close modal

For small Cα we find perfect scaling collapse and the distribution is a Gaussian (black dashed line). We attribute the Gaussian distribution to small Cα occurring at late times t (see Fig. 1) when diffusive dynamics is approached. For intermediate and large Cα we find a non-Gaussian, i.e., non-trivial P(fs,r), and nevertheless very good data collapse. To check quantitatively whether the slight tw dependence in Fig. 11 is systematic, we focus on the location of most discrepancy, the maximum. Fig. 12 shows the maximum value Pmaxα of the distribution, as a function of tw. Consistent with the results above, we find that for small tw scaling is not valid, but for tw ≳ 0.1 ns Pmaxα is approximately independent of tw. We obtain similar results for q > 2.7 Å−1 and also for α = Si.

FIG. 12.

To quantify how well scaling is satisfied in Fig. 11, we show here the maximum PmaxO as a function of the waiting time tw. Here 〈Nr〉 = 5.1 and q = 2.7 Å−1. We find similar results for q > 2.7 Å−1.

FIG. 12.

To quantify how well scaling is satisfied in Fig. 11, we show here the maximum PmaxO as a function of the waiting time tw. Here 〈Nr〉 = 5.1 and q = 2.7 Å−1. We find similar results for q > 2.7 Å−1.

Close modal

To probe the dependence of P(fs,rα) on particle type α, we show in the inset of Fig. 13 the comparison of P(fs,rα) for α = Si, O, all. We conclude that P(fs,rα) does depend on α. In Secs. III and IV we had found that despite different dynamics of silicon and oxygen atoms, their scaling z(tw, t) gives rise to a common aging clock. This common clock allows us to analyze Si and O together (α = all). To test this common aging clock for the case of P(fs,rα), we therefore show in Fig. 13P(fs,rall) and find indeed a data collapse for different tw.

FIG. 13.

Distribution of the local incoherent intermediate scattering function for α = all, i.e., when both silicon and oxygen atoms are analyzed. The color coding for the waiting times is the same as in Fig. 11. The dashed black lines correspond to Gaussian fits. For comparison, the inset shows P(fs,rα) for different α with fixed waiting time tw = 16.7 ns.

FIG. 13.

Distribution of the local incoherent intermediate scattering function for α = all, i.e., when both silicon and oxygen atoms are analyzed. The color coding for the waiting times is the same as in Fig. 11. The dashed black lines correspond to Gaussian fits. For comparison, the inset shows P(fs,rα) for different α with fixed waiting time tw = 16.7 ns.

Close modal

So far we have shown P(fs,rα) for various fixed Cα, for q = 2.7 Å−1, 〈Nr〉 ≈ 5, with either α = O or α = all. When also q and 〈Nr〉 are varied, we find that the predicted data collapse occurs as long as tw ≳ 0.1 ns, q ≥ 2.7 Å−1, and 〈Nr〉 ≈ 5. Next we look at cases when scaling fails. Fig. 14 shows P(fs,rα) for 〈Nr〉 ≈ 5 as before, but now q = 1.7 Å−1. For CO = 0.3 and CO = 0.5 we find that the data collapse for different tw is much worse than before. This is quantified in the inset of Fig. 14, which shows the systematic tw dependence of PmaxO. The inset also shows that in the Gaussian case of CO = 0.1, the scaling does work even for q = 1.7 Å−1. Furthermore, when the sub-box size is chosen to be much larger, scaling does not occur even if q ≥ 2.7 Å−1, as shown in Fig. 15 for the case of 〈Nr〉 ≈ 42. A similarly worse data collapse for larger 〈Nr〉 has been found for a fragile glass former (Fig. 2 of Ref. 13).

FIG. 14.

P(fs,rO) for q = 1.7 Å−1. Scaling is less good when q is too small. The inset shows the maximum value, PmaxO, as a function of tw. Even for tw > 0.1, PmaxO is tw-dependent, which indicates that there is a breakdown of scaling.

FIG. 14.

P(fs,rO) for q = 1.7 Å−1. Scaling is less good when q is too small. The inset shows the maximum value, PmaxO, as a function of tw. Even for tw > 0.1, PmaxO is tw-dependent, which indicates that there is a breakdown of scaling.

Close modal
FIG. 15.

P(fs,rO) for sufficiently large q = 2.7 Å−1 but with a large sub-box size, 〈Nr〉 = 41.8, leading to tw-dependent distributions. To quantify the breakdown of scaling the inset shows the peak value PmaxO as a function of the waiting time tw.

FIG. 15.

P(fs,rO) for sufficiently large q = 2.7 Å−1 but with a large sub-box size, 〈Nr〉 = 41.8, leading to tw-dependent distributions. To quantify the breakdown of scaling the inset shows the peak value PmaxO as a function of the waiting time tw.

Close modal

Thus scaling breaks down for intermediate timescales, corresponding to intermediate Cα, when the regions probed by the local incoherent intermediate scattering function fs,r become too large. Those probed regions can become larger either directly because the coarse graining region is chosen to be larger, or indirectly, because q is chosen to be smaller, thus allowing longer displacements to contribute significantly to fs,r.

One effect that could contribute to the imperfect collapse is that the measured probability distribution could be influenced to some degree by spatial correlation effects, and these effects could vary as the system ages and the typical size of the correlated regions grows. In the case of fragile glasses,12 it was found that the width of P(fs,r) grew with tw at constant Cα. In that case it was argued that P(fs,r) was narrowed by averaging of fs,r over more than one correlated region, but this narrowing became weaker as the size of the correlated regions grew with tw. In our case, however, we notice in the insets of Figs. 14 and 15 that the direction of the effect is not always the same: the distributions widen with increasing tw for CO = 0.3, but they narrow with increasing tw for CO = 0.5. Therefore, for SiO2, although this narrowing effect could in principle play some role, there must also be other effects at play.

In what follows we address the question of why scaling fails and how to adjust the analysis to recover the data collapse even for larger length scales (q = 1.7 Å−1 and 〈Nr〉 ≈ 40). To gain this insight, first a closer look at the details of the analysis is necessary. A crucial point is how we choose the time pairs, i.e., tw and tw + t. We illustrate in Fig. 16 how this is done. First we choose a unique t = tfix — the same for all 200 independent simulation runs — by demanding that the global incoherent scattering function Cα(tw, tw + t) (thick line in Fig. 16) take a certain value at t = tfix, such as Cα(tw, tw + tfix) = 0.5.

FIG. 16.

Ccα for single simulation runs c = 127, 57, 7, and 121 (thin lines) and the averaged incoherent intermediate scattering function Cα (thick line). In this example α = O.

FIG. 16.

Ccα for single simulation runs c = 127, 57, 7, and 121 (thin lines) and the averaged incoherent intermediate scattering function Cα (thick line). In this example α = O.

Close modal

If we now look at individual runs, the intermediate scattering function computed for each run c is Ccα=fsc where 〈…〉c corresponds to an average over q vectors of fixed magnitude but not over simulation runs. Since the system simulated in each independent run contains only 336 particles, it is not large enough for Ccα to be self-averaging: in Fig. 16, the values of Ccα(tw,tw+tfix) for four individual runs are shown with circles, and they differ dramatically from each other and from the value of the fully averaged Cα(tw, tw + tfix). In other words, choosing a unique value t = tfix is equivalent to choosing very different values for Ccα for each simulation run c. Consequently, the distributions P(fs,r(tw, tw + t, q)) at the same times are necessarily very different for different runs c.

In Fig. 17, single simulation run distributions Pc(fs,rO) for four independent runs are shown as thin black lines, and it is clear that the variation of Pc(fs,rO) between runs is very large: there is a nontrivial distribution of distributions. For comparison, the average over runs P(fs,rO) is shown in the same figure with a thick black line. All of the black lines in the figure correspond to tw = 0.49 ns. In the same figure, results are shown for another waiting time, tw = 33.0 ns, as magenta/grey lines: the thin lines corresponding to Pc(fs,rO) for individual runs c, and the thick line corresponding to P(fs,rO). Also, for tw = 33.0 ns the distributions Pc(fs,rO) vary greatly from simulation run to simulation run. Yet, the set of possible Pc(fs,rO) seems to be the same for the two waiting times. We find that for obtaining a particular shape of Pc(fs,rα), the key variable is Ccα(tfixO) (marked by circles in Fig. 16). For example, in Fig. 17 the two thin dashed lines correspond to two different simulation runs, c = 26 for tw = 33.0 ns and c = 57 for tw = 0.49 ns, chosen so that in both cases CcO(tfixO)0.3. This makes the two distributions close enough that their differences are of the order of their statistical error. The figure also shows that the same procedure is successful for obtaining other pairs of nearly identical distributions for CcO(tfixO)0.1,0.5,0.7 (for more details see endnote80).

FIG. 17.

Distributions of local incoherent intermediate scattering functions. Thick lines correspond to the average distribution P(fs,r) and thin lines are examples for individual simulation run distributions Pc(fs,r). Colors (black and magenta) indicate the waiting times.

FIG. 17.

Distributions of local incoherent intermediate scattering functions. Thick lines correspond to the average distribution P(fs,r) and thin lines are examples for individual simulation run distributions Pc(fs,r). Colors (black and magenta) indicate the waiting times.

Close modal

This leads us directly to finding a way to improve the scaling even for longer length scales. We no longer use a unique time tfix which is the same for all simulation runs c. Instead, we specify a fixed value Cmesoα of the correlation, and we define, for each simulation run c, a time tmesoα(c) such that

Cmesoα=fsα(tw,tw+tmesoα(c),q)c,
(13)

as shown in Fig. 16. Thus, for a specified Cmesoα, and for each run c = 1, 2, …, 200, we determine tmesoα(c), the corresponding Pc(fs,rα(tw,tw+tmesoα(c),q)), and then we average all of the individual distributions Pc(fs,rα) to obtain P(fs,rα) (see endnote81).

The resulting average distribution P(fs,rα) is shown in Fig. 18 for q = 1.7 Å−1 and 〈Nr〉 = 5.2. The evolution of the maximum PmaxO with waiting time is shown in the inset. The comparison with Fig. 14 confirms that the scaling is drastically improved by fixing Cmesoα instead of Cα. To test the limits of this improved scaling, we show in Fig. 19 the distributions P(fs,r) for q = 1.7 Å−1, 〈Nr〉 = 42, and α = all. This is the most unfavorable case, with the lowest wavevector we have considered, the largest coarse graining region, and including both O and Si atoms (which additionally tests whether the aging clock is the same for both particle types). Even for this most unfavorable case, we find almost perfect scaling. We thus conclude that Cmesoα is the appropriate scaling quantity.

FIG. 18.

Distribution of the local incoherent intermediate scattering function. P(fs,rα) is an average over single simulation run distributions Pc(fs,rα(tw,tmesoα(c),q)) where tmesoα(c) is chosen such that Cmesoα is fixed. The inset shows the maximum value as a function of waiting time tw. For comparison with Fig. 14 we show the case of q = 1.7 Å−1 and 〈Nr〉 = 5.2 for oxygen atoms, but here we fixed Cmesoα instead of Cα.

FIG. 18.

Distribution of the local incoherent intermediate scattering function. P(fs,rα) is an average over single simulation run distributions Pc(fs,rα(tw,tmesoα(c),q)) where tmesoα(c) is chosen such that Cmesoα is fixed. The inset shows the maximum value as a function of waiting time tw. For comparison with Fig. 14 we show the case of q = 1.7 Å−1 and 〈Nr〉 = 5.2 for oxygen atoms, but here we fixed Cmesoα instead of Cα.

Close modal
FIG. 19.

Probability distributions of local incoherent intermediate scattering functions, with α = all (i.e., both O and Si atoms are included), for small q = 1.7 Å−1 and large 〈Nr〉 = 42, for different waiting times. The color coding for the waiting times is the same as in Fig. 11. Despite this being the most unfavorable case for a successful collapse, we find that choosing Cmesoα to be constant results in an almost perfect data collapse among results for different waiting times.

FIG. 19.

Probability distributions of local incoherent intermediate scattering functions, with α = all (i.e., both O and Si atoms are included), for small q = 1.7 Å−1 and large 〈Nr〉 = 42, for different waiting times. The color coding for the waiting times is the same as in Fig. 11. Despite this being the most unfavorable case for a successful collapse, we find that choosing Cmesoα to be constant results in an almost perfect data collapse among results for different waiting times.

Close modal

In summary, we have performed molecular dynamics simulations of the strong glass former SiO2 to investigate the scaling of dynamical heterogeneities in this system. We have quenched the system from an initial high temperature Ti to a final temperature Tf below the mode-coupling critical temperature Tc and observed the out of equilibrium dynamics as a function of the waiting time tw, the time elapsed since the temperature quench. In particular, we have investigated the global incoherent intermediate scattering function Cα(tw, tw + t, q), the dynamic susceptibility χ4α(tw,tw+t,q), and the distribution P(fs,rα(tw,tw+t,q)) of the local incoherent intermediate scattering function, where q corresponds to the wave vector magnitude and α specifies the particle type.

We have found that for sufficiently long waiting times tw, and when probing small enough regions in the system, the dependence on (tw, t) of χ4α(tw,tw+t,q) and of P(fs,rα(tw,tw+t,q)) is governed by Cα(tw, tw + t, q), up to a tw-dependent scale factor in χ4. This is consistent with predictions for spin glasses and similar to previous results for fragile glass formers. We thus conclude that the behavior of dynamical heterogeneity in the aging regime of glassy systems shows a remarkable degree of universality. A similarity of the behavior for strong and fragile glass formers had previously been shown for the microscopic dynamics of single particle jumps,67 but here we have shown that it extends to the scaling of dynamical heterogeneities. We would like to note that this similarity of the dynamics of SiO2 and the comparison fragile glass formers is rather surprising, since SiO2 is a network of directional bonds, whereas the comparison fragile glass formers are sphere-type models. As such, the common scaling of dynamic heterogeneities seems to give us essential insight into relaxation processes. We leave it for future work to reveal what exactly leads to this similar dynamics.

Furthermore we have studied directly the influence of the particle type α on the dynamics. We have found that Cα(tw, tw + t, q) = C(z(tw, t), q, α), where z(tw, t) is independent of q and α. Thus z(tw, t) plays the role of a “common aging clock” that determines the slow aging behavior of the relaxation for both Si and O atoms. By combining this statement with the fact, discussed above, that the aging of χ4α and of P(fs,rα) is controlled by Cα, it follows that the aging of χ4α and of P(fs,rα) for both Si and O atoms should be controlled by that same unique aging clock. Our results show that this prediction is indeed satisfied. In summary, we have found that both the average and the fluctuations of the slow aging dynamics are controlled by a unique aging clock, which is independent of the wavevector q and is the same for O and Si atoms.

When fluctuations are probed over larger regions, either by taking q ≤ 1.7 Å−1, or by considering larger coarse graining regions containing around 40 particles, new phenomena emerge, presumably due at least in part to the fact that the probed regions contain more than one correlation volume. In particular, the scaling of P(fs,r) discussed above no longer holds in its initial form. It is clear that the probability distributions Pc(fs,r) obtained from the independent runs c = 1, …, 200 vary dramatically from run to run if the time interval is kept the same across runs. In other words, there is a nontrivial distribution of distributions. This spread of distributions becomes more pronounced the larger the spread of Ccα. Since the fluctuations in Ccα are given by χ4α and since χmaxα is increasing with decreasing wave vector, we find a particularly large spread of distributions Pc(fs,r) for the smallest wave vector. The average over Pc(fs,r) is not independent of tw. This is equivalent to the statement that, in a very large system, a measurement of P(fs,r) over a mesoscopic region containing a few hundred particles is not self-averaging, and that a new significant intermediate lengthscale emerges. It is, however, possible to recover an excellent collapse of probability distributions for different waiting times tw if instead of averaging probability distributions Pc(fs,r) from different runs (or mesoscopic regions) at constant time interval t = tfix, one averages probability distributions taken at constant mesoscopic intermediate scattering function Cmesoα. We speculate that this importance of Cmesoα is system independent because the limitations of scaling with Cα have been found also for fragile glass formers.12,13 The testing of the analysis with Cmesoα for other systems remains to be done. Furthermore it might be interesting to test whether a significantly larger system and several orders more independent simulation runs would even recover self-averaging. We leave it for future work to gain further insight into the importance and origin of this intermediate length scale.

We thank A. Parsaeian for preliminary work. K.V.L. and H.E.C. thank A. Zippelius and the Institute of Theoretical Physics, University of Göttingen, for hospitality and financial support. C.H.G. was supported by NSF REU Grant No. PHY-1156964. This work was supported in part by the Deutsche Forschungsgemeinschaft via Nos. SFB 602 and FOR1394, by DOE under Grant No. DE-FG02-06ER46300, and by Ohio University. Numerical simulations were carried out at Bucknell University and Ohio University. Part of this work was performed at the Aspen Center for Physics, which is supported by National Science Foundation Grant No. PHY-1066293.

From Eq. (10) it follows that all the (t, tw)-dependence of χˆ enters only via z. Since z is independent of q, it follows that a parametric plot of χ4α/χmaxα for q = q2 as a function of χ4α/χmaxα for q = q1, with q2q1, should show data collapse. In Fig. 20, we show a parametric plot of this kind, for q2 = 2.7 Å−1, q1 = 1.7 Å−1, and α = O, and find that indeed there is data collapse within the error bars. We also find similar results for other values of q2, q1, and α (not shown).

FIG. 20.

Parametric plot of χ4O/χmaxO(tw,tw+t,q=2.7Å1) versus χ4O/χmaxO(tw,tw+t,q=1.7Å1). We find data collapse among the results for different waiting times.

FIG. 20.

Parametric plot of χ4O/χmaxO(tw,tw+t,q=2.7Å1) versus χ4O/χmaxO(tw,tw+t,q=1.7Å1). We find data collapse among the results for different waiting times.

Close modal

For the dependence on particle type, α, Eq. (10) predicts that there should also be data collapse in a parametric plot of (χ4O/χmaxO) versus (χ4Si/χmaxSi). This data collapse is confirmed with Fig. 21 for q = 2.7 Å−1. We find equally good collapse for all other investigated q values. We emphasize that the data collapse shown in Figs. 20 and 21 is non-trivial, in the sense that the data are not along the diagonal, implying that the shape of χ4α/χmaxα(1Cα) does depend both on α and q. However, z(tw, t) is independent of q and α and therefore the (tw, t)-dependence of χ4 is uniquely specified via Cα, which is a function of z(tw, t).

FIG. 21.

To test the dependence on particle type α, this figure shows the parametric plot of χ4O/χmaxO as a function of χ4Si/χmaxSi, for q = 2.7 Å−1. The results for different waiting times collapse with each other.

FIG. 21.

To test the dependence on particle type α, this figure shows the parametric plot of χ4O/χmaxO as a function of χ4Si/χmaxSi, for q = 2.7 Å−1. The results for different waiting times collapse with each other.

Close modal
1.
K.
Binder
and
W.
Kob
,
Glassy Materials and Disordered Solids—An Introduction to Their Statistical Mechanics
(
World Scientific
,
Singapore
,
2005
).
2.
L.
Berthier
and
G.
Biroli
,
Rev. Mod. Phys.
83
,
587
(
2011
).
4.
G. L.
Hunter
and
E. R.
Weeks
,
Rep. Prog. Phys.
75
,
066501
(
2012
).
5.
Slow Relaxations and Nonequilibrium Dynamics in Condensed Matter
, edited by
J.-L.
Barrat
,
M.
Feigelman
,
J.
Kurchan
, and
J.
Dalibard
(
Springer
,
2003
).
6.
Dynamical Heterogeneities in Glasses, Colloids and Granular Media
, edited by
L.
Berthier
,
G.
Biroli
,
J.-P.
Bouchaud
,
L.
Cipelleti
, and
W.
van Saarloos
(
Oxford University Press
,
2011
).
8.
H.
Sillescu
,
J. Non-Cryst. Solids
243
,
81
(
1999
).
9.
H. E.
Castillo
,
C.
Chamon
,
L. F.
Cugliandolo
, and
M. P.
Kennett
,
Phys. Rev. Lett.
88
,
237201
(
2002
).
10.
H. E.
Castillo
,
C.
Chamon
,
L. F.
Cugliandolo
,
J. L.
Iguain
, and
M. P.
Kennett
,
Phys. Rev. B
68
,
134442
(
2003
).
11.
C.
Chamon
,
P.
Charbonneau
,
L. F.
Cugliandolo
,
D. R.
Reichman
, and
M.
Sellitto
,
J. Chem. Phys.
121
,
10120
(
2004
).
12.
H. E.
Castillo
and
A.
Parsaeian
,
Nat. Phys.
3
,
26
(
2007
).
13.
A.
Parsaeian
and
H. E.
Castillo
,
Phys. Rev. Lett.
102
,
055704
(
2009
).
14.
S.
Golde
,
M.
Franke
, and
H. J.
Schoepe
,
AIP Conf. Proc.
1518
,
304
(
2013
).
15.
L.
Berthier
,
G.
Biroli
,
J.-P.
Bouchaud
,
W.
Kob
,
K.
Miyazaki
, and
D. R.
Reichman
,
J. Chem. Phys.
126
,
184503
(
2007
).
16.
L.
Berthier
,
G.
Biroli
,
J.-P.
Bouchaud
,
W.
Kob
,
K.
Miyazaki
, and
D. R.
Reichman
,
J. Chem. Phys.
126
,
184504
(
2007
).
17.
L.
Berthier
,
Phys. Rev. E
76
,
011507
(
2007
).
18.
M.
Vogel
and
S. C.
Glotzer
,
Phys. Rev. E
70
,
061504
(
2004
).
19.
G.
Parisi
,
J. Phys. Chem. B
103
,
4128
(
1999
).
20.
A.
Parsaeian
and
H. E.
Castillo
,
Phys. Rev. E
78
,
060105(R)
(
2008
).
21.
A.
Parsaeian
and
H. E.
Castillo
, e-print arXiv:0811.3190v1 [cond-mat.dis-nn] (2008).
22.
H.
Oukris
and
N. E.
Israeloff
,
Nat. Phys.
6
,
135
(
2010
).
23.
C.
Maggi
,
R. D.
Leonardo
,
G.
Ruocco
, and
J. C.
Dyre
,
Phys. Rev. Lett.
109
,
097401
(
2012
).
24.
A.
Smessaert
and
J.
Rottler
,
Phys. Rev. E
88
,
022314
(
2013
).
25.
B. S.
Gupta
and
S. P.
Das
,
Phys. Rev. E
90
,
012137
(
2014
).
26.
C.
Chamon
,
M. P.
Kennett
,
H. E.
Castillo
, and
L. F.
Cugliandolo
,
Phys. Rev. Lett.
89
,
217201
(
2002
).
27.
C.
Chamon
and
L. F.
Cugliandolo
,
J. Stat. Mech.
2007
,
P07022
.
28.
H. E.
Castillo
,
Phys. Rev. B
78
,
214430
(
2008
).
29.
G. A.
Mavimbela
and
H. E.
Castillo
,
J. Stat. Mech.
2011
,
P05017
.
30.
K. E.
Avila
,
H. E.
Castillo
, and
A.
Parsaeian
,
Phys. Rev. Lett.
107
,
265702
(
2011
).
31.
K. E.
Avila
,
H. E.
Castillo
, and
A.
Parsaeian
,
Phys. Rev. E
88
,
042311
(
2013
).
32.
G. A.
Mavimbela
,
H. E.
Castillo
, and
A.
Parsaeian
, e-print arXiv:1210.1249v2 [cond-mat.dis-nn] (2013).
33.
B. W. H.
van Beest
,
G. J.
Kramer
, and
R. A.
van Santen
,
Phys. Rev. Lett.
64
,
1955
(
1990
).
34.
K.
Vollmayr
,
W.
Kob
, and
K.
Binder
,
Phys. Rev. B
54
,
15808
(
1996
).
35.
J.
Horbach
and
W.
Kob
,
Phys. Rev. B
60
,
3169
(
1999
).
36.
J.
Badro
,
D. M.
Teter
,
R. T.
Downs
,
P.
Gillet
,
R. J.
Hemley
, and
J.-L.
Barrat
,
Phys. Rev. B
56
,
5797
(
1997
).
37.
S. N.
Taraskin
and
S. R.
Elliott
,
Phys. Rev. B
59
,
8572
(
1999
).
38.
E.
Lascaris
,
M.
Hemmati
,
S. V.
Buldyrev
,
H. E.
Stanley
, and
C. A.
Angell
,
J. Chem. Phys.
142
,
104506
(
2015
).
39.
E.
Lascaris
,
M.
Hemmati
,
S. V.
Buldyrev
,
H. E.
Stanley
, and
C. A.
Angell
,
J. Chem. Phys.
140
,
224502
(
2014
).
40.
C.
Rajappa
,
S. B.
Sringeri
,
Y.
Subramanian
, and
J.
Gopalakrishnan
,
J. Chem. Phys.
140
,
244512
(
2014
).
41.
M. R.
Farrow
and
M. I. J.
Probert
,
J. Chem. Phys.
135
,
044508
(
2011
).
42.
I.
Saika-Voivod
,
F.
Sciortino
,
T.
Grande
, and
P. H.
Poole
,
Phys. Rev. E
70
,
061507
(
2004
).
43.
I.
Saika-Voivod
,
P. H.
Poole
, and
F.
Sciortino
,
Nature
412
,
514
(
2001
).
44.
J.
Badro
,
P.
Gillet
, and
J.-L.
Barrat
,
Europhys. Lett.
42
,
643
(
1998
).
45.
J.-L.
Barrat
,
J.
Badro
, and
P.
Gillet
,
Mol. Simul.
20
,
17
(
1997
).
46.
A.
Saksaengwijit
,
J.
Reinisch
, and
A.
Heuer
,
Phys. Rev. Lett.
93
,
235701
(
2004
).
47.
J.
Reinisch
and
A.
Heuer
,
Phys. Rev. Lett.
95
,
155502
(
2005
).
48.
A.
Saksaengwijit
and
A.
Heuer
,
Phys. Rev. E
73
,
061503
(
2006
).
49.
A.
Saksaengwijit
and
A.
Heuer
,
J. Phys.: Condens. Matter
19
,
205143
(
2007
).
50.
J.
Reinisch
and
A.
Heuer
,
J. Phys. Chem. B
110
,
19044
(
2006
).
51.
K.
Vollmayr-Lee
and
A.
Zippelius
,
Phys. Rev. E
88
,
052145
(
2013
).
52.
L. T.
Vinh
,
N. V.
Huy
, and
P. K.
Hung
,
Int. J. Mod. Phys. B
26
,
1250117
(
2012
).
53.
P.
Scheidler
,
W.
Kob
,
A.
Latz
,
J.
Horbach
, and
K.
Binder
,
Phys. Rev. B
63
,
104204
(
2001
).
54.
S. N.
Taraskin
and
S. R.
Elliott
,
Phys. Rev. B
56
,
8605
(
1997
).
55.
S. N.
Taraskin
and
S. R.
Elliott
,
Physica B
316
,
81
(
2002
).
56.
T.
Uchino
,
J. D.
Harrop
,
S. N.
Taraskin
, and
S. R.
Elliott
,
Phys. Rev. B
71
,
014202
(
2005
).
57.
F.
Leonforte
,
J. Non-Cryst. Solids
357
,
552
(
2011
).
58.
T.
Kawasaki
,
K.
Kim
, and
A.
Onuki
,
J. Chem. Phys.
140
,
184502
(
2014
).
59.
T.
Kawasaki
and
A.
Onuki
,
J. Chem. Phys.
138
,
12A514
(
2013
).
60.
M.
Vogel
and
S. C.
Glotzer
,
Phys. Rev. Lett.
92
,
255901
(
2004
).
61.
M. N. J.
Bergroth
,
M.
Vogel
, and
S. C.
Glotzer
,
J. Phys. Chem. B
109
,
6748
(
2005
).
62.
63.
P. K.
Hung
,
N. T. T.
Ha
, and
N. V.
Hong
,
Eur. Phys. J. E
36
,
60
(
2013
).
64.
J.
Helfferich
,
K.
Vollmayr-Lee
,
F.
Ziebert
,
H.
Meyer
, and
J.
Baschnagel
,
Europhys. Lett.
109
,
36004
(
2015
).
65.
66.
K.
Vollmayr-Lee
,
J. A.
Roman
, and
J.
Horbach
,
Phys. Rev. E
81
,
061203
(
2010
).
67.
K.
Vollmayr-Lee
,
R.
Bjorkquist
, and
L. M.
Chambers
,
Phys. Rev. Lett.
110
,
017801
(
2013
).
68.
This is not a complete list of BKS-simulations. For further work please see references therein.
69.

At high temperatures the self-diffusion constants follow a power law, as predicted by mode coupling theory, with critical temperature Tc = 3330 K. Below Tc the dynamics exhibits Arrhenius behavior.35 When the system is quenched with a sufficiently slow cooling rate a density maximum occurs at about 4800 K.34,82

70.
B.
Ruta
,
G.
Baldi
,
G.
Monaco
, and
Y.
Chushkin
,
J. Chem. Phys.
138
,
054508
(
2013
).
71.
F. A.
de Melo Marques
,
R.
Angelini
,
E.
Zaccarelli
,
B.
Farago
,
B.
Ruta
,
G.
Ruocco
, and
B.
Ruzicka
,
Soft Matter
11
,
466
(
2015
).
72.
A.
Saksaengwijit
and
A.
Heuer
,
Phys. Rev. E
74
,
051502
(
2006
).
73.
W.
Kob
and
J.-L.
Barrat
,
Eur. Phys. J. B
13
,
319
(
2000
).
74.

In the ensemble used in our simulations the numbers of both O and Si atoms are kept constant, and our definition of dynamic susceptibility describes fluctuations in that particular ensemble. The dynamic susceptibility in a different ensemble could have a different value.15 For example, a possible alternative would be an ensemble that allows particle number fluctuations, thus yielding an extra contribution that would increase the value of the dynamic susceptibility. We do not pursue such alternative in the present work.

75.
K. E.
Avila
,
H. E.
Castillo
,
A.
Fiege
,
K.
Vollmayr-Lee
, and
A.
Zippelius
,
Phys. Rev. Lett.
113
,
025701
(
2014
).
76.
J. C.
Dyre
,
J. Chem. Phys.
143
,
114507
(
2015
).
77.

To preserve scaling, the details of smoothing data need care. We used the smoothed data of Fig. 5 to determine χmaxα which is identical to the values of Fig. 6. We then determined for each tw the unsmoothed data χ4α/χmaxα as function of (1 − Cα). The thus obtained unsmoothed data were then smoothed via linear binning of (1 − Cα).

78.

The criterion for a particle j to belong to a sub-box Br requires the particle to be at a position rj(tw) within the sub-box at time tw. The particle need not belong to Br at the later time tw + t.

79.

A desired Cα value can be achieved only up to a certain accuracy because during any simulation run configurations are saved only at certain discrete times. Of the available configurations, we choose the one for which Cα is the closest to Cfix. Results are reported only for cases when Cα is within 1% of the chosen Cfix value.

80.

To be precise, for tw = 0.49 ns the lowest possible CcO is 0.165, which was chosen here. For all other cases shown in Fig. 17 the target value of CcO is achieved with 1% accuracy or better.

81.

In our analysis we only use simulation runs c for which there exist times such that Eq. (13) is satisfied within an accuracy of 1% or better.

82.
J. M. D.
Lane
,
Phys. Rev. E
92
,
012320
(
2015
).