# The nature of non-phononic excitations in disordered systems

Walter Schirmacher, Matteo Paoluzzi, Felix Cosmin Mocanu, Dmytro

Khomenko, Grzegorz Szamel, Francesco Zamponi, Giancarlo Ruocco

► To cite this version:

Walter Schirmacher, Matteo Paoluzzi, Felix Cosmin Mocanu, Dmytro Khomenko, Grzegorz Szamel, et al.. The nature of non-phononic excitations in disordered systems. Nature Communications, 2024, 15 (1), pp.3107. 10.1038/s41467-024-46981-7. hal-04591871

HAL Id: hal-04591871

https://u-paris.hal.science/hal-04591871v1

Submitted on 24 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The nature of non-phononic excitations in disordered systems

Received: 6 October 2023

Accepted: 18 March 2024

Published online: 10 April 2024

<figure><img src="image_3.png" /></figure>

Check for updates

Walter Schirmacher Ⓔ<fnref n="1" />,²✉, Matteo Paoluzzi Ⓔ³,⁴,<fnref n="5" />, Felix Cosmin Mocanu Ⓔ⁶,⁷, Dmytro Khomenko Ⓔ<fnref n="5" />, Grzegorz Szamel Ⓔ⁸, Francesco Zamponi Ⓔ<fnref n="5" />,⁷ & Giancarlo Ruocco Ⓔ²,<fnref n="5" />✉

The frequency scaling exponent of low-frequency excitations in microscopically small glasses, which do not allow for the existence of waves (phonons), has been in the focus of the recent literature. The density of states $g(\omega)$ of these modes obeys an $\omega^s$ scaling, where the exponent $s$, ranging between 2 and 5, depends on the quenching protocol. The origin of these findings remains controversial. Here we show, using heterogeneous-elasticity theory, that in a marginally-stable glass sample $g(\omega)$ follows a Debye-like scaling ($s = 2$), and the associated excitations (type-I) are of random-matrix type. Further, using a generalisation of the theory, we demonstrate that in more stable samples, other, (type-II) excitations prevail, which are non-irrotational oscillations, associated with local frozen-in stresses. The corresponding frequency scaling exponent $s$ is governed by the statistics of small values of the stresses and, therefore, depends on the details of the interaction potential.

Understanding the nature of vibrational states in glasses is crucial for gaining insight into their mechanical and thermal properties<fnref n="1" />. Correspondingly a large amount of experimental<fnref n="1" />⁻¹⁰ theoretical¹¹⁻²² and simulation work²³⁻³⁰ has been undergone in the last 50 years. A paradigm for the anomalous vibrational features of glasses is the so-called boson peak (BP), which is an enhancement of the vibrational density of states (DOS) $g(\omega)$ over Debye's $g(\omega) - \omega^2$ law, where $\omega = 2\pi v$ is the angular frequency, and $v$ is the frequency. Such an enhancement is observed in experimental studies of macroscopically large samples. The nature of the boson-peak anomaly has been debated controversially³¹,³² (see³³ for a listing and a discussion of various proposed models for the the boson peak), but – as we feel – in the light of heterogeneous-elasticity theory (HET)¹⁹,²⁰,³⁴ it became clear that the boson peak marks a cross-over of the vibrational state's nature from occasionally scattered plane waves ("phonons") at low frequencies, as considered by Debye³⁵, to random-matrix type states in the BP region, as demonstrated first by

Schirmacher et al.¹⁴. To avoid confusion, from now on we will call these random-matrix type states 'type-I non-phononic' states.

On the other hand, in simulations of samples of reduced size, a different type of low-frequency modes was detected. We call these modes "type-II" non-phononic modes.

As early as 1991 Laird and Schober²³ found, in a simulation of a glass made up of 500 particles, low-frequency states that appeared to be localized, as estimated from the participation ratio. They coined the term "quasi-localized" excitations, because in a larger system these excitations must hybridize with the waves, like local oscillatory defects in crystals³⁶. Such non-phononic modes, which become visible at low frequency in small samples have more recently attracted a lot of attention³⁷,³⁸,³⁹⁻⁴¹. Many of the reported type-II non-phononic excitations have been found to exhibit a DOS $g(\omega) - \omega^s$ scaling with $s = 4$. Recently, in a detailed MD study of low-frequency non-phononic excitations, a continuous change of the DOS exponent from $s = 4$ to the Debye-like

<footnote n="1">¹Institut für Physik, Staudinger Weg 7, Universität Mainz, D-55099 Mainz, Germany. ²Center for Life Nano Science @Sapienza, Istituto Italiano di Tecnologia, 291 Viale Regina Elena, I-00161 Roma, Italy. ³Istituto per le Applicazioni del Calcolo del Consiglio Nazionale delle Ricerche, Via Pietro Castellino 111, 80131 Napoli, NA, Italy. ⁴Departament de Física de la Matèria Condensada, Universitat de Barcelona, Carrer de Martí i Franquès 1, 08028 Barcelona, Spain.</footnote>

<footnote n="5">⁵Dipartimento di Fisica, Universita' di Roma "La Sapienza", P'le Aldo Moro 5, I-00185 Roma, Italy. ⁶Dept. of Materials, Univ. of Oxford, Parks Road, Oxford OX13PH, UK. ⁷Laboratoire de Physique de l'Ecole Normale Supérieure, ENS, Université PSL, CNRS, Sorbonne Université, Université Paris-Diderot, Sorbonne Paris Cité, Paris, France. ⁸Dept. of Chemistry, Colorado State University, Fort Collins, CO 80523, USA. ✉e-mail: walter.schirmacher@uni-mainz.de; giancarlo.ruocco@iit.it</footnote>
value s = 2 has been reported, depending on the quenching protocol$^{42–44}$. In one approach$^{42}$, a fraction of the particles was fixed in space during the quenching process (pinned particles). With increasing fraction of pinned particles s was found to increase continuously from 2 to 4 (and even above). In another approach$^{43,44}$, the glass at T=0 was produced by quenching from a well equilibrated liquid temperature T (the “parental temperature”). Upon increasing T from the (numerical) dynamical arrest temperature T_d, the exponent of the DOS in the resulting glass was reported to decrease continuously from 4 to 2, up to T roughly twice the value of T_d, and to remain Debye-like for higher parental temperatures. It is worth to emphasise that the Debye-like DOS with s = 2 found in$^{43,44}$ has no relation with wave-like excitations, which, in the simulated small systems, can only exist at high enough frequencies.

As pointed out by Paoluzzi et al.⁴⁴, the observed s = 2 may be related to the vicinity of a marginal stability transition$^{22,45,46}$. At this transition, unstable modes would appear if the glass temperature would be slightly increased. In fact, it has been demonstrated within mean-field spin glass theory$^{22,47,48}$ that in the case of marginal stability the DOS obeys an $\omega^2$ law, which is due to a Gaussian-Orthogonal Ensemble (GOE) type or Marchenko-Pastur-type random-matrix statistics (and not due to waves). In models with spatially fluctuating force constants$^{14,16,18,49}$ and in the HET model (spatially fluctuating elastic constants)$^{19,20,28,50}$, marginal stability appears if the amount of negative force or elastic constants approaches a critical threshold value. If this amount is increased beyond this threshold, again, unstable modes appear in the spectrum (with negative eigenvalues). Such unstable spectra, calculated with HET theory, in fact, have been recently used to model the instantaneous spectrum of liquids⁵¹.

Given the observed scenario, one can argue that quenching from a high parental temperature T' the system may reach a situation of marginal stability, in which the type-I non-phononic excitations extend to zero frequency and s takes the value of 2, as predicted by mean field theory and HET for the case of marginal stability. Quenching from a lower T' might enable the liquid to accomodate in a more comfortable situation, in which the type-I modes are confined to the region above a finite frequency, which would be the BP frequency in macroscopically large systems. In small samples, which do not allow for the existence of waves at small frequency, there is now room for the appearance of type-II non-phononic modes. Evaluating the spectrum of these type-II modes, and explaining their origin, is the main scope of the present paper.

Before characterizing the type-II non-phononic excitations, we show that HET predicts a DOS $g(\omega) - \omega^2$ at the marginally stable limit. We then derive a generalization of heterogeneous elasticity theory, by using the continuum limit of the system's Hessian. We unveil the nature of the type-II non-phononic excitations, being irrotational, vortex-like displacement fields, associated with local, frozen-in stresses.

We find a direct relationship between the statistics of the local stresses, governed by small values of the first derivative of the interaction potential, and the DOS of the type-II non-phononic modes. From this relationship it follows that the value s = 4 often found in numerical simulation⁴¹ is the consequence of the cubic smoothing (tapering) of the potential at its cutoff. We tested this with numerical simulations, demonstrating that the value of s is changed if the tapering function is altered. A further consequence of the relation between the internal stresses and the spectrum of the type-II spectrum is a scaling of s = 5 for potentials with a minimum, such as the Lennard-Jones (LJ) potential. This scaling is expected to be generic for systems with both attractive and repulsive interactions, and has been observed recently in simulations of small disordered LJ systems$^{52,53}$.

## Results

### Type-I nonphononic modes and Heterogeneous-Elasticity Theory (HET)

Heterogeneous-elasticity theory (a derivation from the microscopic Hessian and a brief description of the main steps of the theory are given

in paragraph M1 of the Methods section) is based on the assumption of a spatially fluctuating local shear modulus $G(r) = G_0 + \Delta G(r)$. Here $G_0$ is the average of the shear modulus, and the fluctuations $\Delta G(r)$ are supposed to be short-range correlated (see Methods M1). The bulk modulus $K$ is supposed to be uniform. Such fluctuations can be derived by a coarse-graining procedure (see$^{54,55}$ and Methods M1) from the Hessian matrix of the glass, i.e. the second-order, harmonic Taylor coefficients of the total energy. The statistics of the fluctuations has been verified in a simulation of a soft-sphere glass$^{28}$. Using a mean-field theory derived by field-theoretical techniques (self-consistent Born approximation, SCBA), the mesoscopic spatially fluctuating part of the shear modulus is transformed into a complex, frequency dependent self energy $\Sigma(z)$, which is the central quantity for the discussion of the influence of the disorder on the spectrum. Here $z = \lambda + i0 \equiv \omega^2 + i0$ is the spectral parameter. The imaginary part $\Sigma'(\lambda)$ is proportional to $\Gamma(\omega)/\omega$, where $\Gamma(\omega)$ is the sound attenuation coefficient. It is worth to clarify once more that the HET is fully harmonic, the system dynamics is described by the Hessian, no anharmonic terms are present and the low-frequency sound attenuation is due to Rayleigh scattering of waves from the disordered structure.

It has been shown in$^{20}$ that at low frequencies $\Sigma'(\lambda)$ just adds to the spectrum $\rho(\lambda) \equiv g(\omega)/2\omega$ (where $\lambda = \omega^2$ is understood) that in absence of the self energy would be the Debye phonon spectrum. So by definition $\Sigma'(\lambda)$ describes essentially the non-phononic part of the vibrational spectrum.

The disorder-induced frequency dependence of $\Sigma(z)$ is governed by a dimensionless disorder parameter, $\gamma$, which is proportional to the variance of the shear modulus fluctuations:

$$ \gamma = A \frac{\langle (\Delta G)^2 \rangle}{G_0^2}, \qquad (1) $$

where A is a dimensionless factor of order unity. Within the mean-field HET, on changing $\gamma$, there is a sharp crossover from stability (no negative eigenvalues $\lambda$) to instability (presence of negative $\lambda$). Near the instability, which takes place at the critical value $\gamma_c$, the self energy can be represented as$^{50}$:

$$ \Sigma(z) = \Sigma_c \left\{ 1 + \frac{2}{\gamma_c} \sqrt{\gamma_c - \gamma[1 + zG(z)]} \right\} \qquad (2) $$

where $\Sigma_c$ is $\Sigma(z=0)$ at criticality and $G(z)$ is a linear combination of the longitudinal and transverse Green's function (see$^{50}$ and Methods M1). Near $\lambda = 0$, $G(z)$ can be represented as

$$ G(z) = G_0 + iG_1 z^{(d-2)/2}, \qquad (3) $$

where the imaginary part $G''$ gives rise to the Debye spectrum:

$$ \begin{aligned} \rho_{\varepsilon}(\lambda = \omega^2) &= \frac{1}{2\omega} g_{\varepsilon}(\omega) \sim G''(\lambda) \\ &\sim \lambda^{(d-2)/2} = \omega^{d-2}, \end{aligned} \qquad (4) $$

which is the Debye law

$$ g_{\varepsilon}(\omega) \sim \omega^{d-1}. \qquad (5) $$

The subscript $\varepsilon$ of the spectral function $\rho_{\varepsilon}(\lambda)$ and the DOS $g_{\varepsilon}(\omega)$ indicates that these excitations correspond to the HET theory, which is formulated in terms of the strain tensor $\hat{\varepsilon}$.

For small frequencies, $\Sigma'(\lambda)$ can be expanded with respect to $G''(\lambda)$, leading to Rayleigh scattering $\Sigma'(\lambda) - \lambda d^{220,21,56}$. For larger frequencies ($\lambda > \lambda_c - \gamma_c - \gamma$) the square-root in Eq. (2) produces its own imaginary part, leading to a shoulder in the spectrum. This shoulder appears as a maximum in the “reduced DOS” $g(\omega)/\omega^2$ and has been called ‘boson’
peak' for historical reasons³⁴. So the boson peak is the crossover from a Debye-spectrum to a random-matrix spectrum.

If we deal with finite-size samples (N particles), phonons do not exist below the first resonance frequency $\lambda_0 = \omega_0^2 - N^{-2/d}$. In the absence of phonons, $G_0$ has still a finite value, but the low-frequency imaginary part $G''(\lambda)$ is gone. Instead of a boson peak Eq. (2) predicts now a gap in the spectrum. Beyond the gap edge (situated at $\lambda_c = \omega_c^2$) $\Sigma'(\lambda)$ increases as $\sqrt{\lambda - \lambda_c}$. This feature is shared by the mean-field theory of spin glasses²²,⁴⁷,⁴⁸, associated with a Marchenko-Pastur or GOE-type random-matrix spectrum. As stated above, we call these disorder-dominated random-matrix modes "type-I non-phononic excitations". The eigenvalues of these excitations obey random-matrix statistics of the GOE type, as will be demonstrated below.

The type-I modes obviously dominate the low-frequency spectrum of a marginally stable system with $\lambda_c - \gamma_c - \gamma = 0$, leading to $\rho_c(\lambda) - \lambda^{1/2}$ for $\lambda \to 0$. Because this law (which converts to a $g_c(\omega) - \omega^2$ law) is observed for quenching from a high parental temperature $T'$, we conjecture that for such high parental temperatures a marginally stable system is produced by the quenching procedure. This is rather plausible, because the quenching procedure forces the final glass to have only positive eigenvalues of the Hessian, but still keeps the original liquid structure almost unaltered. Thus, the obtained glass inherent structure lies "high" in the Potential Energy Landscape (PEL). This is different in the case of quenching from a low parental temperature, where the initial liquid structure is equilibrated into a lower energy region of the PEL. The distinction between low and high parental temperature, which is fundamental for understanding the properties of simulated glasses⁴²,⁴³, is irrelevant in the case of real glasses: the quenching of a real system is anyway appreciably slower (by many orders of magnitude) than in the case of simulated glasses, which means that during the quenching the system continuously equilibrates at lower and lower temperature. Thus, except for the case of more complex disordered materials like rubbers, gels, foams or granular materials⁵⁷, an atomic or molecular glass is always somewhat away from marginality.

## Type-II non-phononic excitations and generalized heterogeneous-elasticity theory (GHET)

We now turn to the type-II non-phononic excitations, which we conjecture to appear in small samples (no low-frequency phonons) in a more stable situation $\gamma_c - \gamma \gg 0$ reached by quenching from low parental temperatures (no type-I non-phononic excitations at low frequency). In this case, almost all the numerical simulations, independently of the specific interaction potentials, indicate $s = 4^{41}$, although for Lennard-Jones-type potentials the situation is less clear⁴⁰,⁵²,⁵³,⁵⁸⁻⁶⁰. To clarify the origin and the nature of these type-II excitations, which are predicted neither by the HET as described before, nor by the standard mean-field spin glass theory (see, however, Bouchbinder et al.⁶¹), we have developed a generalized version of the HET (let us call it Generalized HET, i.e. GHET) where the displacement field is described by two, instead than by a single, variables. These variables are derived in terms of a modified continuum description, in which, in addition to the usual strain $\tilde{e}(r, t)$ tensor field³⁴,⁶², a second-non-irrotational-vector field $\eta(r, t)$ (vorticity) is necessary for a correct description of the dynamics. The two fields are coupled. While the features of the strains are basically controlled by the space dependence of the fluctuating elastic constants and give rise to the type-I non-phononic modes (predicted by standard HET), the vorticities are associated with spatially fluctuating local stresses $\bar{\sigma}_{ij}^{\alpha\beta}$, which can be written as (see Methods M2)

$$ \sigma_{ij}^{\alpha\beta} = \frac{1}{\Omega_Z} r_{ij}^{\alpha} r_{ij}^{\beta} \phi'(r_{ij}) / r_{ij} = \frac{1}{\Omega_Z} \tilde{e}_{ij}^{\alpha} \tilde{e}_{ij}^{\beta} r_{ij} \phi'(r_{ij}) \quad (6) $$

Here $\phi'(r)$ is the first derivative of the pair potential, $r_{ij} = r_i - r_j$ is the vector between the positions $r_i$ of a pair of interacting particles,

$\hat{e}_{ij} = r_{ij}/r_{ij}$ the corresponding unit vector, $r_{ij} = |r_{ij}|$ is the interparticle distance, and $\Omega_Z$ is a small volume around the center-of gravity vector $R_{ij} = \frac{1}{2}(r_i + r_j)$

As shown in the Methods M2 section, we are able to demonstrate that the dynamics of $\eta(r, t)$ is similar to that of a set of local oscillators, coupled to the strain field. To the best of our knowledge, such vortex-like harmonic vibrational excitations have not yet been considered.

It has been pointed out by Alexander⁵⁷,⁶³ that the term in the potential energy, which involves the local stresses, violates local rotation invariance. In fact, the presence of frozen-in local stresses is the reason for the existence of the type-II excitations. The local stresses also provide the coupling between the type-II excitations and the strain field.

In order to be more specific, but, as well, avoid the introduction of too complicated technicalities, we use a simplified version of GHET (Methods M2), where the non-irrotational fields and the local stresses are treated as scalars. The additional contribution to the self energy is therefore (see Methods M2)

$$ \Sigma_{\epsilon\eta}(z) \sim \frac{(\sigma_{ij}^{\alpha\beta})^2}{-z\zeta + \sigma_{ij}^{\alpha\beta}} = \int d\sigma \mathcal{P}(\sigma) \frac{\sigma^2}{-z\zeta + \sigma} \quad (7) $$

Here $\sigma$ (without indices) denotes any component of the stress tensor $\sigma_{ij}^{\alpha\beta}$, and $\zeta$ is a local moment-of-inertia density (see Methods M2). The overbar indicates the average over the local stresses, and $\mathcal{P}(\sigma)$ is their distribution density. The local stresses may assume both positive and negative values. Because the local stresses include those due to the image particles of the periodic boundary conditions, even for a strictly repulsive potential the diagonal elements $\sigma^{aa}$ may be positive and negative, see Methods M2 and §17 of⁵⁷. Only positive values of $\sigma$ contribute to the spectrum.

The corresponding contribution to the spectrum $\rho_{\epsilon\eta}(\lambda)$ is

$$ \rho_{\epsilon\eta}(\lambda) \sim \Sigma''_{\epsilon\eta}(\lambda) \sim \sigma^2 \mathcal{P}(\sigma) |_{\sigma=\lambda\zeta} \quad (8) $$

Equation (8) represents one of the main results of the present paper. Once the distribution of the local stresses is known, we have an explicit expression for the contribution of the type-II non-phononic modes to the DOS. If, for example, $\mathcal{P}(\sigma) \sim \sigma^{-1/2}$, from Eq. (8) one obtains $\rho_{\epsilon\eta}(\lambda) - \lambda^{3/2}$, which leads to $g_{\epsilon\eta}(\omega) \sim \omega^4$.

The distribution of the quantity $\sigma_{ij}^{\alpha\beta}$ can be easily derived (Methods M3) assuming an isotropic system, and from the knowledge of the particles' radial pair distribution function $g_2(r)$. Indeed, the stress distribution $\mathcal{P}(\sigma)$ is related to the particle distance distribution via (Methods, M3)

$$ \mathcal{P}(\sigma) = r^{d-1} g_2(r) \left| \frac{d\sigma}{dr} \right|^{-1}, \quad (9) $$

where $\sigma(r) = r\phi'(r)/2\Omega_c$.

How to rationalise the (almost) always observed $s = 4$? In the numerical simulations, in order to speed up the calculations and to optimally deal with the periodic boundary conditions, a cut-off radius is introduced, i.e. the interaction potentials is forced to be zero for distances larger than a given interparticle distance $r_c$. Furthermore, to guarantee the stability of the dynamics, and the numerical conservation of the total energy, the interaction potential is adjusted - using a proper tapering function - in such a way to have at least the first two derivatives continuous at $r_c$. Besides details specific to each simulation, the general form of the employed interaction potential is:

$$ \phi_{DM}(r) = (\phi(r) - \phi(r_c)) T_m(r/r_c), \quad (10) $$

where $T_m(x)$ ($m \ge 2$) is the tapering function, which is zero for $x \ge 1$, and whose first $m$ derivatives are also zero at $x = 1$. This guarantees that the
first *m* derivatives of the potential are continuous at *r* = *r*<sub>c</sub>. In other words, the interaction potentials employed in the numerical simulations always vanish at *r*<sub>c</sub> with a power law (*r*<sub>c</sub> − *r*)<sup>(*m* + 1)</sup>. In most simulations the value *m* = 2 is taken, which is the minimum value to guarantee that the first two derivatives of the potential are continuous at *r* = *r*<sub>c</sub>.

Because the DOS of the type-II non-phononic modes is highly sensitive to the shape of the interaction potentials close to the zeros of its first derivative, one can expect that the tapering function has a large effect on the DOS. Indeed (Methods M3), we find that, for any inverse-power law potential, in presence of a tapering function one has

$$P(\sigma) \sim \sigma^{-1+\frac{1}{m}}, \qquad (11)$$

and consequently

$$\begin{align}
\rho_{\epsilon\eta}(\lambda) &\sim \lambda^{1+\frac{1}{m}} \tag{12} \\
g_{\epsilon\eta}(\omega) &\sim \omega^{3+\frac{2}{m}}. \nonumber
\end{align}$$

Thus, in this case, the DOS of the type-II non-phononic modes is controlled by the tapering function, not by the details of the glass or of the interaction potentials. (The almost) universally used *m* = 2 gives for the type-II low-frequency DOS the scaling *g*<sub>*εη*</sub>(*ω*) ~ *ω*<sup>4</sup>, i.e. the *s* = 4 scaling of the low-frequency DOS, which has been observed in many recent small-sample simulations<sup>41</sup>.

In the case of a potential with a smooth, tapered cutoff, the scale of the relevant values of *σ* is determined, among other constants, by *g*<sub>2</sub>(*r*<sub>c</sub>). The latter quantity actually depends on the thermodynamic state of the system and on the quenching protocol, and this observation could explain the findings that the intensity of the type-II non-phononic mode spectrum does actually depend on the glass structure: see Fig. 2 in Ref. 64 and Fig. 10 in Ref. 60.

A final consideration concerns the rather interesting case of potentials with repulsive and attractive parts like the frequently used Lennard-Jones (LJ) potential. In this case the first derivative of the potential – relevant to the stresses–vanishes linearly at the potential minimum. Because *g*<sub>2</sub>(*r*) has a maximum there, the corresponding small stresses will dominate over the tapering-induced ones. At the potential minimum the stress goes to zero, but the second derivative does not, which corresponds to the *m* = 1 tapering, leading to a DOS *g*<sub>*εη*</sub>(*ω*) ~ *ω*<sup>5</sup>. In fact, such a scaling has been recently observed in simulations of small LJ systems<sup>52,53</sup>. These authors thus conjectured that the *s* = 5 scaling is generic to LJ systems. In our view this is true and not only applies to the LJ case but also to all potentials with attractive and repulsive contributions.

In other simulations of small LJ systems the spectrum does not clearly show a simple power law<sup>58–60</sup>. Here we may assume that in addition to the small stresses due to the potential minimum, contributions from the smooth cutoff and/or type-I excitations may be involved, depending on the quenching protocol. We finally mention a simulation of an LJ system<sup>40</sup>, which was large enough to include Debye phonons. In this simulation the use of a hard cutoff at *r*<sub>c</sub> (no tapering, “*m* = 0”) was compared with *m* = 1 tapering. Interestingly only the in *m* = 1 case non-phononic modes on top of the Debye phonons were observed.

## Resulting scenario for the role of nonphononic excitations in glasses

In a nutshell, the scenario described by the HET and by the GHET, distinguishes (i) between (*i*<sub>A</sub>) large, macroscopic, glasses as those investigated experimentally, where the acoustic waves (phonons) extend to zero frequency, and (*i*<sub>B</sub>) small systems investigated numerically by molecular dynamics simulations and similar techniques where the phonons exist only above a certain lowest resonance frequency; (ii) between (ii*<sub>A</sub>*) “quickly” quenched glasses, such as those obtained by quenching from a high parental temperature in numerical simulations

(here the glass is at criticality *y* ≈ *γ*<sub>c</sub>), and (ii*<sub>B</sub>*) “slowly” quenched glasses, obtained by quenching from a well equilibrated low parental temperature, or in the real world, where the quenching rate is by far slower than in any numerical simulations (here the glass is deeply in the stable region *γ*<sub>c</sub> ≫ *y*). Within the present description (in terms of HET and GHET) in all these cases a rather different scenario is realized.

In Fig. 1 we summarize, in a cartoon-like fashion, what is expected from our theory.

In Fig. 1a the case of a macroscopically large system far from criticality is considered. This is the situation (almost) always found in real experiments. In this figure the reduced DOS *g*(*ω*)/*ω*<sup>2</sup> is sketched (blue line), for a 3d system of infinite size (*N* = ∞), far from criticality (*y* ≪ *γ*<sub>c</sub>). The high-frequency side is characterized by Anderson-localized modes (also called “locons” following the notation of Allen et al.<sup>65</sup>). The rest of the spectrum is populated by extended modes. This remaining region is divided into (i) a low-frequency regime, where acoustical, only occasionally scattered waves (“phonons”), as considered by Debye, exist, and (ii) an intermediate regime, where the modes are strongly scattered by the disordered structure. In this strongly-scattering regime the intensity of the waves obeys a diffusion equation akin to light in milky glass<sup>66</sup>. This is why Allen et al.<sup>65</sup> call these excitations “diffusions”.

It has been suggested<sup>14,22</sup> that the states in this region obey the statistics of random-matrix eigenstates.

In Fig. 1b the reduced density of states, *g*(*ω*)/*ω*<sup>2</sup> is depicted for *γ* approaching its critical value *γ*<sub>c</sub>. The blue line represents, as in panel a, the case *γ* ≪ *γ*<sub>c</sub>. The violet, magenta and red lines correspond to increasing disorder, approaching *γ* → *γ*<sub>c</sub>. The BP onset shifts to lower and lower frequencies, eventually reaching zero frequency at *γ* = *γ*<sub>c</sub>. In parallel, the random-matrix type modes, here named ‘type-I non-phononic’ modes, cover more and more the low-frequency region.

Finally, in Fig. 1c we present the case obtained in numerical simulations of very small systems, where a gap opens in the phonon spectrum. The gap edge is located at the lowest resonance frequency *ω*<sub>o</sub>. This frequency corresponds to the frequency of the transverse phonons in the glass at a wavevector *k* = 2π/*L*, where *L* is the box size, which gives *ω*<sub>o</sub> = 2πν<sub>T</sub>/*L*, where ν<sub>T</sub> is the transverse sound velocity. As sketched in the Figure, the frequency *ω*<sub>o</sub> provides the upper limit for the visibility of the type-II spectrum. At higher frequencies these modes hybridize with the waves and probably can no more be distinguished from them.

If the system is quenched from a high parental temperature, i.e. the resulting inherent structure lies high in the potential energy landscape, *γ* ≈ *γ*<sub>c</sub>, the gap in the phonon spectrum exists no more, i.e the non-phononic type-I modes (random matrix eigenstates) extend towards *ω* = 0. The low-frequency DOS in this case scales as *g*(*ω*) ~ *ω*<sup>2</sup> like a *d* = 3 Debye spectrum (this case is not sketched in the figure).

If, on the contrary (as reported in of Fig. 1c), the quenching is performed starting from a well equilibrated (supercooled) liquid configuration at low parental temperature, the inherent structure reaches a low value in the potential energy landscape, and the glass is far from criticality. At low frequencies we have both a gap in the phonon spectrum (because of the small system size) and a gap in the non-phononic type-I modes (because *γ* ≪ *γ*<sub>c</sub>). Within this gap now emerge the otherwise not visible non-phononic type-II modes.

It is worth to note that the non-phononic type-II modes can be accessed in a very specific situation (small system size) that can be reached only in numerical simulations. As mentioned in the introduction, these modes have been the subject of a large amount of work in the last years, almost all of them reporting a density of states *g*(*ω*) ~ *ω*<sup>s</sup> with *s* = 4. The GHET predicts that this *s* value depends on the details of the potential. It turns out that the value of *s* depends on the tapering function used to smooth the potential around the cutoff in the molecular dynamic simulations, such that the first *m*th derivatives are continuous. For *m* = 2, which is almost always used in numerical
<figure><img src="image_4.png" /></figure>

simulations, the prediction is $s=4$, as observed in the literature. More generally, $s=3+2/m$. As for a correct numerical determination of the Hessian (and thus of the dynamics) $m$ must be $m \ge 2$, $s$ cannot be $>4$. It can be reduced to its minimum value, $s=3$, if a $C^\infty$-function is used as tapering function.

There exists a large body of simulations with a potential, which quadratically becomes zero at a radius $r_c$ (see Ref. 67 for references). Such potentials have been used in jamming studies$^{68,69}$. The Hessian of this potential is discontinuous at $r=r_c$. The Hessian, however, on which our analytical theory is based, and from which the DOS is obtained by diagonalization, plays a pivotal role. Therefore we

**Fig. 1 | Reduced density of states in different situations.** a Sketch of the reduced DOS, $g(\omega)/\omega^2$ (blue line), for a 3d system of infinite size ($N=\infty$), far from criticality ($\gamma \ll \gamma_c$). The high-frequency side is characterized by Anderson-like localized modes (also called 'Locons' following the Allen-Feldman$^{65}$ notation), while the rest of the spectrum is populated by extended modes. These are divided into a low-frequency region where phonon(-like) modes exist (briefly 'Phonons') and an intermediate region where the modes are strongly scattered by the disordered structure, and their dynamics is wave diffusion ('Diffusons'). According to the HET, this region starts at the onset of the BP and the modes in this region are of random-matrix type, i.e. similar to eigenstates of random matrices. b Sketch of the reduced density of states, $g(\omega)/\omega^2$, for a $d=3$ system of infinite size ($N=\infty$), at different level of criticality. The blue line represents, as in (a) the case $\gamma \ll \gamma_c$, while the violet, magenta and red lines correspond to cases with $\gamma \to \gamma_c$. The BP onset shifts to lower and lower frequency, eventually reaching zero frequency at $\gamma=\gamma_c$. In parallel, the random-matrix type modes, here named 'type I non-phononic' modes, cover more and more the low-frequency region. c Sketch of the reduced density of states, $g(\omega)/\omega^2$ (violet line), for a $d=3$ system, with finite size ($N$), and far from criticality ($\gamma \ll \gamma_c$). Due to the finite size of the system, the boundary conditions impose a lowest $k$ value, and consequently a lowest phonon frequency $\omega_o$. Below this value, no phonons are allowed. If-as in the sketch-the system is far from criticality, a gap opens at small frequencies (no phonons and no non-phononic type-I modes are allowed in the low frequency region). On the contrary, different numerical simulation studies indicate that in this condition the low-frequency part of the spectrum shows an $\omega^4$ DOS. We have called these excitations "non-phononic type-II modes", and we have shown that their origin can be explained via the generalized HET (see text).

decide not to consider the cases with $m<2$, where some matrix elements of the Hessian are discontinuous at the truncation point.

Finally, as pointed out above, in the case of a potential with repulsive and attractive parts, which has a minimum at the nearest-neighbour distance-like the Lennard-Jones potential - the GHET predicts that the dominant type-II spectrum scales with $s=5$, in agreement with recent simulations$^{52,53}$.

The GHET, therefore, is in agreement with results previously reported in the literature. The universally observed $s=4$ is due to a specific choice of tapering made in the numerical simulations, whereas $s=5$ is generic for potentials with a minimum.

## Numerical simulations

As the GHET gives us precise predictions on the value of $s$ for different tapering shapes, and links this value to the (distribution of the) local stress tensor, we have performed extensive numerical simulations of a model glass to give strength to the theory. We performed the simulations with different tapering functions (specifically, a $C^2$ and a $C^\infty$ function), and measured both the stress distributions and the non-phononic spectra. Both turn out to be in agreement with GHET.

As a microscopic model of glass, we consider a polydisperse mixture of short-range repulsive particles$^{70}$. The diameters $a_i$, with $i=1, \dots, N$ the particle's label, are drawn from a power law distribution $P(a)$ with $\langle a \rangle = \int_{a_{min}}^{a_{max}} da P(a)a=1$, with $a_{min}=0.73$, $a_{max}=1.62$, and $P(a)=N a^{-3}$, with $N$ the proper normalization constant$^{70-76}$. Particles are arranged into a box of side $L$ with periodic boundary conditions. Indicating with $\mathbf{r}_1, \dots, \mathbf{r}_N$ a configuration of the system, the mechanical energy is

$$V(\mathbf{r}_1, \dots, \mathbf{r}_N) = \frac{1}{2} \sum_{i \neq j} \phi(\mathbf{r}_{ij}), \qquad (13)$$

with the pair potential $\phi(\mathbf{r}_{ij})$

$$\begin{align}
\phi(\mathbf{r}_{ij}) &= \left(\frac{a_{ij}}{r_c}\right)^n \varphi_m(x_{ij}) \tag{14} \\
\varphi_m(x) &= [\varphi(x) - \varphi(1)] T_m(x) \nonumber \\
\varphi(x) &= x^{-n}. \nonumber
\end{align}$$
with $x_{ij} = r_{ij}/r_c$, $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$ and $r_c$ is the cutoff distance. In writing (14) we have introduced $a_{ij} \equiv \frac{a_i + a_j}{2}$ and we indicate with $T_m(x)$ the cutoff function (tapering function) which is different from zero only in the interval $x \in [0, 1]$. $T_m(x)$ has the property that it varies as $(1-x)^m$ for $x \to 1$. Consequently the potential goes to zero for $r \to r_c$ as

$$ \phi(r_{ij}) \sim (r_{ij} - r_c)^{m+1} \qquad (15) $$

This guarantees that the first $m$ derivatives of the potential vanish continuously at the cutoff $r_c$. We specialize our numerical study to the case $m=2$ and $m=\infty$ adopting the following functions

$$ T_2(x) = (1 - x^2)^2, \qquad (16) $$

and

$$ T_{\infty}(x) = \exp \left[ \frac{x}{2(x-1)} \right]. \qquad (17) $$

The cutoff distance $r_c$ is chosen to be larger than the largest particle's diameter $\sigma_{max} = 1.62$ (we set $r_c = 1.955$).

The simulated systems are small enough to allow a direct inspection of the non-phononic modes, be it Type-I or Type-II. Actually, we measure the lowest phonon resonance at $\omega_o \approx 1.5$, and modes are found down to a frequency ten times smaller.

To check the character of these modes, i.e. whether they are extended or localized, we use the fact that in a disordered system eigenvalues cannot be degenerate, which leads to level repulsion. In the theory of random matrices<sup>77</sup> this phenomenon is investigated using the statistics of the eigenvalue differences. For delocalized states this statistics is governed by that of the Gaussian Orthogonal Ensemble (GOE), for which the distribution of small eigenvalue differences increases linearly, whereas for localized states one expects a Poissonian (exponentially decaying) distribution.

For evaluating the GOE spectral statistics we use the method reported in<sup>78</sup> where the distribution of the variable $r$ (the ratio of two differences between adjacent eigenvalues) is considered. The variable $r$ is defined as $r = (\lambda_{p+1} - \lambda_p)/(\lambda_{p+2} - \lambda_{p+1})$ ($p = d-1...[dN-2]$), where $\lambda_p$ are the eigenvalues ordered in ascending order.

Atas et al.<sup>78</sup> have shown that for the GOE case the statistics is well obeyed by their surmise

$$ P(r) = \frac{27}{8} \frac{r + r^2}{(1 + r + r^2)^{5/2}}. \qquad (18) $$

Figure 2 shows the distribution $P(r)$ for three different frequency regions: lowest 50 eigenvalues (panel A), middle frequency range (B) and highest 50 eigenvalues of a $d=3$ and $N=1600$ system. The distributions have been averaged over 900 independent samples. The highest frequency modes, as expected, are localized (the dotted line is an exponential fit to the data), while both middle and low frequency modes show their extended character: they are very well represented by the Atas surmise, Eq. (18), as demonstrated by the parameterless full line, which represents the data very well. It is interesting to note that the level distance statistics is universal, i.e. does not depend on the type of non-phononic excitations: type-I in samples prepared from high $T'$ (blue symbols), type-II in samples prepared from low $T'$ (red symbols).

This result confirms that all low-frequency non-phononic modes, type-I and type-II are extended and not localized.

We now turn to the evaluation of the distribution of the local frozen-in stresses. Liquid theory (see Methods M3) predicts the behaviour of the stress distribution $P(\sigma)$ at low $\sigma$ value. Specifically, Eq. (11) predicts a power law $P(\sigma) - \sigma^\alpha$ with exponent $\alpha = \frac{1}{m} - 1$, that is $\alpha = -1/2$ for $m=2$ and $\alpha = -1$ for $m=\infty$, which are the two simulated cases.

<figure><img src="image_5.png" /><figcaption>Fig. 2 | Eigenvalue statistics. Distribution of the quantity $r = (\lambda_{p+1} - \lambda_p)/(\lambda_{p+2} - \lambda_{p+1})$ ($p = d-1...dN-2$). Blue symbols: High parental temperature, red symbols: low parental temperature; Panel A: low-frequency modes (lowest 50 modes); panel B: middle frequency regime; panel C: High-frequency (highest 50) modes. The highest-frequency modes, as expected, are localized (the full line is an exponential fit to the data), while both middle- and low-frequency modes show their extended character, as they are very well represented (full lines) by the Gaussian-Orthogonal-Ensemble surmise, Eq. (18) of Atas et al.<sup>78</sup>. The data have been obtained by evaluating 900 samples.</figcaption></figure>

Figure 3 reports the distributions $P(\sigma)$ in a double-logarithmic scale for the investigated $m=2$ (black) and $m=\infty$ (red) cases. The dashed lines represent the expected slopes. The present results clarify how the local stress distributions in the numerical simulation are controlled by the need of introducing a tapering for giving continuity to the interaction potential and its first two derivatives.

As last step we evaluate the integrated density of states (The data quality of this function is much better than that of $g(\omega)$, evaluated by a binning procedure.). from the numerical simulations:

$$ \mathcal{F}(\omega) = \int_{0}^{\omega} d\omega' g(\omega'). \qquad (19) $$

It is clear that for a DOS with $g(\omega) \sim \omega^s$ ($s > -1$) the integrated DOS behaves as $\mathcal{F}(\omega) \sim \omega^{s+1}$. The results are reported in Fig. 4. The filled symbols correspond to the “usual” simulations: the tapering function is chosen to ensure the continuity up to the second derivative ($m=2$, Eq. (16)) of the interaction potential, while the open symbols correspond to $m=\infty$ with tapering function (17). In the top panel (A), the quenching is performed starting at a high parental temperature $T'$, i.e. the
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Stress distributions.** Distribution of the (modulus of) the off-diagonal stress components. (A) blue dots, case *m* = 2. The bin size is 10⁻⁴. The dashed line indicates the slope 1/2, as expected according to Eq. (11) for *m* = 2. B red dots, case *m* = ∞. The bin size is 10⁻⁵, smaller than in the previous case as for *m* = ∞ there are much more small stress values. In this case the dashed line indicates the slope 1, again expected according to Eq. (11) for *m* = ∞.</figcaption></figure>

equilibration of the liquid is at $T \gg T_d$. The system is then rapidly quenched, i.e. it does not have enough time to equilibrate the glass structure, which, therefore, is close to marginality. Consequently, as already noticed in literature<sup>4,23</sup> and as predicted by HET, the DOS follows a $\omega^2$, with $s=2$, scaling, which is independent of the potential tapering.

If, on the contrary, the parental temperature $T$ is lowered, and approaches $T_d$, the system is far from marginality, a gap opens in the non-phononic type-I spectrum, and the type-II modes now dominate the low-frequency spectrum. This spectrum is predicted by GHET to be affected by the way the potential is tapered. For better transparency we collected the GHET predictions for different values of the tapering index $m$ in Table 1. Indeed, as reported in the bottom panel (B), we observe a different low-frequency behavior for the two different tapering versions: for the $m=2$ tapering $s=4$ is reached at low frequencies, and for the $m=\infty$ tapering $s=3$, as predicted by GHET for the DOS of the type-II modes.

In other words: The main result of Fig. 4 is the *insensitivity* of the DOS on the tapering exponent $m$ for high parental temperature $T$, as opposed to a *sensitivity* on $m$ for low $T$. The observed low-frequency exponents compare well with the predictions of HET and GHET: $s=2$ for high $T$, signifying type-I excitations at marginality, and $s=4, 3$ for low $T$, corresponding to type-II excitations, which depend on the stress distributions for tapered potentials with $m=2, \infty$.

## Discussion

We have used the heterogeneous elasticity theory (HET) and its generalization (GHET) in order to clarify the behaviour of harmonic vibrational excitations in disordered systems.

In the present study we clarify the interplay between system stability (controlled by the disorder parameter $\gamma$, i.e. the normalized variance of the elastic constant heterogeneity) and the system size in determining the DOS shape at low frequency.

As predicted by HET, in large macroscopic systems (as real glasses are), the disorder parameter controls the frequency position $\omega_{BP}$ of the boson peak (BP), and $\omega_{BP}$ marks the frequency border between the phonon realm at low frequency and that of the non-phononic type-I excitations at $\omega > \omega_{BP}$. If the disorder is increased, i.e. the stability is decreased, the system approaches the marginal state. Here the BP

<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Dependence of the spectra for different tapering versions.** Numerical results for the integrated DOS $F(\omega) = \int_0^\omega d\omega' g(\omega')$ of a 'small' (N=1000) system obtained for a high ($T \gg T_d$, panel A) and low ($T \simeq T_d$, B) parental temperature, where $T_d$ is the temperature of dynamical arrest. The full symbols report the 'usual' tapering case ($m=2$), while the open ones correspond to $m = \infty$, Eq. (17), introduced in this work, in order to test the GHET. For high parental temperature (A) the low-frequency slope of the DOS $s=2$ for both tapering versions, corresponding to the marginal type-I spectrum, while at low parental temperature (B) where a gap opens in the type-I-mode spectrum, the DOS is due to the type-II excitations, which exhibit an $m$ dependent low-frequency scaling: $s=4$ for $m=2$ and $s=3$ for $m = \infty$.</figcaption></figure>

<table><caption>**Table 1 | Spectral exponents for the different tapering schemes**</caption><thead><tr><th>Quantity</th><th>Scaling</th><th>m=1</th><th>m=2</th><th>m=∞</th></tr></thead><tbody><tr><td>P(σ)</td><td>~σ<sup>1+1/m</sup></td><td>σ<sup>0</sup></td><td>σ<sup>-1/2</sup></td><td>σ<sup>-1</sup></td></tr><tr><td>ρ<sub>En</sub>(λ)</td><td>~λ<sup>1+1/m</sup></td><td>λ<sup>2</sup></td><td>λ<sup>3/2</sup></td><td>λ<sup>1</sup></td></tr><tr><td>g<sub>En</sub>(ω)</td><td>~ω<sup>3+2/m</sup></td><td>ω<sup>5</sup></td><td>ω<sup>4</sup></td><td>ω<sup>3</sup></td></tr><tr><td>F<sub>En</sub>(ω)</td><td>~ω<sup>4+2/m</sup></td><td>ω<sup>6</sup></td><td>ω<sup>5</sup></td><td>ω<sup>4</sup></td></tr></tbody></table>

We show the results for $P(\sigma)$, the level density $\rho(\lambda)$, the DOS $g(\omega)$ and the integrated DOS $F(\omega) = \int_0^\omega d\omega' g(\omega')$ for the case of a potential with cutoff and tapering function of order $m$.

moves to zero frequency, the phonons disappear, and the type-I modes dominate the whole frequency region.

When the system size becomes small enough, a new frequency scale enters into the game: the lowest phonon resonance, $\omega_o$, which fixes the lower frequency boundary of waves allowed by the system ($\omega > \omega_o$). If the system is marginally stable, or close to marginal stability, below $\omega_o$ the type-I excitations dominate the scene. Their DOS is – like that of the phonons in $d=3$ – proportional to $\omega^2$ in any dimension $d$, although with a different prefactor. If, on the contrary, the glass is highly stable ($\gamma \ll \gamma_c$) a gap opens in the type-I non-phononic excitations spectrum: the type-I excitations only exists for $\omega > \omega_c$. It is in these conditions (stable small glass) that the type II non-phononic excitations can show up below both $\omega_c$ and $\omega_o$.

In other words—if the overwhelming phonon (Propagon) modes, whose DOS is $g(\omega) \approx \omega^{d-1}$, are removed from the low-frequency realm
by numerically investigating small systems – in addition to the type-I modes situated above the BP frequency – a new kind of excitations emerges from the background. It is worth to note that type-II non-phononic excitations are always present, but they are hidden by the presence of waves or by the presence of the type-I modes. When the waves are “removed” (small systems) and the type-I modes are removed as well (well equilibrated systems), the type-II non-phononic modes become visible.

At this point we note that that the local stresses, which give rise to the type-II modes involve a length scale of a few interatomic distance, whereas the heterogeneous elasticity, responsible for the type-I modes involves the mesoscopic length scale of the coarse-graining procedure. To investigate the difference of the two types of non-phononic excitations in more detail is an interesting task for future research. Similarly, it will be very interesting to investigate, how the local stresses and the associated type-II modes enter into the vibrational spectrum of macroscopic samples, and, in particular, how they may influence the boson peak.

The depicted scenario for small samples allows to clarify recent numerical findings, where a $\omega^s$ DOS has been observed, with s depending on the specific parental temperature: s = 4 for minima deep in the potential energy landscape (type-II modes in well equilibrated small system), and s → 2 on increasing the potential energy of the minimum (type-I modes in badly equilibrated small system).

After having clarified the scenario of the excitations existing in glasses, we have investigated the nature of the type-II non-phononic ones. By performing a proper continuum limit for the harmonic dynamics of the glass, we found that in addition to the strain field considered in heterogeneous-elasticity theory (HET), one has to take the non-irrotational part of the atomic displacement pattern into account. The dynamics of these local vorticities are governed by the statistics of the local frozen-in stresses. The local vorticities act as local oscillators, hybridized with the waves and the type-I modes.

The generalization of HET (GHET) predicts that the scaling of the DOS of type-II non-phononic modes at low frequency is related to the distribution of the local stresses. The low-frequency scaling of the type-II excitations is predicted to be non-universal and related via the stress distributions to the statistics of the interatomic forces, giving a one-to-one relation between the potential statistics and the spectrum. These predictions are verified in our soft-sphere simulation. Both GHET and the simulation show that in simulations type-II modes are governed by the type of the smooth cutoff (tapering), employed in the simulations, i.e. they are not an intrinsic glass property.

For potentials with a minimum near the first-nearest-neighbor potential GHET predicts a scaling with s = 5, which is instead a generic property of the system.

At the end we would like to emphasize that the present theory, along with the previous numerical results, obtained by diagonalizing the Hessian matrix (e.g. 38,38,39,41-44,67), is inherently harmonic. Anharmonic effects, which certainly contribute to experimentally measured spectra (e.g. 4-9) and spectra obtained by evaluating correlation functions, using molecular-dynamics simulations (e.g. 25-28), are not taken into account. The harmonic spectra pertain to zero temperature, while the anharmonic effects vanish in this limit. However, the experimentally measured non-phononic spectra in the THz range (boson peak) are reportedly dominated by the harmonic interaction 10,65.

## Methods

### Heterogeneous elasticity theory

Here we review heterogeneous-elasticity theory (HET) 19,28,34, which accounts for the type-I nonphononic excitations.

Heterogeneous elasticity theory is just ordinary elasticity theory with a spatially fluctuating shear modulus $G(\mathbf{R})$. The equations of

motion for the displacement fields $u(\mathbf{R}, t)$ are

$$ \rho_m \ddot{u}^\alpha(\mathbf{R}, t) = \sum_{\beta} \frac{\partial}{\partial \beta} \sum_{\gamma \delta} C^{\alpha\gamma\delta\beta}(\mathbf{R}) \varepsilon^{\beta\gamma}(\mathbf{R}, t) \quad (20) $$

Within HET one discards the non-isotropic terms and assumes that the bulk modulus $K$ does not exhibit spatial fluctuations$^{28,34}$. This gives

$$ \rho_m \ddot{u}^\alpha(\mathbf{R}, t) = K \frac{\partial}{\partial x_\alpha} \mathrm{Tr}[\varepsilon(\mathbf{R}, t)] + \sum_{\beta} \frac{\partial}{\partial x_\beta} 2G(\mathbf{R}) \varepsilon^{\alpha\beta}(\mathbf{R}, t) \quad (21) $$

where we have defined the traceless strain tensor

$$ \hat{\varepsilon}^{\alpha\beta}(\mathbf{R}, t) \doteq \varepsilon^{\alpha\beta}(\mathbf{R}, t) - \frac{1}{3} \mathrm{Tr}[\varepsilon(\mathbf{R}, t)] \quad (22) $$

Fourier-transformed into the frequency regime, Eq. (21) takes the form

$$ \begin{aligned} & -\omega^2 \rho_m u^\alpha(\mathbf{R}, \omega) \\ &= K \frac{\partial}{\partial x_\alpha} \mathrm{Tr}[\varepsilon(\mathbf{R}, \omega)] + \sum_{\beta} \frac{\partial}{\partial x_\beta} 2G(\mathbf{R}) \hat{\varepsilon}^{\alpha\beta}(\mathbf{R}, \omega) \end{aligned} \quad (23) $$

The stochastic Helmholtz equation (23) can be solved with the help of a mean-field theory, the self-consistent Born approximation$^{19,20,28,34}$ (SCBA). The SCBA arises as a saddle point within a replica-field theoretic treatment. The solution amounts to replacing the true medium with spatially fluctuating shear modulus $G(\mathbf{r}) = G_0 - \Delta G(\mathbf{r})$ by a complex, frequency-dependent one $G(z)$ with $z = \omega^2 + i0$. The input is the mean value $G_0 = \langle G(\mathbf{r}) \rangle$ and the variance $\langle (\Delta G)^2 \rangle$, and the output is $G(z) = G_0 - \Sigma(z)$, where $\Sigma(z)$ is the self energy.

The SCBA equations for $\Sigma(z)$ are formulated in terms of the longitudinal ($L$) and transverse ($T$) Green's functions, which are given by

$$ G_{(L,T)}(k, z) = \frac{1}{-z + k^2 v_{(L,T)}^2(z)} \quad (24) $$

As usual one has to add an infinitesimal positive imaginary part to the spectral variable $\lambda = \omega^2$: $z \doteq \lambda + i0$. The quantities $v_L(z)$, $v_T(z)$ are effective, complex, frequency-dependent sound velocities, given by

$$ v_L^2(z) = \frac{1}{\rho_m} \left[ K + \frac{4}{3} (G_0 - \Sigma(z)) \right] \quad (25) $$

$$ v_T^2(z) = \frac{1}{\rho_m} [G_0 - \Sigma(z)] \quad (26) $$

Here $K$ the bulk modulus that we assume to be uniform along the sample. This expression holds in $d=3$, in other spatial dimensions the factor $4/3$ in the RHS of Eq. (25) will be different. The quantity $M = K + 4/3G_0$ is the mean longitudinal modulus, thus we can define $M(z) = K + 4/3(G_0 - \Sigma(z))$ as the generalised longitudinal modulus. The self consistent HET-SCBA equation for the self energy $\Sigma(z)$ turns out to be$^{79}$:

$$ \Sigma(z) = \int \frac{d^3 \mathbf{k}}{(2\pi)^3} C(k) k^2 \left( \frac{2}{3} G_L(k, z) + G_T(k, z) \right) \quad (27) $$

Here $C(k)$ is the Fourier transform of the correlation function $C(r)$

$$ C(r) = \langle \Delta G(r_0 + \mathbf{r}) \Delta G(r_0) \rangle = \langle (\Delta G)^2 \rangle f(r) \quad (28) $$
of the fluctuations $\Delta G(r) = G(r) - \langle G \rangle$ of the elastic shear moduli. We assume that these correlations are short-ranged, i.e. $f(k=0)$ is supposed to be finite. Because the inverse correlation length $\xi^{-1}$ acts effectively as an ultraviolet cutoff, we schematically use

$$f(k) = f_0 \theta(k_\xi - k), \quad (29)$$

where $\theta(x)$ is the Heaviside step function and $k_\xi$ is inversely proportional to the correlation length. From the condition

$$1 = f(r=0) = \int \frac{d^3 k}{(2\pi)^3} f(k) \quad (30)$$

we obtain

$$f_0 = \frac{6\pi^2}{k_\xi^3} \quad (31)$$

We now introduce dimensionless variables: $q = k/k_\xi, \tilde{z} = z^2 \rho_m / G_0 k_\xi^2$, $\tilde{\Sigma} = \Sigma / G_0$, $\tilde{K} = K / G_0$, $\tilde{M}(\tilde{z}) = \tilde{K} + \frac{4}{3}[1 - \tilde{\Sigma}(\tilde{z})]$, $\tilde{G}_{(L,T)}(q, \tilde{z}) = G_{(L,T)}(k, \tilde{z})$. We also define the dimensionless disorder parameter:

$$\gamma \doteq \frac{1}{G_0^2} \langle (\Delta G)^2 \rangle. \quad (32)$$

In terms of these quantities the self-consistent equation (27) takes the form

$$\begin{aligned} \tilde{\Sigma}(\tilde{z}) &= 3\gamma \int_0^1 dq \, q^4 \left( \frac{2}{3} \tilde{G}_L(q, \tilde{z}) + \tilde{G}_T(q, \tilde{z}) \right) \\ &= \gamma \left( \frac{2}{3} \frac{1}{\tilde{M}(\tilde{z})} [1 + \tilde{z} \tilde{G}_L(\tilde{z})] + \frac{1}{1 - \tilde{\Sigma}(\tilde{z})} [1 + \tilde{z} \tilde{G}_T(\tilde{z})] \right) \end{aligned} \quad (33)$$

with

$$\tilde{G}_{(L,T)}(\tilde{z}) \doteq 3 \int_0^1 dq \, q^2 \tilde{G}_{(L,T)}(q, \tilde{z}) \quad (34)$$

The DOS $g_\epsilon(\omega)$ and the level density $\rho_\epsilon(\tilde{z})$ are given by

$$g_\epsilon(\omega) = 2\omega \rho_\epsilon(\tilde{z}) = \frac{2\omega}{3\pi} \text{Im}\{\tilde{G}_L(z) + 2\tilde{G}_T(z)\}. \quad (35)$$

We introduced the subscript $\epsilon$ in order to distinguish the HET-type-I spectra from those due to the type-II excitations, which we call $g_{en}(\omega)$ and $\rho_{en}(\lambda)$.

We now multiply by the factor $1 - \tilde{\Sigma}(\tilde{z})$ and define a function $\tilde{B}(\tilde{z}) = \frac{2}{3}[1 - \tilde{\Sigma}(\tilde{z})]/\tilde{M}(\tilde{z})$ and obtain

$$\tilde{\Sigma}(\tilde{z}) [1 - \tilde{\Sigma}(\tilde{z})] = \left( [1 + \tilde{B}(\tilde{z})] + \tilde{z} \tilde{G}_0(\tilde{z}) \right) \quad (36)$$

with

$$\tilde{G}_0(\tilde{z}) \doteq \gamma (\tilde{G}_T(\tilde{z}) + B(\tilde{z}) \tilde{G}_L(\tilde{z})) \quad (37)$$

Because $\tilde{B}(\tilde{z})$ depends only weakly on $\tilde{z}$, for small $\tilde{z}$ we may set $\tilde{z}$ in this function equal to zero. This gives

$$\tilde{\Sigma}(\tilde{z}) [1 - \tilde{\Sigma}(\tilde{z})] = \left( \underbrace{\gamma [1 + \tilde{B}(0)]}_{\tilde{\gamma}} + \tilde{z} \tilde{G}_0(\tilde{z}) \right) \quad (38)$$

we may now solve this quadratic equation to obtain

$$\tilde{\Sigma}(\tilde{z}) = \tilde{\Sigma}_c(0) + \sqrt{\tilde{\gamma}_c - \tilde{\gamma} - \tilde{z} \tilde{G}_0(\tilde{z})} \quad (39)$$

with $\tilde{\gamma}_c = \frac{1}{4}$ and $\tilde{\Sigma}_c(0) = \tilde{\Sigma}(0)|_{\tilde{\gamma}=\tilde{\gamma}_c} = \frac{1}{2}$. For $\tilde{\gamma} > \tilde{\gamma}_c$ obviously $\Sigma(0)$ becomes imaginary, which introduces a finite DOS for $\tilde{z} \le 0$, i.e. an instability. This instability is due to the presence of too many negative local shear moduli, which in SCBA have been assumed to obey Gaussian statistics.

In very small samples the integrals in Eqs. (27) and (34) have to be replaced by discrete sums over wavevectors $k = (k_x, k_y, k_z)$ with $k_\alpha = 2\pi n_\alpha / L$, $n_\alpha \in \mathbb{Z}$. For small frequencies the resonance conditions $\tilde{\omega} = kv_{(L,T)}(0)$ in the denominator of the Green's functions cannot be met (the sample is too small for carrying long-wavelength waves), and the Green's functions $G_{(L,T)}$, and hence $G_0$ become real quantities. This, in turn leads to a gap in the DOS for $\tilde{\gamma} < \tilde{\gamma}_c$:

$$\rho_{\epsilon}(\tilde{\lambda}) \begin{cases} = 0 & \text{for } \tilde{\lambda} < \tilde{\lambda}_c \\ \sim \sqrt{\tilde{\lambda} - \tilde{\lambda}_c} & \text{for } \tilde{\lambda} \ge \tilde{\lambda}_c \end{cases} \quad (40)$$

with $\tilde{\lambda}_c = [\tilde{\gamma}_c - \tilde{\gamma}]/G_0(0)$. At marginal stability $\tilde{\lambda}_c = 0$ this leads to $\rho(\tilde{\lambda}) \sim \tilde{\lambda}^{1/2}$ and consequently to a DOS $g(\omega) - \omega^2$.

In macroscopically large, stable samples the gap is filled by the Debye phonons. The onset of the imaginary part of $\tilde{\Sigma}(z)$ at $\tilde{\lambda}_c$ corresponds to the boson peak$^{19,20,28,34}$, and the states above $\tilde{\lambda}_c$ are irregular random-matrix states.

## Derivation of Generalized Heterogeneous-Elasticity Theory (GHET)

We confine our treatment to three dimensions; the generalization to two dimensions is straightforward.

We start with the harmonic expansion of the total energy

$$\begin{aligned} E &= T + V(\mathbf{r}_1, \dots, \mathbf{r}_N) \\ &= T + V_0 - \sum_i \mathbf{f}_i \cdot \mathbf{u}_i + \frac{1}{2} \sum_{ij} \mathbf{u}_i \cdot \hat{\mathbf{H}}_{ij} \cdot \mathbf{u}_j \end{aligned} \quad (41)$$

Here $\mathbf{u}_i$ are infinitesimal displacements from the equilibrium position $\mathbf{r}_i$, the $\mathbf{f}_i$ are local forces, defined by

$$f_i^\alpha = - \frac{\partial}{\partial \mathbf{r}_i^\alpha} V(\mathbf{r}_1, \dots, \mathbf{r}_N), \quad (42)$$

and $\hat{\mathbf{H}}_{ij}$ is the Hessian matrix

$$H_{ij}^{\alpha\beta} = \frac{\partial}{\partial \mathbf{r}_i^{\alpha}} \frac{\partial}{\partial \mathbf{r}_j^{\beta}} V(\mathbf{r}_1, \dots, \mathbf{r}_N) \quad (43)$$

At equilibrium the local forces $\mathbf{f}_i$ are zero, so we discard them from the discussion. We now assume that $V(\mathbf{r}_1, \dots, \mathbf{r}_N)$ can be represented as a sum over pairs with a pair potential $\phi(r)$

$$V(\mathbf{r}_1, \dots, \mathbf{r}_N) = \sum_{(i,j)} \phi(r_{ij}) \quad (44)$$

where $\sum_{(i,j)} = \frac{1}{2} \sum_{ij} \phi$ denotes a sum over pairs $(i,j)$, and $r_{ij} = |\mathbf{r}_i - \mathbf{r}_j|$. We further assume that $\phi(r)$ has a finite range $r_c$.
In terms of the pair potential the Hessian $\tilde{H}_{ij}$ can be represented

by

vector

$$H_{ij}^{\alpha\beta} = -K_{ij}^{\alpha\beta}(1 - \delta_{ij}) + \left(\sum_{l \neq i} K_{il}^{\alpha\beta}\right) \delta_{ij} \quad (45)$$

with the force constants

$$K_{ij}^{\alpha\beta} = \left[ \phi''(r_{ij}) - \frac{1}{r_{ij}}\phi'(r_{ij}) \right] \frac{r_{ij}^{\alpha} r_{ij}^{\beta}}{r_{ij}^{2}} + \frac{1}{r_{ij}}\phi'(r_{ij})\delta_{\alpha\beta} \quad (46) \\ \doteq \phi^{(1)}(r_{ij}) r_{ij}^{\alpha} r_{ij}^{\beta} + \phi^{(2)}(r_{ij}) \delta_{\alpha\beta}$$

with implicit definition of the functions $\phi^{(1)}(r_{ij})$ and $\phi^{(2)}(r_{ij})$.

With these definitions the harmonic part of the potential energy can be represented by

$$V_H(\mathbf{r}_1, \dots, \mathbf{r}_N) = \sum_{\langle i,j \rangle} V_{ij} \quad (47)$$

with

$$V_{ij} = \frac{1}{2} \sum_{\alpha\beta} u_{ij}^{\alpha} K_{ij}^{\alpha\beta} u_{ij}^{\beta} \\ = \frac{1}{2} \left( \phi^{(1)}(r_{ij}) \sum_{\alpha\beta} r_{ij}^{\alpha} r_{ij}^{\beta} u_{ij}^{\alpha} u_{ij}^{\beta} + \phi^{(2)}(r_{ij}) u_{ij}^{2} \right) \quad (48) \\ \doteq V_{ij}^{(1)} + V_{ij}^{(2)}.$$

Here $u_{ij}^{\alpha} = u_{i}^{\alpha} - u_{j}^{\alpha}$.

We now go over to a continuum description by introducing difference and center-of-mass variables $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$ and $\mathbf{R}_{ij} = \frac{1}{2}[\mathbf{r}_i + \mathbf{r}_j]$ and interpreting $\mathbf{R}_{ij} \doteq \mathbf{R}$ as the local vector of the continuum theory. We make a Taylor expansion of $\mathbf{u}(\mathbf{r}_i)$ around $\mathbf{R}$

$$\begin{aligned} \mathbf{u}(\mathbf{r}_i) &= \mathbf{u}(\mathbf{R}) + ([\mathbf{r}_i - \mathbf{R}] \cdot \nabla) \mathbf{u}(\mathbf{R}) \\ &= \mathbf{u}(\mathbf{R}) + \frac{1}{2} (\mathbf{r}_{ij} \cdot \nabla) \mathbf{u}(\mathbf{R}) \end{aligned} \quad (49)$$

It follows

$$\frac{1}{2} [\mathbf{u}(\mathbf{r}_i) + \mathbf{u}(\mathbf{r}_j)] = \mathbf{u}(\mathbf{R}) \quad (50)$$

and

$$u^{\alpha}(\mathbf{r}_i) - u^{\alpha}(\mathbf{r}_j) = \sum_{\gamma} r_{ij}^{\gamma} u_{\alpha|\gamma}(\mathbf{R}) = u_{ij}^{\alpha} \quad (51)$$

with abbreviation $u_{\alpha|\gamma} \doteq \partial_{\gamma} u^{\alpha}$.

We now introduce symmetrized and antisymmetrized spatial derivatives

$$\begin{aligned} \varepsilon^{\alpha\gamma}(\mathbf{R}) &= \frac{1}{2} [u_{\gamma|\alpha}(\mathbf{R}) + u_{\alpha|\gamma}(\mathbf{R})] \\ \eta^{\alpha\gamma}(\mathbf{R}) &= \frac{1}{2} [u_{\gamma|\alpha}(\mathbf{R}) - u_{\alpha|\gamma}(\mathbf{R})] \end{aligned} \quad (52)$$

or

$$\begin{aligned} u_{\alpha|\gamma}(\mathbf{R}) &= \varepsilon^{\alpha\gamma}(\mathbf{R}) - \eta^{\alpha\gamma}(\mathbf{R}) \\ u_{\gamma|\alpha}(\mathbf{R}) &= \varepsilon^{\alpha\gamma}(\mathbf{R}) + \eta^{\alpha\gamma}(\mathbf{R}) \end{aligned} \quad (53)$$

Here $\varepsilon^{\alpha\gamma}(\mathbf{R})$ are local strains and $\eta^{\alpha\gamma}(\mathbf{R})$ are local vorticities. Because this tensor has only three independent entries, one can define a vorticity

$$\boldsymbol{\eta}(\boldsymbol{R}) = \frac{1}{2} [\nabla \times \mathbf{u}(\boldsymbol{R})] = \begin{pmatrix} \eta^{yz}(\boldsymbol{R}) \\ \eta^{zx}(\boldsymbol{R}) \\ \eta^{xy}(\boldsymbol{R}) \end{pmatrix} \quad (54)$$

We define now the symmetric part of the local (Born-type) elastic constants as

$$B^{\alpha\beta\gamma\delta}(\mathbf{R}) \doteq \frac{1}{\Omega_Z} \phi^{(1)}(r_{ij}) r_{ij}^{\alpha} r_{ij}^{\beta} r_{ij}^{\gamma} r_{ij}^{\delta} \bigg|_{\mathbf{R}=\mathbf{R}_{ij}} \quad (55)$$

and the stress tensor as

$$\sigma^{\gamma\delta}(\mathbf{R}) \doteq \frac{1}{\Omega_Z} \phi^{(2)}(r_{ij}) r_{ij}^{\gamma} r_{ij}^{\delta} \bigg|_{\mathbf{R}=\mathbf{R}_{ij}} \quad (56)$$

Here $1/\Omega_Z N(Z-1)/2\Omega$ is the average number of pairs within the interaction range per volume, where $\Omega$ is the total volume and $Z$ is the average number of neighbors within the interaction range.

With these definitions we can define a local potential energy density as

$$\nu(\mathbf{R}) = \nu^{(1)}(\mathbf{R}) + \nu^{(2)}(\mathbf{R}) \quad (57)$$

with

$$\begin{aligned} \nu^{(1)}(\mathbf{R}) &\doteq \frac{1}{\Omega_Z} V_{ij}^{(1)} \bigg|_{\mathbf{R}=\mathbf{R}_{ij}} \\ &= \frac{1}{2} \sum_{\alpha\beta\gamma\delta} B^{\alpha\beta\gamma\delta}(\mathbf{R}) \varepsilon^{\alpha\gamma}(\mathbf{R}) \varepsilon^{\beta\delta}(\mathbf{R}) \end{aligned} \quad (58)$$

$$\begin{aligned} \nu^{(2)}(\mathbf{R}) &\doteq \frac{1}{\Omega_Z} V_{ij}^{(2)}(\mathbf{R}) \bigg|_{\mathbf{R}=\mathbf{R}_{ij}} \\ &= \frac{1}{2} \sum_{\alpha\gamma\delta} \sigma^{\gamma\delta}(\mathbf{R}) u_{\alpha|\gamma}(\mathbf{R}) u_{\alpha|\delta}(\mathbf{R}) \end{aligned} \quad (59)$$

As noted by Alexander$^{57,63}$, the term $\nu^{(2)}$ violates local rotation invariance. Therefore in the literature$^{55,62}$ a symmetrization procedure has been applied to this term and the symmetrized non-irrotational term is then incorporated into the symmetric potential energy density, that we define $\tilde{\nu}^{(1)}$, in the following way:

$$\tilde{\nu}^{(1)}(\mathbf{R}) \doteq \frac{1}{2} \sum_{\alpha\beta\gamma\delta} C_{ij}^{\alpha\beta\gamma\delta}(\mathbf{R}) \varepsilon^{\alpha\gamma}(\mathbf{R}) \varepsilon^{\beta\delta}(\mathbf{R}) \quad (60)$$

with

$$\begin{aligned} C_{ij}^{\alpha\beta\gamma\delta} &\doteq B_{ij}^{\alpha\beta\gamma\delta} + \frac{1}{4} (\sigma^{\gamma\delta}\delta_{\alpha\beta} + \sigma^{\alpha\delta}\delta_{\gamma\beta} + \sigma^{\gamma\delta}\delta_{\alpha\delta} + \sigma^{\alpha\beta}\delta_{\gamma\delta}) \\ &\doteq B_{ij}^{\alpha\beta\gamma\delta} + \text{Sy}\{\sigma^{\gamma\delta}\delta_{\alpha\beta}\} \end{aligned} \quad (61)$$

Here we omitted the arguments $(\mathbf{R})$ for brevity, which we do from now on. The remaining part of the irrotational term is then

$$\begin{aligned} \tilde{\nu}^{(2)} &\doteq \frac{1}{2} \sum_{\alpha\gamma\delta} \sigma^{\gamma\delta} \eta^{\alpha\gamma} (\eta^{\alpha\delta} - \varepsilon^{\alpha\delta}) \\ &\doteq \tilde{\nu}_{\eta\eta}^{(2)} + \tilde{\nu}_{\eta\varepsilon}^{(2)} \end{aligned} \quad (62)$$
In terms of the vorticity vector defined in Eq. (54) the two non-irrotational parts of the potential-energy density take the form

$$ \tilde{V}_{[\eta\eta]ij}^{(2)} = \frac{1}{2} \left( \mathrm{Tr}\{\sigma_{ij}\}\eta^2 - \sum_{\gamma\mu} \sigma_{ij}^{\gamma\delta} \eta^{\gamma} \eta^{\delta} \right) \quad (63) $$

$$ \tilde{V}_{[\eta\epsilon]ij}^{(2)} = \frac{1}{2} \boldsymbol{\tau} \cdot \boldsymbol{\eta} \quad (64) $$

with the coupling vector

$$ \boldsymbol{\tau} = \begin{pmatrix} \sum_{\gamma} (\sigma^{\gamma\gamma} \epsilon^{\gamma z} - \epsilon^{\gamma\gamma} \sigma^{\gamma z}) \\ \sum_{\gamma} (\sigma^{\gamma z} \epsilon^{\gamma x} - \epsilon^{\gamma z} \sigma^{\gamma x}) \\ \sum_{\gamma} (\sigma^{\gamma x} \epsilon^{\gamma y} - \epsilon^{\gamma x} \sigma^{\gamma y}) \end{pmatrix} \quad (65) $$

It is now obvious that the vorticity field $\eta(\mathbf{R})$ is exclusively associated with the local stresses $\vec{\sigma}$. Because these stresses are spatially bounded defects in the glass$^{57,80}$, we now give them the “defect label” $\ell$, and call the associated vorticity fields $\eta_{\ell}(\mathbf{R})$. We also treat local stresses as traceless, and their sign to be independent of the sign of the potential derivative $\phi'(\mathbf{r})$. The latter property is due to the external pressure exerted to the system by the boundary condition.

The Lagrangian density in terms of the two vector variables $\mathbf{u}(\mathbf{R}, t)$ and $\eta_{\ell}(\mathbf{R}, t)$ is given by

$$ \mathcal{L}_{\text{GHET}} = \mathcal{L}_{\text{HET}} + \Delta\mathcal{L}_{\text{GHET}} \quad (66) $$

with

$$ \mathcal{L}_{\text{HET}} = \frac{\rho}{2} [\dot{u}(\mathbf{R})]^2 - V^{(1)} \quad (67) $$

where $\rho$ is the mass density, and

$$ \Delta\mathcal{L}_{\text{GHET}} = \sum_{\ell} (T_{[\eta]\ell} - V_{[\eta\eta]\ell} - V_{[\epsilon\eta]\ell}) \quad (68) $$

with the kinetic-energy density of the vorticities

$$ T_{[\eta]\ell} = \frac{1}{2} \zeta (\boldsymbol{n}_{\ell})^2, \quad (69) $$

and the additional terms of the potential-energy density

$$ \begin{aligned} \tilde{V}_{[\eta\eta]\ell} &= \sum_{\alpha\gamma\delta} \sigma_{\ell}^{\gamma\delta} \eta_{\ell}^{\alpha\gamma} \eta_{\ell}^{\alpha\delta} \\ &= \frac{1}{2} \left( \underbrace{\mathrm{Tr}\{\sigma_{\ell}\}}_{=0} \boldsymbol{n}^2 - \sum_{\gamma\delta} \sigma_{\ell}^{\gamma\delta} \eta_{\gamma} \eta_{\delta} \right) \end{aligned} \quad (70) $$

and

$$ \tilde{V}_{[\epsilon\eta]\ell} = \sum_{\alpha\gamma\delta} \sigma_{\ell}^{\gamma\delta} \eta_{\ell}^{\alpha\gamma} \epsilon^{\alpha\delta} = \frac{1}{2} \boldsymbol{\tau}_{\ell} \cdot \boldsymbol{n}_{\ell} \quad (71) $$

with the coupling vector $\boldsymbol{\tau}_{\ell}$ defined in Eq. (65).

In the following it will be of use to write the components of this vector as linear combination of the strain components

$$ \boldsymbol{\tau}_{\ell}^{\prime\prime} = \sum_{\alpha\beta} \boldsymbol{t}_{\ell}^{(\nu),\alpha\beta} \boldsymbol{\epsilon}^{\alpha\beta} \quad (72) $$

with

$$ \boldsymbol{t}_{\ell}^{(\nu),\alpha\beta} = \frac{\partial}{\partial\boldsymbol{\epsilon}^{\alpha\beta}} \boldsymbol{\tau}_{\ell}^{\prime\prime} \quad (73) $$

The coefficients $\boldsymbol{t}_{\ell}^{(\nu),\alpha\beta}$ are just combinations of local stresses. They are listed explicitly in Table 2.

The equation of motion (in frequency space with $z = \omega^2 + i\epsilon$), for $\mathbf{u}(\mathbf{R}, z)$ is

$$ -\rho_m z \mathbf{u}^{\alpha}(\mathbf{R}, z) = \sum_{\beta} \frac{\partial}{\partial x_{\beta}} \left( \sum_{\gamma\delta} C^{\alpha\gamma\delta\beta}(\mathbf{R}) \epsilon^{\gamma\delta}(\mathbf{R}, z) \right. \\ \left. + \sum_{\nu} \sum_{\ell} s_{\ell}^{(\nu),\alpha\beta} \eta_{\ell}^{\nu}(\mathbf{R}, z) \right) \quad (74) $$

with

$$ s_{\ell}^{(\nu),\alpha\beta} = \frac{1}{2} \begin{cases} t_{\ell}^{(\nu),\alpha\alpha} & \alpha = \beta \\ \frac{1}{2} t_{\ell}^{(\nu),\alpha\beta} & \alpha \neq \beta \end{cases} \quad (75) $$

For $\eta(\mathbf{R}, z)$ we have

$$ -\zeta z \eta_{\ell}^{\nu}(\mathbf{R}, z) = \sum_{\mu} \sigma_{\ell}^{\nu\mu} \eta_{\ell}^{\mu}(\mathbf{R}, z) - \tau_{\ell}^{\mu}(\mathbf{R}, z) \quad (76) $$

Defining the non-coupled Green matrix of $\eta(\mathbf{R}, z)$ as

$$ [G_{\ell}(z)^{-1}]^{\nu\mu} = \sigma^{\nu\mu} + \zeta z \delta_{\nu\mu}, \quad (77) $$

we can solve for $\eta_{\ell}^{\nu}$:

$$ \begin{aligned} \eta_{\ell}^{\nu}(\mathbf{R}, z) &= \sum_{\mu} G_{\ell}^{\nu\mu}(z) \tau_{\ell}^{\nu}(\mathbf{R}, z) \\ &= \sum_{\mu} G_{\ell}^{\nu\mu}(z) \sum_{\gamma\delta} t_{\ell}^{(\nu),\gamma\delta} \epsilon^{\gamma\delta} \end{aligned} \quad (78) $$

Inserting Eq. (78) into (74) we obtain

$$ \begin{aligned} &-\rho_m z \mathbf{u}^{\alpha}(\mathbf{R}, z) \\ &= \sum_{\beta} \frac{\partial}{\partial x_{\beta}} \sum_{\gamma\delta} (C^{\alpha\gamma\delta\beta}(\mathbf{R}) + \Delta C^{\alpha\gamma\delta\beta}(z)) \epsilon^{\gamma\delta}(\mathbf{R}, z), \end{aligned} \quad (79) $$

<table><caption>Table 2 | List of the stress contributions $\boldsymbol{t}_{\ell}^{(\nu),\alpha\beta}$</caption><tr><td>t<sup>(x),xx</sup> = 0</td><td>t<sup>(y),xx</sup> = -σ<sup>xz</sup></td><td>t<sup>(z),xx</sup> = σ<sup>yy</sup></td></tr><tr><td>t<sup>(x),yy</sup> = σ<sup>yz</sup></td><td>t<sup>(y),yy</sup> = 0</td><td>t<sup>(z),yy</sup> = -σ<sup>xy</sup></td></tr><tr><td>t<sup>(x),zz</sup> = -σ<sup>yz</sup></td><td>t<sup>(y),zz</sup> = σ<sup>xz</sup></td><td>t<sup>(z),zz</sup> = 0</td></tr><tr><td>t<sup>(x),yz</sup> = σ<sup>yy</sup> - σ<sup>zz</sup></td><td>t<sup>(y),yz</sup> = -σ<sup>xy</sup></td><td>t<sup>(z),yz</sup> = -σ<sup>xz</sup></td></tr><tr><td>t<sup>(x),xz</sup> = -σ<sup>xy</sup></td><td>t<sup>(y),xz</sup> = σ<sup>xx</sup> - σ<sup>zz</sup></td><td>t<sup>(z),xz</sup> = -σ<sup>yz</sup></td></tr><tr><td>t<sup>(x),xy</sup> = σ<sup>xz</sup></td><td>t<sup>(y),xy</sup> = -σ<sup>yz</sup></td><td>t<sup>(z),xy</sup> = σ<sup>yy</sup> - σ<sup>xx</sup></td></tr></table>

For brevity we omitted the subscript $\ell$ labelling the local regions.
where the additional contributions to the elastic coefficients are given by

$$ \Delta C^{\alpha\gamma\delta\beta}(z) = \sum_{\ell} \sum_{\ell'} s_{\ell}^{(\nu),\alpha\beta} G_{\ell}^{\prime\prime\mu}(z) t_{\ell'}^{(\mu)\nu\delta} \quad (80) $$

We now simplify the treatment as follows:

Because all patches around a local stress $\sigma_{\ell}$ involve a vorticity pattern, which is just described by a single variable, namely $\eta_{\ell}^{(z)} \equiv \eta_{\ell}$, where $\bar{z}$ indicates the z axis of the local coordinate system pointing into the direction of $\eta_{\ell}$, we work in terms of scalar variables $\eta_{\ell}$ with a Green's function

$$ G_{\ell}(z) = \frac{1}{\bar{z} - \sigma_{\ell}} \quad (81) $$

where $\sigma_{\ell}$ is now any of the local stress components. Because the stress-induced correction to the elastic coefficients, given by Eq. (80) is just a quadratic form of the local stresses with the local Green's functions as coefficients we define – in the spirit of the isotropic approximation of the HET, Eq. (21)– the additional contribution to the frequency-dependent elastic self energy as

$$ \Delta C^{\alpha\gamma\delta\beta}(z) \rightarrow \Sigma_{\epsilon\eta}(z) (\delta_{\alpha\gamma}\delta_{\delta\beta} + \delta_{\alpha\delta}\delta_{\gamma\beta}) \quad (82) $$

with

$$ \Sigma_{\epsilon\eta}(z) = \sum_{\ell} \sigma_{\ell}^2 G_{\ell}(z) \quad (83) $$

where $\sigma_{\ell}$ is, again, any matrix element of the local stresses. The simplified version of the generalized heterogenous-elasticity (GHET) is now obtained by replacing

$$ \Sigma(z) \rightarrow \Sigma_{\text{HET}}(z) + \Sigma_{\epsilon\eta}(z) $$

The contribution to the shear-modulus self energy takes the form

$$ \Sigma_{\epsilon\eta}(z) \sim \frac{\overline{\sigma_{\ell}^2}}{-z\bar{\zeta} + \sigma_{\ell}} = \int d\sigma \mathcal{P}(\sigma) \frac{\sigma^2}{-z\bar{\zeta} + \sigma} \quad (84) $$

The overbar indicates the average over the local stresses and $\mathcal{P}(\sigma)$ is their distribution density.

The corresponding contribution to the spectrum $\rho_{\epsilon\eta}(\lambda)$ is given by

$$ \rho_{\epsilon\eta}(\lambda) \sim \Sigma''_{\epsilon\eta}(\lambda) \sim \sigma^2 \mathcal{P}(\sigma)|_{\sigma=\lambda\bar{\zeta}} \quad (85) $$

### Evaluation of $\mathcal{P}(\sigma)$

In this appendix we discuss the general properties of the distribution of the local stress and what is expected for the low stress value behaviour of this distribution.

**Basic.** We are dealing with the local stress:

$$ \sigma_{ij}^{\alpha\beta} \doteq r_{ij}^{\alpha} r_{ij}^{\beta} \frac{\phi'(r_{ij})}{r_{ij}} \quad (86) $$

or, by introducing the unit vector $\hat{\mathbf{r}}=\mathbf{r}/|\mathbf{r}|$

$$ \sigma_{ij}^{\alpha\beta} = \hat{e}_{ij}^{\alpha} \hat{e}_{ij}^{\beta} r_{ij} \phi'(r_{ij}) \quad (87) $$

For sake of simplicity, let's omit the subscript $ij$ and define $\psi(r) = r\phi'(r)$:

$$ \sigma^{\alpha\beta} = \hat{e}^{\alpha} \hat{e}^{\beta} \psi(r) \quad (88) $$

As discussed in paragraph M2, we are interested in the traceless stress tensor:

$$ \tilde{\sigma}^{\alpha\beta} = \left( \hat{e}^{\alpha} \hat{e}^{\beta} - \frac{1}{d} \delta_{\alpha\beta} \right) \psi(r) \quad (89) $$

The angular and radial part of the stress for a given couple of particle enter multiplicatively. Given the cylindrical symmetry of the problem (the angular part only depends on the orientation of the axis defined by the two particles) we expect that the angular part plays no role in the distribution if the system is isotropic. Indeed, we can make a rotation of the matrix $M^{\alpha\beta} = \hat{e}^{\alpha} \hat{e}^{\beta}$ so to orient the vector $\hat{e}$ along the first coordinate. Being the matrix $M^{\alpha\beta}$ dyadic, its eigenvalues are all zero but one, which value is equal to one, and the corresponding eigenvector is parallel to $\hat{e}$. After rotation, the stress matrix in any dimension is thus:

$$ \bar{\bar{\sigma}} = \begin{pmatrix} 1 & 0 & \cdots & 0 \\ 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & 0 \end{pmatrix} \psi(r) $$

and its traceless countepart:

$$ \bar{\bar{\sigma}} = \begin{pmatrix} \frac{d-1}{d} & 0 & \cdots & 0 \\ 0 & -\frac{1}{d} & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & -\frac{1}{d} \end{pmatrix} \psi(r) $$

As expected, the angular part is deterministic, and the distribution of the stress is controlled by the distribution of $\psi(r)$. The sign of the elements of the matrix $\bar{\bar{\sigma}}$ however can be both positive and negative. As for the inverse law potential investigated here $\psi(r)$ is always negative, the traceless stress matrix has $d-1$ positive and 1 negative elements. As discussed, only positive stress values give rise to non-phononic type II modes. Importantly, in any dimension, positive stress values do exist.

**Stress distribution.** The quantity that defines, for each couple of particles, the absolute value of the stress is $|\psi(r)|$. As it represents the value of the stress, let's call it 'σ' (a scalar, not to be confused with the stress tensor). That is

$$ \sigma \doteq |\psi(r)| = |r\phi'(r)| \quad (90) $$

Let's now search for the distribution of $\sigma$, let's say $\mathcal{P}(\sigma)$.

The distribution of the $r_{ij}$ value is governed by the particle pair distribution function $g(r)$. Specifically in dimension $d=2$:

$$ P(r) = 2\pi\rho r^2 g(r) \quad (91) $$

(where $\rho=N/\Omega$ is the particle density) and in $d=3$:

$$ P(r) = 4\pi\rho r^2 g(r) \quad (92) $$

and, using

$$ \mathcal{P}(\sigma)|d\sigma| = P(r)|dr| \quad (93) $$
we get

$$ \mathcal{P}(\sigma) \sim r^{d-1} g(r) \left| \frac{d\sigma}{dr} \right|^{-1} \qquad (94) $$

**The case of a smooth (tapered) cutoff.** In real numerical studies, the interaction potential is “shifted” and “tapered” in such a way to be zero at a finite r value (at the cut-off value $r_c$) and to have continuity up to its m-th derivative at $r_c$. Its expression becomes:

$$ \begin{aligned} \phi(r) &= a(r^{-n} - r_c^{-n}) (r_c^2 - r^2)^m \quad r<r_c \\ \phi(r) &\sim (r_c - r)^{m+1} \quad r<r_c \\ \phi(r) &= 0 \quad r>r_c \end{aligned} \qquad (95) $$

The value of m is (almost) always chosen to be m = 2, as the continuity of the first derivative of the potential (the force) is needed for the correct energy conservation during the dynamics, and the continuity of the second derivative guarantees a good estimation of the Hessian.

By defining $x=r/r_c$ we can rewrite the previous equation as:

$$ \begin{aligned} \phi(x) &= a'(x^{-n} - 1) (1 - x^2)^m \quad x<1 \\ \phi(x) &= 0 \quad x>1 \end{aligned} \qquad (96) $$

Obviously $\sigma=|r\phi'(r)|=|x\phi'(x)|$.

Being interested in the r-region where the stress is small, we have now to consider the region around $r \approx r_c$, i.e. $x \approx 1$. In this case, it's worth to define $\epsilon=1-x^2$ and to expand the potential around $\epsilon=0$ ($\epsilon>0$ hereafter):

$$ \phi(\epsilon) = a' \left( (1 - \epsilon)^{-n/2} - 1 \right) \epsilon^m \sim \epsilon^{m+1} \qquad (97) $$

Now, remembering we are in the limit of small $\epsilon$:

$$ \begin{aligned} \sigma &= |x^2 \phi'(x)| = (1 - \epsilon) \left| \frac{d\phi(\epsilon)}{dx} \right| \sim \left| \frac{d\epsilon^{m+1}}{dx} \right| \\ &= \frac{d\epsilon^{m+1}}{d\epsilon} \left| \frac{d\epsilon}{dx} \right| \sim \epsilon^m \left| \frac{d\epsilon}{dx} \right| = \epsilon^m \left| \frac{d(1 - x^2)}{dx} \right| \\ &= \epsilon^m x \sim \epsilon^m \end{aligned} \qquad (98) $$

Thus

$$ \sigma = \epsilon^m \qquad (99) $$

and

$$ \epsilon = \sigma^{\frac{1}{m}} \qquad (100) $$

Now

$$ \frac{d\sigma}{dr} = \frac{d\sigma}{d\epsilon} \frac{d\epsilon}{dr} \sim \epsilon^{m-1}(-r) \approx \epsilon^{m-1}(-r_c) \qquad (101) $$

thus

$$ \left| \frac{d\sigma}{dr} \right| \sim \epsilon^{m-1} \sim \sigma^{\frac{m-1}{m}} \qquad (102) $$

Finally, using Eq. (94), we have

$$ \mathcal{P}(\sigma) \sim r_c^{d-1} g_2(r_c) \sigma^{\frac{1}{m}} \qquad (103) $$

We observe that the behaviour of $\mathcal{P}(\sigma)$ for $\sigma \to 0$, does not depend on $n$. Rather it only depends on the value of $m$: it is the presence of a cut-off, and the continuity of $\phi(r)$ and of its first $m$ derivatives at the cut-off position that determine the stress distribution. Further, we

observe that for $m=2$ we have $\mathcal{P}(\sigma) \sim \sigma^{-\frac{1}{m}}$, which, in turn, as we will see in the next paragraphs, implies $\rho(\omega) \sim \omega^4$.

**The case of a potential with a minimum.** If an interatomic potential $\phi(r)$ has a minimum, say, at $r=r_0$ like in the case of a Lennard-Jones (LJ) potential, then the first derivative passes linearly through zero at $r=r_0$. Obviously, this corresponds to the case of a tapered potential with $m=1$. Correspondingly we have

$$ \mathcal{P}(\sigma) \sim r_0^{d-1} g_2(r_0) \sigma^0 = \text{const.} \qquad (104) $$

It is plausible that the overall distribution of the small stresses is dominated by this contribution because $g_2(r)$ is maximal at $r_0$.

**From $\mathcal{P}(\sigma)$ to $g(\omega)$**. As discussed in the main text and in M2, there are a direct and an indirect contributions of the type-II excitations to the spectrum, $\rho_{\epsilon}(\lambda)$ and $\rho_{\epsilon}(\lambda)$, given by $(\lambda = \omega^2)$

$$ \rho_{\epsilon}(\lambda) = \frac{1}{2\omega} g(\omega) \sim \sigma^2 P(\sigma)|_{\sigma=\lambda\zeta} \qquad (105) $$

From this, we obtain the following relations for the spectra

### Potential with cutoff and tapering

$$ \begin{aligned} \rho_{\epsilon\eta}(\lambda) &\sim \lambda^{\eta+1} \\ g_{\epsilon\eta}(\omega) &\sim \omega^{\eta+3} \end{aligned} \qquad (106) $$

### Potential with minimum at $r=r_0$

$$ \begin{aligned} \rho_{\epsilon\eta}(\lambda) &\sim \lambda^2 \\ g_{\epsilon\eta}(\omega) &\sim \omega^5 \end{aligned} \qquad (107) $$

### Simulation details

We performed numerical simulations using the MC-swap algorithm of a d-dimensional ($d=3$ here) system composed of $N=10^d$ particles in a box (with periodic boundary conditions) of side $L=N^{1/d}$ so that the number density is $\rho=N/L^d=1$.

To find the inherent glass structures, we consider Monte Carlo (MC) dynamics combined with swap moves for equilibrating the system in a wide range of temperatures. We compute stable glass configurations from different parent temperatures by minimizing the configurational energy $\Phi$ and thus obtaining the corresponding inherent structures through the Limited-memory Broyden-Fletcher-Goldfarb-Shanno (LBFGS) algorithm$^{81}$. The spectrum of the harmonic oscillations around the inherent configurations has been obtained by computing all the $dN$ eigenvalues of the hessian matrix using Python NumPy linear algebra functions$^{82}$. We produced data for two tapering functions $m=2$ and $m=\infty$. The potential is given in Eq. (14) of the main text.

Following early works on MC-swap$^{70,71}$, at each MC step we perform a swap move instead of a standard MC step with probability $p_{swap} = 0.2$. The displacement vector $\Delta r$ of the standard MC move has modulus $|\Delta r| = 0.2$ so that the acceptance probability falls into the interval 30% (low temperature) to 60% at high temperatures. We consider runs of $N_t = 1.1 \times 10^6$ MC steps. The power of the repulsive potential is $n=12$.

### Data availability

All data shown in Figs. 2–4 are available on request from the authors.

### References

1. Pan, Z. et al. Disorder classification of the vibrational spectra of modern glasses. Phys. Rev. B **104**, 134106 (2021).
2. Leadbetter, A. J. Inelastic cold neutron scattering from different forms of silica. *J. Chem. Phys.* **51**, 779 (1969).

3. Zeller, R. C. & Pohl, R. O. Thermal conductivity and specific heat of noncrystalline solids. *Phys. Rev. B* **4**, 2029 (1971).

4. Buchenau, U., Nücker, N. & Dianoux, A. J. Neutron scattering study of the low-frequency vibrations in vitreous silica. *Phys. Rev. Lett.* **53**, 2316 (1984).

5. Wuttke, J., Petry, W., Coddens, G. & Fujara, F. Fast dynamics of glass-forming glycerol. *Phys. Rev. E* **52**, 4026 (1995).

6. Foret, M., Courtens, E., Vacher, R. & Suck, J.-B. Scattering investigation of acoustic localization in fused silica. *Phys. Rev. Lett.* **77**, 3831 (1996).

7. Sette, F., Krisch, M. H., Masciovecchio, C., Ruocco, G. & Monaco, G. Dynamics of glasses and glass-forming liquids studied by inelastic x-ray scattering. *Science* **280**, 1550 (1998).

8. Chumakov, A. I. et al. Collective nature of the boson peak and universal transboson dynamics of glasses. *Phys. Rev. Lett.* **92**, 245508 (2004).

9. Monaco, G. & Giordano, V. M. Breakdown of the debye approximation for the acoustic modes with nanometric wavelengths in glasses. *Proc. Nat. Acad. Sci.* **106**, 3659 (2009).

10. Baldi, G., Giordano, V. M. & Monaco, G. Elastic anomalies at terahertz frequencies and excess density of vibrational states in silica glass. *Phys. Rev. B* **83**, 174203 (2011).

11. Jäckle, J. Low-frequency Raman scattering in glasses, in *Amorphous Solids: Low-Temperature Properties* (ed. Phillips, W. A.) 135 (Springer-Verlag, Berlin, 1981).

12. Karpov, V. G., Klinger, M. I. & Ignat'ev, F. N. Theory of the low-temperature anomalies in the thermal properties of amorphous structures. *Z. Eksp. Teor. Fiz* **84**, 760 (1983).

13. Buchenau, U., Galperin, Y. M., Gurevich, V. L. & Schober, H. R. Anharmonic potentials and vibrational localization in glasses. *Phys. Rev. B* **43**, 5039 (1991).

14. Schirmacher, W., Diezemann, G. & Ganter, C. Harmonic vibrational excitations in disordered solids and the boson peak. *Phys. Rev. Lett.* **81**, 136 (1998).

15. Götze, W. & Mayr, M. R. Evolution of vibrational excitations in glassy systems. *Phys. Rev. E* **61**, 587 (1999).

16. Taraskin, S. N., Loh, Y. L., Natarajan, G. & Elliott, S. R. Origin of the boson peak in systems with lattice disorder. *Phys. Rev. Lett.* **86**, 1255 (2002).

17. Gurevich, V. L., Parshin, D. A. & Schober, H. R. Anharmonicity, vibrational instability, and the boson peak in glasses. *Phys. Rev. B* **67**, 094203 (2003).

18. Grigera, T. S., Martín-Mayor, V., Parisi, G. & Verrocchio, P. Phonon interpretation of the boson peak in supercooled liquids. *Nature* **422**, 289 (2003).

19. Schirmacher, W. Thermal conductivity of glassy materials and the boson peak. *Europhys. Lett.* **73**, 892 (2006).

20. Schirmacher, W., Ruocco, G. & Scopigno, T. Acoustic attenuation in glasses and its relation with the Boson Peak. *Phys. Rev. Lett.* **98**, 025501 (2007).

21. Ganter, C. & Schirmacher, W. Rayleigh scattering, long-time tails and the harmonic spectrum of topologically disordered systems. *Phys. Rev. B* **82**, 094205 (2010).

22. Franz, S., Parisi, G., Urbani, P. & Zamponi, F. Universal spectrum of normal modes in low-temperature glasses. *Proc. Nat. Acad. Sci.* **112**, 14539 (2015).

23. Laird, B. B. & Schober, H. R. Localized low-frequency vibrational modes in a simple model glass. *Phys. Rev. Lett.* **66**, 636 (1991).

24. Léonforte, F., Tanguy, A., Wittmer, J. P. & Barrat, J.-L. Inhomogeneous elastic response of silica glass. *Phys. Rev. Lett.* **97**, 055501 (2006).

25. Shintani, H. & Tanaka, H. Universal link between the boson peak and transverse phonons in glass. *Nat. Mater.* **7**, 870 (2008).

26. Monaco, G. & Mossa, S. Anomalous properties of the acoustic excitations in glasses on the mesoscopic length scale. *Proc. Nat. Acad. Sci.* **106**, 16907 (2009).

27. Mayr, S. G. Relaxation kinetics and mechanical stability of metallic glasses and supercooled melts. *Phys. Rev. B* **79**, 060201 (2009).

28. Marruzzo, A., Schirmacher, W., Fratalocchi, A. & Ruocco, G. Heterogeneous shear elasticity of glasses: the origin of the boson peak. *Sci. Rep.* **3**, 1 (2013).

29. Derlet, P. M., Maass, R. & Löffler, J. F. The boson peak of model glass systems and its relation to atomic structure. *Eur. Phys. J. B* **85**, 148 (2012).

30. Mizuno, H., Mossa, H. & Barrat, J.-L. Elastic heterogeneity, vibrational states, and thermal conductivity across an amorphisation transition. *Europhys. Lett.* **104**, 56001 (2013).

31. Chumakov, A. I. et al. Equivalence of the boson peak in glasses to the transverse acoustic van Hove singularity in crystals. *Phys. Rev. Lett.* **106**, 225501 (2011).

32. Hu, Y.-C. & Tanaka, H. Origin of the boson peak in amorphous solid. *Nat. Phys.* **18**, 669 (2022).

33. Schirmacher, W. and Ruocco, G. Vibrational excitations in disordered solids, in *Encyclopedia of Condensed-Matter Physics* (ed Chakraborty, T.) 298 (Elsevier, 2023).

34. Schirmacher, W., Scopigno, T. & Ruocco, G. Theory of vibrational anomalies in glasses. *J. Noncryst. Sol.* **407**, 133 (2014).

35. Debye, P. Zur Theorie der spezifischen Wärme. *Ann. Phys.* (Leipzig) **39**, 789 (1912).

36. Economou, E. N. Green's Functions in Quantum Physics (Springer-Verlag, Heidelberg, 1979).

37. Baity-Jesi, M., Martín-Mayor, V., Parisi, G. & Perez-Gaviro, S. Soft modes, localization, and two-level systems in spin glasses. *Phys. Rev. Lett.* **115**, 267205 (2015).

38. Lerner, E., Düring, G. & Bouchbinder, E. Statistics and properties of low-frequency vibrational modes in structural glasses. *Phys. Rev. Lett.* **117**, 035501 (2016).

39. Lerner, E. & Bouchbinder, E. Effect of instantaneous and continuous quenches on the density of vibrational modes in model glasses. *Phys. Rev. E* **96**, 020104 (2017).

40. Shimada, M., Mizuno, H. & Ikeda, A. Anomalous vibrational properties in the continuum limit of glasses. *Phys. Rev. E* **97**, 022609 (2018).

41. Lerner, E. & Bouchbinder, E. Low-energy quasilocalized excitations in structural glasses. *J. Chem. Phys.* **155**, 200901 (2021).

42. Angelani, L., Paoluzzi, M., Parisi, G. & Ruocco, G. Probing the non-debye low-frequency excitations in glasses through random pinning. *Proc. Nat. Acad. Sci.* **115**, 8700 (2018).

43. Paoluzzi, M., Angelani, L., Parisi, G. & Ruocco, G. Relatiion between heterogeneous frozen regions in supercooled liquids and non-debye spectrum in the corresponding glasses. *Phys. Rev. Lett* **123**, 155502 (2019).

44. Paoluzzi, M., Angelani, L., Parisi, G. & Ruocco, G. Probing the debye spectrum in glasses using small system sizes. *Phys. Rev. Res.* **2**, 043248 (2021).

45. Franz, S. & Parisi, G. The simplest model of jamming. *J. Phys. A: Math. Theor.* **49**, 145001 (2016).

46. Franz, S., Parisi, G., Sevelev, M., Urbani, P. & Zamponi, F. Universality of the sat-unsat (jamming) threshold in non-convex continuous constraint satisfaction problems. *SciPost Phys.* **2**, 019 (2017).

47. Cugliandolo, L. F. & Kurchan, J. Analytical solution of the off-equilibrium dynamics of a long-range spin-glass model. *Phys. Rev. Lett.* **71**, 173 (1993).

48. Folena, G., Franz, S., Ricci-Tersenghi, F. & Ruocco, G. Rethinking mean-field glassy dynamics and its relation with the energy
landscape: The surprising case of the spherical mixed p-spin model. *Phys. Rev. X* **10**, 031045 (2020).

49. Schirmacher, W., Folli, V., Ganter, C. & Ruocco, G. Self-consistent euclidean random matrix theory. *J. Phys. A: Math. Theor.* **52**, 464002 (2019).

50. Marruzzo, A., Schirmacher, W., Köhler, S., Fratalocchi, A. & Ruocco, G. Vibrational anomalies and marginal stability of glasses. *Eur. Phys. J. spec. Topics* **216**, 83 (2013).

51. Schirmacher, W., Bryk, T. & Ruocco, G. Modeling the instantaneous normal mode spectra of liquids as that of unstable elastic media. *Proc. Nat. Acad. Sci.* **119**, e2119288119 (2022).

52. Krishnan, V. V., Karmakar, S. & Ramola, K. Singularities in hessian element distributions of amorphous media. *Phys. Rev. Res.* **2**, 042025(R) (2020).

53. Krishnan, V. V., Ramola, K. & Karmakar, S. Universal non-debye low-frequency vibrations in sheared amorphous solids. *Soft Matter* **18**, 3395 (2022).

54. Lutsko, J. F. Stress and elastic constants in anisotropic solids: Molecular dynamics techniques. *J. Appl. Phys.* **64**, 1152 (1988).

55. Lutsko, J. F. Generalized expressions for the calculation of elastic constants by computer simulation. *J. Appl. Phys.* **65**, 2991 (1989).

56. Strutt (Lord Rayleigh), J. W. On the scattering of light by small particles. *Philos. Mag.* **41**, 247 (1871).

57. Alexander, S. Amorphous solids: their structure, lattice dynamics and elasticity. *Phys. Reports* **296**, 65 (1998).

58. Wang, L., Fu, L. & Nie, Y. Density of states below the first sound mode in 3d glasses. *J. Chem. Phys.* **157**, 074502 (2022).

59. Wang, L., Szamel, G. & Flenner, E. Scaling of the non-phononic spectrum of two-dimensional glasses. *J. Chem. Phys.* **158**, 126101 (2023).

60. Mocanu, F. C. et al. Microscopic observation of two-level systems in a metallic glass model. *J. Chem. Phys.* **158**, 014501 (2023).

61. Bouchbinder, E., Lerner, E., Rainone, C., Urbani, P. & Zamponi, F. Low-frequency vibrational spectrum of mean-field disordered systems. *Phys. Rev. B* **103**, 174202 (2021).

62. Ashcroft, N. W. & Mermin, D. *Solid State Physics* (Harcourt College Publishers, 1976).

63. Alexander, S. Is the elastic energy of amorphous materials rotationally invariant? *J. Physique* **45**, 1939 (1984).

64. Rainone, C., Bouchbinder, E. & Lerner, E. Pinching a glass reveals key properties of its soft spots. *Proc. Nat. Acad. Sci. USA* **117**, 5228 (2020).

65. Allen, P. B., Feldman, J. L. & Fabian, J. Diffusons, locons and pro-pagons: character of atomic vibrations in amorphous si. *Philos. Mag.* **79**, 1715 (1999).

66. Ishimaru, A. *Wave Propagation and Scattering in Random Media* (Academic Press, 1978).

67. Mizuno, H., Shiba, H. & Ikeda, A. Continuum limit of the vibrational properties of amorphous solids. *Proc. Nat. Acad. Sci. USA* **114**, E9767 (2017).

68. O'Hern, C. S., Silbert, L. E., Liu, A. J. & Nagel, S. R. Jamming at zero temperature and zero applied stress: The epitome of disorder. *Phys. Rev. E* **68**, 011306 (2003).

69. Parisi, G. & Zamponi, F. Mean-field theory of hard sphere glasses and jamming. *Rev. Mod. Phys.* **82**, 789 (2010).

70. Ninarello, A., Berthier, L. & Coslovich, D. Models and algorithms for the next generation of glass transition studies. *Phys. Rev. X* **7**, 021039 (2017).

71. Grigera, T. S. & Parisi, G. Fast Monte Carlo algorithm for super-cooled soft spheres. *Phys. Rev. E* **63**, 045102 (2001).

72. Wang, L. et al. Low-frequency vibrational modes of stable glasses. *Nat. Commun.* **10**, 26 (2019).

73. Khomenko, D., Scalliet, C., Berthier, L., Reichman, D. R. & Zamponi, F. Depletion of two-level systems in ultrastable computer-generated glasses. *Phys. Rev. Lett.* **124**, 225901 (2020).

74. Liao, Q. & Berthier, L. Hierarchical landscape of hard disk glasses. *Phys. Rev. X* **9**, 011049 (2019).

75. Scalliet, C. & Berthier, L. Rejuvenation and memory effects in a structural glass. *Phys. Rev. Lett.* **122**, 255502 (2019).

76. Scalliet, C., Berthier, L. & Zamponi, F. Nature of excitations and defects in structural glasses. *Nat. Commun.* **10**, 1 (2019).

77. Mehta, M. L. *Random Matrices* (Academic, 1967).

78. Atas, Y. Y., Bogomolny, E., Giraud, O. & Roux, G. The distribution of the ratio of consecutive level spacings in random-matrix ensembles. *Phys. Rev. Lett.* **110**, 084101 (2013).

79. Schirmacher, W. et al. Acoustic attenuation in glasses and its relation with the boson peak. *Phys. Stat. Sol.* **5**, 862 (2008).

80. Vitek, V. & Egami, T. Atomic level stresses in solids and liquids. *Phys. Stat. Sol.* **144**, 145 (1987).

81. Bonnans, J.-F., Gilbert, J. C., Lemaréchal, C. & Sagastizábal, C. A. *Numerical Optimization: Theoretical and Practical Aspects* (Springer Science & Business Media, 2006).

82. McKinney, W. *Python for Data Analysis: Data Wrangling with Pandas, NumPy, and IPython* (O'Reilly Media, Inc, 2012).

## Acknowledgements

We are grateful to Hideyuki Mizuno and Vishnu Krishnan for discussions concerning small Lennard-Jones systems. M.P. has received funding from the European Union's Horizon 2020 research and innovation programme under the MSCA grant agreement No. 801370 and by the Secretary of Universities and Research of the Government of Catalonia through Beatriu de Pinós program Grant No. BP 00088 (2018). G.S. acknowledges the support of NSF Grant No. CHE 2154241.

## Author contributions

G.R. and W.S. set up the project and worked out the theory. M.P., F.C.M., and D.K. performed the simulations, G.S. and F.Z. supported the analysis, G.R. and W.S. wrote the manuscript with contributions from all authors.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at
https://doi.org/10.1038/s41467-024-46981-7.

**Correspondence** and requests for materials should be addressed to Walter Schirmacher or Giancarlo Ruocco.

**Peer review information** *Nature Communications* thanks Yuanchao Hu, Yuliang Jin and Jie Zhang for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at
http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
**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/.

© The Author(s) 2024