Effect of instantaneous and continuous quenches on the density of vibrational modes in model glasses

TL;DR

This study examines how the method of quenching affects the vibrational mode density in model glasses, revealing a transition from a $eta ext{≈}3$ to a universal $eta ext{= }4$ as quench rates slow down.

Contribution

It demonstrates that the vibrational mode density exponent depends on quench rate and parent temperature, connecting instantaneous quenches to realistic slow cooling.

Findings

Instantaneous quenches from high T yield $eta ext{≈}3$.

As parent temperature approaches glass transition, $eta$ approaches 4.

Slower, physically realistic quenches produce a universal $eta=4$.

Abstract

Computational studies of supercooled liquids often focus on various analyses of their "underlying inherent states" --- the glassy configurations at zero temperature obtained by an infinitely-fast (instantaneous) quench from equilibrium supercooled states. Similar protocols are also regularly employed in investigations of the unjamming transition at which the rigidity of decompressed soft-sphere packings is lost. Here we investigate the statistics and localization properties of low-frequency vibrational modes of glassy configurations obtained by such instantaneous quenches. We show that the density of vibrational modes grows as $\omega^\beta$ with $\beta$ depending on the parent temperature $T_{0}$ from which the glassy configurations were instantaneously quenched. For quenches from high temperature liquid states we find $\beta\!\approx\!3$, whereas $\beta$ appears to approach the previously-observed value $\beta\!=\!4$ as $T_0$ approaches the glass transition temperature. We discuss the consistency of our findings with the theoretical framework of the Soft Potential Model, and contrast them with similar measurements performed on configurations obtained by continuous quenches at finite cooling rates. Our results suggest that any physical quench at rates sufficiently slower than the inverse vibrational timescale --- including all physically-realistic quenching rates of molecular or atomistic glasses --- would result in a glass whose density of vibrational modes is universally characterized by $\beta\!=\!4$.

Equations

$$\varphi(r_{ij})=\left\{\begin{array}[]{ccc}\varepsilon\left[\left({\textstyle\frac{\lambda_{ij}}{r_{ij}}}\right)^{n}+\sum\limits_{\ell=0}^{q}c_{2\ell}\left({\textstyle\frac{r_{ij}}{\lambda_{ij}}}\right)^{2\ell}\right]&,&{\textstyle\frac{r_{ij}}{\lambda_{ij}}}\leq x_{c}\\ 0&,&{\textstyle\frac{r_{ij}}{\lambda_{ij}}}>x_{c}\end{array}\right.,$$

$$c_{2\ell }=\frac{(-1{)}^{\ell +1}}{(2q-2\ell )!!(2\ell )!!}\frac{(n+2q)!!}{(n-2)!!(n+2\ell )}{x}_c^{-(n+2\ell )}.$$

$$C_{\text{mbox}Ber}\equiv \sqrt{1+\frac{\delta t}{{\tau }_{\text{mbox}Ber}}\frac{T-\tilde{T}(t)}{\tilde{T}(t)}},$$

$$\frac{\tilde{T}(t)-T(t)}{T(t)}\approx \frac{2{\tau }_{\text{mbox}Ber}\dot{T}}{T(t)}.$$

$$n_g\sim L^{\text{mathchar}22d}{\int }_0^{L^{-1}}D(\omega )d\omega \sim L^{\text{mathchar}22d-\beta -1},$$

Full text

Introduction.– Instantaneous quenches of high-temperature configurations into their so-called “underlying inherent states” are a prevalent practice in computational studies of disordered materials Schröder et al. (2000); Heuer (2008); Karmakar et al. (2012); Widmer-Cooper et al. (2008); Hentschel et al. (2011); O’Hern et al. (2003). One conspicuous example of this methodology is found in the large body of numerical work dedicated to the unjamming scenario (see e.g. Liu and Nagel (2010) and references therein), in which simple models of soft repulsive spheres are regularly employed. In these studies, packings of soft spheres at zero temperature are conventionally generated by instantaneous quenches from some random, high-energy disordered states, and later subjected to various structural analyses Wyart et al. (2005); van Hecke (2010); Goodrich et al. (2014); Mizuno et al. (2017) and/or perturbations Lerner et al. (2014); Xu et al. (2017).

A similar methodology is also extensively utilized in computational investigations of the glass transition Schröder et al. (2000); Heuer (2008); Karmakar et al. (2012); Widmer-Cooper et al. (2008), whose structural origin remains a highly-debated topic in condensed-matter physics Debenedetti and Stillinger (2001); Cavagna (2009). Instantaneous quenches that map an equilibrium configuration to a zero-temperature glassy state were first put forward by Stillinger and Weber Stillinger and Weber (1983), and subsequently utilized by many others Schröder et al. (2000); Heuer (2008); Karmakar et al. (2012); Widmer-Cooper et al. (2008), with the general assumption that the structural properties of the inherent states are indicative in some quantitative way of the dynamics of the supercooled configurations from which they are mapped.

A clear advantage of analyzing the structure of glassy inherent states over equilibrium configurations is the ability to cleanly and quickly extract structural observables while avoiding the difficulties that stem from thermal-fluctuations-induced noise, and from the broad spectrum of relaxation times that characterizes these systems. Instantaneous quenches are considered to be unrealistic idealizations of the physical cooling process by which glasses are formed. However, it is regularly assumed that generic properties of the resulting glasses remain unaffected by such protocols. This uncontrolled assumption overlooks the potential physical artifacts involved in performing instantaneous quenches in computational studies of the structural properties of glassy materials.

In this Rapid Communication we question the common practice of investigating glassy states that were instantaneously-quenched from high-temperature configurations, and subsequently deducing conclusions about generic glasses formed via physically-realistic protocols. We focus in particular on the statistical and structural properties of low-frequency vibrational modes measured in ensembles of inherent states created by an instantaneous quench of configurations equilibrated at various parent temperatures $T_{0}$. Recent studies of several structural glass forming models Lerner et al. (2016); Gartner and Lerner (2016); Mizuno et al. (2017) (see additional comments about the relation between Lerner et al. (2016) and the present work in Lerner et al. (2017)) identified a population of quasilocalized low-frequency glassy vibrational modes whose density $D(\omega)$ grows from vanishing frequencies $\omega\!\to\!0$ as $D(\omega)\!\sim\!\omega^{4}$. These modes are either measured below the lowest Goldstone modes’ frequency Lerner et al. (2016); Gartner and Lerner (2016), or identified by classifying vibrational modes according to their degree of localization Mizuno et al. (2017). Similar findings for a three-dimensional Heisenberg spin glass in a random field were put forward in Baity-Jesi et al. (2015). Here we show that the functional form of the density of low-frequency vibrational modes can be affected by instantaneous quenches, and, under some conditions, displays deviations from the $\omega^{4}$ law. We contrast our findings with measurements performed in ensembles of inherent states created by a continuous quench at a cooling rate $\dot{T}$ through our model systems’ glass transition temperature. We find that the $\omega^{4}$ law is robust to very rapid but not overdamped quenches, suggesting that inertia plays a key role in the self-organizational processes that occur as systems tumble down the multi-dimensional potential energy landscape during their quench into a glassy solid.

Models and methods.— Here we briefly review the models and methods used in this work; a detailed description of our model, methods and preparation protocols can be found in the Supplemental Material (SM) SM . We employ a binary mixture of point-like particles in three dimensions (3D) that interact via a purely repulsive inverse-power-law potential. In what follows physical observables (temperatures, frequencies, lengths, etc.) are understood as expressed in terms of the relevant microscopic units as defined in the SM. For visualization purposes alone we rescale frequency axes by a scale $\omega_{0}$, see figure captions. We chose to simulate systems of $N\!=\!2000$ particles for which the linear size of the box is slightly larger than the localization length of low-frequency glassy modes (estimated in our model at about 10 particle diameters Lerner et al. (2016)), but still small enough such that Goldstone modes are sufficiently suppressed, allowing ample exposure of vibrational modes that occur below the lowest Goldstone mode frequency Lerner et al. (2016).

Ensembles of inherent states were created by collecting a large number of independent equilibrium configurations from each parent temperature $T_{0}$, and evolving each one of these forward in time under fully overdamped dynamics $\dot{\vec{x}}\!\propto\!-\frac{\partial U}{\partial\vec{x}}$ until convergence, where $\vec{x}$ denotes particles’ coordinates and $U$ the potential energy. We have also created ensembles of continuously-quenched glasses, starting from independent equilibrium configurations at $T\!=\!1.00$, followed by a quench at a prescribed quench rate $\dot{T}$. Each of the constructed ensembles consists of 10,000 glassy samples, which ensures statistical convergence, see SM for further details.

Results.— Our model system exhibits the conventional phenomenology of computer glass forming models. In Fig. 1 we demonstrate the slowing down in the relaxational dynamics upon supercooling of our model by monitoring the stress autocorrelation function $c(t)\!\equiv\!N\langle\sigma(t)\sigma(0)\rangle$ measured at various equilibrium runs at temperatures $T$. Here $\sigma\!\equiv\!\frac{1}{V}\frac{\partial U}{\partial\gamma}$, $V$ is the volume of the simulation cell, $U$ is the potential energy and $\gamma$ is a simple shear strain. The angular brackets denote an average over the time-translationally-invariant signals of the stress from our equilibrium simulations. The inset of Fig. 1 shows the relaxation time $\tau_{\alpha}$ vs. $1/T$; relaxation times are estimated via $c(\tau_{\alpha})\!=\!1$, as indicated by the horizontal dashed line. The computer glass transition temperature of our model is estimated at $T_{g}\!\approx\!0.5$, where the relaxation time $\tau_{\alpha}(T_{g})\!\approx\!10^{5}$.

We next turn to the investigation of the statistics and properties of vibrational modes in the different ensembles of instantaneously quenched glasses. Each such ensemble was obtained by an instantaneous quench of independent configurations that were equilibrated at some parent temperature $T_{0}$. In Fig. 2 we show the low-frequency tails of the density of vibrational modes $D(\omega)$ measured in all the ensembles of glassy samples that were instantaneously quenched from parent temperatures as indicated by the figure caption. We find that $D(\omega)\!\sim\!\omega^{\beta}$ with $3<\beta\!\leq\!4$ for all ensembles, and $\beta\!\to\!4$ as $T_{0}\!\to\!T_{g}$. These data demonstrate that it is not only that the high $T_{0}$ inherent states possess more soft glassy vibrational modes, but that the actual functional form of the vibrational modes’ distribution function depends explicitly on $T_{0}$, at least up to the vicinity of the accessible equilibrium temperatures using conventional simulation methods. We emphasize at this point that our goal is not to accurately estimate the precise numerical value of the scaling exponents that characterize the density of vibrational modes. Our aim is rather to identify trends in the observed exponents upon systematically varying the preparation protocol of the glassy samples.

In Fig. 3 we show the means and the 2nd-9th deciles of the participation ratio $e$ binned over frequency $\omega$. The participation ratio of a vibrational mode $\hat{\Psi}$, defined as e\!\equiv\!\big{(}N\sum_{i}(\hat{\Psi}_{i}\!\cdot\!\hat{\Psi}_{i})^{2}\big{)}^{-1}, is a simple measure of the degree of localization of a mode: the more localized a mode is, the smaller its participation ratio is expected to be. In Lerner et al. (2016); Gartner and Lerner (2016) it has been shown that the participation ratio of low-frequency glassy modes scales as $N^{-1}$, indicating that they are quasilocalized foo . Fig. 3 here shows that the degree of localization of low-frequency glassy modes increases for deeper supercooling, consistently with the findings of Lerner et al. (2016) that show a decrease in the participation ratio of low-frequency glassy modes for slower cooling rates. The data indicate that the transition of the mean participation ratio from the Goldstone modes’ value to the low-frequency plateau (shown clearly for a much larger data set in Gartner and Lerner (2016)) is faster in ensembles created by an instantaneous quench from deeply supercooled solids. Furthermore, the rapid crossover in the localization properties of modes with increasing frequencies suggests that the exponent $\beta$ can only be read off $D(\omega)$ below frequencies that are roughly a third of the lowest Goldstone mode frequency, see e.g. the data in Fig. 4 below. We note that the crossover from quasilocalized, glassy modes at low frequencies to the first Goldstone modes is broader for smaller $N$; this can be seen, for instance, in Fig. 1 of Ref. Gartner and Lerner (2016). In the SM we show that increasing the system size does not, however, appear to have a substantial effect on our results, which reinforces the statement that the crossover broadening in our systems of $N\!=\!2000$ does not effect our conclusions.

It is natural to contrast our results for instantaneously quenched glasses with similar measurements in glassy samples formed by a continuous quench into solids at finite quench rates. In Fig. 4 we show the density of vibrational modes of systems quenched at rates $\dot{T}$ as described in the legend. Each such quench was preformed from initial equilibrium configurations at temperature $T\!=\!1.00$ (see SM for details). We find $\beta\!=\!4$ at rates $\dot{T}\!<\!10^{-2}$, as shown by the continuous lines. For higher rates, $\beta$ appears to decrease, and for an infinitely-fast quench from $T\!=\!1.00$ we find $\beta\!\approx\!3.3$ as indicated by the dashed line.

In Fig. 5 we plot the average potential energy per particle of the ensembles of instantaneously-quenched and continuously-quenched glassy samples. Interestingly, we find that the mean energy per particle of glasses quenched at the highest continuous rate for which $\beta\!=\!4$ is observed ($\dot{T}\!=\!10^{-3}$, see Fig. 4) is the same as for instantaneously quenched samples from the parent temperature $T_{0}\!=\!0.60$, up to less than a percent. However, in the latter ensemble we clearly find $\beta\!<\!4$, see dash-dotted line in Fig. 2. This observation of two ensembles with the same inherent state energies but different $\beta$ indicates that inertia that is present during the continuous quenches, but absent in the instantaneous quenches, plays an important role in the self-organizational processes that determine the fine details of the microstructure of the resulting glasses.

In order to explore the implications of our results for realistic glasses, we cast our reported observables into dimensionless numbers. We start with forming an atomistic timescale by considering the shear-wave speed $c_{s}\!\equiv\!\sqrt{\mu/\rho}\!\approx\!4$, with an athermal shear modulus Lutsko (1989) $\mu\!\approx\!15$ and a mass density $\rho\!=\!0.82$, and dividing it by an atomistic length $a_{0}\!\approx\!1.0$ to find $c_{s}/a_{0}\!\approx\!4$. We next take $T_{g}\!\approx\!0.5$ as a characteristic temperature scale, such that a dimensionless quench rate is formed as $\frac{\dot{T}a_{0}}{T_{g}c_{s}}$. Our findings suggest that for dimensionless quench rates lower than a crossover value $10^{-3}$, or alternatively, for quench rates $\dot{T}\!\lesssim\!10^{-3}T_{g}c_{s}/a_{0}$, the density of vibrational modes exhibits the $\omega^{4}$ law. To compare to physical glasses, e.g. metallic glasses, we take $T_{g}\!\approx\!500$K, $c_{s}\!\approx\!10^{3}\mbox{m/sec}$, and $a_{0}\!\approx\!10^{-9}\mbox{m}$ Wang et al. (2004), from which we conclude that glasses quenched at rates $\dot{T}\!\lesssim\!10^{11}\mbox{K/sec}$ would exhibit the $\omega^{4}$ law. The fastest rates that these materials can be quenched are typically on the order of $10^{7}\mbox{K/sec}$ (for quasi-2D ribbons; for bulk glasses the fastest cooling rates are slower), some 4 decades slower than our estimated crossover rate. This comparison essentially implies that any laboratory glass formed by quenching a melt would follow the $\omega^{4}$ law.

Summary and discussion.– In this Rapid Communication we have shown that the low-frequency tails of the density of vibrational modes of computer glasses created by an instantaneous quench have qualitatively different features compared to glasses created by a continuous quench. Our results suggest that the presence of inertia is important for the structural relaxation that occurs during quenches which leads to more stable glassy structures with less low-frequency vibrational modes. This suggestion is consistent with the results of Salerno et al. Salerno et al. (2012), who showed that upon reducing the inertia in the microscopic dynamics of sheared model glasses, the nature of avalanches of plastic activity, which depends in turn on the abundance of soft glassy modes, can change dramatically. Similar findings were reported in Karimi et al. (2017). If indeed the presence of inertia in the microscopic dynamics is key in determining the low-frequency spectra, it would be of interest to observe whether overdamped glasses such as emulsions or foams, or computer glasses generated in simulations that employ Brownian dynamics, exhibit observable qualitative differences in their spectra compared to their inertial counterparts.

Our results call for caution when attempts are made to establish general conclusions about glassy solids from studies of model glasses that are created by instantaneous quenches from high temperature liquid states. For instance, it is common practice in studies of the unjamming point to create packings of soft spheres by instantaneous quenches. While the qualitative features of the scaling of most mechanical observables with respect to the distance to the unjamming point do not seem to depend on the protocol with which packings are generated, our results suggests that the density of vibrational modes of those packings might not be representative of the spectra of glasses created by physical quenches.

It is interesting to attempt to relate our findings to the predictions of the Soft Potential Model Buchenau et al. (1991); Ilyin et al. (1987); Karpov and Parshin (1985); Gurarie and Chalker (2003). This theoretical framework assumes that a glass can be decomposed into small subsystems, each possessing a quasilocalized soft glassy mode. Focusing on such a typical subsystem, and assuming that particles are displaced a distance $s$ along the soft mode associated with that subsystem, this framework suggests that if the energy in the vicinity of $s\!=\!0$ satisfies $U(s)\!\geq\!U(0)$, then $D(\omega)$ is expected to grow as $\omega^{4}$. However, relaxing this constraint results in a different prediction, namely that $D(\omega)\!\sim\!\omega^{3}$ Gurarie and Chalker (2003). The condition that the energy only grows in the vicinity of $s\!=\!0$ can be viewed as a stability condition; in instantaneously quenched glasses the overdamped nature of the quench makes it possible to form barely-stable glasses that would possess local soft potentials $U(s)$ that have deeper minima at $s\!\neq\!0$ compared to $U(0)$, e.g. asymmetric double well potentials Karpov and Parshin (1985). Creating such unstable structures in slowly quenched glasses is much less likely. According to the discussed framework, one may hypothesize that $\beta\!=\!3$ should be observed in glasses created by an instantaneous quench; we indeed find $\beta$ very close to 3 in samples that were instantaneously quenched from very high temperatures, see Fig. 2.

In this work we followed the simple approach of Lerner et al. (2016) and investigated the density of vibrational modes in small, three dimensional model glasses, in which Goldstone modes are sufficiently suppressed to expose a population of quasilocalized soft glassy vibrational modes. This approach is, however, still limited in terms of the range of soft glassy modes’ frequencies that can be probed, due to hybridizations with extended Goldstone modes at higher frequencies, as can be seen in Fig. 3. It is therefore of interest to investigate these issues using frameworks that overcome the issue of hybridization with Goldstone modes e.g. Gartner and Lerner (2016); Zylberg et al. (2017); Brito et al. (2013); Baity-Jesi et al. (2015); Wijtmans and Manning (2017), allowing one to probe the density of quasilocalized excitations up to higher frequencies.

A key question to be addressed in future research is whether extremely slow quench rates can result in glasses with $\beta\!>\!4$. Recent developments Berthier et al. (2016); Ninarello et al. (2017) in the computational research of structural glasses allow one to equilibrate a particular model glass well below what is possible using conventional molecular dynamics or Monte Carlo methods. The new methodology introduced in Berthier et al. (2016); Ninarello et al. (2017) will be certainly useful in addressing this question.

Acknowledgments.– We thank Gustavo Düring, Eric DeGiuli, and Matthieu Wyart for fruitful discussions.

I Models and numerical methods

I.1 Model definitions

We employed a 50:50 binary mixture of ‘large’ and ‘small’ particles of equal mass $m$ enclosed in a cubic three dimensional box of linear size $L$, interacting via a radially-symmetric purely repulsive inverse power-law pairwise potential of the form

[TABLE]

where $r_{ij}$ is the distance between the $i^{\mbox{\tiny th}}$ and $j^{\mbox{\tiny th}}$ particles, $\varepsilon$ is a microscopic energy scale, and $x_{c}$ is the dimensionless distance for which $\varphi$ vanishes continuously up to $q$ derivatives. Distances are measured in terms of the interaction lengthscale $\lambda$ between two ‘small’ particles, and the rest are chosen to be $\lambda_{ij}\!=\!1.18\lambda$ for one ‘small’ and one ‘large’ particle, and $\lambda_{ij}\!=\!1.4\lambda$ for two ‘large’ particles. The coefficients $c_{2\ell}$ are given by

[TABLE]

We chose the parameters $x_{c}\!=\!1.48,n\!=\!10$, and $q\!=\!3$. The density was set to be $N/L^{3}\!=\!0.82\lambda^{-3}$, with the total number of particles $N\!=\!2000$ used in the majority of the numerical simulations performed (see additional discussion about system size effects below). Time is expressed in terms of $\tau_{0}\!\equiv\!\sqrt{m\lambda^{2}}/\varepsilon$, temperature in terms of $\varepsilon/k_{B}$ with $k_{B}$ the Boltzmann constant, quench rates in terms of $\varepsilon/(k_{B}\tau_{0})$, stresses in terms of $\varepsilon/\lambda^{3}$, and vibrational frequencies in terms of $\tau_{0}^{-1}$. In plots of the density of vibrational modes here and in the main text, we rescaled the $x$-axis by an arbitrary scale $\omega_{0}$ as reported in the figure caption, for visualization purposes.

I.2 Thermostating

Temperature was controlled using a Berendsen thermostat scheme berendsen1 , which amounts to multiplying the momentum vector of each particle at every integration step by a factor $C_{\mbox{\tiny Ber}}$ calculated as

[TABLE]

where $\tilde{T}(t)\!\equiv\!\frac{m}{3N}\sum_{i}v_{i}^{2}(t)$ is the instantaneous temperature with $v_{i}$ the magnitude of the velocity of the $i^{\mbox{\tiny th}}$ particle, and $\delta t$ is the numerical integration step, chosen to be $0.005\tau_{0}$ for $T\!\leq\!1.0\varepsilon/k_{B}$ and $0.001\tau_{0}$ for $T\!>\!1.0\varepsilon/k_{B}$. This thermostating scheme requires chosing a time parameter $\tau_{\mbox{\tiny Ber}}$ which controls the rate at which heat is injected into or removed from the system; chosing a large value of $\tau_{\mbox{\tiny Ber}}$ is preferable in order to minimize the intervention of the thermostat with the purely Newtonian dynamics of the simulation. In particular, the dynamics reduces to Newtonian in the limit $\tau_{\mbox{\tiny Ber}}\!\to\!\infty$. In our equilibrium runs we chose $\tau_{\mbox{\tiny Ber}}\!=\!10.0\tau_{0}$, and made sure that systems were equilibrated during several $\tau_{\mbox{\tiny Ber}}$ before collecting statistics. This is only relevant for high temperature runs, as for lower temperatures the $\alpha$-relaxation time $\tau_{\alpha}$ (see main text for definition and measurements) becomes much larger than $\tau_{\mbox{\tiny Ber}}$.

I.3 Instantaneous and finite-rate quenches

We employed two schemes for generating glassy samples. In the first scheme, we equilibrated statistically independent systems at various parent temperatures $T_{0}$, and then evolved time using fully overdamped dynamics. During the overdamped dynamics we calculate a characteristic interaction force scale $\bar{f}\equiv\sqrt{\sum_{\alpha}f_{\alpha}^{2}/N}$ and a characteristic net force scale $\bar{F}\equiv\sqrt{\sum_{i}|\vec{F}_{i}|^{2}/N}$, where $\alpha$ labels a pair of interacting particles, $f_{\alpha}\equiv-\frac{\partial\varphi}{\partial r_{\alpha}}$ is the force exerted between the $\alpha^{\mbox{\tiny th}}$ pair, $\vec{F}_{i}\equiv-\frac{\partial U}{\partial\vec{x}_{i}}$ is the net force experienced by the $i^{\mbox{\tiny th}}$ particle. The system was deemed a glass at mechanical equilibrium once the ratio $\bar{F}/\bar{f}$ dropped below $10^{-10}$. The integration step used for these overdamped runs was $\delta t/10$, where $\delta t$ was chosen as reported above.

For the second scheme we started from high temperature equilibrium liquid states at $T\!=\!1.0\varepsilon/k_{B}$, which is roughly $2T_{g}$ for our model. We then used the Berendsen thermostat to cool down our high temperature liquid states at a perscribed rate $\dot{T}$, until the final temperature of $0.05\varepsilon/k_{B}\!\approx\!T_{g}/10$ is reached. Finally, overdamped dynamics were employed as described above to remove the remaining heat.

Chosing the time parameter $\tau_{\mbox{\tiny Ber}}$

During the continuous quenches as described above the instantaneous temperature $\tilde{T}(t)$ lags behind (or strictly speaking, above) the target temperature $T(t)$, which, in this protocol, is also a function of time. We empirically observe that the lag closely follows

[TABLE]

as shown in Fig. S1 above. This means that given a quench rate $\dot{T}$, the time parameter $\tau_{\mbox{\tiny Ber}}$ must be chosen to be small enough for the instantanous temperature of the system to closely follow the target temperature. Following this constraint, we chose the time parameter $\tau_{\mbox{\tiny Ber}}$ for the continuous quenches such that the relative deviation of the instantaneous temperature from the target temperature at $T(t)\!=\!T_{g}$ remains smaller than 3%. This translates to $\tau_{\mbox{\tiny Ber}}\!=\!0.5\tau_{0}$ for $\dot{T}\!=\!10^{-2}\varepsilon/(k_{B}\tau_{0})$, and $\tau_{\mbox{\tiny Ber}}\!=\!4.0\tau_{0}$ for $\dot{T}\!=\!10^{-3}\varepsilon/(k_{B}\tau_{0})$. For the slower cooling rate of $\dot{T}\!=\!10^{-4}$ we chose $\tau_{\mbox{\tiny Ber}}\!=\!10.0\tau_{0}$ which is itself sufficiently larger than the microscopic timescale $\tau_{0}$ to have any observable effects.

To check for effects of the time parameter $\tau_{\mbox{\tiny Ber}}$ on our results, we have carried out independent continous quench runs using different values of $\tau_{\mbox{\tiny Ber}}$. In Fig. S2 we show data for the low-frequency tails of the density of vibrational modes (see details below) obtained for different time parameters $\tau_{\mbox{\tiny Ber}}$ and different quench rates $\dot{T}$. We find that our results are largely insensitive to variations of the time parameter $\tau_{\mbox{\tiny Ber}}$; following the notation of the main text, i.e. $D(\omega)\!\sim\!\omega^{\beta}$, we find $\beta\!=\!4$ for $\dot{T}\!=\!10^{-3}$, and $\beta$ slightly smaller than 4 for $\dot{T}\!=\!10^{-2}$. We clearly cannot rule out that $\beta\!\to\!4$ as $\omega\!\to\!0$. However, these data demonstrate that our conclusions are not biased by our choice of the time parameter $\tau_{\mbox{\tiny Ber}}$.

I.4 Spectra calculations

Normal mode analyses were carried out using the numerical analysis software MATLAB matlab1 . We calculated the first (lowest) 100, 200, and 400 modes for each glassy sample of systems of size $N\!=\!2000,4000$ and $10,000$, respectively.

II Statistical convergence

In this section we provide evidence indicating that our data sets are sufficiently large and therefore statistically converged. We test in particular the instantaneous quench ensemble with parent temperature $T_{0}\!=\!2.00$, and independently generate an additional data set which is 5 times larger, i.e. we instantaneously quenched 50,000 independent equilibrium configurations, using the methods described above. The results for the density of vibrational modes are displayed in Fig. S3. While the larger ensemble reveals lower frequency modes, we do not see any signs of a crossover to a different scaling at lower frequencies, and conclude therefore that our usual ensembles of 10,000 glassy samples are sufficiently large to allow for reasonable statistical convergence.

III Finite size effects

In this section we present and discuss data for the density of vibrational modes calculated in a variety of system sizes, in order to assess to what degree our observations are effected by the finite-sizes of our simulations. We focus on the instantaneous quench protocol, as reported in Fig. 2 of the main text. Before presenting and discussing our data, we remind the reader that in modes_prl1 it was shown that the possibility to observe the non-Debye low-frequency tail of the density of vibrational modes relies on delicately tuning the system sizes studied; on one hand, the spatial structure of soft glassy modes in generic models of structural glasses —such as the one studied here— are characterized by a localization length footnote1 which is of the order of 10 particle sizes, as shown in modes_prl1 . The system sizes considered must be large enough to accommodate this localization length. On the other hand, the number of vibrational modes with frequencies smaller than the lowest frequency phonon vanishes with increasing the system size. This statement can be made more quantitative; to this aim, we assume that the density of vibrational modes grows as $D(\omega)\!\sim\!\omega^{\beta}$. For a given preparation protocol, the number of glassy vibrational modes $n_{g}$ that appear below the lowest frequency phonon in a single sample of linear size $L$ follows

[TABLE]

where ${\,\mathchar 22\relax\mkern-12.0mud}$ denotes the spatial dimension. Since $\beta\!\geq\!3$ is observed in all cases, we conclude that the number of modes observed with frequencies lower than the lowest frequency phonon vanishes at least as $L^{-1}$. For protocols that generate the $\omega^{4}$ law (see main text), it vanishes as $L^{-2}$. In any event, we conclude that observing the $\omega^{\beta}$ tail of the density of vibrational modes in a statistically robust manner in larger systems requires increasingly larger ensembles of glassy samples to be generated, which quickly becomes computationally challanging.

For the reasons discussed above, we focused on the analysis of systems of $N\!=\!2000,4000$ and $10,000$ particles. Each ensemble is generated by first producing 10,000 independent equilibrium configurations at the parent temperatures $T_{0}\!=\!2.00,0.70$ and $0.60$, and performing an instantaneous quench of each of these independent equilibrium configurations as described in the previous Section. We calculated the first (lowest) 100, 200, and 400 modes for each member of our ensembles of systems of size $N\!=\!2000,4000$ and $10,000$, respectively, and obtained the distributions $D(\omega)$ of vibrational modes.

Our results are presented in Fig. S4. Our data for systems of size $N\!=\!2000$ and $N\!=\!4000$ do not show any systematic trend, and the respective slopes of $D(\omega)$ at low frequencies do not seem to depend on system size. For the larger systems of $N\!=\!10,000$ the intrusion of the lowest frequency phonon begins to ‘pull’ the distribution upwards, which leads to what could be interpreted as a steeper exponent $\beta$. We reiterate here that we do not aim at accurately determining the numeric value of the exponent $\beta$, but rather identify the trends that it follows as the preparation protocol of our glassy samples is systematically varied. We assert that, within the statistical limitations presented by our data, and within the window of system sizes in which the distributions can be robustly observed, finite size effects do not appear to be significant.