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

TL;DR

This study uses molecular dynamics simulations to explore how elastic heterogeneity influences vibrational wave propagation in glasses at finite temperatures, revealing complex transport behaviors.

Contribution

It provides a detailed analysis of vibrational excitations in glasses, linking elastic heterogeneity to transport properties through simulation and theoretical comparison.

Findings

Elastic heterogeneity significantly affects sound dispersion and attenuation.

Vibrational transport properties are strongly temperature-dependent.

Simulation results align with heterogeneous elastic theory predictions.

Abstract

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.

Tables

Table 1. Table 1: Values of the HET parameters G 0 subscript 𝐺 0 G_{0} , K 0 subscript 𝐾 0 K_{0} , and γ 𝛾 \gamma , 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 γ 𝛾 \gamma are estimated at T = 10 − 3 𝑇 superscript 10 3 T=10^{-3} . Note that all values of γ 𝛾 \gamma determined from the MD configurations exceed the limit of validity of HET, γ c = 0.226 subscript 𝛾 c 0.226 \gamma_{\text{c}}=0.226 .

	HET fitting	Fully local approach	Affine components	Frozen matrix approach
$G_0$	$25.0$	$14.0$	$35.6$	$26.3$
$K_0$	$60.0$	$59.7$	$60.2$	$60.0$
$\gamma$	$0.221$ (${\gamma }_{\text{c}}=0.226$)	$3.49$	$0.388$	$0.820$

Equations

$$S_{\text{T}}(q,\omega)=\frac{1}{2\piup N}\left(\frac{q}{\omega}\right)^{2}\int\mathrm{d}t\big{\langle}\mathbf{j}_{\text{T}}(q,t)\cdot\mathbf{j}^{\dagger}_{\text{T}}(q,0)\big{\rangle}e^{\mathrm{i}\omega t},$$

$$I_{\text{B}}(q,\omega )\propto \frac{\Gamma (q){\Omega }^2(q)}{{\left[{\omega }^2-{\Omega }^2(q)\right]}^2+{\omega }^2{\Gamma }^2(q)}.$$

$$C_{G^m}(r)=\frac{{⟨\left[G^m(\mathbf{r})-G\right]\left[G^m(0)-G\right]⟩}_0}{{⟨{\left[G^m(0)-G\right]}^2⟩}_0}.$$

$$C_{\delta \tilde{u}}(r)=\frac{⟨\delta {\tilde{\mathbf{u}}}_i({\mathbf{r}}_i)\cdot \delta {\tilde{\mathbf{u}}}_j({\mathbf{r}}_j)⟩}{⟨\delta {\tilde{\mathbf{u}}}_i({\mathbf{r}}_i)\cdot \delta {\tilde{\mathbf{u}}}_i({\mathbf{r}}_i)⟩},$$

$$\hat{g}(\omega )=\frac{3}{{\omega }_{\text{D}}^3}+\frac{4}{\text{piup}{q}_{\text{D}}^2{c}^2(\omega )}\left[\frac{\Gamma (\omega )}{{\omega }^2+{\Gamma }^2(\omega )}\right].$$

$$\Sigma (\omega ,T)={\Sigma }_{\text{harm}}(\omega ,T)+{\Sigma }_{\text{anh}}(\omega ,T),$$

$${\Sigma }_{\text{harm}}(\omega ,T)=\gamma (T){\left(\frac{G_0}{\rho }\right)}^2\frac{1}{N}\sum _{|\mathbf{q}|<q_{\text{D}}}\left[{\chi }_{\text{T}}(q,\omega )+\frac{2}{3}{\chi }_{\text{L}}(q,\omega )\right],$$

$${\chi }_{\text{T}}(q,\omega )=\frac{q^2}{-{\omega }^2+q^2\left[\frac{G_0}{\rho }-\Sigma (\omega ,T)\right]},{\chi }_{\text{L}}(q,\omega )=\frac{q^2}{-{\omega }^2+q^2\left[\frac{K_0}{\rho }+\frac{4}{3}\left(\frac{G_0}{\rho }-\Sigma (\omega ,T)\right)\right]},$$

$$c(\omega )=\sqrt{\frac{G_0}{\rho }-\text{Re}\left[\Sigma (\omega ,T)\right]},\Gamma (\omega )=\frac{\omega \text{Im}\left[\Sigma (\omega ,T)\right]}{c^2(\omega )},$$

$$g(\omega )=\frac{2\omega }{3\text{piup}}\text{Im}\left\{\frac{1}{N}\sum _{|\mathbf{q}|<q_{\text{D}}}\frac{1}{q^2}\left[2{\chi }_{\text{T}}(q,\omega )+{\chi }_{\text{L}}(q,\omega )\right]\right\},$$

Full text

Key words: quasiparticles and collective excitations, amorphous materials, molecular dynamics

PACS: 63.50.+x, 43.20.+g

• åëàñòèчíó íåîäíîðiäíiñòü. Áóëî ïîêàçàíî, ùî ïîøèðåííÿ çâóêîâèõ õâèëü çàçíàє ñèëüíîãî âïëèâó âiä íåîäíîðiäíèõ ìåõàíiчíèõ îñîáëèâîñòåé ìàòåðiàëiâ, ÿêi äîäàþòüñÿ äî ñòàíäàðòíèõ àíãàðìîíiчíèõ âçàєìîäié, ðåçóëüòóþчè â íàäçâèчàéíî ñêëàäíèõ ÿâèùàõ ïåðåíîñó ïðè ñêiíчåíèõ òåìïåðàòóðàõ. Òóò ìè äîñëiäæóєìî öi ïðîáëåìè, ìàþчè âèïàäîê ïðîñòî¿ àòîìàðíî¿ ñêëîôîðìóþчî¿ ñèñòåìè, êîìï’þòåðíèì ìîäåëþâàííÿì ìåòîäîì ìîëåêóëÿðíî¿ äèíàìiêè. Çîêðåìà, ìè ôîêóñóєìñÿ íà ïîïåðåчíèõ êîìïîíåíòàõ êîëèâíèõ çáóäæåíü чåðåç äèíàìiчíi ñòðóêòóðíi ôàêòîðè, i õàðàêòåðèçóєìî òåìïåðàòóðíó çàëåæíiñòü äèñïåðñi¿ çâóêó òà çàãàñàííÿ â øèðîêié îáëàñòi чàñòîò. Ìè ïîäàєìî ïîâíó êàðòèíó òîãî, ÿê åëàñòèчíà íåîäíîðiäíiñòü âèçíàчàє ïåðåíîñ êîëèâíèõ çáóäæåíü, òàêîæ íà îñíîâi ïðÿìîãî ïîðiâíÿííÿ чèñåëüíèõ ðåçóëüòàòiâ ç ïåðåäáàчåííÿì òåîði¿ íåîäíîðiäíî¿ åëàñòèчíîñòi.

Ключовi слова: êâàçiчàñòèíêè i êîëåêòèâíi çáóäæåííÿ, àìîðôíi ìàòåðiàëè, ìîëåêóëÿðíà äèíàìiêà

1 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 $\propto\Omega^{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/\tau=\Gamma$, from a Rayleigh-like scattering, $\Gamma\propto\Omega^{4}$, to a disorder-induced broadening, $\Gamma\propto\Omega^{2}$, moving from low to high frequency. We also found that the crossover frequency, $\Omega_{\text{co}}$, is not far from the Ioffe-Regel frequency, $\Omega_{\text{IR}}$, implying that sound waves start to loose their plane wave character at $\Omega_{\text{co}}$. In addition, in the same $\Omega$-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 observations still presents problematic aspects. One can, for instance, reasonably expect that vibrations interact with a spatially varying field associated to the disorder, undergoing substantial damping which adds to that associated to anharmonic couplings. On the one hand, this poses the problem of establishing the relative impact of the two mechanisms on the departure from the macroscopic limit, on approaching the nanoscale. On the other hand, it forces to identify the true nature of the interacting field.

Recently, evidences have accumulated (including some of our own works) which point to a mechanical origin of sound damping, in the form of substantial elastic heterogeneities [22]. Simulation [23, 24, 25] and experimental [26, 27] works have indeed demonstrated that disordered solids exhibit inhomogeneous mechanical response at the nano-scale, i.e., elastic moduli do not simply assume the hydrodynamic (continuum) values but rather they fluctuate around the hydrodynamic values, with a non-negligible distribution width. These elastic heterogeneities generate in turn non-affine deformations [28, 29], which add to the applied affine deformation field inducing a significant reduction in elastic moduli [30, 31]. Following the non-affine deformation, particles turn out to be displaced in a correlated manner, characterized by a typical mesoscopic correlation length-scale [32, 33]. It is natural to expect that interaction with this non-affine displacement amounts to profound modifications of the sound waves transport. In recent simulation works [34, 35, 36], we have provided strong evidences of this direct correlation between sound waves features and the heterogeneous mechanical properties. Based on these ideas, W. Schirmacher and co-workers [37, 38, 39] have developed a heterogeneous elasticity theory (HET), capable of reproducing numerous of the above features, including sound softening, Rayleigh-like scattering, and the Boson peak.

Differently to recent insightful theoretical works [40, 41, 42, 43, 44, 45, 46, 47, 48] which have focused on disordered systems at zero temperature, therefore disregarding the effect of anharmonicities, here we investigate the sound properties of a standard atomic glass, at variable, finite temperatures. In particular, we have determined the transverse component of the dynamic structure factor, by simulating extremely large glassy samples, and provided sound dispersions and attenuation rates in an extended frequency range. Building on our previous work, in particular [49], and adopting the heterogeneous elastic theory as a guide-line for the discussion, we highlight clearly and in a very direct way, the temperature dependent correlation between the sound transport features and the elastic heterogeneity modifications.

2 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 $\sigma$, interacting via a pairwise Lennard-Jones (LJ) potential, $V(r)=4\epsilon[(\sigma/r)^{12}-(\sigma/r)^{6}]$, with $r=r_{ij}$ the distance between particles $i$ and $j$. The potential $V(r)$ is cut-off and shifted at $r_{\text{c}}=2.5\sigma$. 111The discontinuity of the derivative $V^{\prime}(r)$ at $r_{\text{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. $m$, $\sigma$, $\epsilon$, and $\tau=(m\sigma^{2}/\epsilon)^{1/2}$ are the units of mass, length, energy, and time, respectively, while temperature is in units of $\epsilon/k_{\text{B}}$ ($k_{\text{B}}$ is the Boltzmann constant). If we consider the case of argon as a reference, temperature is in units of $\epsilon/k_{\text{B}}=125.2$ K, while length and time scales are expressed in units of $\sigma=3.405$ Å  and $\tau=2.11$ ps. We have used periodic boundary conditions in cubic simulation boxes of size $L$, such that the number density is $\hat{\rho}=N/L^{3}=1.015$. At $\hat{\rho}$, the melting and glass-transition temperatures are $T_{\text{m}}\simeq 1.0$ and $T_{\text{g}}\simeq 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 ($NVE$) ensemble at the temperature $T=2.0$ in the liquid phase, followed by a fast quench with a rate $\mathrm{d}T/\mathrm{d}t\approx 400$ down to $T=10^{-3}$, well below the glass-transition temperature $T_{\text{g}}$. Next, the systems were heated to $T=10^{-2}$, $10^{-1}$, $2\times 10^{-1}$, still below $T_{\text{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}\simeq 3$. Following thermalization, we performed the production runs for a ($N$-dependent) total time sufficient to obtain the desired $\omega$-resolution, always well below the smallest calculated line width. We have used the velocity Verlet algorithm for the integration of the equations of motion, with a time step $\delta t=5\times 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 $\mathbf{q}$ and frequency $\omega$,

[TABLE]

with the transverse current vector $\mathbf{j}_{\text{T}}(q,t)=\sum_{i=1}^{N}\left\{\mathbf{v}_{i}(t)-\left[\mathbf{v}_{i}(t)\cdot\widehat{\mathbf{q}}\right]\widehat{\mathbf{q}}\right\}\exp\{\mathrm{i}\mathbf{q}\cdot\mathbf{r}_{i}(t)\}$. Here, $q=|\mathbf{q}|$, $\hat{\mathbf{q}}=\mathbf{q}/q$, and $\langle\rangle$ is the thermodynamic average. Note that the simulation box with $N=1000188$ accommodates vibrational modes with $q=2\piup/L\simeq 0.06$. Since the first diffraction peak in the static structure factor, $S(q)$, is at $q_{m}\simeq 7$ (corresponding to an average nearest-neighbors distance $\simeq 2\piup/q_{m}=0.9$), some of our spectra refer to $q$-values $\approx 10^{2}$ times smaller than the border of the pseudo-Brillouin zone at $q\simeq q_{m}/2$.

The vibrational density of states. In the discussion we will also consider the vibrational density of state (vDOS), $g(\omega)$ [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(\omega)$ was then calculated by populating an histogram with the square-root of the obtained (positive) eigenvalues. Similar to $S_{\text{T}}(q,\omega)$, we have considered system sizes up to $N=256000$, in order to adequately sample the $g(\omega)$ on all the relevant $\omega$-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\gg 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\simeq 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], where $G^{m}$ is expressed as $G^{m}=G_{\text{A}}^{m}-G_{\text{N}}^{m}$. Here, $G_{\text{A}}^{m}$ encodes the elastic response to affine deformations, where the particles follow the applied affine strain field. By contrast, the non-affine modulus, $G_{\text{N}}^{m}$, 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_{\text{c}}$. Indeed, in [51] we demonstrated that the local modulus values are modified by non-linearities associated to discontinuities in the forces at $r_{\text{c}}$.

3 Results

3.1 Characterizing the transverse dynamical structure factor

The transverse dynamical structure factor. In figure 1, we show a selected set of data for $S_{\text{T}}(q,\omega)$, calculated from equation (2.1) for $N=1000188$, at the lowest wave number, $q=0.06$. The indicated values of $T$, all below $T_{\text{g}}\simeq 0.4$, increase from top to bottom. The $S_{\text{T}}(q,\omega)$ spectra are characterized by two symmetric Brillouin peaks, in addition to the elastic feature at $\omega=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],

[TABLE]

Here, the parameters $\Omega(q)$ and $\Gamma(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 $\Omega(q)$ we can obtain the transverse sound phase velocity, $c(q)=\Omega(q)/q$, shown in figure 2 (a)–(d) (symbols) as a function of the corresponding $\Omega$, 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 $\Omega>1$. Next, as already observed in [17], the velocity shows a minimum before reaching a region of softening, where it increases by decreasing $\Omega$. Eventually, the sound velocity saturates at the correct acoustic value, $c(q\rightarrow 0)\rightarrow v_{T}=\sqrt{G/\rho}\simeq 3.65$ (dashed line), where $G$ is the shear modulus, and $\rho=m\hat{\rho}$ 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\times 10^{-1}\simeq T_{\text{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_{\text{L}}^{2}(\omega)\simeq K/\rho+(4/3)c_{\text{T}}^{2}(\omega)$, where we have discarded the frequency dependence of $K(\omega)$, with the replacement of the zero-frequency limit, $K=K(\omega=0)\simeq 59$. This implies that the knowledge of one of the two velocities allows one to infer the $\omega$-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 [60, 45]. 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 $\Gamma$ as extracted from equation (3.1), as a function of the corresponding $\Omega$. 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 $\propto\Omega^{2}$ [61], and a Rayleigh-like scattering contribution, $\propto\Omega^{4}$, at lower frequencies [62]. As already noticed, the crossover frequency $\Omega_{\text{co}}\simeq 1$, is not far from the calculated BP frequency, $\omega_{\text{BP}}\simeq 2$, which is compatible with the Ioffe-Regel frequency, $\Omega_{\text{IR}}\simeq\omega_{\text{BP}}$, as also discussed in [17] 222In [45] we have shown that the crossover occurs at the frequency $\omega_{0}$ where the vDOS converges to the Debye limit. For the present LJ system, $\omega_{0}\simeq 1$ [17], confirming $\Omega_{\text{co}}\simeq\omega_{0}$.. 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 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 $\Omega_{\text{co2}}\simeq 0.6$, between the Rayleigh region and an unexpected low-frequency $\propto\Omega^{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 ($\Omega_{\text{co2}}\simeq\Omega_{\text{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 $\Omega$-range. This is indeed what we observe at $T=10^{-1}$ in (g), where the Rayleigh contribution is taken over by the fractal $\sim\Omega^{3/2}$ behaviour at low frequency. In addition, a quadratic $T$-independent contribution is still visible at high $\Omega$. Finally, at the highest $T=2\times 10^{-1}\simeq T_{\text{g}}/2$ in panel (h), the intensity is fully $T$-dependent in the entire $\Omega$-range, with an Akhiezer-like $\Omega^{2}$ scaling prevailing in a large frequency range, while a remnant $\Omega^{3/2}$ behaviour can be only recovered at low $\Omega$. 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 $\Omega^{3/2}$ fractional scaling described above and predicted in the context of the HET [64, 63], 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.

3.2 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}_{\text{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 $\gamma(T)=\hat{\rho}w^{3}\delta G^{2}/G^{2}$, involving the ratio of $\delta 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, $\gamma_{\text{A}}(T)=\hat{\rho}w^{3}\delta G_{\text{A}}^{2}/G_{\text{A}}^{2}$ (blue squares), with again a $T$-dependence analogous to that of $\gamma(T)$. Note that $\gamma(T)$ ($\gamma_{\text{A}}(T)$) plays a crucial role in the HET [37, 38, 39, 63, 64], where the parameter $\delta\gamma=\gamma_{\text{c}}-\gamma$ (inset of figure 3 (d), with $\gamma_{\text{c}}=0.226$) controls the approach to the mechanical instability at $\gamma_{\text{c}}$ and induces the fractal frequency dependence of broadening, $\Gamma\propto\Omega^{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 $\Omega=\Omega_{\text{co}}\simeq 1$, where the scattering mechanism is modified from $\Gamma\propto\Omega^{4}$ to $\propto\Omega^{2}$. As we already mentioned, $\Omega_{\text{co}}\simeq 1$ is not far from the Ioffe-Regel frequency $\Omega_{\text{IR}}\simeq 2$ [17]. As a consequence, at $\Omega_{\text{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 $\xi_{\text{co}}$ (or a wave number $q_{\text{co}}=2\piup/\xi_{\text{co}}$), with $\xi_{\text{co}}=2\piup v_{T}/\Omega_{\text{co}}\simeq 23$ ($q_{\text{co}}=\Omega_{\text{co}}/v_{T}\simeq 0.27$) [32]. We stress that both $\Omega_{\text{co}}$ and the length-scale $\xi_{\text{co}}$ are $T$-independent and related to the $T=0$ structural properties only. Here, we attempt to characterize $\xi_{\text{co}}$ in terms of the local shear modulus distributions.

In figure 4 (a), we show the spatial correlation function of $G^{m}$, defined as

[TABLE]

Here, we have explicitly represented $G^{m}(\mathbf{r})$ as a function of the position $\mathbf{r}$, and the terms $[G^{m}(\mathbf{r})-G]$ are the fluctuations around the average (macroscopic) value, $G$. $\left<\right>_{\mathbf{0}}$ denotes the spatial average over $\mathbf{r}=\mathbf{0}$ and $r=\left|\mathbf{r}\right|$. 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 333More precisely, fluctuations of the affine elastic moduli, which are also randomly distributed in space [57], drive the non-affine displacement field., 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 $\xi_{\text{na}}$. In [32], it has been estimated that $\xi_{\text{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 $\xi_{\text{na}}\simeq\xi_{\text{co}}$, extracted from the total transverse sound broadening, a point highlighted in the experimental work of [69].

To verify this point, we have extracted $\xi_{\text{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 ($\gamma=5\times 10^{-4}$), and next relax the system to the closest local potential energy minimum. Throughout the relaxation process, we keep track of the total non-affine displacement, $\delta{\mathbf{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, $\delta\tilde{\mathbf{u}}_{i}=\delta{\mathbf{u}}_{i}-\sum_{k;\ \omega_{k}<\omega_{0}}\left[\delta\mathbf{u}_{i}\cdot\mathbf{e}^{k}_{i}\right]\mathbf{e}^{k}_{i}$, with $\omega_{k}$ and $\mathbf{e}^{k}_{i}$ the eigenfrequencies and eigenvectors extracted by a normal mode analysis of the Hessian matrix [17]. Note that we set $\omega_{0}=1$ where the $T=0$ vDOS converges to the Debye limit [17].

We represent the $\delta\tilde{\mathbf{u}}_{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_{\delta\tilde{u}}(r)$,

[TABLE]

with $r=\left|\mathbf{r}_{i}-\mathbf{r}_{j}\right|$, and $\left<\right>$ denotes the average over all pairs of particles, $i$ and $j$. [Note that our $C_{\delta\tilde{u}}(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_{\delta\tilde{u}}(r)$ decays on length scales $r=\xi_{\text{na}}\simeq 20$, consistently with [32] and very close to the $\xi_{\text{co}}\simeq 23$ associated to $\Omega_{\text{co}}$ as discussed above.

This observation ultimately allows us to rationalize the value of the crossover frequency $\Omega_{\text{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 $\lambda<\xi_{\text{co}}$, while for $\lambda>\xi_{\text{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 ($\Omega<\Omega_{\text{co}}$) and high ($\Omega>\Omega_{\text{co}}$) frequencies, respectively; iv) These distinct effects induce the observed distinct scattering scalings, $\Gamma\propto\Omega^{4}$ and $\propto\Omega^{2}$, together with the associated cross-over at $\Omega_{\text{co}}$.

3.3 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 $\hat{g}(\omega)=g(\omega)/\omega^{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 $\Gamma$ to the BP features,

[TABLE]

Here, $q_{\text{D}}=(6\piup^{2}\hat{\rho})^{1/3}\simeq 3.92$, $c_{\text{D}}=\left[\left(v_{\text{L}}^{-3}+2v_{\text{T}}^{-3}\right)/3\right]^{-1/3}\simeq 4.13$, and $\omega_{\text{D}}=q_{\text{D}}c_{\text{D}}\simeq 16.19$ are the the Debye wave-number, velocity and frequency, respectively, with $v_{\text{L}}=\sqrt{(K+4G/3)/\rho}\simeq 8.71$ and $v_{\text{T}}=\sqrt{G/\rho}\simeq 3.65$ the longitudinal and transverse sound speeds. The term $3/\omega^{3}_{\text{D}}$ is the Debye limit. Hence, $\hat{g}(\omega)$ can be determined by equation (3.4) if the sound speed, $c(\omega)$, and broadening, $\Gamma(\omega)$, 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 $\hat{g}(\omega)$ 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 $\omega_{\text{BP}}\simeq 2$. The situation is similar at $T=10^{-2}$ in (b), although anharmonicities already start to alter the small-$\omega$ behaviour, a modification which is completed at the two highest temperatures, $T=10^{-1}$ (c) and $2\times 10^{-1}$ (d). Now, the data decrease by increasing frequency, following an $\omega^{-1/2}$ power-law, before directly matching the BP intensity. The simultaneous presence of the $\omega^{-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:

[TABLE]

where $\chi_{\text{T,L}}(q,\omega)$ 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 $\Sigma_{\text{anh}}(\omega,T)\equiv\text{i}\nu\omega T$ [70, 71]. Equations (3.5) contain four free parameters: the shear ($G_{0}$) and bulk ($K_{0}$) moduli, the disorder parameter $\gamma=\hat{\rho}w^{3}\delta G_{0}^{2}/G_{0}^{2}$, and the anharmonic coupling $\nu$. Since $\Sigma$ is directly related to sound velocity, broadening and vDOS [37, 38, 39], as

[TABLE]

(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). We have extracted the optimized ($T$-independent) values $G_{0}=25$, $K_{0}/G_{0}=2.4$, and $\nu\rho/G_{0}=3.0$ 444For the case of $T=10^{-3}$, we set $\nu=0$ because we did not recognize any anharmonic contribution to the simulation data of figure 2., while the variable strength of the elastic heterogeneity is encoded in the distance from the elastic instability, $\delta\gamma(T)=\gamma_{\text{c}}-\gamma(T)$, with $\gamma_{\text{c}}=0.226$ [inset of figure 3 (d)].

In figures 2 and 5, we plot by solid lines the predictions for $c(\Omega)$, $\Gamma(\Omega)$, and $\hat{g}(\omega)$ 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(\Omega)$, the crossover in $\Gamma$ from the $\propto\Omega^{4}$ to the $\propto\Omega^{2}$ scaling, and the BP intensity and position in $\hat{g}(\omega)$. While in general the intensity of the broadening is considerably underestimated, the theory is capable of pin-pointing the crossover frequency, $\Omega_{\text{co}}\approx 1$, and thus the length scale, $\xi_{\text{co}}\approx 20$. Note that these values are directly determined by the parameter $\gamma$, 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(\Omega)$, the increase of anharmonic damping with temperature, and the fractal scaling laws $\propto\Omega^{3/2}$ in $\Gamma$ and $\propto\omega^{-1/2}$ in $\hat{g}(\omega)$ at the highest $T=2\times 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 $\gamma(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\simeq 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_{\text{A}}$ and $K_{\text{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 elasticity on length scales of order $w$, which, therefore, must be absorbed in the values of $G_{0}$ and $\gamma$. 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 $\gamma(T)$ of figure 3 (d), also reported in table 1. These MD values exceed in all cases the limit of validity of the HET, $\gamma_{\text{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.

4 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(\Omega,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, $\Gamma(\Omega,T)$, analyzing in depth the evolution of the sometimes elusive Rayleigh-like scattering, $\Gamma\sim\Omega^{4}$, in the intermediate $\Omega$-regime, and the disorder-controlled channel, $\Gamma\propto\Omega^{2}$, in the high-$\Omega$ region. We have also provided a complete description of the anharmonic channel, highlighting a fractal-like frequency scaling $\Gamma\propto\Omega^{3/2}$ at low frequencies.

In addition, we have demonstrated that the heterogeneous elasticity theory developed by W. Schirmacher and co-workers [37, 38, 39, 64, 70, 71], is capable of acceptably reproducing the simulation data when fed with realistic nanoscale mechanics information. At the lowest $T$, the theory reproduces the sound softening as well as the crossover from $\Gamma\propto\Omega^{4}$ to $\propto\Omega^{2}$, both of which are disordered effects. The crossover frequency $\Omega_{\text{co}}\simeq 1$ and length $\xi_{\text{co}}=2\piup v_{\text{T}}/\Omega_{\text{co}}\simeq 23$ are intimately related to the length $\xi_{\text{na}}$ associated to the non-affine displacement field, which is indeed induced by randomly fluctuating local elastic moduli. At the higher $T$, where the anharmonicities couple with the elastic heterogeneities, HET indeed reproduces a complicated sound propagation behavior, including disappearance of the sound softening, an increase of anharmonic damping with temperature, and the fractal-like scaling laws in $\Gamma(\Omega)$ and $\hat{g}(\omega)$. Remarkably, the theory predicts that the fractal scaling is related to the approach to the instability $\gamma_{\text{c}}$ [64, 63], a picture reasonably supported by our simulation data. Our results are an additional evidence that only developments which integrate local mechanical features of materials and a full treatment of the anharmonic couplings will be capable of providing the complete picture for sound waves propagation in disordered solids.

Acknowledgements

We thank Giancarlo Ruocco for useful discussions, and Walter Schirmacher for feedback in the early stage of this work. H. M. is supported by JSPS KAKENHI Grant Numbers 19K14670 and the Asahi Glass Foundation. S. M. is supported by ANR-18-CE30-0019 (HEATFLOW).