content_image/413854.jpg
system was quenched to a low temperature T$_{f}$ < T$_{c}$ , fol- lowed by the production runs. Wahlen and Rieger [26] analyze time-dependent correlation functions at different waiting times t$_{w}$ and Berthier [25] and Scala et al. [27] study the generalized fluctuation dissipation relation and the energy landscape [27]. All three studies [25–27] in- vestigate the early stages of the aging dynamics, i.e. the dynamics was explored on time scales that were much smaller than the equilibrium relaxation time τ$_{eq}$ at the temperature T$_{f}$ . 

In this work, we also consider quenches in BKS sil- ica from a high temperature T$_{i}$ to a low temperature T$_{f}$ . Different from previous simulation studies, we aim at elu- cidating the full transient dynamics at T$_{f}$ from the initial state at t$_{w}$ = 0 to equilibrium. To this end, temperatures T$_{f}$ are chosen such that the system can be fully equili- brated on the typical time span of the MD simulation. Note that the considered temperatures T$_{f}$ ∈ { 2500 K, 2750 K, 3000 K, 3250 K } are below the critical MCT tem- peratures T$_{c}$ . Thus, we have access to the full aging dy- namics in the experimentally relevant Arrhenius temper- ature regime that we have mentioned above. 

The analysis of time-dependent density correlation functions C$_{q}$ ( t$_{w}$, t$_{w}$ + t ) (with q the wavenumber) and the mean square displacement Δ r $^{2}$( t$_{w}$, t$_{w}$ + t ) reveal three different regimes of waiting times t$_{w}$ : In the case of C$_{q}$ ( t$_{w}$, t$_{w}$ + t ) (and similarly for Δ r $^{2}$( t$_{w}$, t$_{w}$ + t )) at early t$_{w}$ , a rapid decay to zero is seen, without forming a plateau at intermediate times. Then, for larger values of t$_{w}$ a plateau is formed. The height of this plateau grows with waiting time and becomes more pronounced, before in the final regime the plateau height is indepen- dent of t$_{w}$ and T$_{i}$ . In the latter regime, time superposition holds, i.e. by scaling time with a decay time t Cq r the C$_{q}$ for the different values of t$_{w}$ fall onto a master curve at a given wavenumber q . This behavior is very similar to that found for the KALJ mixture. However, it is differ- ent from the behavior predicted by mean-field spin glass models and the activated dynamics scaling, as proposed by Wahlen and Rieger. Thus, these results suggest that the aging dynamics in silica, the prototype of a glass- former with a tetrahedral network structure, is very sim- ilar to that of simple glass-formers with a closed-packed hard-sphere-like structure. We find a difference between the KALJ mixture and SiO$_{2}$, however, in the parametric plot of C ′ $_{q}$( C$_{q}$ ). For SiO$_{2}$ C ′ $_{q}$( C$_{q}$ ) shows a data collapse for different sufficiently large t$_{w}$ and thus C$_{q}$ = C ( q, z ( t$_{w}$, t )) whereas this data collapse does not hold as well in the case of the KALJ mixture. 

The rest of the paper is organized as follows: In the next Sec. we give the details of the BKS potential and the simulation. Then, we present the results in Sec. III, before we summarize in Sec. IV. Appendix A describes the implementation of the Nos´e-Hoover thermostat used in our simulation. 

The interactions between the particles are modeled by the BKS potential [40] which has been used frequently and has proven to be reliable for the study of the dy- namics of amorphous silica [25–39] The functional form of the BKS potential is given by a sum of a Coulomb term, an exponential and a van der Waals term. Thus the potential between particles i and j , a distance r$_{ij}$ apart, is given by 

where e is the charge of an electron and the con- stants A$_{ij}$ , B$_{ij}$ and C$_{ij}$ are A$_{SiSi}$ = 0 . 0 eV, A$_{SiO}$ = 18003 . 7572eV, A$_{OO}$ = 1388 . 7730eV, B$_{SiSi}$ = 0 . 0 ˚ A − $^{1}$, B$_{SiO}$ = 4 . 87318 ˚ A − $^{1}$, B$_{OO}$ = 2 . 76000 ˚ A − $^{1}$, C$_{SiSi}$ = 0 . 0 eV ˚ A − $^{6}$, C$_{SiO}$ = 133 . 5381 eV ˚ A − 6 and C$_{OO}$ = 175 . 0000eV ˚ A − 6 [40]. The partial charges q$_{i}$ are q$_{Si}$ = 2 . 4 and q$_{O}$ = − 1 . 2 and e 2 is given by 1602 . 19 / (4 π 8 . 8542)eV ˚ A. 

The Coulombic part of the interaction was computed by using the Ewald method [41, 42] with a constant αL = 6 . 3452, where L is the size of the cubic box, and by using all q -vectors with | q | ≤ 6 · 2 π/L . [51]. We ensure that the Ewald term in real space is also differ- entiable at the cutoff by smoothing similarly to Eq. (3) in [43] with r$_{c}$ = 8 ˚ A and d = 0 . 05 ˚ A $^{2}$. To increase com- putation speed the non-Coulombic contribution to the potential was truncated, smoothed and shifted at a dis- tance of 5.5 ˚ A. Note that this truncation is not negligible since it affects the pressure of the system. In Ref. [43] further slight variations on the potential are described in detail [52]. In order to minimize surface effects pe- riodic boundary conditions were used. The masses of the Si and O atoms were 28.086 u and 15.9994 u, re- spectively. The number of particles was 336, of which 112 were silica atoms and 224 were oxygen atoms. For all simulation runs the size of the cubic box was fixed to L = 16 . 920468 ˚ A which corresponds to a density of ρ = 2 . 323 g/cm $^{3}$, a value that is very close to the one of real silica glass, ρ = 2 . 2 g/cm 3 [44]. 

We investigated the aging dynamics for systems which were quenched from a high temperature T$_{i}$ to a low tem- perature T$_{f}$ . To increase the statistics, for each ( T$_{i}$, T$_{f}$ ) 20 independent simulation runs were performed. To obtain 20 independent configurations we carried out molecular dynamics (MD) simulations using the velocity Verlet al- gorithm with a time step of 1 . 6 fs at 6000 K. The tem- perature was kept constant at 6000 K with a stochastic heat bath by replacing the velocities of all particles by new velocities drawn from the corresponding Boltzmann distribution every 150 time steps. Independent configura- tions were at least 3 . 27 ns apart. Each of these configura- tions undergoes the following sequence of simulation runs (see also Fig. 1). After fully equilibrating the samples at 

## II. MODEL AND DETAILS OF THE SIMULATION 

<img src='content_image/413854.jpg'>

