Low-frequency vibrational density of states of ordinary and ultra-stable glasses

A remarkable feature of disordered solids distinct from crystals is the violation of the Debye scaling law of the low-frequency vibrational density of states. Because the low-frequency vibration is responsible for many properties of solids, it is crucial to elucidate it for disordered solids. Numerous recent studies have suggested power-law scalings of the low-frequency vibrational density of states, but the scaling exponent is currently under intensive debate. Here, by classifying disordered solids into stable and unstable ones, we find two distinct and robust scaling exponents for non-phononic modes at low frequencies. Using the competition of these two scalings, we clarify the variation of the scaling exponent and hence reconcile the debate. Via the study of both ordinary and ultra-stable glasses, our work reveals a comprehensive picture of the low-frequency vibration of disordered solids and sheds light on the low-frequency vibrational features of ultra-stable glasses on approaching the ideal glass.

A remarkable feature of disordered solids distinct from crystals is the violation of the Debye scaling law of the low-frequency vibrational density of states.Because the low-frequency vibration is responsible for many properties of solids, it is crucial to elucidate it for disordered solids.Numerous recent studies have suggested power-law scalings of the low-frequency vibrational density of states, but the scaling exponent is currently under intensive debate.Here, by classifying disordered solids into stable and unstable ones, we find two distinct and robust scaling exponents for non-phononic modes at low frequencies.Using the competition of these two scalings, we clarify the variation of the scaling exponent and hence reconcile the debate.Via the study of both ordinary and ultra-stable glasses, our work reveals a comprehensive picture of the low-frequency vibration of disordered solids and sheds light on the low-frequency vibrational features of ultra-stable glasses on approaching the ideal glass.
Low-temperature properties of solids, such as specific heat and thermal conductivity, are closely related to the excitation of low-frequency vibrational states.For crystals, it is well-established that the vibrational states, i.e., phonons, form a low-frequency vibrational density of states (VDOS) following the Debye scaling law: D(ω) ~ωd−1 , where ω is the frequency and d is the spatial dimension, resulting in the T d scaling of the specific heat at low temperatures T 1 .The thermal conductivity is believed to be governed by the specific heat, phonon mean free path, and sound velocity.In crystals, because the phonon mean free path and sound velocity remain approximately constant in temperature, the thermal conductivity follows the low-temperature scaling of the specific heat 1 .
However, we face great challenges when dealing with disordered solids such as glasses.The low-temperature scalings of the specific heat and thermal conductivity are no longer T d2-4 .When T < 1K, the specific heat is linearly scaled with T [2][3][4] , which is attributed largely to the existence of two-level systems instead of the VDOS 5,6 .It is also believed that the two-level systems change the mean free path, causing anomalous behaviors of the thermal conductivity.At higher temperatures, the VDOS matters.The disordered structure of glasses causes the coexistence of phonon-like and non-phononic modes at low frequencies [7][8][9][10][11][12][13][14] , so the VDOS is at least a superposition of the Debye scaling and that of the non-phononic modes.The excess non-phononic modes form a peak in D(ω)/ω d−1 , defined as the boson peak 11,12,15 .It has been shown that the boson peak may be correlated with the simultaneity of the peak in c p /T 3 , with c p being the constant-pressure specific heat, and the plateau in the thermal conductivity at the boson peak temperature (~10K for typical glasses such as vitreous silica) 4 .Both simulations and experimental measurements such as the neutron scattering and X-ray, have significantly advanced our understanding of the constituent modes of the boson peak [11][12][13][14][15][16] .However, what the VDOS of non-phononic modes looks like below the boson peak frequency is still an unsettled issue 7,8,[17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] , which is crucial to understanding the thermal properties in the 1-10K temperature regime.In addition to the thermal properties, the anomalous low-frequency non-phononic modes have been successfully applied to understand various other properties of disordered solids, e.g., mechanical failure [32][33][34][35][36][37] , glass transition 13,38,39 , and heterogeneous dynamics of glass-forming liquids [40][41][42] .
Note that generic glasses lie at local minima of the complex energy landscape 56 , whose stabilities can vary a lot from each other.One can tell that the variation of α mentioned above is more or less related to the stability.However, the local minima with various degrees of stability were always mixed up to calculate the VDOS in previous studies.Moreover, probably limited by the development of experimental techniques, as far as we know, there have been no direct experimental measurements of α for molecular glasses.Therefore, the examination of α has heavily relied on simulations.In most of previous simulations, the VDOS was calculated for systems with periodic boundary conditions, whose shapes were not allowed to change.However, it has been shown that some glasses that are stable under periodic boundary conditions may be unstable under certain deformations 57,58 .Apparently, the effects of such deformation stability on the VDOS were completely overlooked.
Here, we systematically study the low-frequency VDOS for both ordinary and ultra-stable model glasses quenched from different parent temperatures T p .Remarkably different from previous approaches, we divide all glasses into two categories: stable ones, which can resist any infinitesimal deformations, and unstable ones, which are unstable subject to some infinitesimal deformations, and calculate their VDOSs separately.The VDOSs for stable and unstable solids depart from each other below a crossover frequency ω d , where they have different scaling exponents.For unstable solids, α = α u ≈ 3.3, independent of system size and spatial dimension.For stable solids, α = α s ≈ 5.5 and 6.5 in 2D (d = 2) and 3D (d = 3), respectively, which does not vary with system size either.The superposition of these two VDOSs results in the VDOS studied in previous approaches.This explains the variation of α under various circumstances.Moreover, we observe the emergence of an ω 4 scaling right above the ω α s one when the system size of stable solids increases for both ordinary and ultra-stable glasses in 3D.Interestingly, our results suggest that the number of non-phononic modes forming the ω α s and ω 4 scalings decays with the decrease of T p , possibly vanishing at a sufficiently low T p .Therefore, our study may shed light on the perspective of the vibrational features of the ideal glass.

Results
In this work, we mainly show results for systems composed of polydisperse soft particles interacting via the inverse-power-law (IPL) potential (see Methods for details), which have been widely used to study the glass transition 8,26,28,[59][60][61][62] .In Supplementary Fig. 1 of the Supplementary Information and a parallel study, we also show consistent results for Lennard-Jones and harmonic potentials, suggesting the generality of our findings.We obtain the zero-temperature glasses by instantaneously quenching liquids equilibrated at the parent temperature T p .It is well-known that the stability of quenched glasses increases with the decrease of T p when T p is lower than the onset temperature T on , i.e., the crossover temperature from Arrhenius to super-Arrhenius dynamics 56 .We will first study glasses obtained from a given T p and discuss the T p dependence afterward.

VDOSs for stable and unstable solids
In most of the previous simulations, the normal modes of vibration were obtained from the diagonalization of the normal Hessian matrix, with the elements being the second derivatives of the potential energy with respect to particle coordinates.No boundary deformation was taken into account in such an approach.A glass was treated as a stable one if all nontrivial eigenvalues of the normal Hessian matrix were positive.However, this cannot guarantee that the glass is stable subject to boundary deformations.If we introduce the d(d + 1)/2 degrees of freedom corresponding to the boundary deformations (shear and compression) and construct the extended Hessian matrix (see Methods), the matrix of some glasses may have negative eigenvalues, indicating that the glasses are unstable under some deformations.We thus define these glasses as unstable glasses.On the other hand, the glasses whose extended Hessian matrix has no negative eigenvalues are defined as stable glasses.Note that the extended Hessian matrix is only used to classify all glasses into stable and unstable ones, and VDOS is still calculated from the normal Hessian matrix.Here, we denote D s (ω), D u (ω), and D(ω) as the VDOSs of stable, unstable, and all glasses, respectively.
Figure 1a, b compares D(ω), D s (ω), and D u (ω) in 2D and 3D, respectively.They collapse above a crossover frequency ω d , and depart from each other otherwise.Both low-frequency tails of D s (ω) and D u (ω) display a clear power-law scaling behavior, D s ðωÞ ∼ ω α s and D u ðωÞ ∼ ω α u .Beyond that, D u (ω) forms a valley bottomed at ω d , while D s (ω) still monotonically increases and transits to ω d .However, α s and α u are apparently different.In 2D and 3D, α s = 5.5 ± 0.2 and 6.5 ± 0.2, respectively.In contrast, α u = 3.3 ± 0.1 in both 2D and 3D.This α u value is close to the α ≈ 3 arguments of the fold instability model 36 .Note that α ≈ 3 is obtained based on the approximation that the distribution of the stress distance to instabilities is constant 36 , which may fluctuate if the distribution is not strictly flat.The fold instability model is raised for glasses with weak stability and is prone to rearrangement upon deformations and does not rely on spatial dimension.This agreement thus proposes a plausible physical origin of the scaling behavior of D u (ω).
By definition, the low-frequency part of D(ω) should be the superposition of D s (ω) and D u (ω): where f s is the fraction of stable glasses, and A s and A u are prefactors of D s (ω) and D u (ω), respectively.In Fig. 1a, b, we compare the simulated D(ω) with the prediction by Eq. (1) (dashed line).They are in excellent agreement at low frequencies.
As done in previous studies, the low-frequency part of D(ω) can be fitted with ω α .Figure 1a, b indicates that α should be between α u and α s , if we perform the fitting.Now, the excellent agreement between D(ω) and Eq. ( 1) provides another interpretation of the α value at the lowfrequency tail.If the values of α s and α u are definite, the α value is jointly determined by f s , A s , and A u , which may change with parameters such as system size and parent temperature.We are thus able to understand why α was reported to vary under some circumstances [20][21][22][23][24] .Moreover, at low enough frequencies, D u (ω) dominates.This may be the reason why lower values of α were always observed when rather low-frequency regimes were accessed 24,27,28 .
Figure 1c, d compares the participation ratio, P s (ω) and P u (ω), of stable and unstable solids.A mode with a lower participation ratio is more localized.We can see that, below the first Goldstone (phononlike) mode, the low-frequency modes forming the ω α s and ω α u scalings have the lowest participation ratios and are thus most quasi-localized on average.However, the degrees of quasi-localization of stable and unstable solids are similar, only that unstable solids extend to lower frequencies.
Figure 2a, b visualizes the structures of the modes lying in the ω α u and ω α s scaling regimes.They both exhibit the typical feature of quasilocalized modes with localized regions hybridizing with the planewave-like background.For the unstable solid in Fig. 2a, we show in Fig. 2c its unstable mode of the extended Hessian matrix, whose eigenvalue is negative.It looks almost identical to the mode in Fig. 2a.The dot product of the two normalized modes in Fig. 2a, c is 0.997.Note that when a disordered solid approaches the fold instability under load such as shear and compression, its lowest-frequency mode is responsible for the instability, whose frequency decays to zero following a power law while its structure remains unchanged 36 .This type of mode contributes to the ω 3 behavior predicted by the fold instability argument 36 .Therefore, the perfect agreement between the lowestfrequency mode of the unstable solid and the unstable mode of the extended Hessian matrix is the evidence supporting our argument that α u ≈ 3.3 originates from fold instabilities.Figure 2d  compression, which is the typical form of the boundary deformation of unstable modes.

System size dependence
Recently, it was reported that the value of α for D(ω) increased with the growth of system size for ordinary glasses quenched from high parent temperatures 20,21,24 .As shown in Fig. 3a, for 2D systems, α indeed grows from 3.4 to 4 when the system size N changes from 256 to 4096.Interestingly, Fig. 3b, c shows that α s and α u remain constant in N.However, both A s and A u grow with N.Meanwhile, f s increases when N increases, which can be fitted well with 1 − f s ~N −1.4 , as illustrated in Fig. 3d.Therefore, the system size dependence of α in 2D directly reflects the competition among f s , A s , and A u .
Figure 3 e-h indicates that similar system size evolution happens in 3D.When system size increases, α gradually increases.Again, α s and α u are insensitive to the change in system size, while A s , A u , and f s grow with N.However, the comparison between Fig. 3b, f demonstrates a seeming difference between 2D and 3D.In 2D, the ω α s scaling extends all the way to the crossover frequency ω d , above which D s (ω) and D u (ω) collapse.In 3D, the ω α s scaling is deviated above another crossover frequency ω s < ω d .Figure 3f shows that, when system size increases, ω s decreases so that the frequency regime for the ω α s scaling to survive is suppressed.
In addition to ω d and ω s , there is another characteristic frequency ω p < ω d of the first peak in D u (ω).In Fig. 3d, h, we show the system size dependence of these three characteristic frequencies.In both 2D and 3D, ω d is approximately scaled with N −1/d .Since D s (ω) and D u (ω) deviate below ω d , if such system size dependence persists on approaching the thermodynamic limit, we would expect that the ω α s and ω α u scalings tend to disappear so that D s (ω) and D u (ω) eventually become identical to D(ω).Note that the Goldstone modes have the same system size dependence.It may be plausible to ask whether ω d is associated with some inherent properties of disordered solids such as the elastic moduli, which contribute to the Goldstone modes.Figure 3d, h shows that ω p is approximately scaled with N −0.55 in both 2D and 3D.As seen from Fig. 3h, ω s (N) in 3D roughly agrees with ω p (N).At the current stage, we are not able to confirm whether there are any physical origins of these characteristic frequencies and hope to leave them to future investigations.
In Fig. 4, we collapse the low-frequency parts of D s (ω) and D u (ω) for different system sizes by plotting N Àν s D s ðωÞ and N Àν u D u ðωÞ against ωN ν s and ωN ν u , respectively.These scalings conserve the integrals of the VDOSs.Our best data collapse gives ν s ≈ 0.21 for both 2D and 3D and ν u ≈ 0.35 and 0.28 for 2D and 3D, respectively.The scaling collapse indicates that A s ∼ N ðα s + 1Þν s and A u ∼ N ðα u + 1Þν u , respectively.
Seen from Fig. 3f, the ω > ω s part of D s (ω) in 3D shows the trend to converge to a master curve when system size increases.Figure 3f also indicates that D s (ω) reaches the maximum at ω = ω * ≈ 5, above which D s (ω) is plateau-like and gradually decreases.In Fig. 5a, we focus on ω < ω * with more system sizes.There seem to be three consecutive frequency regimes with different scalings: (i) ω α s when ω < ω s , (ii) ω α 1 when ω s < ω < ω 0 , and (iii) ω α 2 when ω 0 < ω < ω * .Unlike the size-independent α s , α 1 and α 2 evolve with N. For the largest system sizes studied here, we can observe the emergence of α 1 ≈ 4 and α 2 ≈ 1.5.Right below the plateau of the VDOS (ω 0 < ω < ω * ), mean-field theories predict an ω 2 behavior due to marginal stability 63,64 .For the system sizes studied here, α 2 slightly varies with system size.Although we are not able to exclude the possibility that α 2 could approach 2 in sufficiently large systems, α 2 ≈ 1.5 observed here is still apparently lower than the mean-field value.It thus remains a question whether α 2 is meaningful and related to marginal stability.Recent studies suggest that quasi-localized modes below ω 0 could form the ω 4 scaling.Here, we see this scaling right above ω s .However, whether this scaling is real or is just a crossover still needs to be examined in sufficiently large systems with good statistics.Note that, even if ω 4 could be real, our results suggest that it does not generally exist.As shown in Fig. 3b, in 2D, there is no sign for the ω 4 behavior to emerge in D s (ω).
Assuming that the system size evolution of ω s is still valid in much larger systems, we can expect that there is always a contribution of the ω α s scaling below ω s , as long as the system size is finite.Therefore, the low-frequency tail of D(ω) is always jointly determined by D s (ω) and D u (ω) according to Eq. (1).e-h Results in 3D.In (h), we also show the system size depends of the frequency ω s below which the ω α s scaling exists.The parent temperature T p is approximately the onset temperature T on for both 2D and 3D systems.The dashed lines show the power-law scalings.

Parent temperature dependence
It was also reported that α grew when the parent temperature T p decreased 22,23 .To understand this T p dependence, we study the VDOSs of glasses quenched from various T p ranging from above the onset temperature T on to near the glass transition temperature T g .Figure 6 compares D(ω), D s (ω), and D u (ω) near the three representative temperatures, T on , T mc (mode-coupling temperature), and T g , for both 2D and 3D glasses.Figure 6a, e shows that α evolves roughly from 3.4 to 4 when T p decreases from T on to T g in both 2D and 3D, consistent with previous studies 22,24 and similar to the evolution with system size.The difference is that the low-frequency part of D(ω) apparently decays with the decrease of T p . Figure 6b, c (f and g) indicates that α s and α u also remain constant in T p in 2D (3D).Unlike the system size dependence, A u is insensitive to the change of T p .Therefore, the evolution of D(ω) at low frequencies is jointly determined by A s and f s .When T p decreases, Fig. 6b, f shows that A s decreases, while f s increases, as shown in Fig. 6d, h.
For 3D ultra-stable glasses quenched from T p ≈ T g , Fig. 5b shows that the three frequency regimes in D s (ω) discussed above also emerge.Compared to the high T p case, the intermediate ω 4 scaling seems to be more pronounced.For example, for smaller systems with N ≤ 1000, there is no apparent ω 4 scaling in Fig. 5a, but we can already see it in Fig. 5b.Again, the authenticity of the ω 4 scaling needs to be verified in sufficiently large systems, which is however still absent in 2D ultra-stable glasses.
Figure 7 directly displays the T p dependence of A s .For 3D systems, we also show the prefactor A 4 of the ω 4 scaling above ω s .Both A s and A 4 keep decreasing with the decrease of T p .Similar T p dependence was reported for the prefactor of the ω 4 scaling of D(ω) 8,31 .At the current stage, it is difficult to obtain reliable results at much lower T p .If such T p dependence persists at even lower T p , the number of low-frequency non-phononic modes significantly decreases and could be expected to vanish at low enough T p .If this is the case, such low-temperature ultrastable glasses will only have phonon-like modes at low frequencies.Although structurally disordered, evaluated by conventional criteria, the glasses could behave like crystals at long wavelengths.In fact, there was experimental evidence of the low-temperature Debye scaling for ultra-stable glasses [65][66][67] , supporting that the non-phononic mode contribution can be negligible if the glass reaches the highest stability.It is thus interesting to figure out whether such ultra-stable glasses are prototypes of ideal glasses and under what temperatures they could exist.

Discussion
By classifying disordered solids into stable and unstable ones, we find that their VDOSs, D s (ω) and D u (ω), depart from each other when ω < ω d , with the low-frequency tails following distinct scaling laws, D s ðωÞ ∼ ω α s and D u ðωÞ ∼ ω α u , respectively.The robustness of the values of α s and α u is verified by the solids with different sizes and quenched from different parent temperatures.Using this classification, we can understand the variation of the scaling exponent α reported previously.For finite-size disordered solids, it is due to the existence of unstable solids and the competition among the fraction of (un)stable solids and prefactors of the two scalings.Because unstable solids are inevitable in confined systems, our study can explain a recent experimental observation of α ≈ 3 in a confined quasi-2D nanosystem 68 .
We also find that when system size increases, ω d decreases so that the two scalings are pushed to lower frequencies.For stable solids in 3D, the ω α s scaling only exists below another crossover frequency ω s < ω d , which also decreases with the increase of system size.When ω > ω s , our results show the trend of the emergence of the ω 4 scaling in the largest systems studied and ultra-stable glasses.At the current stage, we cannot confirm the authenticity of the ω 4 scaling, which requires the verification of sufficiently large systems in future studies.For stable solids in 2D, we do not see the emergence of the ω 4 scaling.Moreover, the prefactors of the ω α s and ω 4 scalings both decrease with the decrease of parent temperature, implying the possible existence of ultra-stable glasses with only crystal-like low-frequency vibrations at low enough temperatures.Such glasses may act as prototypes of the ideal glass.
In this work, we are focused on generic glasses, which are constrained well above isostaticity 43,44,69,70 .Marginally jammed solids near isostaticity are less stable than generic glasses concerned here.It is thus interesting to know whether and to what extent our findings here are applicable to marginally jammed solids.There are mean-field theories proposing the α = 2 scaling of the VDOS 63,64 .The competition between different theoretical frameworks may complicate the vibrational features of marginally jammed solids moving away from the jamming transition 43,44 .We leave these discussions to a separate study.

Simulation model
Our systems contain N polydisperse particles in a simulation cell with side length L and periodic boundary conditions in all directions.All particles have the same mass m.Particles i and j interact via the IPL  potential: when their separation r ij ≤ 1.25σ ij , and zero otherwise.The coefficients c 0 , c 2 , and c 4 ensure the continuity of the potential up to the second derivative at the cutoff.The particle diameter σ is extracted from a continuous distribution P(σ) = Aσ −3 , where A is the normalization factor and σ ∈ [σ m , σ M ] with σ m /σ M = 0.4492.To enhance the glass-forming ability, we adopt a non-additive mixing rule to determine σ ij in Eq. ( 2): where ϵ measures the degree of non-additivity.We choose ϵ = 0.2 to achieve a better performance 59 .We set the average particle diameter σ, particle mass m, and the Boltzmann constant k B to be 1.The number density ρ = N/L d is 1.01 and 1.0 for 2D and 3D, respectively.
We use an efficient swap Monte Carlo algorithm 59 to prepare wellequilibrated liquids at parent temperatures T p .The onset, mode-coupling, and glass transition temperatures for our IPL model systems are T on ≈ 0.25(0.20),T mc ≈ 0.123(0.108),and T g ≈ 0.082(0.072) in 2D (3D), respectively 8,60 .After equilibration at the parent temperature T p , the liquids are rapidly quenched to zero temperature to obtain the zero-temperature glasses (inherent structures) via the fast inertial relaxation engine algorithm 71 .

Vibrational quantities
We consider two types of Hessian matrix.The normal Hessian matrix is defined as where R = (r 1 , r 2 , …, r N ) with r i (i = 1, 2, …, N) being the location of particle i.The normal Hessian matrix does not take any boundary deformation into account.In comparison, the extended Hessian matrix with (dN + n ex ) × (dN + n ex ) dimensions is 58 where n ex = d(d + 1)/2 is the extra degrees of freedom of the system and R = ðr 1 ,r 2 , . . .,r N , ϵ 1 , ϵ 2 , . . .,ϵ n ex Þ with ϵ i (i = 1, 2, …, n ex ) being the strain of the i − th deformation.The strains ϵ i are upper triangular elements of the d × d strain tensor where β j (j = 1, 2, …, d) denotes the Cartesian coordinates.These n ex degrees of freedom involve boundary deformations, including compression (expansion) and shear.More details and the stability analysis using the extended Hessian matrix can be found in Ref. 58.The normal modes of vibration are obtained by diagonalizing the matrix using the Intel Math Kernel Library 72 .If the extended Hessian matrix has negative eigenvalues, the system is unstable to certain boundary  deformations 58 .Otherwise, the system is stable to compression and shear in arbitrary directions.The participation ratio of a normal mode j is calculated as where e j i is the polarization vector of particle i in mode j.In the calculation of the VDOSs and the participation ratio, we exclude some lowest-frequency localized modes caused by rattler-like particles, as explained in Supplementary Fig. 2 of the Supplementary Information.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.

Fig. 1 |
Fig. 1 | Comparison of VDOS and participation ratio of stable, unstable, and all glasses.a VDOSs of 2D systems with N = 256 and T p = 0.12.b VDOSs of 3D systems with N = 1000 and T p = 0.18.The solid lines are power-law fittings to D s (ω) and D u (ω) at low frequencies.The red dashed lines are results from Eq. (1).They are in excellent agreement with the simulated D(ω) at low frequencies.c and d show the participation ratio of stable and unstable solids for the same systems in (a) and (b), respectively.The vertical dashed lines mark the frequency of the first Goldstone mode.

Fig. 2 |
Figure1c, d compares the participation ratio, P s (ω) and P u (ω), of stable and unstable solids.A mode with a lower participation ratio is more localized.We can see that, below the first Goldstone (phononlike) mode, the low-frequency modes forming the ω α s and ω α u scalings have the lowest participation ratios and are thus most quasi-localized on average.However, the degrees of quasi-localization of stable and unstable solids are similar, only that unstable solids extend to lower frequencies.Figure2a, b visualizes the structures of the modes lying in the ω α u and ω α s scaling regimes.They both exhibit the typical feature of quasilocalized modes with localized regions hybridizing with the planewave-like background.For the unstable solid in Fig.2a, we show in Fig.2cits unstable mode of the extended Hessian matrix, whose eigenvalue is negative.It looks almost identical to the mode in Fig.2a.The dot product of the two normalized modes in Fig.2a, c is 0.997.Note that when a disordered solid approaches the fold instability under load such as shear and compression, its lowest-frequency mode is responsible for the instability, whose frequency decays to zero following a power law while its structure remains unchanged36 .This type of mode contributes to the ω 3 behavior predicted by the fold instability argument36 .Therefore, the perfect agreement between the lowestfrequency mode of the unstable solid and the unstable mode of the extended Hessian matrix is the evidence supporting our argument that α u ≈ 3.3 originates from fold instabilities.Figure2dillustrates how the boundary deforms associated with the unstable mode of the extended Hessian matrix shown in Fig.2c.It involves both shear and

Fig. 3 |
Fig. 3 | System size dependence of the VDOSs.a-d VDOSs of all stable and unstable glasses, D(ω), D s (ω), and D u (ω), in 2D and system size evolution of the fraction of stable glasses f s , the frequency ω p of the first peak in D u (ω), and the frequency ω d below which D u (ω) and D s (ω) depart from each other, respectively.

Fig. 4 |
Fig. 4 | Scaling collapse of the low-frequency parts of the VDOSs for different system sizes.The VDOSs of stable and unstable glasses, D s (ω) and D u (ω), collapse at low frequencies, when D s ðωÞN Àν s and D u ðωÞN Àν u are plotted against ωN ν s and ωN ν u , respectively.Results of 2D glasses are shown in (a) and (b), while (c) and (d) show results of 3D glasses.Here ν s = 0.21 for both 2D and 3D; ν u = 0.35 and 0.28 for 2D and 3D, respectively.The solid lines are power-law fittings to the collapsed curves.

Fig. 6 |
Fig. 6 | Parent temperature dependence of the VDOSs.a-d VDOSs of all stable, and unstable glasses, D(ω), D s (ω), and D u (ω), and parent temperature evolution of the fraction of stable glasses f s for N = 256 systems in 2D.e-h Results for N = 1000 systems in 3D.The vertical dashed lines in (d) and (h) locate T on , T mc , and T g , respectively.The dashed lines in the other panels show the power-law scalings.

Fig. 7 |
Fig. 7 | Parent temperature dependence of prefactors of the VDOS of stable glasses.a Prefactor A s versus T p in 2D.b Prefactors A s and A 4 versus T p in 3D.The vertical dashed lines locate T on , T mc , and T g , respectively.