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 SiO_{2} 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 *t*_{w} 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*(*f*_{s,r}) of the local incoherent intermediate scattering function *f*_{s,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*(*f*_{s,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.

## I. INTRODUCTION

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 *t*_{w} — 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*(*f*_{s,r}) of the local incoherent intermediate scattering function *f*_{s,r}. In previous work, there have been relatively few studies on *P*(*f*_{s,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 SiO_{2} specifically see Refs. 15–18) There have been fewer studies for the aging dynamics, i.e., the dependence of *χ*_{4} on the waiting time *t*_{w}, which we discuss in this paper.^{13,19–25}

Predictions for the scaling of *P*(*f*_{s,r}) with respect to *t*_{w} 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 *f*_{s,r}(*t*_{w}, *t*_{w} + *t*), which depends on the position **r**, the waiting time *t*_{w}, and the time interval *t*. The Goldstone modes correspond to space dependent shifts of the time variable *t* → *ϕ*_{r}(*t*) such that

Here *C*(*t*_{w}, *t*_{w} + *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 *t*_{w} and the time interval *t* essentially only through the global two-time correlation function *C*(*t*_{w}, *t*_{w} + *t*).^{9–11} Thus it is expected that the probability distribution *P*(*f*_{s,r}(*t*_{w}, *t*_{w} + *t*)) should collapse for different waiting times *t*_{w}, for (*t*_{w}, *t*) pairs chosen such that *C*(*t*_{w}, *t*_{w} + *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 glasses^{12,13,20,21,25,30,31} and experimental results for a polymer glass^{22} find this predicted scaling to hold. The numerical simulations also show that *χ*_{4}(*t*_{w}, *t*_{w} + *t*) is a product of two factors: a waiting-time dependent scale that grows with *t*_{w} and a scaling function that depends on (*t*_{w}, *t*) only through the value of *C*(*t*_{w}, *t*_{w} + *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 SiO_{2} which is a strong glass former. The van Beest-Kramer-van Santen (BKS) potential^{33} which we use has not only been shown to be an excellent model for real silica^{34–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} structure^{34,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 SiO_{2} 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*(*f*_{s,r}) uncovers the importance of intermediate length scales and the scaling of *C*, *χ*_{4}, and *P*(*f*_{s,r}) all indicates a common aging clock which is the same for Si and O atoms.

## II. MODEL AND SIMULATION DETAILS

To model amorphous SiO_{2}, we used the BKS potential.^{33} We carried out molecular dynamics (MD) simulations with *N*_{Si} = 112 silica atoms and *N*_{O} = 224 oxygen atoms, at a constant volume $V=16.920\u2009468\xc53$ which corresponds to a density *ρ* = 2.323 g/cm^{3}. 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 *T*_{i} = 5000 K for 3.27 ns, followed by an instantaneous quench to lower temperature *T*_{f} = 2500 K, i.e., below *T _{c}* = 3330 K.

^{69}Unique to our simulations is that we applied the Nosé-Hoover temperature bath at

*T*

_{f}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

*T*

_{f}stays constant and is similar to

*T*

_{f}(

*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

*T*

_{f}.

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 *T*_{f} 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).

## III. GLOBAL INCOHERENT INTERMEDIATE SCATTERING FUNCTION

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

with

Here **r**_{j}(*t*) is the position of particle *j* at time *t*, *t*_{w} 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 $\u2026$ 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*(*f*_{s,r}).

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

Instead of the commonly used choice $Ccut\alpha =1/e$, we adjust to the varying plateau height of *C*^{α}(*q*) for different *q*. Hence for $Ccut\alpha (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*).

q
. | Si . | O . | All . |
---|---|---|---|

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 |

q
. | Si . | O . | All . |
---|---|---|---|

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 *t*_{w}. As in Ref. 66, different regimes for *t*_{w} can be identified from Figs. 1 and 2. In Ref. 66 it was found that for small waiting times, *t*_{w} ≲ 0.1 ns, *C*^{α}(*t*_{w}, *t*_{w} + *t*) does not form a plateau, for intermediate *t*_{w} a plateau is formed, and time superposition applies, and for sufficiently large waiting times *C*^{α} becomes *t*_{w}-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 *t*_{w} as a change in the *t*_{w} dependence of the relaxation times *τ*^{α}(*t*_{w}) from non-power law to power law behavior. Fig. 2 also shows that $\tau q\alpha (tw)$ does not reach a plateau, i.e., the waiting times are not long enough to reach equilibrium.

Let us next investigate further the dynamics for intermediate waiting times *t*_{w} ≳ 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 *τ*(*t*_{w}) have been found experimentally for a metallic glass^{70} 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 Heuer^{72} 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*^{α}(*q*_{2}) versus *C*^{α}(*q*_{1}) for various *t*_{w}. Whereas for their system, a binary Lennard-Jones system, they found no data collapse,^{73} for our system, SiO_{2}, it was found in Ref. 66 that data collapse indeed happens. This indicates that $C\alpha (tw,tw+t,q)=C\alpha (z\u0303(tw,t,\alpha ),q)$.^{66} In other words, for each particle type *α*, there is a unique *q*-independent aging clock represented by $z\u0303(tw,t,\alpha )$. 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 *t*_{w} ≥ 0.14 ns. When we include in a parametric plot of *C*^{α}(*q*_{2}) versus *C*^{α}(*q*_{1}) also too short waiting times *t*_{w} < 0.14 ns, we obtain no data collapse for SiO_{2} 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 SiO_{2}.

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\u0303(tw,t,\alpha )$ is independent of *α*. The question is therefore whether it is true that

To answer this question, we investigate directly the *α*-dependence via a parametric plot of *C*^{O}(*t*_{w}, *t*_{w} + *t*, *q*) versus *C*^{Si}(*t*_{w}, *t*_{w} + *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 (*t*_{w}, *t*)-dependence is solely governed by one function *z*(*t*_{w}, *t*), i.e., an “inner aging clock” which is not only *q*-independent but also the same for different particle types.

This may appear surprising at first, since, as mentioned before, it is known that in SiO_{2} 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 (*t*_{w}, *t*_{w} + *t*)-pairs for fixed *C*^{α}. For example, the blue circles of the bottom curve were obtained by finding for each *t*_{w} the corresponding *t*_{w} + *t* for which *C*^{Si}(*t*_{w}, *t*_{w} + *t*, *q* = 2.7 Å^{−1}) = 0.575 with 1% accuracy. For equilibrium dynamics, this curve would be trivial: it would be the set of (*t*_{w}, *t*_{w} + *t*)-pairs for a certain constant value of *t*. For aging dynamics, the curve is non-trivial: *t* changes as *t*_{w} changes. Nevertheless, we obtain the same non-trivial curve for O atoms (blue triangles) for fixed *C*^{O} = 0.463. Similarly we obtain identical curves for Si atoms (circles) and oxygen atoms (triangles) for different choices of *C*^{Si} and *C*^{O} (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 *C*^{O} being lower than the corresponding constant value of *C*^{Si} on the same curve.

## IV. DYNAMIC SUSCEPTIBILITY

In this section, we study the dynamic susceptibility, $\chi 4\alpha $, which is a four-point correlation function. $\chi 4\alpha $ quantifies thermal fluctuations of the incoherent intermediate scattering function. Following the notation of Berthier,^{17} $\chi 4\alpha $ is defined^{74} to be

where $fs\alpha (tw,tw+t,q)$, *N*_{α}, and $\u2026$ are as defined at the beginning of Sec. III. To obtain error bars for $\chi 4\alpha $, we divide the 200 independent simulation runs into 20 subsets each of 10 independent simulation runs and compute a value $\chi 4(\alpha ,i)$ for each subset *i*, with *i* = 1, …, 20. The error bar of $\chi 4\alpha $ is the standard deviation of the mean over those 20 independent $\chi 4(\alpha ,i)$ values. As in the case of *C*^{α}, we further smooth the results for $\chi 4\alpha $ 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\xaf,\chi 4\xaf,\Delta \chi 4\xaf)$.

Fig. 5 shows the resulting dynamic susceptibility for oxygen atoms and *q* = 2.7 Å^{−1}. Since $\chi 4\alpha $ 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 $\chi max\alpha $ at an intermediate time $tmax\alpha $. 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 $\chi 4\alpha $ on waiting time *t*_{w}. To quantify the dependence of this maximum on *t*_{w}, *q*, and *α*, we show the peak height $\chi max\alpha $ and the peak position $tmax\alpha $ as functions of *t*_{w} in Figs. 6 and 7, respectively. Similar to previous results in fragile glass formers,^{13,20,21,23–25} $\chi max\alpha $ and $tmax\alpha $ increase with increasing *t*_{w}. It is possible that $\chi max\alpha $ reaches a plateau for large *t*_{w}, but the noise in the results is too large to allow for any definite conclusions to be drawn. As above, the dependence of $tmax\alpha (tw,q)$ on *t*_{w} is consistent with the presence of two regimes: *t*_{w} ≲ 0.1 ns and *t*_{w} ≳ 0.1 ns. As in the case of $\tau q\alpha (tw)$, in the longer time regime, the time scale $tmax\alpha (tw,q)$ has a power law dependence on *t*_{w}, with exponents *μ*_{max} which are again independent of particle type, as shown in the inset of Fig. 7.

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 find^{13,20,21,25} a scaling behavior of the (*t*_{w}, *t*) dependence of the dynamic susceptibility given by

A similar more general result follows from Ref. 76. Without loss of generality, we choose $\chi 40(tw,q,\alpha )$ in Eq. (7) to be the maximum height $\chi max\alpha $. To test the validity of Eq. (7), we plot in Fig. 8 $\chi 4\alpha /\chi max\alpha $ as a function of $1\u2212C\alpha $ for the case of oxygen atoms (*α* = O) and *q* = 2.7 Å^{−1} (for details see endnote^{77}). 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

Fig. 9 shows the resulting (1 − *C*^{α})_{cross} as a function of *t*_{w} 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 *t*_{w} appears to be present, which may be evidence for the existence of small corrections to Eq. (7).

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

The dependence on time *t* enters into Eq. (9) only via *C* and therein only through *z*(*t*_{w}, *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

where $\chi \u02c6(z,q,\alpha )$ is a *t*_{w}-independent function in the sense that all (*t*, *t*_{w})-dependences enter only via *z*. As a consequence of *z* being independent of *q*, we obtain data collapse for a parametric plot of $\chi 4\alpha /\chi max\alpha $ for *q* = *q*_{2} as a function of $\chi 4\alpha /\chi max\alpha $ for *q* = *q*_{1}, with *q*_{2}≠*q*_{1}, and similarly as a consequence of *z* being independent of *α* we obtain data collapse for a parametric plot of $(\chi 4O/\chi maxO)$ versus $(\chi 4Si/\chi 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 $\chi 4all/\chi maxall\u2009(1\u2212Call)$ for different *t*_{w}, which is shown in Fig. 10 for *q* = 2.7 Å^{−1} and quantified for all investigated *q* values via (1 − *C*^{all})_{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}

## V. DISTRIBUTION OF LOCAL INCOHERENT INTERMEDIATE SCATTERING FUNCTION

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\alpha (tw,tw+t,q)$. In this section, we present results on the probability distribution for the *local* coarse grained intermediate scattering function

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

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

The theory therefore predicts that if (*t*_{w}, *t*) pairs are chosen such that *C*^{α}(*t*_{w}, *t*_{w} + *t*, *q*) is fixed, the corresponding $P(fs,r\alpha )$ should be *t*_{w} independent. This data collapse has been confirmed for spin glasses^{9,11} and for fragile structural glass formers.^{12,13,21,30} In this section, we investigate the scaling of $P(fs,r\alpha (tw,tw+t,q))$ for our SiO_{2} 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 *B*_{r}. 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 *B*_{r} is then an integer *b* times this very small sub-box, i.e., *L* × *b*/*M*. The average number of particles in *B*_{r} is therefore 〈*N*_{r}〉 = *N*_{α}(*b*/*M*)^{3}. We present in this paper results for 〈*N*_{r}〉 ≈ 5 and 〈*N*_{r}〉 ≈ 40 as listed in Table II. The distribution $P(fs,r\alpha (tw,tw+t,q))$ is then obtained for a specific set of (*t*_{w}, *t*, *q*, *α*) and for a specific simulation run via measurements of $fs,r\alpha $ for all *M*^{3} possible *B*_{r} and for all **q**-vectors of magnitude *q*. To obtain the *M*^{3} possible sub-boxes *B*_{r}, the sub-box is shifted in the three directions and periodic boundary conditions were used.

α
. | M
. | b
. | 〈N_{r}〉
. | M
. | b
. | 〈N_{r}〉
. |
---|---|---|---|---|---|---|

Si | 14 | 5 | 5.1 | 14 | 10 | 40.8 |

O | 7 | 2 | 5.22 | 7 | 4 | 41.8 |

All | 8 | 2 | 5.25 | 8 | 4 | 42.0 |

α
. | M
. | b
. | 〈N_{r}〉
. | M
. | b
. | 〈N_{r}〉
. |
---|---|---|---|---|---|---|

Si | 14 | 5 | 5.1 | 14 | 10 | 40.8 |

O | 7 | 2 | 5.22 | 7 | 4 | 41.8 |

All | 8 | 2 | 5.25 | 8 | 4 | 42.0 |

To test whether the (*t*_{w}, *t*) dependence of $P(fs,r\alpha )$ is governed by *C*^{α}(*t*_{w}, *t*_{w} + *t*, *q*), we use the same approach as in previous work.^{12,13,21,30} We choose a fixed value *C*_{fix} of the global *C*^{α}, and for each waiting time *t*_{w} we determine the time *t*_{fix} such that *C*^{α}(*t*_{w}, *t*_{w} + *t*_{fix}, *q*) = *C*_{fix} (see endnote^{79}). For these time pairs (*t*_{w}, *t*_{fix}) we then determine for each independent simulation run *c* = 1, 2, …, 200, a distribution $Pc(fs,r\alpha (tw,tw+tfix,q))$, using for Eq. (11) the positions $ri(tw)c$ and $ri(tw+tfix)c$. Please note that the parameters *t*_{w}, *t*_{w} + *t*_{fix}, and *q* are common to all simulation runs. For each bin of the distribution, we obtain the mean $Pc(fs,r\alpha (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 〈*N*_{r}〉 = 5.1.

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*(*f*_{s,r}), and nevertheless very good data collapse. To check quantitatively whether the slight *t*_{w} dependence in Fig. 11 is systematic, we focus on the location of most discrepancy, the maximum. Fig. 12 shows the maximum value $Pmax\alpha $ of the distribution, as a function of *t*_{w}. Consistent with the results above, we find that for small *t*_{w} scaling is not valid, but for *t*_{w} ≳ 0.1 ns $Pmax\alpha $ is approximately independent of *t*_{w}. We obtain similar results for *q* > 2.7 Å^{−1} and also for *α* = Si.

To probe the dependence of $P(fs,r\alpha )$ on particle type *α*, we show in the inset of Fig. 13 the comparison of $P(fs,r\alpha )$ for *α* = Si, O, all. We conclude that $P(fs,r\alpha )$ does depend on *α*. In Secs. III and IV we had found that despite different dynamics of silicon and oxygen atoms, their scaling *z*(*t*_{w}, *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\alpha )$, we therefore show in Fig. 13 $P(fs,rall)$ and find indeed a data collapse for different *t*_{w}.

So far we have shown $P(fs,r\alpha )$ for various fixed *C*^{α}, for *q* = 2.7 Å^{−1}, 〈*N*_{r}〉 ≈ 5, with either *α* = O or *α* = all. When also *q* and 〈*N*_{r}〉 are varied, we find that the predicted data collapse occurs as long as *t*_{w} ≳ 0.1 ns, *q* ≥ 2.7 Å^{−1}, and 〈*N*_{r}〉 ≈ 5. Next we look at cases when scaling fails. Fig. 14 shows $P(fs,r\alpha )$ for 〈*N*_{r}〉 ≈ 5 as before, but now *q* = 1.7 Å^{−1}. For *C*^{O} = 0.3 and *C*^{O} = 0.5 we find that the data collapse for different *t*_{w} is much worse than before. This is quantified in the inset of Fig. 14, which shows the systematic *t*_{w} dependence of $PmaxO$. The inset also shows that in the Gaussian case of *C*^{O} = 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 〈*N*_{r}〉 ≈ 42. A similarly worse data collapse for larger 〈*N*_{r}〉 has been found for a fragile glass former (Fig. 2 of Ref. 13).

Thus scaling breaks down for intermediate timescales, corresponding to intermediate *C*^{α}, when the regions probed by the local incoherent intermediate scattering function *f*_{s,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 *f*_{s,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*(*f*_{s,r}) grew with *t*_{w} at constant *C*^{α}. In that case it was argued that *P*(*f*_{s,r}) was narrowed by averaging of *f*_{s,r} over more than one correlated region, but this narrowing became weaker as the size of the correlated regions grew with *t*_{w}. 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 *t*_{w} for *C ^{O}* = 0.3, but they

*narrow*with increasing

*t*

_{w}for

*C*= 0.5. Therefore, for SiO

^{O}_{2}, 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 〈*N*_{r}〉 ≈ 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., *t*_{w} and *t*_{w} + *t*. We illustrate in Fig. 16 how this is done. First we choose a unique *t* = *t*_{fix} — the same for all 200 independent simulation runs — by demanding that the global incoherent scattering function *C*^{α}(*t*_{w}, *t*_{w} + *t*) (thick line in Fig. 16) take a certain value at *t* = *t*_{fix}, such as *C*^{α}(*t*_{w}, *t*_{w} + *t*_{fix}) = 0.5.

If we now look at individual runs, the intermediate scattering function computed for each run *c* is $Cc\alpha =\u3008fs\u3009c$ 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\alpha $ to be self-averaging: in Fig. 16, the values of $Cc\alpha (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*^{α}(*t*_{w}, *t*_{w} + *t*_{fix}). In other words, choosing a unique value *t* = *t*_{fix} is equivalent to choosing very different values for $Cc\alpha $ for each simulation run *c*. Consequently, the distributions *P*(*f*_{s,r}(*t*_{w}, *t*_{w} + *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 *t*_{w} = 0.49 ns. In the same figure, results are shown for another waiting time, *t*_{w} = 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 *t*_{w} = 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\alpha )$, the key variable is $Cc\alpha (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 *t*_{w} = 33.0 ns and *c* = 57 for *t*_{w} = 0.49 ns, chosen so that in both cases $CcO(tfixO)\u22480.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)\u22480.1,0.5,0.7$ (for more details see endnote^{80}).

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

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

The resulting average distribution $P(fs,r\alpha )$ is shown in Fig. 18 for *q* = 1.7 Å^{−1} and 〈*N*_{r}〉 = 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\alpha $ instead of *C*^{α}. To test the limits of this improved scaling, we show in Fig. 19 the distributions *P*(*f*_{s,r}) for *q* = 1.7 Å^{−1}, 〈*N*_{r}〉 = 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\alpha $ is the appropriate scaling quantity.

## VI. CONCLUSIONS

In summary, we have performed molecular dynamics simulations of the strong glass former SiO_{2} to investigate the scaling of dynamical heterogeneities in this system. We have quenched the system from an initial high temperature *T _{i}* to a final temperature

*T*below the mode-coupling critical temperature

_{f}*T*and observed the out of equilibrium dynamics as a function of the waiting time

_{c}*t*

_{w}, the time elapsed since the temperature quench. In particular, we have investigated the global incoherent intermediate scattering function

*C*

^{α}(

*t*

_{w},

*t*

_{w}+

*t*,

*q*), the dynamic susceptibility $\chi 4\alpha (tw,tw+t,q)$, and the distribution $P(fs,r\alpha (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 *t*_{w}, and when probing small enough regions in the system, the dependence on (*t*_{w}, *t*) of $\chi 4\alpha (tw,tw+t,q)$ and of $P(fs,r\alpha (tw,tw+t,q))$ is governed by *C*^{α}(*t*_{w}, *t*_{w} + *t*, *q*), up to a *t*_{w}-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 SiO_{2} and the comparison fragile glass formers is rather surprising, since SiO_{2} 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*^{α}(*t*_{w}, *t*_{w} + *t*, *q*) = *C*(*z*(*t*_{w}, *t*), *q*, *α*), where *z*(*t*_{w}, *t*) is independent of *q* and *α*. Thus *z*(*t*_{w}, *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 $\chi 4\alpha $ and of $P(fs,r\alpha )$ is controlled by *C*^{α}, it follows that the aging of $\chi 4\alpha $ and of $P(fs,r\alpha )$ 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*(*f*_{s,r}) discussed above no longer holds in its initial form. It is clear that the probability distributions *P _{c}*(

*f*

_{s,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\alpha $. Since the fluctuations in $Cc\alpha $ are given by $\chi 4\alpha $ and since $\chi max\alpha $ is increasing with decreasing wave vector, we find a particularly large spread of distributions

*P*(

_{c}*f*

_{s,r}) for the smallest wave vector. The average over

*P*(

_{c}*f*

_{s,r}) is not independent of

*t*

_{w}. This is equivalent to the statement that, in a very large system, a measurement of

*P*(

*f*

_{s,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

*t*

_{w}if instead of averaging probability distributions

*P*(

_{c}*f*

_{s,r}) from different runs (or mesoscopic regions) at constant time interval

*t*=

*t*

_{fix}, one averages probability distributions taken at constant

*mesoscopic*intermediate scattering function $Cmeso\alpha $. We speculate that this importance of $Cmeso\alpha $ 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\alpha $ 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.

## Acknowledgments

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.

### APPENDIX: DYNAMIC SUSCEPTIBILITY DEPENDENCE ON WAVE NUMBER AND ON PARTICLE TYPE

From Eq. (10) it follows that all the (*t*, *t*_{w})-dependence of $\chi \u02c6$ enters only via *z*. Since *z* is independent of *q*, it follows that a parametric plot of $\chi 4\alpha /\chi max\alpha $ for *q* = *q*_{2} as a function of $\chi 4\alpha /\chi max\alpha $ for *q* = *q*_{1}, with *q*_{2}≠*q*_{1}, should show data collapse. In Fig. 20, we show a parametric plot of this kind, for *q*_{2} = 2.7 Å^{−1}, *q*_{1} = 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 *q*_{2}, *q*_{1}, and *α* (not shown).

For the dependence on particle type, *α*, Eq. (10) predicts that there should also be data collapse in a parametric plot of $(\chi 4O/\chi maxO)$ versus $(\chi 4Si/\chi 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 $\chi 4\alpha /\chi max\alpha \u2009(1\u2212C\alpha )$ does depend both on *α* and *q*. However, *z*(*t*_{w}, *t*) is independent of *q* and *α* and therefore the (*t*_{w}, *t*)-dependence of *χ*_{4} is uniquely specified via *C*^{α}, which is a function of *z*(*t*_{w}, *t*).

## REFERENCES

At high temperatures the self-diffusion constants follow a power law, as predicted by mode coupling theory, with critical temperature *T _{c}* = 3330 K. Below

*T*the dynamics exhibits Arrhenius behavior.

_{c}^{35}When the system is quenched with a sufficiently slow cooling rate a density maximum occurs at about 4800 K.

^{34,82}

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.

To preserve scaling, the details of smoothing data need care. We used the smoothed data of Fig. 5 to determine $\chi max\alpha $ which is identical to the values of Fig. 6. We then determined for each *t*_{w} the unsmoothed data $\chi 4\alpha /\chi max\alpha $ as function of (1 − *C*^{α}). The thus obtained unsmoothed data were then smoothed via linear binning of (1 − *C*^{α}).

The criterion for a particle *j* to belong to a sub-box *B*_{r} requires the particle to be at a position **r**_{j}(*t*_{w}) within the sub-box at time *t*_{w}. The particle need not belong to *B*_{r} at the later time *t*_{w} + *t*.

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 *C*_{fix}. Results are reported only for cases when *C*^{α} is within 1% of the chosen *C*_{fix} value.

To be precise, for *t*_{w} = 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.

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.