Front. Phys.Frontiers in PhysicsFront. Phys.2296-424XFrontiers Media S.A.10.3389/fphy.2017.00029PhysicsOriginal ResearchMoment-Preserving Theory of Vibrational Dynamics of Topologically Disordered SystemsFolliViola1RuoccoGiancarlo12*SchirmacherWalter121Center for Life Nano Science, Fondazione Istituto Italiano di TecnologiaRome, Italy2Department of Physics, University of Rome ‘La Sapienza’Rome, Italy
Edited by: Gang Zhang, Institute of High Performance Computing (A*Star), Singapore
Reviewed by: Byungchan Han, DGIST, South Korea; Jiebin Peng, National University of Singapore, Singapore
This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
We investigate a class of simple mass-spring models for the vibrational dynamics of topologically disordered solids. The dynamical matrix of these systems corresponds to the Euclidean-Random-Matrix (ERM) scheme. We show that the self-consistent ERM approximation introduced by Ganter and Schirmacher [1] preserves the first two nontrivial moments of the level density exactly. We further establish a link between these approximations and the fluctuating-elasticity approaches. Using this correspondence we derive and solve a new, simplified mean-field theory for calculating the vibrational spectrum of disordered mass-spring models with topological disorder. We calculate and discuss the level density and the spectral moments for a model in which the force constants obey a Gaussian site-separation dependence. We find fair agreement between the results of the new theory and a numerical simulation of the model. For systems with finite size we find that the moments strongly depend on the number of sites, which poses a caveat for extrapolating finite-system simulations to the infinite-size limit.
glassesdisordered systemsvibrational dynamicsdensity of statestheorySCBAheterogeneous elasticityPacs numbers:63.50.-×61.43.Fs65.60.+a1. Introduction
The vibrational properties of disordered solids at high frequencies, in the THz range, have been subject to great attention both from the experimental [2] and from the theoretical side [3] mainly because of the anomalies observed in the specific heat and in the thermal conductivity of glasses [4, 5]. The origin of an excess of the vibrational density of states (DOS) with respect to the Debye prediction in the THz region (an excess that is often referred to as “boson peak”) is thought to be the common origin of both anomalies. Also of interest, and at the origin of a lively debate [6–10], is the anomalous sound attenuation and the presence of a non-linear sound dispersion observed in the same frequency regime.
There is, nowadays, agreement on the fact that both the boson peak and the frequency dependence of the sound attenuation are the dynamical manifestation of the underlying topological (structural) disorder and that these features are not related to anharmonic interactions, to hopping processes or to other exotic phenomena like van-Hove singularities [11, 12]. The rich phenomenology observed in the dynamics of glasses in the THz frequency region can all be considered of harmonic origin [13, 14]. What is still missing is a deep understanding of the mechanisms that correlate the specificity of the topological disorder with the details of the observed dynamical anomalies.
In this respect, it becomes evident how important is the development of models of topological disorder with tunable properties.
The model approaches present in the literature for the influence of the disorder on the harmonic vibrational spectrum of solids can be grouped into models with waves in a spatially inhomogeneous media and mass-spring models. The wave approaches can again be grouped into a discription in terms of defect models [15–21] and fluctuating-elasticity approaches [9, 22–26].
In the defect models additional irregularities like atoms with heavy masses [15–17] or quasi-local harmonic oscillators [18–21], which couple to the Debye phonons, are introduced. The defect and coupling parameters are subject to a random distribution. In the fluctuating elasticity approach [9, 22–25, 27, 28] it is assumed that the elastic constants, (in particular the shear modulus) fluctuate around a mean value at different locations in space in a spatially correlated or uncorrelated fashion. Within this latter framework the disorder can be shown to lead to the anomalous increase of the DOS (and hence specific heat) with respect to the Debye prediction and to a contribution to the sound attenuation. It has been pointed out [27, 28] that in principle the defect description can be subsumed under the fluctuating-elasticity approach via a coarse-graining procedure.
In the disordered mass-spring models (which we shall address in the present contribution) one considers a set of mass points (mass M), which are fixed in space according to a given distribution and are connected by harmonic springs with spring constants Kij. The equation of motion for the (scalar) displacements ui from the rest positions is
Md2dt2ui(t)=-∑j≠iKij[ui(t)-uj(t)]
The disorder may be introduced by assuming a distribution of the Kij on an otherwise ordered lattice (force-constant disorder) or placing the sites ri randomly in space and then making the force constant depend on the site separations rij = |ri − rj|, i.e., Kij = K(rij) (topological disorder).
There is an extended literature on such disordered mass-spring models [1, 29–42], which is due to the fact that they are considered to be relevant to the high-frequency vibrational dynamics of amorphous solids and liquids. They exhibit a wealth of anomalous phenomena including a boson peak [29, 32, 34–36, 38, 39, 41] and localization [33, 35, 42]. There is another point, which makes the investigation of models, which obey (Equation 1), interesting. If one replaces the double time derivative by a single one (and takes M = 1), one obtains the equation of motion of a random walker in a disordered system. The ui(t) then take the meaning of the probability to find the walker at site i. It has been pointed out [1, 32, 43] that all anomalous features observed in the vibrational system (Equation 1) are also observed in the random-walk system: For example Rayleigh scattering in the vibrational system corresponds to an algebraic long-time tail in the velocity autocorrelation function in the random-walk system [1]. The boson peak corresponds to the onset of a strong frequency dependence of the conductivity in the random-walk system [32, 43].
The so-called Euclidean-Random-Matrix (ERM) approach introduced by Grigera et al. [39] is a systematic analytic approach to calculate the averaged vibrational spectrum and dynamical structure factor of systems govened by Equation (1). This is achieved by a high-frequency and high-density expansion, which is organized with the help of a diagram technique [1, 39, 40]
The ERM formalism allows for formulating self-consistent effective-medium approximations for the averaged Green's function corresponding to Equation (1). Such effective-medium approximations (“self-consistent ERM (SCERM) approximations,” which include the Rayleigh-scattering property1 have been introduced by Ganter and Schirmacher [1].
In view of the fact that the ERM formalism enables to directly link the structure of a disordered solid to its vibrational spectrum it is desirable to establish a connection between this theory and the phenomenological heterogeneous elasticity theory. This, among others, will be done in the present paper. We start our investigation of the ERM mass-and spring model by showing that the self-consistent ERM approximations preserve the first two nontrivial moments (1st and 2nd) of the averaged level density. We then establish a link between the SCERM approximations and the self-consistent Born approximation (SCBA) of the heterogeneous elasticity theory. By means of this correspondence we devise a simplified SCERM theory (ERM-SCBA), which allows for solving the equations by a few iterations instead solving a set of integral equations.
In the second part of the paper we calculate numerically and in ERM-SCBA the level density and the first two nontrivial moments for a model with a Gaussian r dependence of the force constants. The moment calculations are done both for an infinite and a finite system. For the latter we find a strong N dependence of the moments, which is relevant for simulations in which only finite systems can be dealt with.
2. Formalism2.1. Mass-spring model
We consider a mass-spring model with N mass points distributed at random in a volume V = L3 = N/ρ. We shall also consider the limit N → ∞ with fixed density ρ.
Equation (1) can be re-written as
d2dt2ui=-∑jDijuj
with the dynamical matrix D, having off-diagonal elements Dij = −K(rij)/M and diagonal elements Dii=-∑i≠jDij=∑i≠jK(rij)/M. D has N eigenvalues λp=ωp2, where ωp are the eigenfrequencies.
2.2. Density of states and spectral moments
The density of states g(ω) and the density of eigenvalues (“levels”) ρ(λ) are given by
ρ(λ)=12ωg(ω)=1N〈∑pδ(λ−λp)〉=〈1NπIm{tr[s𝟙+D]−1}〉
where s = −ω2 − iϵ, 𝟙 is the unit matrix and 〈… 〉 denotes a configurational average. Again λ = ω2 denotes the eigenvalues of the dynamical matrix. Due to the sum rule ∑iDij=0 (translation invariance) the dynamical matrix has always a “trivial” eigenvalue λ = 0 and, more important, there exists always an interval 0 < λ < λ*, in which a Debye law is valid, i.e., g(ω) ∝ ω2 ↔ ρ(λ) ∝ λ1/2.
In the present paper we thoughout assume a random distribution of points with density P(r1…ri…rN)=V-N. The n-th spectral moment of the level density is defined as
M(n)=∫dλλnρ(λ).
With the present definition of ρ(λ), the zeroth moment turns out to be equal to one, therefore the moments are automatically “normalized”. Let's also define the “central” moments:
Mc(n)=∫dλ(λ-λ̄)nρ(λ)
with λ̄≐M(1).
As stated above, the dynamical matrix has always the “trivial” eigenvalue, λ = 0. This has a non-negligible effect on the higher moments (specifically, it effects Mc(n)). Therefore, we define a “non trivial” density of states, ρ^(λ), where the vanishing eigenvalue is excluded from the sum:
ρ^(λ)=1N-1∑p|λp≠0δ(λ-λp)
Let's call M^(n) and M^c(n) the corresponding moments.
As far as the “hat” moments, we note that all the moments with n > 0 are unaffected (add or subtract λpn, with λp = 0), but the normalization ((N − 1) instead of N) is different. The “hat” moments are thus obtained by multiplying the “non-hat” ones with N/(N − 1).
The moments M(n) can be directly related to the dynamical matrix via:
M(n)=∫dλλnρ(λ)=1N〈∑pλpn〉=1N〈tr{Dn}〉
We represent the moments in terms of the dimensionless function f(r) defined as K(rij) = K0f(rij/σ) = K0fij, where σ is the characteristic decay length and their average
In the N → ∞ limit it is useful to define the dimensionless constants and Γ=1ρσ3, f0=∫d3rf(r), and f0(2)=∫d3rf(r)2, where the integral is over the entire space.
In this limit the M^(2) moments become equal to the ones without hat and we have
M∞(1)=K0f0MΓM∞(2)=2K0f0(2)MΓ+[M∞(1)]2Mc,∞(2)=2K0f0(2)MΓ2.3. Euclidean random matrix (ERM) formalism
To make contact to the ERM formalism [1, 39, 40] we consider the Fourier transformed force constants
t(q)≐ρMK~(q)=K0f0MΓf~(σq)=t0f~(σq)
with
f~(k)=f~(k)=1f0∫d3reikrf(r)
so that f~(k=0)=1 and t0=t(k=0)=M∞(1).
The quantity of interest is the q dependent averaged Green's function
where Σ(q, s) is the self-energy function and, as before, s = −λ − ϵ = −ω2 − ϵ. The level density is obtained from the Green's function via
ρ(λ)=1πIm{G∞(s)}
with G∞(s)=lim|q|→∞G(q,s).
The self energy can (and has been [1, 39, 40]) calculated with increasing powers of the inverse density ρ. The lowest-order result is Σ(0)(q, s) = 0, so that in the high-density limit
G(0)(q,s)=1s+t0−t(k)
One can make contact between this high-density limit and the traditional liquid-phonon theory of Hubbard and Beeby [44] by identifying the “liquid dispersion” with Ω=t0-t(k) and the Einstein frequency with ΩE=t0.
The lowest-order nontrivial contribution to the self energy is [1, 39, 40]
Σ(1)(q,s)=∑p[t(p-q)-t(p)]2G(0)(p,s)
Here ∑p denotes 1ρ(2π)3∫d3p.
It has been suggested [39, 40] that one might obtain a suitable self-consistent (i.e., non-perturbative) approximation by replacing G0(p, s) in Equation (17) by the full Green's function G(p, s). However, as shown by Ganter and Schirmacher [1], this self-consistent scheme violates the requirement that the sound attenuation coefficient Γ(ω)∝lim|q|→0ωIm{Σ(q,s)}/q2 should vary ∝ω4 (Rayleigh scattering). A self-consistent scheme, which preserves this property has been proposed [1] by introducing auxiliary quantities g(q, s) and σ(q, s).
Within the self-consistent ERM (SCERM) approximation the self energy is given by
Σ(q,s)=∑p[t(q-p)-t(p)]2g(p,s)
This relation between Σ(q, s) and the auxiliary quantities g(q, s) and σ(q, s) has been called “cactus-1-approximation” in Ganter and Schirmacher [1]. A more complicated relation (“cactus-2”) suggested also by Ganter and Schirmacher [1] shall not be considered here.
The auxiliary quantities σ(q, s) and g(q, s) obey the self-consistent set of equations
g(q,s)=[s+t0-t(q)-σ(q,s)]-1σ(q,s)=∑p[t(q−p)−t(p)][t(q−p)−t(q)]g(p,s)2.3.1. The spectral moments in the ERM formalism
In order to relate the SCERM scheme outlined in the previous paragraph to the spectral moments, we perform a high frequency expansion of the Green's function. We take advantage of the Hilbert-Stieltjes transformation relating G(q, s) to its imaginary part:
G(q,s)=1π∫G″(q,λ)s-λdλ
Performing the |q| → ∞ limit as prescribed in Equation (15), we obtain
In order to expand the self energy it is advisable to reformulate Equation (19), using the identy,
∑ph1(q-p)h2(p)=-∑ph2(q-p)h1(p)
as
Σ(q,s)=∑p[t(p)-t(q-p)]t(p)[g(p,s)+g(q-p,s)]
The |q| → ∞ limit is then
Σ∞(s)=∑pt2(p)[g(p,s)+g∞(s)]
We are interested in the first two nontrivial moments, therefore we need the expansion in 1/s up to the cubic, 1/s3, term. By inspecting Equation (22), we note that it is sufficient to expand Σ∞(s) up to 1/s to have G∞(s) correct up to 1/s3. In Equation (25) the leading contribution is given by the high-frequency limit g(p, s) → 1/s. The leading result for Σ∞(s) is then
Σ∞(s)=2s∑pt2(p)+O(1s2)
Inserting this equation into Equation (22):
G∞(s)=1s−1s2to+1s3[to2+2∑pt2(p)]−…
By comparing this equation with Equation (21) we have:
M(0)=1M(1)=t0M(2)=t02+2∑pt2(p)=t02[1+2∑pf˜2(σp)]
Using the Parseval theorem one can easily show
∑pt2(p)=K0f0(2)MΓ=12Mc,∞(2)
so that the results Equations (11) and (28) are the same.
This means that the SCERM approximation not only preserves the Rayleigh scattering property, as emphasized in Ganter and Schirmacher [1] but also preserves the first three spectral moments. Because we used only the property lim|s|→∞g(q,s)=1/s, any theory for calculating the “internal” Green's function will have this property. The nontrivial task is to come up with such a theory. In the next but one paragraph we show that one can establish a link between the self-consistent ERM approximations [1] and the self-consistent Born approximation of the heterogeneous elasticity theory [9, 22–25]. By this correspondence we shall find a simplified theory for calculating the DOS.
2.3.2. The DOS in the SCERM formalism
For calculating the level density we need the large-wavenumber limit of the self energy Σ∞(s), Equation (25). In these expressions appears the large-wavenumber limit of the auxiliary quantities σ∞=limq→∞σ(q,s) and g∞=limq→∞g(q,s). In order to calculate these quantities we take the large wavenumber limit in the self-consistent equations (Equation 19) applying again the convolution identity (Equation 23) to obtain
σ∞=∑pt(p)2g∞g∞=1s+t0−σ∞=1s+t0−∑pt(p)2g∞
Equations (30) and (31) lead to a quadratic equation for g∞ with the solution (“Hubbard Green's function” [15])
g∞(s)≐gHE(s)=2s+t0+[s+t0]2-4∑pt(p)2
Its imaginary part is a half-ellipse (HE), centered around t0 with half width 2∑pt2(p)=2Mc,∞(2).
Combining now Equations (25) and (30) we have
Σ∞(s)=Σ∞(1)(s)+σ∞(s)
with
Σ∞(1)(s)=∑pt(p)2g(p,s)
Combining this with the expression (Equation 22) for G∞(s) we get
G∞(s)=g∞(s-Σ∞(1)(s))=gHE(s-Σ∞(1)(s))
For the level density we then obtain
ρ(λ)=1πIm{gHE(s−Σ∞(1)(s))}
with gHE(s) given by Equation (32) and Σ∞(1)(s) according to Equation (34).
2.3.3. Relation to heterogeneous-elasticity theory and new self-consistent ERM theory
If we consider the self-consistency Equation (19) we realize that the main contributions to the integrals over wavevectors are restricted to |q| < 1/σ. This means that one does not make a big error, if one makes an expansion with respect to the parameter |q|σ, i.e., a hydrodynamic expansion. For t(q) we can write
t(q)=t0-q2c02
where
c02=-12∂2∂q2t(q)|q=0
is the unrenormalized sound velocity. We now define a hydrodynamic self energy
σ1(s)=limq→01q2σ(q,s)
so that the auxiliary Green's function becomes
g(q,s)=1s+q2(c02-σ1(s))
In order to perform the limit of Equation (39) we use the Taylor Formula
where t′(q) and t″(q) are the derivatives with respect to q = |q|. Realizing that only even powers of qp contribute to the p integral in Equation (19) we obtain the following self-consistent equation for σ(q, s)
σ(q,s)=1qt′(q)∑p1pt′(p)[(qp)2+14q2p2]g(q,s)
Taking into account the angle average 〈(qp)2〉=13q2p2 and replacing t′(q)/q by its low-q limit -2c02 we get
σ1(s)=76c02∑pp2[-t′(p)/p]g(p,s)
Introducing the dimensionless wavenumber p~=σp and the dimensionles self energy σ~1(s)=σ1(s)/c02 we obtain
We mention that the prefactors of Γ in Equation (45) are model dependent constants, which do not depend on any physical parameter. For the Gaussian force constants t0σ2/c02=2 and I0=π/2.
Equation (45) is mathematical identical with the version of the SCBA for correlated fluctuating elastic constants [24, 25]. One has to make the replacements C(q) ↔ [−f′(q)/q], where C(q) is the Fourier-transformed correlation function of the elastic moduli, divided by their variance. The disorder parameter γ ∝ Γ = 1/ρσ3 is proportional to the relative variation of the fluctuating elastic modulus of a mass-spring model, evaluated by a coarse-graining procedure, as shown by Ganter and Schirmacher [1]. Because both for a Gaussian as for an exponential f~(p~) the integral in Equation (44) can be done analytically, the solution of this equation is a matter of a few iterations and is much easier and quicker than solving the integral Equation (19).
Taking these considerations into account we now propose the following simplified SCERM scheme, which we call ERM-SCBA: In this approximation we use in formula (Equation 34) for Σ∞(1)(s) the hydrodynamic g(q, s) of Equation (40) i.e.,
Σ(1)(s)=∑pt(p)21s+q2(c02-σ1(s))
with the SCBA self energy σ1(s) given by the SCBA Equation (43)
σ1(s)=76c02∑pp2[-t′(p)/p]1s+q2(c02-σ1(s))
This new self-consistent scheme for calculating approximately the spectrum of a disordered mass-spring model is the most important result of the present paper.
We emphasize that this approximation still meets the requirements of
Preserving the lowest two nontrivial moments,
Giving a Debye DOS for λ < λ*,
Leading to Rayleigh scattering in the same frequency range.
We shall investigate the accuracy of this approximation in the next chapter.
3. The gaussian force constants
In this section we report the results for the case of the Gaussian shape for the force constants K(r)=K0e-12(r/σ)2, so that we have f(r) = e−r2/2, f0=(2π)3/2 and f0(2)=π3/2. For the “vertex function” in Equation (44) we have
-1q~f′~(q~)=f~(q~).3.1. Numerical simulation
In order to compare the analytical expression for the spectral moments summarized in the next paragraph, and, more important, to test the level of accuracy of the different theoretical approximation for the DOS spectral shape, we run a numerical simulation of the dynamics in the specific case of Gaussian shape for K(rij). For different values of the couple of parameters (N, Γ) we randomly chose N points (r̄i) distributed in a box of size L = 1, set up the dynamical matrix with K0=1 and M = 1, diagonalized this matrix using the Jacobi diagonalization algorithm [45], and set up the DOS as histogram of the obtained eigenvalues by repeating ≈1,000 times the random points choosing. The so obtained spectra are compared (see later) with the results of the theoretical approximations, and their spectral moments are derived by numerical integrations.
3.2. Results for the level density
Let us collect the equations we need for calculating the ERM-SCBA density of levels for the Gaussian model:
λ~=λ/t0 and s~=s/t0 are the dimensionless frequency variables and w(z) = e−z2erfc(−iz) is the Faddeeva function [46]. It is clear that Equation (50d) is the one, which has to be iterated.
In Figure 1 we compare the level density of the numerical diagonalization for Γ = 1 with the result of the SCERM calculation of Ganter and Schirmacher [1] and our ERM-SCBA calculation for Γ = 1. We include also the level density corresponding to the half-elliptic DOS. It is seen that both self-consistent approximations give a fair description of the numerical DOS. The ERM-SCBA gives even a more symmetrical result for the “Einstein peak” around t0, compared with the SCERM result. It smoothly interpolates from the “hydrodynamic” DOS to the half-elliptic Einstein peak.
Comparison between the numerically calculated level density (N = 800) for Γ = 1 (red connected symbols) with the SCERM result of Ganter and Schirmacher [1] (blue dash-double dots), our ERM-SCBA result (red dashes) and the half-elliptic DOS (green dash-dots).
In Figure 2 we compare the level density of the numerical simulation with the ERM-SCBA for different values of Γ. The decreasing accuracy of the theoretical results with degreasing Γ might be due to the appearance of finite-size effects (see next paragraph).
Comparison between the numerically calculated level density (N = 200, blue connected symbols) with the results of ERM-SCBA (red lines) for Γ = 0.3, 0.65, 1.5, and 3.5 (from top to bottom).
We conclude that the ERM-SCBA introduced by us is a handy tool to calculate almost analytically the DOS of a spring-mass model. While the ERM-SCBA is not as accurate as the SCERM theory, it is much easier tractable, because here we have to solve only the SCBA Equation (50c), which requires a few numerical iterations, whereas for the SCERM theory one has to solve the set of three-dimensional integral Equation (19).
3.3. Results for the moments3.3.1. The average 〈<italic>f</italic>〉 and 〈<italic>f</italic><sup>2</sup>〉 in the gaussian case
We are now calculating the averages 〈f〉 and 〈f2〉 needed for the moments. We do this calculation for a finite volume V as present in our numerical calculation.
Similarly, 〈f2〉 is obtained by the previous equation with the substitution σ→σ/2:
〈f2〉=[πσLerf(L2σ)]3
Summing up:
〈f〉=(2π)3/2(σL)3erf3(2L4σ)〈f2〉=π3/2(σL)3erf3(L2σ)
In terms of (Γ, N) the averages become:
〈f〉=1N(2π)3/2(1Γ)erf3(24(NΓ)1/3)〈f2〉=1Nπ3/2(1Γ)erf3(12(NΓ)1/3)3.3.2. Check at finite size
We now compare the exact results for the moments in a finite volume, Equations (9), (11), together with Equations (55), (56) with the results of the numerical simulation as described before.
In the next figures the analytical results are compared, in log log scale, with the simulation for N = 50 and 150. The data are reported as function of Γ in the range 10−3–103. Figure 3 reports the comparison of the analytic calculation (black lines) of the first moment M(1) with the results of the numerical simulation (dots). Similarly, Figure 4 reports the comparison of the analytic calculation of the second moment M(2) with the simulation. Finally Figure 5 reports the comparison of the analytic calculation of the second central moments Mc(2) (black) and M~c(2) (blue) with the results of the numerical simulation. The red lines in the figure represent the large N (full) and small Γ (dashed) limits. All the comparisons are favorable.
Comparison between the analytical (full lines, from Equations 9 and 55) and numerical (black dots) results for the first moment M(1) of the spectral density of states ρ(λ). λ = ω2 are the eigenvalues of the dynamical matrix and ω are the characteristic frequencies. The data at the indicated selected N values are reported as a function of the parameter Γ. The N values (N = 50, 100, 150, 200, 300, 500, 1,000, and 10,000 for the analytical case and 50 and 150 for the numerical one) are indicated on the black curves. The red line is the asymptotic, N→∞, limit reported in Equation (58).
Comparison between the analytical (full lines, from Equations 11, 55, and 56) and numerical (black dots) results for the second moment M(2) of the spectral density of states ρ(λ). The data at the indicated selected N values are reported as a function of the parameter Γ. The N values (N = 50, 150, and 500 for the analytical case and 50 and 150 for the numerical one) are indicated on the black curves. The red line is the asymptotic, N→∞, limit reported in Equation (58).
Comparison between the analytical (full lines, from Equations 11, 55, and 56) and numerical (dots) results for the second central moment M(2) (black) and for the “non trivial” second central moment M~(2) (blue) of the spectral density of states ρ(λ). The data at the indicated selected N values are reported as a function of the parameter Γ. The N values (N = 50, 150, and 500 for the analytical case and 50 and 150 for the numerical one) are indicated on the black and blue curves. The red full line is the asymptotic, N→∞, limit reported in Equation (58), it turns out to be the same for M(2) and M~(2). The dashed red lines are the low Γ limit, reported in Equation (60).
3.3.3. Infinite size limit
In the limit of infinite size, L→∞, N→∞, ρ = N/L3 = constant, from Equations (55) and (56)
we get (with erf(x) → 1 for large x):
N→∞{N〈f〉=(2π)3/2(1Γ)N〈f2〉=π3/2(1Γ)
Thus, the moments become (in the infinite size limit M~(n)=M(n)):
In the infinite size limit the moments only depend on Γ, which justify the choice of this parameter instead of ρ or/and σ.
These infinite size equations are reported as full red lines in Figures 3–5.
3.3.4. Large and small Γ limits
Large Γ is similar to large L, thus the asymptotic expansions in Equation (58) are valid at large Γ with the further simplification that the term in Γ−2 can be neglected in the expression for M(2):
These limits are easily verified in the three figures previously reported. The equation for M~c(2) is reported in Figure 5 as dashed red line.
4. Discussion and conclusion
In the paper we have shown that the self-consistent euclidean-matrix (SCERM) approximation of Ganter and Schirmacher [1] preserves the two nontrivial moments of the eigenvalue spectrum of mass-and spring models. This confirms the conclusion that the ERM formalism is a very powerful one and is a good starting point for further work.
We have calculated the exact expressions for the 1st and 2nd spectral moments for finite and infinite samples. We show that in the limit of large interaction parameter σ there is a strong size (N) dependence of the spectral moments, thus posing a caveat on the extrapolation of finite-size molecular dynamics simulations to the thermodynamic, large N, limit.
We have established a link between the self-consistent ERM (SCERM) approximation of Ganter and Schirmacher [1] and the self-consistent Born (SCBA) scheme of the heterogeneous-elasticity theory of vibrational anomalies [9, 23]. By this we constructed a simplified version of the SCERM approximation, the ERM-SCBA, which compares well to our numerical diagonalization of the Gaussian-force-constant model. As the SCBA version of herogeneous-elasticity theory has served [9, 23, 47] to explain the vibrational anomalies related to the boson peak this is—in a nutshell—also the case for the new ERM-SCBA.
A more realistic theory, applicable to amorphous solids and liquids in the high-frequency regime should not be based on a scalar mass-spring model but on the full vectorial equations of motions of a disordered solid, interacting via pair potentials. The force constants, which appear in this theory are the second derivatives of the pair potentials. Furthermore the configurational averages should not be based on the random statistics of points but on the statistics of atoms in an amorphous solid, consistent with the pair potentials. Such a theory has been formulated within the ERM scheme by Ciliberti et al. [48], but unfortunately contained a mistake concerning the summation of diagrams, such that the Rayleigh-scattering property was not obtained [1, 49], and the other results are, therefore, questionable. A vectorial version of the ERM-SCBA for the description of amorphous solids will be worked out by the present authors.
Let us mention a drawback of the present approach: The ERM approach is based on an expansion of the averaged Green's function with respect to the inverse density. This means that the theory looses its application range once the density parameter Γ∝ρ−1 becomes much larger than one. In fact, in this regime the spectrum becomes unstable, i.e., negative values of λ are predicted [39, 40, 48]. However, as the force constants Kij are all positive definite so must be the ensemble of eigenvalues. This artifact can be avoided by using the coherent-potential approximation (CPA) [43] instead of the SCBA. A connection between the ERM and CPA formalism which gives a theory valid for low densities, will be published shortly by the present authors.
Author contributions
All authors listed, have made substantial, direct and intellectual contribution to the work, and approved it for publication.
Conflict of interest statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
ReferencesGanterCSchirmacherW. Rayleigh scattering, long-time tails, and the harmonic spectrum of topologically disordered systems. Phys Rev B (2010) 82:094205. 10.1103/PhysRevB.82.094205SetteFKrischMHMasciovecchioCRuoccoGMonacoG. Dynamics of glasses and glass-forming liquids studied by inelastic X-ray scattering. Science (1998) 280:1550–5. 10.1126/science.280.5369.1550BinderKKobW. Glassy Materials and Disordered Solids: An Introduction. London: World Scientific (2011).ZellerRCPohlRO. Thermal conductivity and specific heat of noncrystalline solids. Phys Rev. B (1971) 4:2029. 10.1103/PhysRevB.4.2029ShintaniHTanakaH. Universal link between the boson peak and transverse phonons in glass. Nat Mater (2008) 7:870–7. 10.1038/nmat229318849975HorbachJKobWBinderK. High frequency sound and the boson peak in amorphous silica. Eur Phys J B (2001) 19:531–43. 10.1007/s100510170299NakayamaT. Boson peak and terahertz frequency dynamics of vitreous silica. Rep Prog Phys. (2002) 65:1195. 10.1088/0034-4885/65/8/203GurevichVLParshinDAShoberHR. Anharmonicity, vibrational instability, and the Boson peak in glasses. Phys Rev. B (2003) 67:094203. 10.1103/PhysRevB.67.094203SchirmacherWRuoccoGScopignoT. Acoustic attenuation in glasses and its relation with the boson peak. Phys Rev Lett. (2007) 98:025501. 10.1103/PhysRevLett.98.02550117358618RuffléBParshinDACourtensEVacherR. Boson peak and its relation to acoustic attenuation in glasses. Phys Rev Lett. (2008) 100:015501. 10.1103/PhysRevLett.100.01550118232782ChumakovAIMonacoGMonacoACrichtonWABosakARfferR. Equivalence of the boson peak in glasses to the transverse acoustic van hove singularity in crystals. Phys Rev Lett. (2011) 106:8519. 10.1103/PhysRevLett.106.22550121702612ChumakovAIMonacoGFontanaABosakAHermannRPBessasD. Role of disorder in the thermodynamics and atomic dynamics of glasses. Phys Rev Lett. (2014) 112:025502. see also the article of A. L. Chumakov and G. Monaco in this volume. 10.1103/PhysRevLett.112.02550224484025RuoccoGFSetteRDi LeonardoDFiorettoMLorenzenMKrischC. Nondynamic origin of the high-frequency acoustic attenuation in glasses. Phys Rev Lett. (1999) 83:5583. 10.1103/PhysRevLett.83.5583RuoccoGSetteFDi LeonardoRMonacoGSampoliMScopignoT. Relaxation processes in harmonic glasses?Phys Rev Lett. (2000) 84:5788–91. 10.1103/PhysRevLett.84.578810991055EconomouEN. Green's Function in Quantum Physics. Heidelberg: Springer-Verlag (1971).BurinALMaksimovLAPolishchukIY. The phonon transport in crystals with heavy defects. Phys B (1995) 210:49–54. 10.1016/0921-4526(94)00909-FPolishchukIYMaksimovLABurinAL. Localization and propagation of phonons in crystals with heavy impurities. Phys. Reports (1997) 288:205–22. 10.1016/S0370-1573(97)00025-2KarpovVGKlingerIIgnat'EvFN. Theory of the low-temperature anomalies in the thermal properties of amorphous structures. Zh Eksp Teor Fiz (1983) 84:760–75.BuchenauUGalperinYMGurevichVLSchoberHR. Anharmonic potentials and vibrational localization in glasses. Phys Rev B (1991) 43:5039–45. 10.1103/PhysRevB.43.50399997881BuchenauUGalperinYMGurevichVLParshinDARamosMASchoberHR. Interaction of soft modes and sound waves in glasses. Phys. Rev. B (1992) 46:2798. 10.1103/PhysRevB.46.279810003968GurevichVLParshinDASchoberHR. Pressure dependence of the boson peak in glasses. Phys Rev B (2005) 71:014209. 10.1103/PhysRevB.71.014209MaurerESchirmacherW. Local oscillators vs. elastic disorder: a comparison of two models for the boson peak. J Low-Temperature Phys. (2004) 137:453–70. 10.1023/B:JOLT.0000049065.04709.3eSchirmacherW. Thermal conductivity of glassy materials and the “boson peak. Europhys Lett. (2006) 73:892. 10.1209/epl/i2005-10471-9SchirmacherWSchmidBTomarasCVilianiGBaldiGRuoccoG. Vibrational excitations in systems with correlated disorder. Phys Status Solidi (c) (2008) 5:862. 10.1002/pssc.200777584SchirmacherWSchmidBTomarasCBaldiGVilianiGRuoccoG. Sound attenuation and anharmonic damping in solids with correlated disorder. Condensed Mater Phys. (2010) 13:23606. 10.5488/CMP.13.23606TomarasCSchirmacherW. High-frequency vibrational density of states of a disordered solid. J Phys Condens Matter (2013) 25:495402. 10.1088/0953-8984/25/49/49540224214818SchirmacherW. The boson peak. Phys Stat Sol. (b) (2013) 250:937. 10.1002/pssb.201248544SchirmacherWScopignoTRuoccoG. Theory of vibrational anomalies in glasses. J Noncryst Sol. (2014) 407:133–40. 10.1016/j.jnoncrysol.2014.09.054SchirmacherWWagenerM. In: DianouxAJPetryWTeixeiraJ, editors. Dynamics of Disordered Materials.Heidelberg: Springer (1989).XuBCStrattRM. Liquid theory for band structure in a liquid. II. P orbitals and phonons. J Chem Phys. (1990) 92:1923. 10.1063/1.458023WuTMLoringRF. Phonons in liquids: a random walk approach. J Chem Phys. (1990) 97:8568. 10.1063/1.463375SchirmacherWWagenerM. Analogies between hopping and phonon propagation in disordered solids. Philos Magn B (1992) 65:861. 10.1080/13642819208204892SchirmacherWWagenerM. Vibrational anomalies and phonon localization in glasses. Sol State Comm. (1993) 86:597–603. 10.1016/0038-1098(93)90147-FKühnRHorstmannU. Random matrix approach to glassy physics: low temperatures and beyond. Phys Rev Lett. (1997) 78:4067. 10.1103/PhysRevLett.78.4067SchirmacherWDiezemannGGanterC. Harmonic vibrational excitations in disordered solids and the “boson peak”. Phys Rev Lett. (1998) 81:136. 10.1103/PhysRevLett.81.136SchirmacherWDiezemannGGanterC. Model calculations for the vibrational anomalies of a disordered Lennard Jones solid. Phys B (2000) 284–288:1147. 10.1016/S0921-4526(99)02550-8Martín-MayorVParisiGVeroccioP. Dynamical structure factor in disordered systems. Phys Rev E (2000) 62:2373. 10.1103/PhysRevE.62.237311088717TaraskinSNElliottSLohYHNataranjanG. Origin of the boson peak in systems with lattice disorder. Phys Ref Lett. (2001) 86:1255–8. 10.1103/PhysRevLett.86.125511178057GrigeraTSMartin-MayorVParisiG. Verrocchio P. Vibrational spectrum of topologically disordered systems. Phys Rev Lett. (2001) 87:085502. 10.1103/PhysRevLett.87.085502Martín-MayorVMzardMParisiGVerrocchioP. The dynamical structure factor in topologically disordered systems. J Chem Phys. (2001) 114:8068. 10.1063/1.1349709GrigeraTSMartín-MayorVParisiGVerrocchioP. Phonon interpretation of the ‘boson peak’ in supercooled liquids. Nature (2003) 422:289–92. 10.1038/nature0147512646916AmirAKrichJJVitelliVOregYImryY. Emergent percolation length and localization in random elastic networks. Phys Rev X (2013) 3:021017. 10.1103/PhysRevX.3.021017KöhlerSRuoccoRSchirmacherW. Coherent potential approximation for diffusion and wave propagation in topologically disordered systems. Phys Rev B (2013) 88:064203. 10.1103/PhysRevB.88.064203HubbardJBeebyJL. Collective motion in liquids. J Phys C (Solid St. Phys.) (1969) 2:556.GolubGHvan der VorstA. Eigenvalue computation in the 20th century. J Comput Appl Math. (2000) 123:35–65. 10.1016/S0377-0427(00)00413-1ElliottSR. Handbook of Mathematical Functions. New York, NY: Wiley–Interscience (1973).MarruzzoASchirmacherWFratalocchiARuoccoG. Heterogeneous shear elasticity of glasses: the origin of the boson peak. Nat Sci Rep. (2013) 3:1407. 10.1038/srep0140723470597CilibertiSGrigeraTSMartin-MayorVParisiGVerrocchioP. Brillouin and boson peaks in glasses from vector Euclidean random matrix theory. J Chem Phys. (2003) 119:8577. 10.1063/1.1610439GrigeraTSMartín-MayorVParisiGUrbaniPVerrocchioP. On the high-density expansion for Euclidean random matrices. J Stat Mech. (2011) 2011:P02015. 10.1088/1742-5468/2011/02/P02015AppendixA.1. Detailed calculation of the momentsA.1.1. First moment
From Equation (7) the first (n = 1) moment is given by the trace of the dynamical matrix:
The sum appearing in the rightmost term in this equation is extended over N(N − 1) terms, all different each others, thus statistically independent, and all arising from the same probability distribution function P(f). Therefore we can write:
M(1)=K0M(N-1)〈f〉
The whole set of first moments, according to their definition, turn out to be:
M(1)=K0M(N-1)〈f〉Mc(1)=0M~(1)=K0MN〈f〉M~c(1)=0A.1.2. Second moment
Following the same steps as in the case of the first moment, we have:
In the first sum of the latter expression there are N(N − 1) terms fik2 (arising from k = k′ ≠ i) and N(N − 1)(N − 2) terms fikfik′ with all the three indexes i, k, k′ different each others. In the second summation there are N(N − 1) terms fij2 (i ≠ j). Therefore:
1In their original papers on the ERM formalism the authors claimed that the model would not exhibit Rayleigh scattering [39, 40]. Ganter and Schirmacher [1] and later Grigera et al. [49] showed that this statement was in error.