Vibrational density of states of amorphous solids with long-ranged power-law correlated disorder in elasticity

A theory of vibrational excitations based on power-law spatial correlations in the elastic constants (or equivalently in the internal stress) is derived, in order to determine the vibrational density of states $D(\omega)$ of disordered solids. The results provide the first prediction of a boson peak in amorphous materials where spatial correlations in the internal stresses (or elastic constants) are of power-law form, as is often the case in experimental systems, leading to logarithmic enhancement of (Rayleigh) phonon attenuation. A logarithmic correction of the form $\sim -\omega^{2}\ln\omega$ is predicted to occur in the plot of the reduced excess DOS for frequencies around the boson peak in 3D. Moreover, the theory provides scaling laws of the density of states in the low-frequency region, including a $\sim\omega^{4}$ regime in 3D, and provides information about how the boson peak intensity depends on the strength of power-law decay of fluctuations in elastic constants or internal stress. Analytical expressions are also derived for the dynamic structure factor for longitudinal excitations, which include a logarithmic correction factor, and numerical calculations are presented supporting the assumptions used in the theory.


I. INTRODUCTION
Understanding the physics of vibrational spectra of disordered systems is a classical topic in condensed matter physics [1][2][3][4]. Glasses and other disordered solids exhibit anomalous features, compared with their crystalline counterparts. Concerning the thermal properties, at few tens of Kelvin, the specific heat of glasses exhibits an excess over the Debye prediction, in the form of a characteristic maximum in the plot of C(T )/T 3 . The peak is ascribed to the presence of an excess of states over the Debye density of states (DOS) ∼ ω 2 , known as the boson peak since its temperature dependence conforms with that of the Bose function, and thus appears to strongly depend on the features of the vibrational modes in the THz frequency [5][6][7][8].
Among previous theories, the heterogeneous elasticity theory (HET) [42,[55][56][57][58] uses a field-theoretical scheme to derive the DOS, by assuming Gaussian uncorrelated * az302@cam.ac.uk spatial fluctuations in the elastic constants of the system [4]. The theory provides a quantitative relation between the boson peak and the Brillouin width (sound attenuation coefficient) Γ, and reproduces the Rayleigh scattering law Γ ∼ ω d+1 . However, following numerical evidence of a logarithmic enhancement correction of the form Γ(k) ∼ −k d+1 ln(k) to the Rayleigh scattering law (with wavenumber k, in d-dimension) [59], it has been shown analytically that long-ranged power-law spatial correlations in elasticity, or equivalently in the internal stresses, are the cause of such enhancement [57].
Previous attempts to derive the logarithmic Rayleigh law using HET with power-law correlations in elasticity by Caroli and Lemaitre [60] were not successful due to two major simplifying approximations used in their theory, namely the assumption of perfectly isotropic wave propagation (with completely decoupled longitudinal and transverse propagators), which leads to a cancellation of terms and to the vanishing of the logarithmic correction. Caroli and Lemaitre's oversimplifying assumption of isotropic wave propagation is at odds with numerical evidence from Ref. [59], which showed that wave propagation in the presence of power-law correlated elasticity is locally anisotropic over relatively large length-scales, leading to at least 5 non-vanishing local elastic constants. In Ref. [57], by finding the rigorous solution to the selfconsistent anistropic wave propagation problem, it was possible to derive the logarithmic Rayleigh scattering law, which is ubiquitously observed in experiments and simulations [61][62][63][64][65][66][67][68], and to show that it is the direct result of the power-law correlation in internal stresses or elastic constants.
Independent evidence supporting the existence of power-law spatially decaying correlations in elasticity have been shown in recent works [69,70]. All these facts point towards the importance of properly accounting for long-ranged power-law elastic correlations in the description of the vibrational properties of disordered systems.
A fundamental unanswered question, therefore, is what impact the underlying long-ranged power-law correlations of elasticity (or internal stresses) may have on the DOS. The answer is presented in this article, where we exploit the successful framework of Ref. [57] for the acoustic attenuation, and apply it to study the properties of the DOS. We reveal that the boson peak picks up a logarithmic correction which is most evident in the excess DOS. We also show how the boson peak sensitively depends on the strength of power-law correlations of elasticity. We also examine the asymptotic scaling behavior of the DOS in the frequency regimes where modes are quasilocalized due to the disorder (hence undergoing diffusivelike propagation instead of ballistic propagation typical of standard phonons, as demonstrated for glasses in earlier works [52]).

II. THEORY FOR LONGITUDINAL EXCITATIONS IN 2D
Numerical simulations, supported by theoretical analysis, and analysis of experimental data, suggest logarithmic enhancement of the Rayleight law in a certain frequency domain [59]. In appendix, we review the model that only predicts the Rayleigh scattering law, such that the mean free path (ω) scales as ω −4 for small ω.
According to the theoretical analysis in Ref. [57], such enhancement to Rayleigh scattering of phonons in amorphous solids, originate from long-range power-law spatial correlations of elastic constants or internal stress. For example, the shear stress tensor, σ(r) = σ 0 + ∆σ(r) is expressed in terms of its mean value plus a random fluctuation, i.e. ∆σ(r) = 0 and ∆σ(r )∆σ(r + r) = B(r) = κ 2 cos(4θ)/(r 2 + ζ 2 ) for some constants κ, ζ. The parameter κ describes the strength of the disorder, while ζ controls the regime of frequency where logarithmic enhancement occurs. The resultant 2D self-consistent equations between self-energy Σ(k, z) and the Green's function G(k, z) for longitudinal waves read [57]: where the prefactor a is a new parameter reflecting the strength of elastic heterogeneity. Determining Σ(k, z) via solving self-consistent equations above, one can compute the Green's function and hence obtain the density of states (DOS), D(ω), via the standard Plemelj identity: From Ref. [57], we can approximateB(k − q) ∼ − ln(bk), which is valid upon assuming the linear (acoustic) dispersion relation between wavenumber q and frequency ω. The parameter b depends on ζ in B(r), in a way such that the larger ζ is, the larger b turns into, thus the lower frequency regime that log-effect emerges. Substituting this into Eq. (2) and re-introducing the parameter a, we solve the self-consistent equation of self-energy Σ(z) in 2D: from which the scaled DOS is obtained as .

III. GENERAL THEORY FOR AMORPHOUS SOLIDS IN 3D
The model in the last section describes purely longitudinal waves in 2D and is served to illustrate the basic functioning of the theoretical framework on an easy example. Now we present the full theory, with the inclusion of transverse waves, in 3D. Making similar assumptions of elastic disorder as in Refs. [55,56], we consider an elastic medium with a mass density m 0 , shear modulus G, bulk modulus K = λ + 2G/3 where λ is the longitudinal Lamé's constant. The elastic constants are related to the longitudinal and transverse local sound velocities as c T ≡ G/m 0 , c 2 L ≡ (K + 4G/3)/m 0 = (λ + 2G)/m 0 , respectively. The Lamé's constant is set to be λ = λ 0 , while the shear modulus includes a random spatial variation G(r) = G 0 [1 + ∆G(r)]. The random function ∆G(r) is supposed to have a long-ranged power-law decay ∆G(r )∆G(r + r) = B (r) ∝ γ 2 /(r 2 + ξ 2 ) 3/2 for some constants γ and ξ, in agreement with recent evidence for glasses and granular materials [69,70]. The explicit form of angular component in B(r) is not relevant to results and is not shown in the last expression. The self-consistent Born approximation for the complex selfenergy Σ(ω), based on the standard replica-trick, leads to the set of self-consistent equations: .
We again assume the linear dispersion relation between k and ω, which is verified by dynamical structure factor calculations in Appendix A. Likewise, the new parameter g, which absorbs the disorder-strength parameter γ, is the prefactor of the self-consistent equation for Σ(k, ω). The Debye length and frequency are given by k −1 [56]. The DOS can be calculated as In Fig. 3, we show the typical reduced DOS D(ω)/ω 2 , i.e. the usual boson peak representation, as well as the reduced excess DOS, D(ω)/ω 2 −1, against the ω/ω D . It can be seen from Fig. 3 (a) that the boson peak frequency decreases sharply upon increasing g, hence upon increasing the degree of elastic disorder, γ. Hence, larger disorder lifts up the boson peak and shifts it to lower frequencies, in accordance with earlier findings from simulations [72]. We also note from Fig. 3(a) that the boson peak drops exponentially with its frequency ω BP upon increasing g, which is a new law found here by our theory. This might be related to the exponential decaying mode in spectra of activation energies in metallic glasses [73]. In Fig. 3(b), obviously, the excess over the Debye level is different form zero only above a certain frequency threshold. The excess DOS turns out to vanish as ω 4 for ω → 0. Upon approaching the boson peak frequency, i.e. ω ω BP , where the DOS D(ω) displays the ω 2 dependence, the excess DOS tends to flatten out. The nature of vibrational eigenmodes varies as ω changes. In particular, the asymptotic behaviour ∼ ω 4 as ω → 0, is consistent with results of previous work using HET and Gaussian disorder in elastic constants in [55], and with numerical evidence in [68,[74][75][76][77]. Remarkably, when the frequency becomes comparable to the boson peak frequency, an additional trend representing a logarithmic correction dependency is observed, as shown in Fig. 3(b). This is the first prediction of the logarithmic correction in the reduced DOS, which is expected based on the generic direct proportionality relation between the excess DOS and the phonon attenuation coefficient highlighted in Refs. [56,68].
We also show how the boson peak changes with different values of b in Fig. 4(a), where a similar monotonic relation as in the purely longitudinal case is observed. In Fig. 4(b), clearly the excess D(ω) − ω 2 is reduced upon increasing b because of the interplay between the prominent ∼ ω 4 behavior at lower ω and the influence of logarithmic enhancement at higher ω.

IV. CONCLUSIONS
In summary, we developed a theory of vibrational excitations in disordered media with long-ranged powerlaw correlated disorder to extract the density of states (DOS). The assumption of power-law correlated disorder in elastic properties (elastic constants or stresses), supported by evidence found in glasses [69] and granular materials [70], has been key to derive the logarithmic enhancement ∼ −k d+1 ln k of Rayleigh scattering in glasses in our previous work [57], by accounting for the anisotropic character of wave propagation in the solid locally. The theory reproduces the boson peak in the DOS along with its dependence on the strength γ (or g), and on the characteristic scale b, of power-law correlated disorder. Importantly, the theory predicts a logarithmic correction ∼ −ω 2 ln ω visible in the reduced excess DOS around the boson peak frequency, predicted here for the first time. The theory also predicts that the boson peak decays exponentially with its frequency ω BP upon increasing the strength of disorder g, which might be related to evidence recently found in metallic glasses [73]. The theory predicts the existence of a ∼ ω 4 regime (in 3D) at low frequency below the boson peak, which can be ascribed to Rayleigh scattering. Similar ω 4 modes have been recently discovered in the non-Debye part of the spectrum, which may be ascribed to localized anharmonic modes [78,79], which have been demonstrated to be universal in recent work [80,81]. It appears that the present theory, which is rather on the continuum level and does not account for anharmonicity, cannot predict those modes, while the predicted ω 4 refers most likely to Rayleigh scattering since we have checked that for a 2D system the scaling is much closer to ω 3 . The logarithmic feature in the reduced DOS predicted by this theory calls for more detailed investigation of experimental data in future analysis. Studying the influences of anharmonicity [71,82], nonaffine elasticity [83][84][85] as well as glass stability [67] on the vibrational excitation modes within the current theoretical model, will be the object of future The self-consistent Born approximation, using the replica trick to evaluate the Green's function of an elastic Lagrangian with quenched Gaussian disorder in the elastic constant, was proposed in Ref. [42]. This leads to a self-consistent relation between the (complex) selfenergy Σ(z) and the 3D Green's function G(z) of the longitudinal waves as: where parameters c 0 , γ correspond to mean of the sound velocity and to the variance of the elastic autocorrelations, respectively. The Debye wavenumber is k D = (6π 2 N/V ) 1/3 for a system with N particles and volume V , so that the frequency z = ω 2 + i0 has a (Debye) cutoff value at ω D = c 0 k D . With standard identification k 2 dk, we can transform the discrete sums into continuous integrals over momentum space: The integral in Eq. (A2) can be calculated analytically, giving Setting k D , c 0 , γ and a proper initial value of self-energy, Σ 0 (z), we can use an iteration scheme to numerically determine Σ(z). Figure 5 shows a typical plot of DOS calculated in this way, where k D = q D = 1 and c 0 = 0.5. For convenience, we let the prefactor (3N γ/2k 3 D ) on the RHS in Eq. (A2) be d.

Appendix B: Dynamical structure factor
According to [57,86], the 3D longitudinal dynamical structure factor S L (k, ω) has the following expression where n(ω) + 1 = [1 − exp(− ω/k B T )] −1 is the Bose factor, c L (ω) is the (generalized) longitudinal sound speed. In the classical limit, ω/k B T → 0, and using the fact that ImΣ(ω) ∼ −ω 2 ln(ωb), we have Taking same parameters as in Fig. 3 in the maintext, we show the longitudinal dynamical structure factor in Fig. 6, where the linear dispersion relation is evident between the peak position and wavevector k. At low ω, Eq. (B2) can be fitted with a damped harmonic oscillator (DHO) model: where Ω(k) corresponds to the excitation frequency and Γ(k) is the width of the Brillouin line (full width at half-maximum of the excitations). This is consistent with the proportionality coefficient between peak position fre-quency and k identifying the longitudinal speed of sound (c L = 1).