Universal disorder-induced broadening of phonon bands: from disordered lattices to glasses

The translational symmetry of solids, either ordered or disordered, gives rise to the existence of low-frequency phonons. In ordered systems, either crystalline solids or isotropic homogeneous continua, some phonons characterized by different wavevectors are degenerate, i.e. they share the exact same frequency ω; in finite-size systems, phonons form a discrete set of bands with nq(ω)-fold degeneracy. Here we focus on understanding how this degeneracy is lifted in the presence of disorder, and its physical implications. Using standard degenerate perturbation theory and simple statistical considerations, we predict the dependence of the disorder-induced frequency width of phonon bands to be Δ ω ∼ σ ω n q / N , where σ is the strength of disorder and N is the total number of particles. This theoretical prediction is supported by extensive numerical calculations for disordered lattices—characterized by topological, mass, stiffness and positional disorder—and for computer glasses, where disorder is self-generated, thus establishing its universal nature. The predicted scaling of phonon band widths leads to the identification of a crossover frequency ω † ∼ L − 2 / ( đ + 2 ) in systems of linear size L in đ > 2 dimensions, where the disorder-induced width of phonon bands becomes comparable to the frequency gap between neighboring bands. Consequently, phonons continuously cover the frequency range ω > ω†, where the notion of discrete phonon bands becomes ill-defined. Two basic applications of the theory are presented; first, we show that the phonon scattering lifetime is proportional to (Δω)−1 for ω < ω†. Second, the theory is applied to the basic physics of glasses, allowing us to determine the range of frequencies in which the recently established universal density of states of non-phononic excitations can be directly probed for different system sizes.


Introduction
The low-frequency spectra of condensed matter is generically populated by Goldstone modes that arise due to broken continuous symmetries [1,2]. In solids, these low-frequency Goldstone modes are long-wavelength phonons (plane-waves) associated with translational symmetry. The distribution of long-wavelength phonons over frequency is of fundamental importance as it determines many thermodynamic and transport properties of solids [3,4]. In ordered systems, either crystalline solids or isotropic homogeneous continua, some phonons characterized by different wavevectors are degenerate, i.e. they share the exact same frequency ω. Moreover, in finite-size systems, phonons are quantized into a discrete set of bands with n q (ω)-fold degeneracy, where q is a phononic band index.
Low-frequency phonons likewise exist in disordered systems, such as disordered crystals/lattices and structural glasses, because they also feature translational symmetry. The presence of disorder, however, is expected to affect the way long-wavelength phonons are distributed over frequency. In particular, the degeneracy of phonons is expected to be lifted by the presence of disorder such that phonon bands broaden, featuring a finite frequency width ∆ω each. In this paper we address the following basic question: how do the phonon band widths ∆ω depend on the strength of disorder σ, on the frequency ω, on the original degeneracy level n q (ω), and on the system size N ? To the best of our knowledge, despite its basic nature and the longtime extensive efforts devoted to understanding phonons in disordered systems [5][6][7][8][9][10][11][12], this question has not been addressed before.
We show that a simple statistical scaling theory predicts ∆ω(σ, ω, n q , N ) ∼ σ ω √ n q / √ N . An immediate implication of the existence of finite band widths ∆ω(σ, ω, n q , N ) is that there exists a crossover frequency ω † (L) ∼ L − 2 d+2 in a system of linear size L ∼ N 1/d in d > 2 spatial dimensions (thed = 2 case is discussed separately), beyond which the width of phonon bands becomes comparable to the frequency gap between neighboring bands. That is, disordered phonons continuously cover the frequency range ω > ω † (L), where the notion of discrete phonon bands becomes ill-defined. The theoretical prediction for ∆ω is supported by extensive numerical calculations for disordered lattices characterized by topological, mass, stiffness and positional disorder, over a large range of system sizes. Quite remarkably, the analytic prediction for ∆ω is shown to be valid also for computer glass-formers, where disorder is self-generated and frustration-induced internal stresses generically emerge, thus establishing its universal nature.
Two basic applications of the developed theory are then presented; first, the theory is applied to phonon dynamics, where the effect of disorder on the scattering lifetime of phonons is considered. In the frequency regime ω < ω † (L), where phonon bands are welldefined, the theory predicts that the lifetime of phonons becomes finite and is proportional to (∆ω) −1 . This prediction is quantitatively verified by dynamic phonon scattering calculations in harmonic disordered lattices.
Second, the theory is applied to a basic problem in glass physics; it has been recently shown that structural glasses feature non-phononic low-frequency vibrational modes that are quasilocalized in space and follow a universal density of states D(ω) ∼ ω 4 [13][14][15][16][17][18][19]. Such nonphononic excitations have been hypothesized for decades to be the origin of various universal low-temperature anomalies in glasses [20][21][22][23][24]. Efforts to observe these glassy excitations of a given frequency ω have been hampered for a long time because their unique quasilocalized nature is destroyed due to hybridization and mixing with extended phonons that share similar frequencies [14,25]. Consequently, quasilocalized glassy excitations can be observed in gaps/holes in the phononic spectrum that exist for ω < ω † (L), explaining the recent observations of [18,19]. Finally, as ω † (L) → 0 in the L → ∞ limit, quasilocalized glassy excitations cannot be directly observed in the thermodynamic limit.

Scaling theory for disorder-induced broadening of phonon bands
The vibrational modes of a solid are determined by the eigenvalue equation M|ψ = ω 2 |ψ , where M ≡ ∂ 2 U ∂x∂x is the Hessian matrix of dimensiondN ×dN (x is the vector ofdN particles' coordinates and U is the potential energy), |ψ are the normalized eigenvectors of dimensiondN , and ω are the vibrational frequencies. Unless noted otherwise, we assume all particles share the same unit mass. There aredN eigenvectors, but there can be less vibrational frequencies due to degeneracy. For crystalline (ordered) solids or for isotropic homogeneous continua, the vibrational modes are phonons [3,4]. Some low-frequency phonons that share the same wavelength/wavenumber are degenerate, i.e. they also share exactly the same frequency ω; in finite-size systems, phonons form a discrete set of bands with n q (ω)fold degeneracy. Our first goal is to understand how this degeneracy is lifted in the presence of disorder.
Let us focus on a single degenerate band of low frequency ω and denote the eigenvectors within the band by |ψ i , with i = 1, 2, . . . , n q , and the rest (whether degenerate or not) by |ψ k , with k = n q + 1, n q + 2, . . . ,dN . We proceed in two steps; first, we use standard degenerate perturbation theory [26] to obtain an equation for the frequency shifts ∆ω due to the presence of disorder characterized by strength σ. Second, we use the law of large numbers and Wigner's semicircle law for the eigenvalues of random matrices [27] to derive the scaling of ∆ω in terms of σ, N , n q , and ω.
We denote by M (0) the Hessian matrix of the system in the absence of disorder. Disorder, which may be realized in many ways (see below), is assumed to be characterized by a typical width σ, which quantifies its strength. For example, if the disorder is extracted from some distribution, then σ is its standard deviation. The presence of disorder modifies the Hessian matrix, which now reads M = M (0) + δM, where the disorder-induced perturbation δM gives rise to shifts ∆ω in the originally degenerate frequency ω. That is, the disorder is expected to lift the degeneracy such that each frequency becomes ω + ∆ω, with n q different shifts ∆ω. The disorder-induced contribution δM is a random matrix in which the variability of the elements is determined by σ.
Within the degenerate band, any linear combination of |ψ m , with coefficients β m , is also an eigenvector The linear combination which is relevant for the perturbation problem, which we denote by |Ψ ≡ nq m=1 β m |ψ m , is not known a priori (there are n q such linear combinations, but we do not introduce another index for the ease of notation). The crux of standard degenerate perturbation theory is that |Ψ , i.e. the coefficients β m , and the frequency shifts ∆ω are simultaneously and self-consistently selected according to [26] M |β = δω 2 |β .
HereM is an n q × n q matrix whose elements are given bỹ |β is an n q -dimensional vector whose components are the coefficients β m and δω 2 ∼ ω∆ω is the leading order correction to the frequency squared, simply obtained from (ω + ∆ω) 2 . The matrix elementsM ij correspond to contracting the disorder-induced perturbation δM of the Hessian matrix with the degenerate normalized modes/eigenvectors |ψ j , from which various scaling properties can be readily derived. First, note that the degenerate phonons |ψ i are extended objects, i.e. theirdN components are of the same order of magnitude, and are normalized, ψ i |ψ j = δ ij . This implies that the components of |ψ i scale as 1/ √ N . Second,M ij is a sum ofdN random numbers characterized by a statistical width σ, which scales as σ √ N for large N , according to the law of large numbers. The last two properties imply thatM ij scales as 1/ Finally, as δM is a second order spatial differential operator contracted with modes of well-defined elastic stiffness ω 2 , we also expectM ij ∼ ω 2 . Consequently, we predict thatM ij ∼ σ ω 2 / √ N . We therefore conclude that the eigenvalue problem defined in Eq. (2) concerns the eigenvalues of an n q × n q random matrix whose elements scale as σ ω 2 / √ N . The dependence on the degeneracy level n q is obtained from Wigner's semicircle law which predicts that the eigenvalues of the random matrixM/ √ n q follow a semicircle distribution in the large n q limit [27]. In particular, this distribution has a compact support, which implies that the eigenvalues scale as √ n q . Combining this result withM ij ∼ σ ω 2 / √ N and δω 2 ∼ ω∆ω, we obtain the final scaling prediction which is a major result of this paper.

A crossover frequency ω † (L)
The theory for the band widths ∆ω developed above makes sense as long as ∆ω does not exceed the gap between neighboring phononic bands. The crossover frequency, denoted hereafter by ω † (L), at which the gap between consecutive bands becomes comparable to the band width, is of special interest; it implies that the spectrum of disordered phonons continuously covers the frequency range ω > ω † (L), with no unoccupied holes. In this frequency range, the notion of phononic bands becomes ill-defined. To calculate the scaling behavior of ω † (L) we need to estimate the gaps in frequency between consecutive bands in the ordered system and compare it to ∆ω. To that aim, we need to express both ω † and ∆ω in terms of a single quantity that characterizes phononic bands. A natural choice would be simply the serial index q of phononic bands; as will be discussed below, the integer q is directly related to the wavenumber of members of the q th phononic band for almost -but not strictly -all integers.
We start by estimating the degeneracy level n q , which is proportional to the number of solutions to the integer sum of squares problem ind-dimensions where the p i 's and q are integers, and by definition q ≥ 0. Equation (5) is nothing but a relation between the wavenumber squared of a phononic band, which is proportional to q, and the squares of its Cartesian components. In dimensionsd ≥ 4, a solution to Eq. (5) exists for any integer q, as implied by Lagrange's four-square theorem [28]. Importantly, in three dimensions a solution of Eq. (5) exists for most integers q (over 83% as q → ∞, where the excluded integers are obtained by Legendre's three-square theorem [29]), therefore ford ≥ 3 one can treat q as a band index for scaling purposes (see some comments in this context aboutd = 2 below). For large q, n q is simply given by the surface of ad-dimensional sphere in reciprocal space, i.e. n q ∼ qd −2 2 ford > 2. The frequency of phonons belonging to the q th band follows the dispersion relation of plane-waves ω(q) ∼ √ q/L, therefore according to Eq. (4) the width of the q th phonon band reads Next, we estimate the frequency gap g(q, L) between the q th and the (q+1) th bands as a function of q. For sufficiently large q, g(q, L) is simply given by dω/dq, i.e. it follows Consequently, the gaps g(q, L) between bands and the band widths ∆ω(q, L) become comparable at band index and the crossover frequency ω † (L) follows as We expect this scaling prediction to hold ford > 2; all of the concepts discussed above are valid also ford = 2, except that a positive integer q is not a proper band index in this case, see additional discussion below (this discussion suggests that Eq. (9) is in fact valid also ford = 2). Some of the physical implications of this crossover frequency, which evidently vanishes in the thermodynamic limit L → ∞, will be discussed later in the paper.

Numerical experiments
We next put our theoretical prediction in Eq. (4) for the broadening of phonon bands in disordered solids to a direct test. We conducted computer experiments on disordered lattices in which quenched stiffness, mass, and topological disorder are introduced, in addition to measurements performed on generic models of computer glasses, in which disorder is emergent. Details about the models employed, numerical methods and numerical analyses are provided in Appendix B.

Disordered lattices
We employ simple lattices of unit masses (unless otherwise stated) connected by relaxed Hookean springs under periodic boundary conditions: square lattices ind = 2 (hereafter denoted as 2D), with springs on the diagonals, and cubic lattices ind = 3 (hereafter denoted as 3D) with the first (short) diagonals connected by springs, as shown in Fig. 1. The connectivities of the square and cubic lattices are z = 8 and z = 18, respectively. We systematically vary the degree of disorder and the system size, and measure the broadening of phonon bands.
We first assign a random stiffness to each of the lattice springs, drawn from a uniform distribution over the interval (1−σ k , 1+σ k ). The widths of phonon bands of the resulting stiffness-disordered lattices for σ k = 10 −3 , 10 −2 and 10 −1 are shown in Fig. 2. We find that ∆ω ∼ σ k ω √ n q / √ N , in perfect agreement with the theoretical prediction in Eq. (4). We next randomly select a fraction f of the springs and remove those from the lattice. The widths of phonon bands of the resulting topologically-disordered lattices for f = 10%, 3% and 1% are shown in Fig. 3. Note that in this case √ f plays the role of the disorder width f of the springs are randomly selected and removed. The outlined, full, and empty symbols correspond to f = 10%, 3% and 1%, respectively. Note that here √ f plays the role of the disorder width, as explained in the text, and hence it is used in the rescaled plots.
because the number of nonzero elements of δM is proportional to f ; this implies that the sum in Eq. (3), which involves ∼ f×dN random numbers characterized by an f -independent statistical width, scales as √ f N and consequently thatM ij scales as 1/ For our last computer experiment on lattices we randomly selected 1% of the lattice nodes, and assign them a mass that is different by a factor c m compared to the mass of the majority of the particles. If masses differ between different degrees of freedom, vibrational frequencies squared are given by eigenvalues of the matrix where i and j are particle indices, m i and m j are their associated masses, and x i is the three-dimensional vector of the Cartesian components of x pertaining to the i th particle. We measured the widths of phonon bands for c m = 0.94, 0.99, 1.03 and 1.10, and the results are displayed in Fig. 4. We find that the theoretical prediction in Eq. (4) is exactly followed, with a prefactor given by |1 − c m | ∼ ∆m, where ∆m is the difference between the minorityspecies and the majority-species masses.

Glass-derived disordered networks
We next build random networks of Hookean springs by utilizing glassy samples of a generic glass-former, see Appendix B for details about the model and the glass preparation protocol. For each glassy sample, we assign a node as the center of each particle, and connect a relaxed Hookean spring with unit stiffness between every pair of particles that interact in the original glass. This procedure leaves us with an ensemble of positionally-disordered networks of Hookean springs, with mean connectivities of z ≈ 6.5 in 2D, and z ≈ 16.5 in 3D. We then measured the phonon-band widths in these random networks. The results are plotted in Fig. 5. We find perfect agreement with the theoretical prediction in Eq. (4), despite that these random networks cannot strictly be considered as perturbed lattices.

Generic computer glass-formers
We end this section with presenting results for phonon band widths measured in a generic computer glass model. We employ a 50:50 binary mixture of point-like particles that interact via a pairwise repulsive inverse-power law interaction. This generic model has been shown [15,30] to reproduce all of the well-known phenomenology of glasses. Glassy samples were prepared by first equilibrating the system in the high temperature liquid phase, and then performing a continuous quench into the glassy phase with a finite quench rate. Further details about the model, analysis and preparation protocol can be found in Appendix B.
Two key features distinguish this generic structural glass model from the lattice models with quenched disorder discussed above; first, disorder here is spontaneously self-generated during glass formation, rather than being controlled externally. Second, frustration-induced internal stresses generically emerge in glasses [31], while they are completely absent (by construction) from the quenched disorder lattice models and the glass-derived disordered networks investigated above. These stresses have been shown to be the origin of several intriguing glassy phenomena, including poor heat transport [32], anomalous scattering of acoustic excitations in the harmonic regime [33][34][35], the existence of quasilocalized glassy modes at the low-frequency end of the vibrational density of states [13][14][15][16][17][18], and the singu- larity of athermal nonlinear elastic moduli [36]. It is therefore not a priori clear whether the presence of these internal stresses alters the scaling behavior of the phonon band widths in structural glasses or not. In Fig. 6 we plot the phonon band widths measured for different system sizes in the model structural glasses in 2D (top panels) and 3D (bottom panels). Unlike in the other disordered solids discussed above, here we find two branches when the phonon band widths of the 2D glassy samples are plotted against the rescaled frequency ω √ n q / √ N ; these correspond to longitudinal (sound) and transverse (shear) phonon bands, where each branch follows the prediction in Eq. (4), further establishing its generality and universality. The clear splitting into two branches is not observed in the other models discussed above; we suspect that it is related to the quite different shear G to bulk K elastic moduli ratio, G/K, in the different models. In the structural glass model (in 2D) we find G/K ≈ 0.15, i.e. G and K are rather well-separated, whereas in the square lattice in 2D G/K = 2/3 and in the cubic lattice in 3D G/K = 3/5, i.e. G and K are of similar magnitude. These differences, while suggestive, do not entirely explain why in the disordered lattices case the prefactors in the ∆ω scaling relations for the longitudinal (sound) and transverse (shear) phonon bands are essentially indistinguishable.
The 3D data plotted in Fig. 6c,d correspond only to the scaling of phonon band widths of transverse (shear) phonons. This is the case for two reasons: first, for small system sizes it is impossible to reliably extract the width of longitudinal (sound) phonon bands since different bands are not well-separated from each other. For larger systems, we cannot access any longitudinal (sound) phonon bands since those occur at frequencies that are too high for us to compute due to computational bottlenecks. Nevertheless, the scaling we observe , that correspond to shear and sound waves, is found to be ≈ 5.7. Note that in panel (c), the raw 3D data of ∆ω vs. ω appear to follow a power-law. This is a coincidence and the power itself corresponds to a significantly stronger-than-linear relation, as highlighted by the added guides to the eye. When these data are plotted against the rescaled frequency ω √ n q / √ N in panel (d), the predicted linear scaling is demonstrated.
for the transverse (shear) phonon band widths in 3D structural glasses convincingly follows the theoretical prediction as well. With this we conclude our numerical experiments, which verified the universal validity of the theoretical prediction in Eq. (4). We next apply the theory to two basic problems in disordered systems.

Application to phonon dynamics: Phonon scattering lifetime
Having established that the disorder-induced broadening of phonon band frequency widths follows the theoretical prediction in Eq. (4), independently of the type of disorder, we next discuss its implications for phonon dynamics. In ordered systems, phonons propagate indefinitely in the harmonic regime, i.e. in the absence of anharmonicity. The presence of disorder is expected to qualitatively change this picture, i.e. we expect pure phonons that are eigenmodes of the original ordered system to feature a finite lifetime even in the harmonic regime. In the frequency range ω < ω † (L), where phononic bands are well-defined, we expect an excited phononic state that belongs to the original degenerate band -which is not an eigenmode of the disordered system -to have significant projections only on members of its own band -which is not degenerate anymore in the presence of disorder. As the spectral width of the band in the presence of disorder is ∆ω, we expect the lifetime to be proportional to (∆ω) −1 .
To test this prediction, we perform numerical scattering simulations on harmonic square lattices with quenched stiffness disorder of strength σ k , as described in Sect. 4 4.1. At time t = 0, we introduce at the lattice positions r i (i is the particle index) velocity vectors of the formu which correspond to pure modes of the perfectly ordered lattice. Here the wavevector is given by k = (2π/L)(p x e x + p y e y ), where p x , p y are integers such that p 2 x + p 2 y = q, cf. Eq. (5), and e x , e y are the corresponding unit vectors. The normalized polarization vector a is defined by a·k = 0 and a·a = 1, i.e. we consider shear plane-waves, which for small p x , p y correspond to the lowest frequency excitations. The 2N -dimensional velocity vectoru(t), composed of the Cartesian velocity vectors {u i (t)} of all of the particles, is calculated for t > 0 using harmonic lattice dynamics and the scattering lifetime is probed through the time-evolution of the velocity auto-correlation function The mean envelopes (i.e. the amplitude of the function, disregarding its oscillatory part) of c(t),c(t), are plotted for various system sizes, wavenumbers, disorder strength and degeneracy levels in Fig. 7a, see caption for details. The theoretical prediction is that the inverse scattering lifetime, i.e. the so-called damping coefficient Γ, scales with ∆ω. Consequently, we expect the velocity auto-correlation function c(t) to collapse onto a master curve once plotted against t ∆ω. This is explicitly demonstrated in Fig. 7b, wherec(t) is plotted against where the phononic dispersion relation ω ∝ √ q/L was used. This result establishes that the phonon scattering lifetime for ω < ω † (L) is indeed determined by (∆ω) −1 .

Application to glass physics: Coexistence of quasilocalized glassy modes and phonons
The theoretical prediction in Eq. (4) for the disorder-induced broadening of phonon band frequency widths has important consequences not only for phonon dynamics, but also for non-phononic vibrational modes in structural glasses, which are of fundamental importance in the physics of glasses. It has been recently discovered [13] that non-phononic vibrational modes populate the low-frequency tails of the vibrational density of states of structural glasses, whose frequencies are distributed according to a universal gapless ω 4 law, independently of microscopic details. It is now broadly accepted that these soft non-phononic excitations are spatially quasilocalized [13,14,17], i.e. they feature a disordered core characterized by a localization length of a few atomic sizes in linear dimension and are accompanied by either a power-law decay (unlike Anderson-localized modes that appear at the high-frequency end of the phononic spectrum and decay exponentially in space [37]), or an extended, wave-like background that spans the entire system [18,38,39]. Furthermore, these soft non-phononic excitations are known to emerge from the presence of frustration-induced internal stresses in the glass [16]. The gapless ω 4 law of these soft quasilocalized glassy modes (QLGMs) had been predicted decades ago [40,41], but only recently has it been firmly established for the first time in simple, finite-size computer glass-formers.
Notwithstanding, one of the central questions left unresolved concerns the survival of QLGMs as a function of the system size L and in particular in the thermodynamic limit L → ∞. To address this important question one needs to understand how phonons cover the frequency axis as a function of L and whether QLGMs retain their quasilocalized nature when they share the same frequencies with phonons. The former was completely addressed in Sect. 2, where we have derived the system-size scaling of the crossover frequency ω † ∼ L −2/5 in 3D, above which phonon bands overlap and merge, leaving the density of vibrational modes free of gaps and holes. In particular, in the thermodynamic limit the crossover frequency vanishes ω † → 0, implying that phonons then fully occupy the entire low-frequency regime. Consequently, the possibility that QLGMs and phonons coexist in the thermodynamic limit in 3D crucially depends on whether QLGMs can share the same frequencies with phonons, without undergoing hybridizations and mixing that destroy their quasilocalized nature. This point is addressed next.
A systematic, large-scale study of hybridization and mixing phenomena of QLGMs and phonons, including its frequency and system-size dependence, is left for a future report; here we show preliminary data indicating that strong hybridizations of QLGMs with phonons of similar frequencies occur, and, consequently, no QLGMs survive within phonon bands. We generated 100 independent glassy samples of N = 10 6 particles in 3D, and calculated the first 40 vibrational modes Ψ with nonzero vibrational frequencies, see Appendix B for details about the model and the calculation. The degree of localization of the calculated vibrational modes is effectively captured by the participation ratio e, defined as where Ψ i is the three-dimensional vector of the Cartesian components of Ψ pertaining to the i th particle. Extended vibrational modes, such as phonons, are characterized by e ∼ O(1). Localized or quasilocalized modes, such as QLGMs, are characterized by e ∼ O(N −1 ). In Fig. 8 we scatter-plot the participation ratio e of the vibrational modes vs. their frequency ω; panel (a) shows the entire calculated frequency range, whereas panels (b) and (c) zoom into the first and second phonon bands, respectively. The horizontal continuous lines mark the characteristic value of the participation ratio of QLGMs, shown in [13,14] to scale as N −1 .
The dashed vertical lines in panels (b) and (c) enclose a frequency window of 4 standard deviations of each phonon band, obtained from the density of vibrational modes as described in Sect. 4 4.3 above.
The presented data indicate that within the frequency intervals covered by phonon bands, QLGMs have a very strong tendency to hybridize with phonons of those bands, which destroys their quasilocalized nature. In the presented ensemble of 100 glassy samples of N = 10 6 , we expect (see Appendix C for details about this estimation) that in the absence of phonons we would find ≈ 3.3 QLGMs within the frequency interval covered by the first phonon band, and ≈ 27 QLGMs within the frequency interval covered by the second phonon band. However, as clearly seen in Figs. 8b,c, we do not find a single QLGM within the frequency intervals covered by both the first and second phonon bands.
The strong hybridizations of QLGMs seen within phonon bands imply that they can only exist at frequencies ω ω † , i.e. below the crossover frequency, where gaps in the spectrum open between subsequent phonon bands. Can QLGMs exist at such low frequencies? The answer to this question was put forward in [13], where it was shown that the characteristic scale ω g at which QLGMs start to appear depends on the linear system size L as ω g ∼ L −3/5 ω † in 3D, and generally as L −d/5 ind dimensions. This scaling is a direct consequence of the universal distribution of QLGMs that varies as ω 4 independently of spatial dimension [17], and of the quasilocalized nature of those modes. Combining the scaling results for ω † and ω g , we conclude that QLGMs and phonons coexist for frequencies ω in the interval Our results suggest that in the thermodynamic limit L → ∞ in 3D all QLGMs cease to exist as harmonic vibrational modes, at odds with the recent claims of [18,19], that argue that these two types of low-frequency modes can be distinguished by virtue of their participation ratio in the thermodynamic limit. While quasilocalized excitations certainly do not disappear due to the presence of phonons with similar frequencies [14,42,43]as they correspond to particularly soft structures within the glass [44] -, their realization as harmonic vibrational modes (i.e. as eigenvectors of the Hessian matrix) with participation ratios of order N −1 becomes impossible due to strong hybridizations with phonons. Consequently, their detection in the thermodynamic limit by conventional techniques is hindered. The important result that harmonic QLGMs do not coexist with phonons in the thermodynamic limit is further strengthened by the existence of nonlinear QLGMs within the frequencies occupied by phonon bands, which is demonstrated in Appendix C using the framework introduced in [14].
In dimensionsd > 3, the coexistence window of QLGMs and phonons follows L −d/5 ω L −2/(2+d) , i.e. it grows in relative terms with increasing dimension. We note, however, that for dimensionsd > 5 this coexistence window extends well-below the longest-wavelength phonons, of frequencies L −1 . This means that ford > 5 QLGMs appear undisturbed by phonons between L −d/5 ω L −1 . Finally, we comment on the behavior expected in 2D structural glasses; under assumptions spelled out and motivated in Appendix D, we expect ω † ∼ L −1/2 in 2D. On the other hand, the onset of glassy modes is expected to follow L −2/5 in 2D (with possible logarithmic corrections, see discussion in [45]), i.e. it is larger than the crossover frequency ω † . This means that in 2D, above some system size, we expect harmonic QLGMs to be unobservable altogether, and no coexistence regime with phonons to exist, as indeed reported in [18].

Summary and discussion
In this work we have theoretically derived the dependence of the disorder-induced frequency widths of phonon bands ∆ω on the strength of disorder, on the number of phonons in a band, on the band's frequency, and on the number of particles in solids. Our result, which is obtained using degenerate perturbation theory and simple statistical considerations, was then validated against computer experiments on simple lattices with quenched stiffness, mass, and topological disorder, and on generic computer models of structural glasses, both in two and three dimensions, and for a wide variety of system sizes. In all cases we found excellent agreement between our theoretical prediction and the numerical data, establishing the universal nature of the theory. Of particular interest is the robustness of our scaling theory to the presence of self-generated disorder and frustration-induced internal stresses in structural glasses, which are known to affect many phonon-related glassy phenomena [32,35].
The derived broadening of phonon band widths with increasing frequency gives rise to the identification of a crossover frequency ω † ∼ L − 2 d+2 in a system of linear size L ind > 2 dimensions, above which the notion of discrete phonon bands becomes ill-defined. Instead, above ω † phonons are expected to cover the entire frequency axis, leaving no holes or gaps. A first implication of these results is that in the frequency range ω < ω † , where phonon bands are well-defined, phonons exhibit a finite scattering lifetime that scales with (∆ω) −1 , as we directly demonstrate by dynamic calculations for harmonic disordered lattices.
We further discussed a key implication of our results on the persistence of harmonic quasilocalized glassy modes (QLGMs) -recently shown to populate the low-frequency tails of the spectrum of 3D structural glasses -in the thermodynamic limit. We presented extensive computer simulation data that suggest that QLGMs loose their quasilocalized nature if they occur within frequency intervals occupied by phonon bands, due to hybridizations and mixing with those phonons. This, in turn, implies that a coexistence frequency window of phonons and harmonic QLGMs opens at intermediate system sizes, explaining the observations of [18,19] that report coexistence of these two types of low-frequency excitations, distinguished by their participation ratio. Our results further indicate that the said coexistence frequency window vanishes in the thermodynamic limit, seriously questioning the claim in [18,19] that coexistence persists in the continuum (i.e. thermodynamic) limit.
Our work opens up various directions for future investigations. For example, the prefactor in the main theoretical result in Eq. (4) might offer a general, dimensionless quantifier of disorder in any condensed matter system featuring Goldstone modes. In the structural glasses analyzed in this work we found χ ∼ O(1) (see Fig. 6), whereas in the positionally-disordered (glass-derived) Hookean-spring networks we found χ ∼ O(10 −1 ) (see Fig. 5). Future investigations should resolve the relative variations observable in χ under different preparation protocols in structural glasses and disorder realizations in disordered lattices/crystals. Moreover, the prediction that the phonon damping coefficient should cross over at ω † from being proportional to ∆ω to following the Rayleigh scattering scaling (or the modified-Rayleigh scattering scaling in glasses [35]) should be systematically tested. Several recent efforts, e.g. [8][9][10][11][12], are devoted to resolving disorder-induced variations in the statistics of phonons of a prescribed wavevector. These studies exclusively focus on frequencies ω ω † , and show that the major effect of disorder is observed at intermediate to high frequencies. While our theoretical approach is not valid above the crossover frequency, where the notion of discrete phonon bands becomes ill-defined, it would be interesting to resolve whether, if at all, our findings bare relevance to this closely related question.
In our theoretical treatment of the phonon band broadening we assumed that the disorder in the system is not spatially correlated; it would be interesting to understand the effect of spatial correlations of structural disorder on our results. In addition, other applications to a broad range of phenomena involving phonons in disordered systems should be systematically explored.
In this Appendix we briefly derive Eqs. (2)-(3), which emerge from standard degenerate perturbation theory. The latter formalism can be found in many Quantum Mechanics textbooks [26] and is briefly repeated here for completeness. We start with n q degenerate eigenvectors |ψ i , with i = 1, 2, . . . , n q , corresponding to an eigenvalue ω 2 , and denote all other eigenvectors by |ψ k , with k = n q + 1, n q + 2, . . . ,dN . The unperturbed Hessian matrix is denoted by M (0) such that any linear combination |Ψ ≡ nq m=1 β m |ψ m is an eigenvector where the relevant combination |Ψ for the perturbation theory is not known a priori. The presence of disorder modifies the Hessian matrix, which now reads M = M (0) +δM, where the disorder-induced perturbation δM gives rise to corrections δω 2 to the eigenvalue ω 2 and |δΨ to the (still unknown) eigenvector |Ψ . Perturbing Eq. (16) accordingly and keeping terms to linear order, we obtain The correction to the eigenvector |δΨ can be expressed in terms of the eigenvectors outside the degenerate band Substituting Eq. (18) in Eq. (17) and contracting with an eigenvector ψ i | within the degenerate band, we obtain nq j=1 which is identical to Eqs. (2)- (3). Note that if instead we contract with an eigenvector ψ l | outside the degenerate band, we can calculate the coefficients α l and hence the correction to the eigenvector |δΨ , though we do not study the latter in this paper.

Appendix B Numerical methods
In this Appendix we describe the procedures we have followed in the analysis of phonon band widths, we describe how the positionally-disordered networks were created, and we describe the model glass-former we employed.

B.1 Phonon band widths calculation
The calculation of vibrational modes was carried out using the linear algebra package ARPACK [46]. In the majority of the calculations performed, extracting phonon widths from the spectra can be done by directly identifying the bands and calculating the standard deviation of the frequencies associated with the members of each band, averaged over several (100 or 200) realizations of the disorder. This is possible since the bands at low frequencies are well-separated, as demonstrated in Fig. 9a. shows the raw histograms of the eigenfrequencies, collected over several independently-quenched glassy samples. Panel (c) displays the calculated Gaussian fits; notice that for these fits we evaluate the histogram for each peak separately, using the same number (12) of bins for each peak. This explains their different heights compared to the peaks of the raw histogram of panel (b). The peak marked by the arrow at ω ≈ 0.116 represents the first sound modes, which explains its minute width.
In our 2D and 3D model glasses a sharp identification of bands is not possible due to the presence of non-phononic modes at low frequencies. We therefore calculated the histogram of the vibrational frequencies, and fitted a Gaussian curve to each peak that was cleanly distinguishable, as demonstrated in Figs. 9b,c. For our 3D glasses, we followed the approach of [18] and calculated histograms of vibrational frequencies after filtering phonons from nonphononic modes by discarding of modes with participation ratios e < 0.03. The standard deviation of the Gaussian fits are reported as the phonon band widths.

B.2 Disordered networks of Hookean springs
To create disordered networks of Hookean springs, we utilized the computer glass model that is described in the next subsection. We set the center of each particle as a node, and connected a relaxed Hookean spring of unit stiffness between every pair of interacting particles of the original glass. The mean connectivities of the resulting spring networks we obtained following this procedure are z ≈ 6.5 in 2D, and z ≈ 16.5 in 3D, i.e. very far from the Maxwell threshold 2d, and, subsequently, far from the unjamming point [47].

B.3 Model glass-former
We employed a 50:50 binary mixture of 'large' and 'small' particles that interact via a purely-repulsive inverse-power-law (∼ r −10 ) pairwise interaction potential. Details of this model system can be found in [15]. In terms of the microscopic units of energy and time τ , the system undergoes a computer glass transition at a temperature T ≈ 0.50 /k B . Systems were initially equilibrated in the high temperature liquid phase at T = 1.0 /k B , before quenching them to zero temperature at a rate ofṪ = 10 −3 /(k B τ ).

Appendix C QLGMs within phonon bands
In this Appendix we first explain how the estimation of the number of QLGMs we would expect to find in a frequency interval ω 1 ≤ ω ≤ ω 2 is made. Recall that frequencies associated with QLGMs are distributed according to a universal ω 4 law. The prefactor of the ω 4 law, c g , can be extracted from the data of Fig. 8 by e.g. counting the total number n QLGM of QLGMs that were measured over n samples glassy samples between the two frequencies ω 1 and ω 2 . That is, we have .
We find c g ≈ 7.6×10 −4 by considering the number of hits between ω 1 = 0.27 and ω 2 = 0.33, see Fig. 8. The number of QLGMs expected to fall between the dashed horizontal lines of Figs. 8b,c is then easily obtained by inverting Eq. (20) in favor of n QLGM , using the previously extracted prefactor c g . We followed this procedure to reach the expectation that we would observe ≈ 3.3 QLGMs within the first phonon band, and ≈ 27 QLGMs within the second phonon band. We also calculated nonlinear QLGMs; following the definitions put forward in [14], nonlinear QLGMs π are solutions to the equation where the symbols :, .
: and :: denote double, triple and quadruple contrations, respectively, M ≡ ∂ 2 U ∂x∂x , and U ≡ ∂ 4 U ∂x∂x∂x∂x are the second and fourth order tensors of spatial derivative of the potential energy U, respectively. For each member of the ensemble of glassy samples of N = 10 6 particles in 3D we calculated a single solution to Eq. (21), following the methods explained in [48]. The frequency of each solution π is defined as ω ≡ √ M : ππ, in analogy with the frequencies of vibrational modes.
Nonlinear QLGMs were shown in [14] to (i) possess very small frequencies, (ii) to be quasilocalized in the same fashion as harmonic QLGMs, and (iii) to be entirely indifferent to the presence of phonons of similar frequencies. For these reasons, they serve as a useful tool to identify QLGMs within the frequency intervals occupied by phonon bands. In in the main text), plotted against their frequency ω. Data were obtained in systems of N = 10 6 particles in 3D; the left and right panels show data for the first and second phonon bands, as seen in Fig. 8. The orange diamonds represent nonlinear QLGMs (see text for definitions and details), whose participation ratio precisely matches that of QLGMs outside of phonon bands, as indicated by the horizontal continuous lines, which are deduced from Fig. 8a. This data set demonstrates that the soft, quasilocalized excitations with frequencies that fall inside the intervals occupied by phonons are still embedded in the microstructure, despite that they cannot assume the form of harmonic vibrational modes. Fig. 10 we plot the same data as in Fig. 8b,c, but this time superimpose the frequencies of nonlinear QLGMs (marked by orange diamonds in the figure) detected in the same glassy samples. Since methods to exhaustively detect all nonlinear QLGMs are still unavailable, the number of nonlinear QLGMs we detected within the frequency intervals occupied by the two first phonon bands falls below our prediction as spelled out above. Notwithstanding, the demonstrated existence of nonlinear QLGMs within these frequency intervals unequivocally shows that (i) QLGMs do not assume the form of harmonic vibrations within phonon band frequency intervals due to hybridizations, and (ii) soft quasilocalised excitations whose frequencies fall within those intervals are certainly embedded in the microstructure.
Appendix D The crossover frequency ω † in 2D In Sect. 3 we derived the scaling of the crossover frequency ω † ∼ L −2/(d+2) for a disordered solid of linear size L ind > 2 spatial dimensions, above which phonon bands overlap and merge, leaving the density of vibrational modes free of gaps. Resolving the L-dependence of ω † in 2D is less straightforward, as a result of the erratic behavior of the degeneracy n q of the q th phonon band, which hinders a scaling estimation of the broadening of phonon band widths ∆ω. In particular, despite that q q =0 n q ∼ q , n q is not a constant function of q in 2D. Instead, n q features two trends with increasing q, illustrated in Fig. 11; first, the envelope of n q (see figure caption for definition) increases with increasing q, as shown in the left panel of Fig. 11. On the other hand, the number of integers q for which n q = 0 grows with increasing q as well [49], many of which occupy consecutive intervals of integers, referred to here as holes. The latter, denoted by h q , are defined for integers q, that feature finite degeneracies n q > 0, as the number of consecutive integers q > q for which n q = 0. In Fig. 11b we plot the envelope of h q , which is seen to grow in a similar manner as the envelope of n q . The appearance of growing holes in n q hinders the clean estimation of the scaling of the frequency gaps between consecutive phonon bands with band index q and system size L, as we have done ford > 2 in Eq. (7) of Sect. 3. The circles represent the envelopes defined as pairs (q, n q ) such that n q > n q for all q < q, while the squares represent running averages binned over q. (b) Intervals h q of consecutive integers q for which n q = 0, defined for integers q that feature n q > 0, see text for further details. The envelopes and running averages are as defined in panel (a).
In Fig. 11 we have also plotted the running averages of the degeneracies n q (squares, left panel) and of the holes h q (squares, right panel), binned over logarithmic intervals of the band index q. The running average of n q is a constant of order one, whereas the running average of h q grows slower than logarithmically. We therefore neglect their weak q-dependence altogether, and assume that n q and h q are constants.
Under the assumptions discussed above, which we cannot fully justify, our estimation of the broadening of phonon band follows verbatim Eqs. (6)-(9), and thus for 2D disorder solids we expect the crossover frequency ω † ∼ L −1/2 . This prediction is in fact consistent with Eq. (9), whend = 2 is used.