Impact of elastic heterogeneity on the propagation of vibrations at finite temperatures in glasses

Some aspects of how sound waves travel through disordered solids are still unclear. Recent work has characterized a feature of disordered solids which seems to influence vibrational excitations at the mesoscales, local elastic heterogeneity. Sound waves propagation has been demonstrated to be strongly affected by inhomogeneous mechanical features of the materials which add to the standard anharmonic couplings, amounting to extremely complex transport properties at finite temperatures. Here, we address these issues for the case of a simple atomic glass former, by Molecular Dynamics computer simulation. In particular, we focus on the transverse components of the vibrational excitations in terms of dynamic structure factors, and characterize the temperature dependence of sound dispersion and attenuation in an extended frequency range. We provide a complete picture of how elastic heterogeneity determines transport of vibrational excitations, also based on a direct comparison of the numerical data with the predictions of the heterogeneous elastic theory.


Introduction
The notion of phonon [1,2] together with the Boltzmann transport equation [3][4][5] efficiently describes vibrational excitations in crystals, ranging from the harmonic approximation to temperatures where anharmonic couplings induce finite life-times of the former. (The Akhiezer effect [6,7], which is strictly valid for low-frequency sound modes and produces a ∝ Ω 2 T dependence for the anharmonic sound attenuation, will be mentioned below.) The situation is comprehensively more complex for disordered solids [8][9][10][11][12], in particular glasses, where the anharmonic effects coexist with those due to the structural disorder. In these conditions, the phonon-based picture is recovered in the continuum limit, at small wave-vectors (frequencies). By contrast, when mechanisms taking place at (nano-)scales comparable to the typical atomic distances are relevant, more sophisticated concepts are needed.
Experimental evidences point to anomalous acoustic excitations in glasses in the THz-GHz regime, including sound softening [13] and Rayleigh-like strong scattering [14][15][16]. Among many other works, our [17] has also provided insight on these matters, by numerical simulations. For instance, we highlighted a clear crossover in the apparent inverse life-time of acoustic-like excitations, 1/τ = Γ, from a Rayleighlike scattering, Γ ∝ Ω 4 , to a disorder-induced broadening, Γ ∝ Ω 2 , moving from low to high frequency. We also found that the crossover frequency, Ω co , is not far from the Ioffe-Regel frequency, Ω IR , implying that sound waves start to loose their plane wave character at Ω co . In addition, in the same Ω-regime, peculiar features such as excess vibrational states, the so-called Boson Peak (BP) [18,19], were observed. Since sound waves in glasses are described as superpositions of different vibrational modes [20], these features are strongly related to the anomalous sound waves propagation. As an example, an universal connection between transverse sound waves and the BP has been proposed [17,21]. Rationalizing these

Methods
The model and simulation details. We have studied by Molecular Dynamics (MD) simulations glassy systems formed by N mono-dispersed point particles, of mass m and diameter σ, interacting via a pairwise Lennard-Jones (LJ) potential, V(r) = 4 [(σ/r) 12 − (σ/r) 6 ], with r = r i j the distance between particles i and j. The potential V(r) is cut-off and shifted at r c = 2.5σ. 1 m, σ, , and τ = (mσ 2 / ) 1/2 are the units of mass, length, energy, and time, respectively, while temperature is in units of /k B (k B is the Boltzmann constant). If we consider the case of argon as a reference, temperature is in units of /k B = 125.2 K, while length and time scales are expressed in units of σ = 3.405 Å and τ = 2.11 ps. We have used periodic boundary conditions in cubic simulation boxes of size L, such that the number density isρ = N/L 3 = 1.015. Atρ, the melting and glass-transition temperatures are T m 1.0 and T g 0.4 [52], respectively. In order to access the relevant small wave vector (q) region relevant here, we have considered N = 4000 to 1000188, corresponding to L = 15.80 to 99. 51. In what follows, we show data pertaining to all values of N together, immediately verifying the absence of any finite size effects.
Initialization runs were conducted in the micro-canonical (NV E) ensemble at the temperature T = 2.0 in the liquid phase, followed by a fast quench with a rate dT/dt ≈ 400 down to T = 10 −3 , well below the glass-transition temperature T g . Next, the systems were heated to T = 10 −2 , 10 −1 , 2 × 10 −1 , still below T g . At all values of T, the quenched glass samples were relaxed for a time (dependent on N) sufficient to eliminate any drift in the total energy, and recover energy fluctuations conforming to the equipartition theorem, with a specific heat per particle c v 3. Following thermalization, we performed the production runs for a (N-dependent) total time sufficient to obtain the desired ω-resolution, always well below the smallest calculated line width. We have used the velocity Verlet algorithm for the integration of the 1The discontinuity of the derivative V (r) at r c can, in general, modify some harmonic vibrational features at T = 0 [50]. Based on previous works [17,51], however, we do not expect notable modifications in sound wave propagation at finite temperatures.

43604-2
equations of motion, with a time step δt = 5 × 10 −3 . We emphasize that no aging effects were observed in the time evolution of total energy throughout the production runs. All simulations have been realized with the large-scale massively parallel MD simulation tool LAMMPS [53].
The dynamical structure factor. Sound waves propagation has been investigated in terms of the transverse component of the dynamical structure factors [17,35] at wave vector q and frequency ω, Here, q = |q|,q = q/q, and is the thermodynamic average. Note that the simulation box with N = 1000188 accommodates vibrational modes with q = 2π/L 0.06. Since the first diffraction peak in the static structure factor, S(q), is at q m 7 (corresponding to an average nearest-neighbors distance 2π/q m = 0.9), some of our spectra refer to q-values ≈ 10 2 times smaller than the border of the pseudo-Brillouin zone at q q m /2.
The vibrational density of states. In the discussion we will also consider the vibrational density of state (vDOS), g(ω) [1,2]. This was determined by numerically diagonalizing the Hessian matrix (second derivatives of the Hamiltonian), calculated from the coordinates of systems completely relaxed in local potential energy minima (inherent structures). The g(ω) was then calculated by populating an histogram with the square-root of the obtained (positive) eigenvalues. Similar to S T (q, ω), we have considered system sizes up to N = 256000, in order to adequately sample the g(ω) on all the relevant ω-range [17].
The spatial distribution of the local shear modulus. In our discussion we build over the concept of elastic heterogeneity, i.e., the existence of probability distributions of local elastic moduli. In systems like our LJ glass, the bulk (K) and shear (G) moduli are such that K G [54], and the shear modulus mostly determines the low-frequency transverse modes behaviour [34][35][36]. We, therefore, focus on the probability distributions of the local shear modulus G m only, which are determined by partitioning the simulation box into an array of cubic domains of linear size w 3.16, identified by an index m. We note that a domain includes on average around 30 particles. For each m, G m was computed by the fluctuation formula [23,55,56] Here, G m A encodes the elastic response to affine deformations, where the particles follow the applied affine strain field. By contrast, the non-affine modulus, G m N , which contributes negatively to the overall modulus, includes the effects associated to additional particle displacements that deviate from the applied affine field. We have followed the same procedure indicated as the "fully local approach" in [57]. We add two observations. First, since we have demonstrated that the G m calculation is insensitive to finite system size effects [57], we have used N = 4000 (L = 15.80). Second, in the calculations we have shifted potential and its first derivative, so that both potential and forces are continuous at r = r c . Indeed, in [51] we demonstrated that the local modulus values are modified by non-linearities associated to discontinuities in the forces at r c .

Characterizing the transverse dynamical structure factor
The transverse dynamical structure factor. In figure 1, we show a selected set of data for S T (q, ω), calculated from equation (2.1) for N = 1000188, at the lowest wave number, q = 0.06. The indicated values of T, all below T g 0.4, increase from top to bottom. The S T (q, ω) spectra are characterized by two symmetric Brillouin peaks, in addition to the elastic feature at ω = 0. As T increases, the Brillouin peaks move towards higher frequencies, with an increasing total intensity and broadening, as expected. We can extract quantitative information from these data by the customary procedure of fitting the points in the spectral region around the Brillouin peaks to the damped harmonic oscillator (DHO) model [58], Here, the parameters Ω(q) and Γ(q) represent the characteristic frequency and inverse life-time (or broadening, full width at half maximum) of the Brillouin excitations, respectively, which we discuss below. In what follows, we will systematically dump the subscript T.

The sound velocity.
From the values of Ω(q) we can obtain the transverse sound phase velocity, c(q) = Ω(q)/q, shown in figure 2 (a)-(d) (symbols) as a function of the corresponding Ω, at the indicated values of T. In (a) we show the data at the lowest temperature, T = 10 −3 , with a velocity increasing with the frequency (the positive dispersion [59]) for Ω > 1. Next, as already observed in [17], the velocity shows a minimum before reaching a region of softening, where it increases by decreasing Ω. Eventually, the sound velocity saturates at the correct acoustic value, where G is the shear modulus, and ρ = mρ is the mass density. The same pattern is still qualitatively present at T = 10 −2 in (b), although the acoustic limit is reached at a slightly higher frequency. At T = 10 −1 in (c), softening has almost completely disappeared, and a direct crossover from the positive dispersion to the zero frequency limit can be recognized. At the highest temperature, T = 2 × 10 −1 T g /2 in (d), the behaviour just discussed is even more evident, with an acoustic limit showing a mild T-dependence.
Note that a very similar scenario holds for the longitudinal modes at all values of T (not shown here). Indeed, in reference [17] we have suggested a common origin for both longitudinal and transverse excitations, based on the isotropic elastic medium equation c 2 L (ω) K/ρ + (4/3)c 2 T (ω), where we have discarded the frequency dependence of K(ω), with the replacement of the zero-frequency limit, K = K(ω = 0) 59. This implies that the knowledge of one of the two velocities allows one to infer the ω-dependence of the other by the above relation. We have checked that this result holds at all values of T, and similar observations were also reported in the simulation works of [45,60]. We note, however, that this is probably not a general conclusion since in [35], for instance, we have demonstrated that, for soft spheres, K does depend on frequency and the above equation is therefore not valid.
The sound broadening. Our results for the transverse sound broadening, of unprecedented quality, are shown in figure 2 (e)-(h), where we plot Γ as extracted from equation (3.1), as a function of the corresponding Ω. The frequency and temperature dependences of these data are very complex [49]. At the lowest T = 10 −3 in (e), a clear crossover occurs between the high-frequency disorder-controlled behaviour ∝ Ω 2 [61], and a Rayleigh-like scattering contribution, ∝ Ω 4 , at lower frequencies [62]. As already noticed, the crossover frequency Ω co 1, is not far from the calculated BP frequency, ω BP 2, which is compatible with the Ioffe-Regel frequency, Ω IR ω BP , as also discussed in [17] 2. We note that, even at this very low T, anharmonic interactions are obviously present and, for instance, still contribute to the thermal conductivity. Their intensity, however, is very low in the investigated frequency range compared to other contributions, while non-negligible effects should be visible at frequencies lower than 2In [45] we have shown that the crossover occurs at the frequency ω 0 where the vDOS converges to the Debye limit. For the present LJ system, ω 0 1 [17], confirming Ω co ω 0 .  our available spectral window. By increasing T, by contrast, we expect the strength of anharmonicities to increase, by eventually entering our accessible frequency range.
This is indeed what we observe in (f), where at T = 10 −2 we detect a second T-dependent crossover, at Ω co2 0.6, between the Rayleigh region and an unexpected low-frequency ∝ Ω 3/2 regime, reminiscent of the fractal-like attenuation of [63,64] (see below). Note that this fractional scaling mechanism is evidently T-dependent, whereas the Rayleigh and disorder-controlled regimes are not. Also, by increasing T, we expect the two crossover frequencies to eventually merge (Ω co2 Ω co ) when the strength of the anharmonic couplings becomes comparable to that associated to the effect of the disorder, and the two mechanisms entangle in the entire Ω-range. This is indeed what we observe at T = 10 −1 in (g), where the Rayleigh contribution is taken over by the fractal ∼ Ω 3/2 behaviour at low frequency. In addition, a quadratic T-independent contribution is still visible at high Ω. Finally, at the highest T = 2 × 10 −1 T g /2 in panel (h), the intensity is fully T-dependent in the entire Ω-range, with an Akhiezer-like Ω 2 scaling prevailing in a large frequency range, while a remnant Ω 3/2 behaviour can be only recovered at low Ω. This scenario is similar to that reported recently in the experimental work of reference [65] for a network glass (sodium silicate), confirming the relevance of our findings for realistic systems.
In [49] we have shown how to rationalize the Ω 3/2 fractional scaling described above and predicted in the context of the HET [63,64], by disentangling the effects due, on the one hand, to the anharmonicities, and the elastic heterogeneity on the other hand. Below, we focus on characterizing the local elastic  heterogeneity at variable T, and linking the latter directly to the spectroscopic parameters in terms of the above theory.

Uncovering the elastic heterogeneity
The local shear modulus distributions. We now discuss our numerical determination of the local shear modulus distributions, recalled in section 2. In figure 3 (a) we show the probability distributions P(G m ) of the local shear modulus G m at the indicated temperatures. In the figure, we also plot P(G m ) at T = 0, which is obtained with the harmonic approximation formulation [66,67]. We confirm that the distributions can be fitted by Gaussian functions at all values of T, and the average values are close to the T = 0 value (red dashed line) and T-independent, as demonstrated in figure 3 (b) (red circles). By contrast, although for T < 10 −1 , the widths of the P(G m ) stay close to the T = 0 value, for higher temperatures, the distributions substantially broaden [figure 3 (c), red circles]. Note that an analogous behaviour is hold by the distributions of the purely affine components, G m A (blue squares). In figure 3 (d) (red circles), we plot our data in a different representation which will be useful below, and show the parameter γ(T) =ρw 3 δG 2 /G 2 , involving the ratio of δG over G, normalized to the CG domain size w. This obviously follows the same behaviour as above, staying close to the T = 0 value (red dashed line) for T < 10 −1 and next increasing rapidly when thermal fluctuations set in and induce substantial broadening. In the figure we also show the affine contribution to the shear moduli, γ A (T) =ρw 3 δG 2 A /G 2 A (blue squares), with again a T-dependence analogous to that of γ(T). Note that γ(T) (γ A (T)) plays a crucial role in the HET [37-39, 63, 64], where the parameter δγ = γ c − γ (inset  figure 3 (d), with γ c = 0.226) controls the approach to the mechanical instability at γ c and induces the fractal frequency dependence of broadening, Γ ∝ Ω 3/2 . We will address this point more in detail in section 3.3.
Spatial correlations and the non-affine displacement field. We have described above a clear crossover of broadening at Ω = Ω co 1, where the scattering mechanism is modified from Γ ∝ Ω 4 to ∝ Ω 2 . As we already mentioned, Ω co 1 is not far from the Ioffe-Regel frequency Ω IR 2 [17]. As a consequence, at Ω co , the character of the vibrational excitations is substantially modified, crossing-over from quasi-plane waves to completely disordered envelopes of vibrational modes. We can associate to the crossover frequency a characteristic length-scale ξ co (or a wave number q co = 2π/ξ co ), with ξ co = 2πv T /Ω co 23 (q co = Ω co /v T 0.27) [32]. We stress that both Ω co and the length-scale ξ co are T-independent and related to the T = 0 structural properties only. Here, we attempt to characterize ξ co in terms of the local shear modulus distributions.
In figure 4 (a), we show the spatial correlation function of G m , defined as Here, we have explicitly represented G m (r) as a function of the position r, and the terms [G m (r) − G] are the fluctuations around the average (macroscopic) value, G. 0 denotes the spatial average over r = 0 and r = |r|. Clearly, the C G m (r) decay on length scales which, although mildly increasing with T at the highest temperatures, always stay very close to the CG domains size w, implying no appreciable spatial correlations. We, therefore, conclude that G m fluctuates in space randomly, without any spatial correlation length, as also demonstrated for a-thermal jammed solids, both numerically [67] and in experiments [68]. We now focus on the non-affine displacement field of particles. As reported in [28], randomly fluctuating local elastic moduli generate a non-affine displacement field 3, which exhibits a vortex-like structure. Interestingly, while the local elastic moduli show no correlations on length scales larger than that associated to the CG procedure, the non-affine displacement field is in contrast characterized by a typical length ξ na . In [32], it has been estimated that ξ na amounts to about 23 particles size (diameter) for a 3-dimensional LJ system similar to the one studied here. It would therefore turn out that ξ na ξ co , extracted from the total transverse sound broadening, a point highlighted in the experimental work of [69].
To verify this point, we have extracted ξ na from the non-affine displacement field of our system, by following the procedure of [32]. This is based on a quasi-static shear deformation, where we first apply a small affine shear strain (γ = 5 × 10 −4 ), and next relax the system to the closest local potential 3More precisely, fluctuations of the affine elastic moduli, which are also randomly distributed in space [57], drive the non-affine displacement field. In the inset we represent by arrows the non-affine displacement field, δũ i , within an arbitrary layer, for a system with N = 64000.

43604-7
energy minimum. Throughout the relaxation process, we keep track of the total non-affine displacement, δu i , for all particles i. As shown in [29], additional elastic waves are excited during the non-affine relaxation, inducing system-spanning correlated motions of particles. It is, therefore, crucial to remove this contribution, by re-defining the non-affine displacement field as, δũ i = δu i − k; ω k <ω 0 δu i · e k i e k i , with ω k and e k i the eigenfrequencies and eigenvectors extracted by a normal mode analysis of the Hessian matrix [17]. Note that we set ω 0 = 1 where the T = 0 vDOS converges to the Debye limit [17].
We represent the δũ i with arrows in the inset of figure 4 (b) where, by visual inspection, we indeed clearly recognize the vortex-like structure. To quantify the representative size of the latter, we have computed the spatial correlation function C δũ (r), with r = r i − r j , and denotes the average over all pairs of particles, i and j. [Note that our C δũ (r) does not show system-size effects, in contrast to the results of [29], since we removed the elastic waves contributions described above]. The data of figure 4 (b) show that C δũ (r) decays on length scales r = ξ na 20, consistently with [32] and very close to the ξ co 23 associated to Ω co as discussed above. This observation ultimately allows us to rationalize the value of the crossover frequency Ω co in terms of the shear modulus heterogeneities as follows: i) The spatial fluctuations of the local shear modulus induce the non-affine character of the displacement field; ii) The latter is modulated in space, with correlations extending to λ < ξ co , while for λ > ξ co correlations are lost; iii) As a consequence, the displacement field interacts with transverse modes differently for long and short length scales, corresponding to low (Ω < Ω co ) and high (Ω > Ω co ) frequencies, respectively; iv) These distinct effects induce the observed distinct scattering scalings, Γ ∝ Ω 4 and ∝ Ω 2 , together with the associated cross-over at Ω co .

Putting everything together: The heterogeneous elasticity theory
Sound broadening and the Boson peak. In reference [17] we demonstrated a possible connection between the sound softening encoded in the pseudo-dispersion curves and the BP, by assuming q as a good parameter for labelling vibrations in glasses and counting the number of acoustic modes in the low-q region. This procedure reproduced the BP feature in the reduced vDOSĝ(ω) = g(ω)/ω 2 . Here, we adopt a different point of view, based on the HET and developed in reference [60] (the so-called generalized Debye model [45]), which directly connects the transverse broadening Γ to the BP features, .
4.13, and ω D = q D c D 16.19 are the the Debye wave-number, velocity and frequency, respectively, with v L = (K + 4G/3)/ρ 8. 71 and v T = G/ρ 3.65 the longitudinal and transverse sound speeds. The term 3/ω 3 D is the Debye limit. Hence,ĝ(ω) can be determined by equation (3.4) if the sound speed, c(ω), and broadening, Γ(ω), are known, including all anharmonic effects at finite T.
In figure 5 (a)-(d), we plot the right-hand side of equation (3.4) at the indicated values of T (closed symbols), together withĝ(ω) at T = 0, obtained by diagonalizing the Hessian matrix [17] (black open circles). The two sets of data in panel (a) are in nice qualitative agreement, with equation (3.4) grasping the Debye limit, increasing quadratically in the Rayleigh range, and matching the BP intensity at ω BP 2. The situation is similar at T = 10 −2 in (b), although anharmonicities already start to alter the smallω behaviour, a modification which is completed at the two highest temperatures, T = 10 −1 (c) and 2 × 10 −1 (d). Now, the data decrease by increasing frequency, following an ω −1/2 power-law, before directly matching the BP intensity. The simultaneous presence of the ω −1/2 dependence and the BP corroborates the predictions of [64], where the theory was modified to include an anharmonic scattering contribution [70,71] responsible for the observed fractal behaviour. Similar data have been reported in the experimental work of [65]. Adjusting the HET to the numerical data. Equation (3.4) (and, therefore, the HET) indeed provides a convincing framework to describe our T-dependent vibrational data. It misses, however, an important aspect: it does not directly (and transparently) involve any quantity related to the local elastic properties, which ultimately control the observed physics. Note that in the theory these are adjustable parameters. In our simulations, by contrast, they can be measured (figure 3) and are directly determined by the interaction potential and by the imposed external thermodynamic conditions. It, therefore, makes sense to fit the HET to our sound data, and next compare the obtained best-fit parameters to the elasticity measurements, providing a complete consistency check of the entire picture.
We have self-consistently solved the momentum conservation equation in the effective medium approximation [37][38][39]. In particular, we followed reference [64], and solved the equations: where χ T,L (q, ω) are the transverse and longitudinal dynamical susceptibilities, respectively. The theory assumes an anharmonic coupling in the form of the Akhiezer-like scattering, amounting to an imaginary term Σ anh (ω, T) ≡ iνωT [70,71]. Equations (3.5) contain four free parameters: the shear (G 0 ) and bulk (K 0 ) moduli, the disorder parameter γ =ρw 3 δG 2 0 /G 2 0 , and the anharmonic coupling ν. Since Σ is directly  related to sound velocity, broadening and vDOS [37][38][39], as (where Re and Im denote the real and imaginary parts, respectively), we can tune the parameters to best reproduce those (MD) data via equations (3.5). In figures 2 and 5, we plot by solid lines the predictions for c(Ω), Γ(Ω), andĝ(ω) for all values of T, based on the procedure above. At the lowest T = 10 −3 , the theory acceptably reproduces the most distinctive features of the vibrational excitations, including the sound softening in c(Ω), the crossover in Γ from the ∝ Ω 4 to the ∝ Ω 2 scaling, and the BP intensity and position inĝ(ω). While in general the intensity of the broadening is considerably underestimated, the theory is capable of pin-pointing the crossover frequency, Ω co ≈ 1, and thus the length scale, ξ co ≈ 20. Note that these values are directly determined by the parameter γ, since we do not explicitly insert those scales into the theory. We also note that the effective medium theory for a-thermal systems [40,41] also strongly underestimates the attenuation, although correctly predicting the crossover frequency and length scales [45].
In addition, the HET also reasonably reproduces the effects due to anharmonicities at finite T visible in our data, accounting for the disappearance of the sound softening in the c(Ω), the increase of anharmonic damping with temperature, and the fractal scaling laws ∝ Ω 3/2 in Γ and ∝ ω −1/2 inĝ(ω) at the highest T = 2 × 10 −1 . In [49] we have been able to isolate the effects related to the anharmonic couplings and to demonstrate a clear correlation between the temperature variation of γ(T) (related to the braoadening the distributions due to an increasing fraction of negative shear stiffnesses on increasing T) and that (linear) of the strength of the Akhiezer-like mechanisms.
We can now provide additional insight by comparing the above best-fit values for G 0 and K 0 to those calculated from the MD configurations. In table 1 we report the values extracted by means of the indicated procedures [57], for a CG scale w 3.16. We have considered: i) the values shown in figure 3 (b), (see section 2); ii) the affine component values [also shown in figure 3 (b)], G A and K A , where the non-affine contributions have been neglected; and iii) the values extracted via the frozen matrix approach [72]. In the latter, the local elastic properties are determined by freezing all atoms located outside a target local region. The frozen part is next constrained to only relax affinely, amounting to a system restricted to only localized non-affine deformations [57].
We immediately note that the optimized value for K 0 is matched by all methods, while completely disregarding the (negative) non-affine component, provides a G 0 which is strongly overestimated, as expected. On the contrary, we find that the fully local approach strongly underestimates G 0 , which is not completely surprising either. Indeed, the HET does not explicitly treat the non-affine contributions to the 4For the case of T = 10 −3 , we set ν = 0 because we did not recognize any anharmonic contribution to the simulation data of figure 2. Table 1. Values of the HET parameters G 0 , K 0 , and γ, as extracted by fitting the theory to the MD data, or determined directly from the MD configurations by the computational methods summarized in [57]. The values of γ are estimated at T = 10 −3 . Note that all values of γ determined from the MD configurations exceed the limit of validity of HET, γ c = 0.226.

HET fitting
Fully local approach Affine components Frozen matrix approach elasticity on length scales of order w, which, therefore, must be absorbed in the values of G 0 and γ. This is in contrast to the fully local approach, which includes all deformations on length scales w, allowing one to dissipate more efficiently the non-affine motions comprised in the CG domain, and substantially decrease the shear modulus. As a consequence, it turns out that the frozen matrix procedure is the only one correctly grasping the mechanisms accounted for by the HET, providing a very similar value for G 0 .
We conclude with an observation on our data for γ(T) of figure 3 (d), also reported in table 1. These MD values exceed in all cases the limit of validity of the HET, γ c = 0.226. This implies that if one simply inserts them in the HET, the theory would become unstable and fail badly. This limitation might originate from over-simplified assumptions, as the absence of spatial fluctuations of the local bulk modulus K m , and could possibly be solved. Already at this level, however, the theory clearly provides a valuable insight, as demonstrated in the discussion above.

Conclusions
In this work we have elucidated the interplay of anharmonic couplings and the heterogeneous mechanical response at the nano-scale in determining transverse sound waves propagation in glasses. By modifying temperature, we have been able to control the relative strengths of these mechanisms. We have provided clear evidence that the sound softening encoded in the phase velocity c(Ω, T), which at very-low T is completely determined by the disorder, is substantially decreased and eventually completely suppressed by anharmonicities. Even more intriguing, we have provided a complete characterization of the frequency dependence of the sound broadening, Γ(Ω, T), analyzing in depth the evolution of the sometimes elusive Rayleigh-like scattering, Γ ∼ Ω 4 , in the intermediate Ω-regime, and the disorder-controlled channel, Γ ∝ Ω 2 , in the high-Ω region. We have also provided a complete description of the anharmonic channel, highlighting a fractal-like frequency scaling Γ ∝ Ω 3/2 at low frequencies.