Low-frequency vibrational modes of stable glasses

Unusual features of the vibrational density of states D(ω) of glasses allow one to rationalize their peculiar low-temperature properties. Simulational studies of D(ω) have been restricted to studying poorly annealed glasses that may not be relevant to experiments. Here we report on D(ω) of zero-temperature glasses with kinetic stabilities ranging from poorly annealed to ultrastable glasses. For all preparations, the low-frequency part of D(ω) splits between extended and quasi-localized modes. Extended modes exhibit a boson peak crossing over to Debye behavior (Dex(ω) ~ ω2) at low-frequency, with a strong correlation between the two regimes. Quasi-localized modes obey Dloc(ω) ~ ω4, irrespective of the stability. The prefactor of this quartic law decreases with increasing stability, and the corresponding modes become more localized and sparser. Our work is the first numerical observation of quasi-localized modes in a regime relevant to experiments, and it establishes a direct connection between glasses’ stability and their soft vibrational modes

A morphous solids exhibit universal low-temperature properties, seen for instance in the heat capacity and thermal conductivity 1 , that differ remarkably from crystal physics. These properties are related to the vibrational density of states D(ω). For a continuous elastic medium in three dimensions, low-frequency excitations are phonons, and the density of states follows D(ω) = A D ω 2 , where A D is given by Debye theory 2 . A well-known universal feature of amorphous solids is an excess in vibrational modes over the Debye prediction that results in a peak in D(ω)/ω 2 at an intermediate frequency, called the boson peak [3][4][5][6] .
More recently, another source of 'excess modes' has been identified in computer simulations of model glasses [7][8][9][10][11][12] . It is composed of quasi-localized low-frequency modes with a density obeying D loc (ω)~ω 4 . Quasi-localized modes are observed at frequencies significantly lower than the boson peak and the link between the two phenomena is not immediate, despite some indications that they may be connected 8,13 . The quartic law was predicted long ago using phenomenological models 14,15 , reanalyzed over the years [16][17][18] , and remains the focus of intense research 19,20 . These predictions differ from two recent mean-field approaches 21,22 , which predict instead a universal non-Debye behavior that is quadratic in all spatial dimensions, also reported numerically 23 . Interest in the low-frequency localized modes extends beyond connections to theoretical models and the boson peak. It was suggested that these modes are correlated with irreversible structural relaxation in the supercooled liquid state 24 , and that the spatial distribution of these soft modes is correlated with rearrangements upon mechanical deformation and plasticity [25][26][27][28] . Localized defects are also central to theoretical descriptions of glass properties at cryogenic temperatures 29,30 .
Recent numerical insights were obtained for glasses that are very different from the ones studied experimentally, since they are prepared with protocols operating on timescales that differ from experimental ones by as many as ten orders of magnitude 31 . It is therefore unknown whether any of the vibrational, thermal, or mechanical properties derived from earlier computational study of the density of states is experimentally relevant. For example, it was reported 9,32 that D loc (ω)~ω β with β ranging from 3 to 4 depending on the glass's stability, with β = 4 for the two most stable simulated glasses created by cooling at a constant rate. It remains unclear, however, whether β would be different for glasses with stability comparable to that of the experimental glasses.
Our main achievement is to extend studies of the vibrational density of states of computer glasses to an experimentally relevant regime of glass stability for the first time. To this end, we build on the recent development of a Monte Carlo method that allows us to equilibrate supercooled liquids down to temperatures below the experimental glass transition [33][34][35] to prepare glasses that cover an unprecedented range of kinetic stability, from extremely poorly annealed systems to ultrastable glasses. We thus match the large gap between previous numerical findings and the experimental regime 36 . Recent studies have shown that that such stable glasses may differ qualitatively from ordinary computer glasses 35,37,38 . For example, qualitatively different yielding behavior of well-annealed glasses compared to that of poorly annealed glasses was reported in ref. 38 . Since rearrangements upon mechanical deformation are correlated with the spatial distribution of soft modes, this result suggested that the density of states could also evolve dramatically with the stability.

Results
System preparation. We prepare glasses by instantaneously quenching supercooled liquids equilibrated at parent temperature T p to T = 0, so that T p uniquely controls the glass stability. We find that the low-frequency part of the vibrational density of states changes considerably when T p varies, thus offering a direct link between soft vibrational modes and kinetic stability. Following earlier work 8,10 , we divide modes into extended and quasilocalized ones. As found for high parent temperature glasses 7-10 , the density of states of the quasi-localized modes follows D loc = A 4 ω 4 , with the same quartic exponent for all glass stabilities. Our work thus establishes the relevance of earlier findings about quasi-localized modes and their effect on the density of states in the experimentally relevant regime of glass stability. In addition, we find that the overall scale A 4 decreases surprisingly rapidly when T p decreases, showing that the density of the quasi-localized modes is highly sensitive to the glass stability. This rapid decrease contrasts with the modest changes found for other structural quantities, such as mechanical moduli, sound speed, and Debye frequency. Quasi-localized modes also become sparser and increasingly localized at low T p , and so the identification of soft localized modes as relevant glassy defects controlling the physics of amorphous solids becomes more convincing near the experimental glass transition. Our results also suggest that ultrastable glasses contain significantly fewer localized excitations than ordinary glasses, which appears consistent with recent experiments [39][40][41] .
We simulate a polydisperse glass forming system in three dimensions, which is a representative glass-forming computer model 33 . We use the swap Monte Carlo algorithm to prepare independent equilibrated configurations at parent temperatures T p ranging from above the onset temperature of slow dynamics T o ≈ 0.200, down to T p = 0.062, which is about 60% of the modecoupling temperature T c ≈ 0.108 (T c marks a crossover to activated dynamics and corresponds typically to the lowest temperature accessed by standard molecular dynamics). Importantly, our lowest T p is lower than the estimated experimental glass temperature T g ≈ 0.072 33 , and no previous computational study has explored such range of glass stability. In addition, we also use a very high parent temperature which we refer to as T p = ∞. We then probe vibrational properties of zerotemperature glasses produced by an instantaneous quench from equilibrated configurations at different T p . The specific simulation details are provided in Methods.
Macroscopic properties. We begin by presenting macroscopic properties of the glasses as a function of the parent temperature T p . The inherent structure energy E IS is directly related to the mobility of the particles 42 , and thus we show E IS in Fig. 1a as an indicator of the increased stability of the glass. E IS deviates from its high-temperature plateau when T p becomes smaller than the onset temperature, and decreases further with decreasing T p 43 . Similarly, the bulk modulus B decreases modestly with decreasing T p (Fig. 1b). By contrast, the shear modulus G in Fig. 1c remains nearly temperature-independent until themode-coupling temperature, which is below the onset temperature, and then the shear modulus increases with decreasing T p . Associated with the increase in the shear modulus is a decrease in the Debye level controlled by the increase of the shear modulus since the transverse speed of sound c t ¼ ffiffiffiffiffiffiffiffi G=ρ p is 2.4-2.6 times smaller than the longitudinal speed of sound c l . The overall relative variations of mechanical moduli and Debye frequency are, however, relatively mild given the broad range of glass stabilities covered in Fig. 1.
Classification of quasi-localized and extended modes. By examining the participation ratio P(ω) as a function of ω at different parent temperatures (see Fig. 2), we observe all the features that characterize the T p -dependence of the density of states. A value of P(ω) = 1 indicates a mode where all the particles participate equally, a value of P(ω) = N −1 indicates a mode where only one particle participates, and a value of P(ω) = 2/3 indicates a plane wave. The sharp peaks in P(ω) at low frequencies are due to the phonon modes, with the first peak corresponding to the first allowed transverse phonon at ω t = c t 2π/L, L being the box length. An increase in ω t indicates an increase in c t ¼ ffiffiffiffiffiffiffiffi G=ρ p . The low-frequency modes can be naturally divided into quasilocalized modes (small P) and extended modes (large P) through an appropriate thresholding procedure 8,10 , this decomposition becoming sharper as L increases and T p decreases. The value P 0 = 0.006 is appropriate, as shown in Fig. 2, but we checked that our results are not qualitatively affected by a reasonable change of P 0 . As T p decreases, phonon modes shift to larger frequencies, as expected from the evolution of the mechanical moduli, whereas quasi-localized modes become increasingly localized and well-separated from the phonons. We also checked that our results hold for small system sizes where allowed phonon modes are shifted to much higher frequencies 7 .
Properties of quasi-localized modes. We examined the density of states for the quasi-localized modes D loc (ω), which are shown in Fig. 3a for a few representative T p . At low frequencies, D loc (ω) = A 4 ω 4 for each parent temperature with a prefactor A 4 that depends on the glass stability. We show the resulting A 4 (T p ) in Fig. 3b. The prefactor A 4 stays nearly constant for high enough T p , but decreases sharply when T p decreases below the modecoupling temperature T c . This observation is robust against changing the system size. The decrease of A 4 at low T p correlates well with the evolution of shear modulus and Debye level in Fig. 3c, d. We note that a study of less stable glasses 32 found an increase in the lowest frequency of quasi-localized modes with decreasing parent temperature, which, under certain assumptions, may be related to the change of A 4 reported here. A major result of our study is that the quartic law governing D loc (ω) is obeyed irrespective of the glass stability, thus extending the validity of previous findings to the experimentally relevant regime.
In Fig. 3c we show the probability distribution for finding a mode with a participation ratio P for the modes with P < P 0 for N = 48,000 particles. With decreasing T p , the distribution becomes narrower and the peak position shifts to smaller P values. We find that the average participation ratio decreases with decreasing T p , which is evident from Fig. 3c. This confirms that these modes become more localized with decreasing parent temperature, which had been observed for less stable glasses 9,13,32 . Since the density of states is a function of the structure of the quenched system, we conclude that subtle local structural changes occur for T p below T c that strongly affect soft vibrational motion in the quenched glass.
To visualize the increasing mode localization, we define a 'softness' 25 for particle i as AðiÞ immobile background. To quantify these observations, we measured the probability distribution of A(i) (Fig. 4c). These distributions show a power-law tail at large A values, PðA i Þ ¼ λðT p ÞA Àα i with α ≈ 3.7. At low T p the tail is well separated from the core of the distribution at small A, and mobile particles with large A are better defined. There is also a pronounced decay of the probability of finding large A values at low T p since λ(0.2)/λ (0.062) ≈ 4.3, which indicates a greater than four fold decrease in the number of soft particles with large vibrational amplitudes. The interpretation of quasi-localized modes as relevant glassy defects controlling mechanical and thermal properties of glasses is therefore more convincing for stable glasses than it is for conventional computer glasses.
Properties of extended modes. Next, we examine the density of states of extended modes, D ex (ω), with a participation ratio greater than P 0 . In Fig. 5a, b we show the reduced density of states D ex (ω)/ω 2 for two parent temperatures. For each temperature, the Debye level is reached at low enough ω and a boson peak is observed at larger frequencies. Using our localization criterion, we find that modes near the boson peak are not localized, but this does not imply that they have a phononic character. The boson peak narrows slightly with decreasing T p . The Debye level, the boson peak location, height, and width all change modestly as T p is varied over the entire range studied. The changes observed in our study agree qualitatively with those found by Grigera et al. 4 .
In Fig. 5c we examine scaling properties of the density of states of extended modes. We rescale ω by the boson peak frequency, ω BP , and plot the rescaled density of states D ex /(A D ω 2 ). We observe an excellent collapse on the low-frequency side of the boson peak. This shows that in this frequency range the reduced frequency dependence has a universal shape, as reported before 44 .
Second, the collapse also shows that the height of the boson peak correlates with the Debye level A D . These results agree with experiments on molecular glass formers [45][46][47] . However, some of the same experiments report that the boson peak position scales as the Debye frequency 45,47 , which is not consistent with our results. We also find that a scaling of ω BP with the bulk modulus suggested in ref. 48 is inconsistent with our results. Note that we study the evolution of the boson peak as a function of the preparation temperature, while experiments sometimes examine the temperature evolution of the boson peak for a given glass preparation. We also note that a correlation between the boson peak and quasi-localized modes has been proposed by studying systems at different pressures around the unjamming transition 49 .
Since the boson peak occurs in a different frequency range than the ω 4 scaling of D loc (ω), it is not clear that there could be a relationship between the boson peak and the low-frequency quasi-localized modes. Simulations close to jamming suggest that A 4 $ ω 4 BP 8 , but we do not find that this relation holds with changing T p . An alternative possibility can be obtained from dimensional analysis, where a characteristic frequency for quasilocalized modes can be defined as A À1=5 4 32 . We find that A À1=5 4 $ ω BP for glasses with T p < T c (Fig. 6), but this relation does not hold for glasses created with T p > T c . We note that ω BP is constant for T p > T c , see the inset to Fig. 6, and only changes for T p < T c . Again we find that T c marks a change in the behavior of D(ω).
Given the relatively small changes in both ω BP and A entire range of parent temperatures studied, it is not clear that a power law is the proper relationship between these quantities and further work is needed to verify it.

Discussion
In summary, we report the first characterization of the vibrational density of states of computer glasses prepared over a range of glass stability that bridges the gap between ordinary simulations and experimental studies. At low-frequency extended and quasilocalized modes coexist, and both types of modes evolve differently when the glass stability is varied. We find a relatively mild temperature dependence of extended modes, with a strong correlation between the Debye level and the boson peak. By contrast, quasi-localized modes evolve more strongly when T p decreases below the mode-coupling temperature, but their density of states is always described by D loc $ A 4 ω 4 . Unexpectedly, the temperature dependence of the prefactor A 4 (T p ) is more interesting than the value of the quartic exponent, which is insensitive to the degree of annealing. The increasing localization of the modes implies that subtle yet significant changes occur in the local structure of the glass that are not reflected in the pair correlation function, which is nearly identical for parent temperatures below T c . Since soft modes have been linked to irreversible relaxation 24 and rearrangements under shear [25][26][27][28] , it follows that the reduction of these soft modes can have significant implications for glassy dynamics. In turn this reduction indicates that there are fewer soft spots, which should increase the strength of the glass. This hypothesis is supported by the observation that the decrease in D loc (ω) mirrors the increase of the shear modulus, and also correlates very well with the evolution of the ductility of the produced glasses 38,50 . Since we can now equilibrate amorphous systems at temperatures low enough so that they do not flow, another perspective would be to analyze the density of states at finite temperatures through the Fourier transform of the velocity autocorrelation function 51 , or by diagonalizing the covariance matrix of displacements 52 . Future studies should examine the difference between these procedures to provide insights into thermal anharmonicities of stable glasses, and more generally into their low-temperature transport properties.

Methods
Simulations. We simulate a polydisperse model glass former of sizes between N = 48,000 and 450,000 particles with equal mass at a number density ρ = 1.0 33 . The interaction between two particles i and j is given by Vðr ij Þ ¼ σ ij r ij 12 þvðr ij Þ when their separation r ij r c ij ¼ 1:25σ ij and zero otherwise. We use , where the coefficients c 0 , c 2 , and c 4 ensure the continuity of V(r ij ) up to the second derivative at the cutoff r c ij . The probability of particle diameters σ is P(σ) = A/σ 3 , where σ∈[0.73,1.63] and we use a non-additive mixing rule, σ ij ¼ σ i þ σ j 2 ð1 À 0:2jσ i À σ j jÞ. For N ≤ 192,000 we use the swap Monte Carlo algorithm to prepare independent equilibrated configurations at parent temperatures T p ranging from above the onset temperature of slow dynamics (T o ≈ 0.200) down to T p = 0.062, which is about 60% of the mode-coupling temperature (T c ≈ 0.108), and is lower than the estimated experimental glass temperature (T g ≈ 0.072) 33 . In addition, we also use a very high parent temperature, which we refer to as T p = ∞. Due to very long equilibration times for systems of N > 192,000 particles we only study systems with N > 192,000 for T p = ∞.
Density of states calculation. Following equilibration at a temperature T p , zerotemperature glasses are produced by instantaneously quenching equilibrium configurations to their inherent structures using the Fast Inertia Relaxation Engine algorithm 53 . We then calculate the modes by diagonalizing the Hessian matrix using Intel Math Kernel Library (https://software.intel.com/en-us/mkl/) and ARPACK (http://www.caam.rice.edu/software/ARPACK/). We calculate all the normal modes for the 48,000 particle systems, but only the low-frequency part of the spectrum in systems with N > 48,000. We characterize the modes through the density of states DðωÞ ¼ 1 3NÀ3 P 3NÀ3 l¼1 δðω À ω l Þ and the participation ratio  $ ω BP for glasses whose T p < T c , which is the parent temperature range where we see an increase in ω BP with decreasing T p , see inset where the vertical dasheddotted, dashed lines mark the positions of T g and T c , respectively NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-018-07978-1 ARTICLE NATURE COMMUNICATIONS | (2019) 10:26 | https://doi.org/10.1038/s41467-018-07978-1 | www.nature.com/naturecommunications Pðω l Þ ¼ P N i¼1 je l;i j 2 À Á 2 N P N i¼1 je l;i j 4 , where e l,i is the polarization vector of particle i in mode l with frequency ω l . For a mode localized to one particle P(ω) = N −1 , and for an ideal plane wave P(ω) = 2/3. The phonon modes occur at discrete frequencies, and care has to be taken in the binning procedure to calculate the density of states of extended modes, D ex (ω). To perform this calculation, we determine the phonon frequencies from the peak positions of the participation ratio versus frequency, and tune the bin size to smooth D ex (ω).To obtain the shear modulus G and the bulk modulus B we use the method described in ref. 54 .

Data availability
All data will be available from the authors upon request.