2 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. 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 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]. 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 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 II. MODEL AND DETAILS OF THE SIMULATION φ ( r$_{ij}$ ) = q$_{i}$q$_{j}$e 2 r$_{ij}$ + A$_{ij}$e − B$_{ij}$r$_{ij}$ − C$_{ij}$ r 6 ij , (1)