The mixing of polarizations in the acoustic excitations of disordered media with local isotropy

An approximate solution of the Dyson equation related to a stochastic Helmholtz equation, which describes the acoustic dynamics of a three-dimensional isotropic random medium with elastic tensor fluctuating in space, is obtained in the framework of the Random Media Theory. The wavevector-dependence of the self-energy is preserved, thus allowing a description of the acoustic dynamics at wavelengths comparable with the size of heterogeneity domains. This in turn permits to quantitatively describe the mixing of longitudinal and transverse dynamics induced by the medium's elastic heterogeneity and occurring at such wavelengths. A functional analysis aimed to attest the mathematical coherence and to define the region of validity in the frequency-wavevector plane of the proposed approximate solution is presented, with particular emphasis dedicated to the case of disorder characterized by an exponential decay of the covariance function.


INTRODUCTION
Most materials we encounter on a daily basis, such as glasses, polycrystalline aggregates, ceramics, composites, geophysical materials, and concrete can be classified as heterogeneous materials, being composed by domains with different physical characteristics. An acoustic wave propagating in a three-dimensional system is characterized by its phase velocity, amplitude and polarization. In a heterogeneous medium the acoustic excitations experience retardation, attenuation (Rayleigh anomalies) and depolarization. Strong attention has been deserved in literature both to the Rayleigh anomalies and to the mixing of polarizations . They have been, however, designed as disjoint phenomena and never been addressed by an analytical theory as related aspects originating from a common root, the disordered nature of the medium. An analytical theory describing the mixing of polarizations of acoustic excitations in disordered systems is, furthermore, so far lacking. One of the challenge in obtaining an unified picture of the above-metioned phenomena is their occurrence on different length-scales. In the so-called Rayleigh region, i.e., for values of wavelength (λ) of elastic excitations much lower than the characteristic size (a) of inhomogeneity domains, the phase velocity of acoustic modes shows a softening with respect to its hydrodynamic value (retardation). It is observed, moreover, a strong increase of the acoustic wave attenuation (Rayleigh scattering), the two quantities being related to each other by Kramers-Kroning relations [11]. The coupling between longitudinal and transverse polarizations is instead maximum beyond the Rayleigh region when λ ∼ a [23]. The basic analytical instrument to describe the ensemble averaged elastodynamic response of a heterogeneous system to an impulsive force is the so-called Dyson equation for the mean field [23][24][25][26][27], introduced in the framework of the Random Media Theory (RMT), or Heterogeneous Elasticity Theory when referring to the specific case of elastic inhomogeneity [1,[7][8][9][10][11]28]. The solution of the Dyson equation, however, can be obtained only under suitable approximations, which have a limited wavelengths range of validity [24][25][26], thus avoiding a unified theoretical description of experimental observations.
The exact solution of the Dyson equation can be formally cast in a Neumann-Liouville series, the so-called perturbative series expansion. Even if in most real cases a direct sum, or even establishing the criteria of convergence of the series, is not possible, it constitutes the general starting point for smoothing methods or approximations [24]. Its truncation to the lowest non-zero order leads to the so-called Born Approximation [9, 23-27, 29, 30]. We propose an approximate solution of the vectorial Dyson equation, which takes into account in an approximate form terms of the Neumann-Liouville series up to the second order, thus introducing corrective terms to the Born Approximation. We will refer to it as to a Generalized Born Approximation (GBA). We first derive an analytical expression for the GBA and state the general conditions for its validity in a given region of wavelengths. In a second stage, we analyze the specific case of an exponential decay of the covariance function of elastic fluctuations. We show how in such a case the GBA can be applied up to wavelengths of the order of the average size of heterogeneity domains. We then calculate in the GBA frame the current spectra related to acoustic dynamics and show that the GBA allows for a description of all the effects that the topological disorder has on acoustic dynamics, including the Rayleigh anomalies and, for the first time, the mixing between longitudinal and transverse polarization.
The rest of the paper is organized as in the following. In section 2, we introduce the GBA, we discuss its physical significance with the support of the Feynman diagram technique and its relationship with the perturbative series expansion. In section 3.1, we describe with mathematical detail the proposed approximation and demonstrate its validity in a proper domain of the frequency (ω) -wavevector (q) plane. In section 3.2, we deal with the specific case of an exponential decay of the covariance function and define in this case the domain of validity of the GBA. In section 3.3 we discuss how the GBA can account for the mixing of polarizations. In section 3.4, we show what the acoustic dynamics properties accounted by the GBA are, allowing a qualitative comparison with existing experimental results. Conclusions are outlined in section 4. Technical details in addition to the main text are reported in two appendices (Appendices A, B).

The Dyson Equation and Its Approximate Solutions in the Random Media Theory
The elastic response of an unbounded and elastic medium to an impulsive force can be obtained as a function of the Green's dyadic by solving the so-called stochastic Helmholtz equation [24,27], Summation over repeated indices is assumed. The second-rank Green's dyadic, G ij (x, x ′ , t), is the response of the system at the spatial point of vectorial coordinate x in the i-th direction at the delay time t to a unit-impulse at the point x ′ in the j-th direction. The function δ kj is a Kronecker delta function, δ(x) and δ(t) are Dirac delta functions. The elastic tensor, C ijkl (x), and the density, ρ(x), of the system are spatially heterogeneous.
. We hypothesize statistical homogeneity, thus C 0 ijkl = <C ijkl (x) > and ρ 0 = <ρ(x) >. The brackets < > denote the ensemble average. The operatorsL 0 (x, t) andL s (x, t) are respectively a deterministic differential operator related to the average, constant in space, elastic tensor and density and a linear stochastic operator accounting for the fluctuating, space-dependent, terms of the same quantities. We assume statistical and local isotropy and express the elastic tensor as a function of the shear modulus, µ(x) = µ 0 + δµ(x), and of the Lamé parameter, λ(x) = λ 0 + δλ(x). Under these hypotheses the operatorsL 0 (x, t) andL s (x, t) are given by [27] We take under exam the case of spatial fluctuations of the elastic tensor, thus the first term in Equation (3) is zero. The quantity physically relevant, related to the dynamic structure factor, which can be accessed, e.g., by Inelastic X-ray Scattering (IXS) or Inelastic Neutron Scattering (INS), is the ensemble averaged Green's dyadic, < G(x, x ′ , t) >. In place of solving the Helmholtz equation (impossible in most cases) and then averaging, one can look for a suitable expression of an effective deterministic operator,D, such that [24] The latter equation is referred to as the Dyson equation. Drawbacks in the definition ofD comes, however, from the fact that the operatorL(x, t) =L 0 (x, t) +L s (x, t) cannot be inverted.
The Dyson equation can be rephrased by setting a formal expression for the average Green's dyadic [27], The integrals are extended to R 3 . G 0 (x, x ′ , t) is the Green's dyadic of the "bare" medium, solution of the deterministic Helmholtz equation related to the operatorL 0 (x, t). Equation (5) is based on the introduction of the so-called mass operator or self-energy, (x, x ′ , t), which embeds all the information related to disorder. The problem of solving the Dyson equation thus translates into finding a suitable expression for the self-energy. The self-energy can be cast in a Neumann-Liouville series by starting from the related stochastic Helmholtz equation, giving rise to the so-called perturbative series expansion [24]. This generally constitutes the starting point for smoothing methods or approximations [24,25].
In the Fourier space Equation (5) becomes where q and ω are respectively the conjugate variables of x and t.

The Born Approximation
Under the hypothesis of statistical homogeneity, truncation of the perturbative series expansion to the lowest non-zero order leads to the so-called Bourret or Born Approximation [9, 23-27, 29, 30]. In the Fourier space it states The operatorL 1 is related to the operatorL s defined in Equation (3) by ensemble averaging and Fourier transforming [27]. Since we only account for fluctuations of the elastic tensor, the operator L 1 can be expressed by introducing the covariance function of the elastic tensor fluctuations,R γ αjlβkiδ (x = x 1 − x 2 ) =< δC γ αjl (x 1 )δC βkiδ (x 2 ) >. Equation (7) becomes [27] B where the wavevector q ′ is the variable of integration. It is q = |q|. The self-energy in the Fourier space can thus be written as a convolution between the "bare" Green's dyadic and the Fourier transform of the covariance function of the elastic tensor fluctuations. Despite simplicity, the Born Approximation imposes rather strong restrictions both on the intensity of the elastic constants fluctuations and on the values of q and ω with respect to a. A necessary condition for the validity of the Born Approximation is indeed to deal with small values of the intensity of elastic fluctuations and of wavevector and frequency [25,29,30]. The conditionǫaq(q 0i ) ≪ 1 shall be met. It is q 0i = ω/c 0 i , where c 0 i is the phase velocity of the acoustic excitations with i-th polarization in the "bare" medium. The parameterǫ 2 is the "disorder parameter" [7,11], i.e., the square of the intensity of spatial fluctuations of elastic constants normalized to their average value, whereas ǫ 2 represents the same unrenormalized quantity.

The Self-consistent Born Approximation
The so-called Self-Consistent Born Approximation (SCBA) [1,7,11,28] or Kraichnan model [29,31] can be derived from the more general mean field theory, the Coherent Potential Approximation (CPA) [32][33][34][35][36], under the hypothesis of small fluctuations [36]. It is, however, not affected by the same small wavevectors and frequencies limitation than the Born Approximation [29,36]. In place of truncating the Neumann-Liouville series in the SCBA frame it is constructed an effective nonlinear deterministic equation defining the average Green's dyadic. Because this latter can be related to a realizable model and it can be exactly solved, the SCBA solution will guarantee certain consistency properties [31]. The self-energy in the SCBA is given by Equation 9 and the Dyson equation, Equation 6, form a set of selfconsistent equations. At the step n = 0 it is < G(q, ω) > n=0 = G 0 (q, ω). Even if the logic behind the two approaches, i.e., the truncation of the Neumann-Liouville series defining the exact solution of the problem leading to the Born Approximation, or the mean-field approach behind the CPA leading to the SCBA, is different, Equation (9) can easily allow a connection between the two: the expression at the first step of the SCBA is the same obtained from the Born approximation. Accordingly we expect for those cases where it is applicable, the SCBA to provide a better approximation than the Born Approximation. Generalizations of the Born Approximation in the framework of the SCBA have attracted interest in several fields of physics [37][38][39][40]. A link between the SCBA and the perturbative series expansion is discussed in section 2.2 by exploiting the Feynman diagram technique. An analytical calculation of the self-energy in the SCBA frame is possible by assuming at each step of the self-consistent calculation q = 0 in the expression of the mass operator [1,7,11,28]. By exploiting this approach the SCBA revealed to correctly describe the Rayleigh anomalies of acoustic waves in a topologically disordered medium [1,7,11]. The SCBA thus revealed to give an answer to important questions such as how does the attenuation and phase velocity vary with the wavevectors in the Rayleigh region. We could also expect that the SCBA in a three-dimensional space can carry information about the polarization properties of the acoustic waves. This kind of study, however, can be hindered by the impossibility to obtain an analytical calculation of the SCBA self-energy for aq ∼ 1, that is at the edge of the Rayleigh region where the strong acoustic wave attenuation starts to slow down and the mixing of polarizations is expected to get in.

The Mixing of Polarizations Beyond the Born Approximation
We introduce an approximate method (GBA) for the calculation of (q, ω). It permits to obtain corrective terms to the Born Approximation in the context of the perturbative series expansion, as discussed in the next section. We discuss in section 3.3 how the GBA permits to describe the mixing of polarizations at the boundary of the Rayleigh region (aq ∼ 1) and in section 3.4 how it allows to describe, together with the mixing of polarizations, also the acoustic anomalies occurring in the Rayleigh region (aq < 1). Similar results cannot be achieved by using the Born Approximation, thus making the two approximations qualitatively different. A sharp increase in the attenuation of the acoustic excitations and a related kink in the phase velocity at aq ∼ 1 are features related to the coupling between longitudinal and transverse dynamics [23]. They can be described by the Born Approximation in the three-dimensional space [23]. By exploiting the Born Approximation, however, we couldn't unravel the presence of a clear "projection" of the transverse into the longitudinal acoustic dynamics, as instead attested by experimental observations in several topologically disordered systems [16][17][18]20]. With "projection" it is meant the occurrence in the longitudinal dynamic structure factor of a peak-like feature centered at frequencies characteristic of the transverse excitations and occurring at sufficiently high wavevectors. We can attribute such a failure to the fact that the Born Approximation has a limited range of validity in the wavevector space, as discussed above. In particular it shifts toward lower values of wavevector for higher values of the disorder parameter. Depending on the value of the disorder paramater, its validity at q ∼ a −1 can thus be questioned. Since most of the phenomenology observed in real systems, including the Rayleigh anomalies, can however be qualitatively grasped even by the Born Approximation [23] corresponding to first order truncation of the Neumann-Liouville series, we choose to take under consideration the next order approximation in the perturbative series expansion. It corresponds to truncate the SCBA to the second order [25]. It not only permits to obtain a qualitative description of the phenomenology but also to fit experimental outputs for a real system [41]. In particular, as we discuss in section 3.3, the GBA permits to describe the "projection" of the transverse dynamics observed in the longitudinal dynamics obtaining results which are qualitatively different from what it is possible to achieve with the Born Approximation.

Basic Considerations
The physical meaning of the Dyson equation as well as of the related approximations can be better understood with the aid of the Feynman diagram technique [25]. The perturbative series can be rephrased as a sum of appropriate infinite subsequences of the same series. The exact series cannot be summed up, but some of the subsequences can [25]. Approximations, among which the Born Approximation, are constructed by summation of one or more of the infinite subsequence extracted from the perturbative series [25]. It is possible to establish a one to one correspondence between the analytic expressions, which we exploit in this text, and the Faynman diagrams. The diagram technique, however, has the advantage to permit to classify the infinite subsequence entering in the perturbative series depending on scattering events. Within this outlook the Feynman diagrams are classified as strongly or weakly connected diagrams [25]. Weakly connected diagrams are those that can  be always divided into strongly connected diagrams. The selfenergy can hence be represented as the hierarchical sum of all the strongly connected diagrams. The topology of different strongly connected diagrams is finally related to different kind of multiple scattering events. The Born Approximation, for example, accounts only for double scattering from the same inhomogeneity of an otherwise freely propagating wave, see Figures 1, 2. It is indeed obtained through the sum of an infinite subsequence of diagrams, which contains one only kind of strongly connected diagram, whose topology describes the above-quoted process. The next order approximation, which will include the next infinite subsequence of diagrams from the exact expansion of the mass-operator, can be obtained with the analytical expression of the mass-operator stated in the Born Approximation (Equation 8) by substituting the "bare" Green's dyadic, G 0 (q, ω), with the approximate expression of the mean Green's dyadic obtained by the Born Approximation itself [25]. This corresponds to the expression obtained by truncation of the SCBA expression to the second iteration step. In terms of diagram technique, it permits the inclusion of Feynman diagrams accounting for a sequence of scattering between two different inhomogeneities [25], see all possible multiple scattering events are, clearly, included. The approximation presented here can thus be view as a scheme for a partial inclusion of contributions from multiple scattering events. On this perspective also the SCBA can be thought as a sum of some of the infinite subsequences composing the perturbative series. Truncation of the SCBA to the third iteration step will include further scattering events not accounted for when the self-energy is obtained by truncation to the second step, and so on. The exact solution of the Dyson equation is unknown.
It is thus not possible to establish what is the error related to the SCBA as well as to its truncation to the second step of the iterative procedure. We can however assume that a necessary condition for the truncation to the n-th step of the SCBA to give an approximate expression of the self-energy is On this ground, Rytov and Kravtsov [25] established the necessary condition for the validity of the Born approximation previously stated. A necessary condition of validity for the proposed approximation can thus be given by the inequality | 3 (q, ω) − 2 (q, ω)| ≪ | 2 (q, ω)|. It is shown in Appendix B that in the domain of the (ω, q) plane where the series representation introduced in section 3.1 approximates the quantity 2 (q, ω) this inequality is satisfied if the magnitude of the remainder function of order one of the series representation of 2 (q, ω) is small enough. It is furthermore shown that in such a domain the necessary condition of validity for the GBA is less stringent than for the Born Approximation.
Depolarization effects in the scattering of electromagnetic waves by an isotropic random medium has been predicted by exploiting a second order representation for the scattered intensity [42]. The scattering of electromagnetic waves by the random media is cast in terms of Green's dyadic and the formal solution of the problem is given in terms of a Neumann iteration series. The n-th order of the scattered intensity is obtained by truncation of the Neumann series and ensemble averaging. Depolarization effects are also observable even in the first-order scattered intensity from an anisotropic random medium [43]. We recall that the optical theorem establish a connection between the self-energy and the intensity operator characterizing the Bethe-Salpter equation, which permits to describe the intensity of the mean field [25]. These results thus emphasize the soundness of our findings.
The input parameters of the theory are the correlation length, a, the disorder parameter,ǫ 2 , and the longitudinal and transverse phase velocity of the "bare" medium, c 0 L(T) .

The Generalized Born Approximation
Under the hypothesis of local isotropy it is convenient to introduce the orthonormal basis defined by the direction of wave propagation,q, and the two orthogonal ones [27]. On this basis all the "bare" Green's dyadic, average Green's dyadic and self-energy are diagonal. The "bare" Green's dyadic becomes with "T" and "L" labeling transverse and longitudinal modes, respectively. The longitudinal and transverse "bare" Green's functions, g 0 L (q, ω) and g 0 T (q, ω) respectively, can be formally written by following a regularization procedure [44] as where . In Equation (11) η is a positive real variable, the symbol p.v. states for the Cauchy principal value and sgn(x) is the sign function of argument x. The retarded solution is selected as required by the causality principle [45]. Furthermore, The GBA address an approximate expression of the self-energy obtained by truncation of the perturbative series expansion to the second order. This is obtained by substituting the "bare" Green's dyadics in Equation (8) (Born Approximation) with the expression of the average Green's dyadic obtained by the Born Approximation [25]. It is thus equivalent to truncate Equation 9 to the second iteration step. We obtain for the diagonal terms of the self-energy, The suffix 1 marks a quantity calculated to the first step of the self-consistent procedure. The repeated indexes kk, ii = L, T. The longitudinal and transverse self-energy are thus respectively composed by two terms accounting for the coupling with longitudinal and transverse dynamics respectively, i.e., L(T) = LL(TT) + LT(TL) .
The expression in curly bracket in Equation (14), < G ii (q, ω) > 1 , is then formally expanded in a Taylor series with respect to ǫ 2 ii (q, ω η ). Theorem I below states that this series is convergent almost everywhere (a.e.) in the domain of the (ω, q) plane where the conditions ǫ 2 we can then possibly find a sub-domain, as specified in Corollary II, where θ (x) being the Heaviside function of argument x and q i Max representing the q-boundary of the domain of convergence of the series representation of < G ii (q, ω) > 1 . Equation (15) provides the expression of the self-energy in the GBA.
The wavevector-dependence of ˜ 1 ii (q, ω) is determined by the covariance function used to statistically describe the elastic heterogeneity of the system, as established in Equation (8) above. We analyze in detail the case of an exponential decay of the covariance function with correlation length a and amplitude ǫ 2 . This choice grounds on simplicity and on the fact that several systems can be described by such a covariance function [46]. We show in section 3.2, in particular, that in this case the domain of validity of the GBA includes the region aq ∼ 1, where the mixing of polarization is expected.
In the following we demonstrate that it exists a domain of the (ω, q) plane, where the function < G(q, ω) > 1 admits a.e. the power series expansion specified in the following Theorem I.
Before to proceede with the proof, we observe that Theorem I can be applied to cases easily realizable by real systems. Because Im[ ˜ 1 ii (q, ω)] is proportional to the attenuation of the acoustic excitations in the random medium calculated to the first order of the perturbative series expansion, in point iii) it is required that such an attenuation is finite for finite values of q and ω. Furthermore, the series a. e. convergence is ensured in a region of wavevectors where the acoustic excitations in the random medium can still be described through a finite and sufficiently small correction with respect to the acoustic excitations in the "bare" medium, being there the self-energy sufficiently small.
From the algebraic equality 1 The remainder function is defined as .
We demonstrate in the following that |R i N | admits an upper bound, i.e. where . The latter quantity exists as a consequence of the hypotheses of Theorem I. It is indeed immediate to recognize that for q, ω ∈ R and ω = 0, if Im[ ˜ 1 ii (q, ω)] is strictly positive, the function is bounded, do not having poles.
To prove the inequality in the third side of Equation (16) we need to show that where for sake of simplicity we renamed the variable η c i as η. We discuss only the caseq 0i < q i Max . The caseq 0i > q i Max can be easily reconducted to the former by noticing that where z is a generic complex variable, z is its complex conjugate and θ = arg(z). Furthermore Frontiers in Physics | www.frontiersin.org 6 October 2018 | Volume 6 | Article 108 The inequality 1 = [sin 2 (x) + cos 2 (x)] 1 2 ≤ | sin(x)| + | cos(x)| has been exploited. We furthermore considered that |cos(x)| = cos(x) · sgn{cos(x)}. The same applies to sin(x). It is thus (20) where θ (η) = arg{q 2 − (q 0i + iη) 2 }. In the framework of a generalization of the Sokhotski-Plemelj theorem [47] due to Fox [48] it is possible to show that [49] where a, b, x 0 and x are real variables: a < x 0 < b, f (x) is a function which admits a complex extension f (z) that is analytic in a region of the complex plane containing the interval [a, b] but not x 0 , R x 0 = R\{x 0 }, γ ± (x 0 ) is a path of the region R x 0 from a to b belonging to the upper (lower) half-plane of the complex plane. The second side of Equation (20) can thus be rephrased as where γ ± (q 0i ) is the contour of the upper (lower) complex halfplane obtained by deformation of the segment [0, q i Max ] around q 0i by an infinitesimal arc of circle of radius φ passing aroundq 0i clockwise (counterclockwise). Furthermore it is θ = Arg(z 2 −q 2 0i ) and θ = Arg(z 2 −q 2 0i ). The symbol # denotes the Hadamard Finite-Part Integrals (or Cauchy Principal Value when N = 0) [47 -49]. We observe that i) θ = −θ ; ii) for z ∈ R it is θ = θ = 0; iii) for z ∈ C it is sgn{sin[(N + 1)θ ]}| z∈γ + (q 0i ) = sgn{sin[(N + 1)θ ]}| z∈γ − (q 0i ) , sgn{cos[(N + 1)θ]}| z∈γ + (q 0i ) = sgn{cos[(N + 1)θ ]}| z∈γ − (q 0i ) . The last passage in Equation (22) follows from (i) the fact that [49] and ii) from the Residue Theorem, where Res (N) [f (z), x 0 ] is the residue of order N of the function f (z) around the pole at z = x 0 enclosed in the closed path of the complex plane γ − (x 0 ) − γ + (x 0 ). It is Using integration by parts and the Stokes' formula one obtains [50] The latter equality ensures the existence of the Hadamard Finite-Part Integral, because the Cauchy p.v. exists as a consequence of the fact that, the N-th order derivative of the function satisfies the Lipschitz condition. It is possible to exploit the Residue Theorem to calculate the integral in Equation (26) and verify that such an integral is equal to zero. In Appendix B, we furthermore prove that The validity of Equation (17) hence follows from Equation (22). From Equation (16), we can finally conclude that if (ω) < 1 and ω = 0, it is lim N→∞ |R i N | = 0. Under these conditions we proved that thus finally proving [51] Theorem I .

Validity of the Generalized Born Approximation
It is worth at this point to provide the expression for the operator L 1 , starting from Equation (8), under the hypothesis of local isotropy in the orthonormal basis defined above. By performing appropriate inner product on the tensorR(q), it is obtained [27], where c(q) is the scalar covariance function of the elastic constants fluctuations, real and positive-defined [46]; L kkii (q q ′ ) is a function of the angleq q ′ between the two versorsq andq ′ , resulting from the inner product [23,27] also accounting for the transverse degeneracy. The assumption of isotropy allows this function to depend only on the angleq q ′ . By making use of spherical coordinates we finally achieve with x = cos(q q ′ ). The function˜ kk (q, ω) is implicitly defined in Equation (30). The validity of the Generalized Born Approximation follows from the following two corollaries of Theorem I.
The function q ′2 c(q, q ′ , x) for x ∈ [−1, 1] and q ′ ∈ [0, q i Max ] is continuous and bounded. Corollary I thus follows immediately from Theorem I.
We recast the function˜ kk (q, ω) in Equation (30) as +R k (q, ω, ǫ 2 ), (32) where the remainder function,R k (q, ω, ǫ 2 ), is . (33) Corollary II. For those values of ǫ 2 , q and ω: |R k (q, ω, ǫ 2 )| ≪ 1, it is˜ The generic term of the series, F n k (q, ω), is implicitly defined in Equation (34). Corollary II follows from Corollary I. We emphasize, furthermore, that the validity of of Corollary II is constraint to the assumption of negligible contribution of the self-energy to the average Green's dyadic if ǫ 2 Corollary II is, however, still valid when this hypothesis is violated but it is possible to show that |R k (q, ω, ǫ 2 )| ≪

The Case of an Exponential Decay of the Covariance Function
We analyze in detail the case of a covariance function cast in the form of an exponential decay function, finding the domain of validity of the GBA in the (ω, q) plane. We furthermore consider only spatial fluctuations of the shear modulus. The results, however, can be easily generalized by including also spatial fluctuations of the Lamé parameter [27].
In a three-dimensional Fourier space the covariance function in this case reads as where d 3 q c(q) = 1. Furthermore it is ǫ 2 =ǫ 2 µ 2 0 = δµ 2 , being δµ the intensity of the spatial fluctuations of the shear modulus per density. From Equation (35) we obtain that in Equation (30) it is c(q, q ′ , x) = a π 2 (aq ′ ) 2 (1+(aq ′ ) 2 +(aq) 2 −2(aq ′ )(aq)x) 2 . We first verify in the following that the hypotheses of Theorem I are verified when the covariance function is given by Equation (35). The validity of Corollary I follows immediately from the continuity of the function in Equation (35). It is finally possible to find a domain of the (ω, q) plane where Corollary II holds. In this domain the GBA can be exploited. It covers a q-range up to ∼ a −1 .
We show that in the case of an exponential decay of the covariance function , as required by the hypotheses of Theorem I, for q, ω ∈ R, i) | ˜ 1 ii (q, ω)| is continuous and bounded for can in such a case be calculated by exploiting the Sokhotski-Plemelj theorem and the Cauchy's Residue Theorem [23,27], finding 4 x 4 (−(aq 0j ) 2 + (aq) 2 x 2 ) 3 +ã 8 (5(aq 0j ) 2 + 6(aq) 2 x 2 ) + a 6 (7(aq 0j ) 4 + 17(aq) 2 (aq 0j ) 2 x 2 +4(aq) 4 x 4 ) +ã 4 (3(aq 0j ) 6 + 16(aq) 2 (aq 0j ) 4 x 2 Frontiers in Physics | www.frontiersin.org 8 October 2018 | Volume 6 | Article 108 +13(aq) 4 (aq 0j ) 2 x 4 + 16(aq) 6 x 6 ) +ã 2 (−3(aq) 2 (aq 0j ) 6 x 2 −(aq) 4 (aq 0j ) 4 x 4 − 5(aq) 6 (aq 0j ) 2 x 6 + 9(aq) 8 We define in the following the domain of the (ω, q) plane where Corollary II holds. Specifics of the mathematical passages are outlined in Supplementary Note 1. We first specify the domain of convergence a.e. in the (ω, q) plane of the series representation of < G ii (q, ω) > 1 . We then show that the magnitude of the remainder function, |R k (q, ω, ǫ 2 )|, is as small as required by Corollary II when aq i Max is large enough. We finally note that such a condition corresponds to deal with small values of ǫ 2 As we infer from Equation (37) and observe in Figure 3, ii (q, ω)|: i) definitively and independently from ω increases by increasing q with a q 2 leading term ; ii) it has a local maximum in q ∼ q 0i . The value of this maximum increases by increasing q 0i , as it is possible to observe in Figure 3, Panels 2. It follows from Theorem I that the condition ǫ 2 discriminate the values of q and ω where the series representation of < G ii (q, ω) > 1 is a.e. convergent. Given the properties of | ˜ 1 ii (q, ω)| stated in points i) and ii) above we infer that this inequality is fulfilled for values of ω and q ∈ [0, q i Max ]: a) ii (q 0i , ω)| increases by increasing q 0i with a q 3 0i leading term, see Equation 37. The frequency values where the condition a) is satisfied are thus q 0i = ω c 0 i ≪ q i Max , as it is possible to infer also by the observation of Figure 4. We then fix a value of frequency where the conditions a) and b) are satisfied and observe that for q ≫ q 0i , i.e., q ∼ q i Max , Im[ ˜ 1 ii (q, ω)] ≪ 1. This is verified in Supplementary Note 1. Consequently, the function has a local maximum at At larger wavevectors this function monotonically decreases by increasing q until to be lower of 1 |q 2 0i −q 2 | . This behavior can be observed in Figure 4, Panels 3. The trends outlined permit finally to asses that for q 0i ≪ q i Max and q ≪ Min {i} [q i Max ] and for values sufficiently large The mathematical passages are shown in detail in Supplementary Note 1. It is thus |R k (q, ω, ǫ 2 )| ≪ 1 when aq i Max is large enough. In this case Corollary II holds and we can exploit the GBA. We infer from point i) and observe in Figure 3, Panels 2, that as smaller it is ǫ 2 c 2 i as larger is aq i Max . It furthermore follows that for a given value of ǫ 2 the larger it isc i , the larger aq i Max . The magnitude of the remainder function remain thus uniquely linked to the . Finally, noting the inverse proportionality between aq i Max and the magnitude of the remainder function, we expect that the domain of the (ω, q) plane where the GBA can be applied includes values of q : aq ∼ 1. A specific case is treated in Figures 4-6, where the quantities shown have been obtained for theory's input parameters which allow to describe the features of longitudinal dynamics for a real system [41]. They are c 0 L = 2.29 meV/nm −1 , c 0 T /c 0 L = 0.53,ǫ 2 = ǫ 2 µ 2 0 = 0.4 and a = 1.1 nm. Figures 4-5 we find aq L Max = 18, aq T Max = 8. For such theory's input parameters the GBA hence holds for q : aq ≪ 8. The mixing of polarization is expected for aq ∼ 1. In this case the proposed approximation can thus be used to describe such a phenomenon.

By inspection of
For the input parameters specified above we calculate the longitudinal and transverse self-energies in the GBA. In sections 3.3 and 3.4 we analyse the features of the acoustic dynamics focusing, in particular, on the mixing of polarizations. To this aim we truncate the series in Equation (34) to the order n = 1, obtaining kk (q, ω) ∼ F 0 k (q, ω) + F 1 k (q, ω). In Supplementary Note 2, we numerically retrieve the value of |F 1 k (q, ω)| and |R k (q, ω, ǫ 2 )| for selected values of wavevector and frequency and show that the former is significantly larger than the latter. We notice that the series in Equation (34) is obtained as the integral of a power series a.e. convergent. It is thus expected that the leading-order terms will be the ones with smaller n. Such an order of approximation is furthermore sufficient to obtain a realistic description of FIGURE 5 | Projection on the (ω, q) plane of the acoustic dynamics, including the Rayleigh anomalies and the mixing of polarizations [41]. In order to facilitate such calculation we furthermore assume ˜ 1 ii (q ′ , ω) ∼ ˜ 1 ii (0, ω). When the approximation ˜ 1 ii (q ′ , ω) ∼ ˜ 1 ii (0, ω) applies, the Hadamard principal value of the integral defining F 1 k (q, ω) can be obtained straightforwardly by exploiting the Residue Theorem. In such a case we can indeed extend the upper integration boundary of the integral to infinity while maintaining unaffected the order of magnitude of the error related to the GBA, as discussed in Supplementary Note 3. It is furthermore shortly discussed in Supplementary Note 4 that as long as the | < 1 2 is fulfilled the dominant contribution to the integral defing F 1 k can be obtained trough the approximation ˜ 1 ii (q ′ , ω) ∼ ˜ 1 ii (0, ω). Figure 5 shows | | for an exponential decay of the covariance function and for the given input parameters. This condition is fulfilled up to frequencies and wavevectors aq L(T) (q 0L(T) ) ∼ 2.
Since the shape of the covariance function makes that the larger contribution to the integral defining F 1 k is for q ′ ∼ q ± a −1 , the approximation ˜ 1 ii (q ′ , ω) ∼ ˜ 1 ii (0, ω) is assumed to give a significant estimation of the integral up to wavevectors of the order of a −1 , where we aim to focus in the present study.  Figure 6, Panels 1, shows the longitudinal and transverse currents (black and red bold lines respectively) obtained by exploiting the GBA for two different values of wavevector: a), aq ≪ 1 and, b), aq ∼ 1. The maximum of the current is normalized to one. The current C L(T) (q, ω) is obtained from the dynamic structure factor S L(T) (q, ω), being C L(T) (q, ω) = ω 2 q 2 S L(T) (q, ω). The dynamic structure factors are related to the average Green functions through the fluctuation-dissipation theorem,

The Mixing of Polarizations in the Born and Generalized Born Approximations
The longitudinal and transverse self-energies defining the average Green functions are obtained from the GBA as described in section 3.2. The longitudinal and transverse currents calculated by exploiting the Born Approximation (dot-dashed lines) are furthermore shown together with the currents of a "bare" medium with phase velocities equal to the average first-order perturbed phase velocities,c L(T) . They are shown respectively as black and red straight lines and referred to be C 0 L(T) (q, ω). The main peak in the longitudinal currents obtained by the GBA points the inelastic excitation centered at the characteristic frequency determined by the longitudinal phase velocity. For the only case aq ∼ 1 it is furthermore observed a low-frequency feature, which can be described as a secondary peak centered at frequencies characteristic of the transverse excitations. This kind of feature has been observed experimentally and by Molecular Dynamics (MD) simulations in several topologically disordered systems [16-18, 20, 52-56]. In the current spectrum related to the Born Approximation we observe a shoulder-like feature at the same frequency, but such a feature is clearly more pronounced when the GBA is used. The endorsement of the fact that the secondary peak is related to the mixing of polarizations comes from the fact that it disappears when the cross term accounting for the coupling with transverse dynamics, LT (q, ω), is removed from the longitudinal self-energy. This is emphasized in Figure 6, Panels 2, where they have been shown the currents obtained by using respectively the full expression of the self-energy, L (q, ω), (full line) and the only longitudinal contribution to the selfenergy, LL (q, ω), (dot-dashed line). Both in the case of the Born Approximation and of the GBA the features observed at the characteristic frequencies of the transverse excitations, which are present when the full expression of the self-energy is taken under account, disappear when the only term related to the longitudinal contribution to the self-energy is considered. This is not the case when we are dealing with the transverse dynamics, as it is possible to infer by observing Figure 6, Panels 3. In this case the secondary peak observed at frequencies higher than the one defined by the transverse phase velocity is in part preserved when the longitudinal contribution to the selfenergy is left out. The occurrence of the secondary peak can be related to the existence of a two-modes regime, observed in a random media whose covariance function can be described by an exponential decay function [23] for wavevectors and frequencies aq(q 0i ) > 1. This behavior can be reproduced also by using the scalar Born Approximation [23,25]. Similar feature is furthermore observable in the longitudinal dynamics for wavevectors higher than the ones considered in Figure 6. It can coexist with the feature related to the mixing of polarization.

Features of the Acoustic Dynamics in the Generalized Born Approximation
The features of the longitudinal acoustic dynamics obtained by GBA up to wavevector of the order of a −1 have been derived from the calculated longitudinal currents by a fitting procedure described in the following and qualitatively compared with experimental finding reported in the literature. The dynamic structure factors related to longitudinal acoustic dynamics in topologically disordered systems have been experimentally characterized by several studies mostly based on IXS or INS measurements in different wavevectors regions. An universal behavior emerged, which can be qualitatively described as i) the presence of the so-called Rayleigh anomalies for values of wavevctors much lower than a −1 , i.e. the phase velocity of the acoustic excitations shows a softening with respect to its macroscopic value while the acoustic mode attenuation is affected by a strong increase and follows roughly a q 4 trend; ii) increase of the phase velocity at higher wavevector values, resulting in a minimum in its wavevector-trend, and crossover from q 4 to q 2 trend of the acoustic modes attenuation [1][2][3][4][5]; (iii) mixing of polarizations for q ∼ a −1 manifesting in the presence of a peak-like feature in the longitudinal current spectra at the characteristic frequencies of transverse excitations [16-18, 20-22, 52]. The GBA can grasp all these characteristics. While the Rayleigh anomalies can be obtained also by exploiting the Born Approximation [8,23] or the SCBA [1,7,11], the mixing of polarizations can be accounted only by the GBA. It is worth, furthermore to observe that at the boundary of the Rayleigh region when depolarization effects begin to affect the acoustic dynamics, the coupling between longitudinal and transverse dynamics, though not manifesting in a clear peaklike feature in the longitudinal currents, can have an impact on the effective experimentally observed attenuation and phase velocity. To obtain a realistic description of the acoustic dynamics also in this wavevectors region it is thus more appropriate to consider a vectorial model in the RMT frame, such as the GBA. The features of the longitudinal acoustic dynamics can be derived from the longitudinal currents obtained by the GBA by fitting the calculated spectra with a fitting model composed by one or two Damped Harmonic Oscillator (DHO) functions, following the same protocol usually used to analyze the experimental IXS or INS data. This approach, also referred to be as spectral function approach, has also been used in the analysis of theoretical results aimed to characterize the acoustic dynamics in random media [23,35]. Because most of the experimental data presented in literature have been analyzed with the above quoted fitting model, the spectral function approach permits a clear connection between the GBA theoretical outputs and the literature results.
where n = 1 in the Rayleigh region q << a −1 and n = 1, 2 in the region q ∼ a −1 . In the transition region, where the feature related to the tranverse dynamics starts to show up, the fitting has been performed by exploiting both 1-or 2-DHO model fit functions. The phase velocity, c(q) is related to the characteristic frequency, , trough the relationships c(q) = (q) q , while the parameter Ŵ is directly related to the acoustic mode attenuation. The wavevector-dependent outcomes of the fitting procedure are diplayed in Figure 7. The theory's input parameters are the same listed above. Both the Rayleigh anomalies and the mixing of polarizations are clearly observed in qualitative agreement with most of experimental outcomes reported in literature. A quantitative comparison with experimental outcomes is reported in Izzo et al. [41].

CONCLUSION
By introducing corrective terms to the Born Approximation we obtained in an analytic form an expression for the self-energy related to the stochastic Helmholtz equation describing the acoustic dynamics in an elestically heterogeneous medium. In the frame of the perturbative series expansion of the Dyson equation the proposed approximation accounts in an approximate form up to the second order term, whereas the Born Approximation stops to the first order. The Feynman diagram technique permits to clarify which multiple scattering events are included in the Generalized Born Approximation. The case of a covariance function given by an exponential decay function is analysed in some detail. In such a case it was proved the validity of the proposed approximation in a domain of the (ω, q) plane of interest in most topologically disordered systems (e.g., glasses). It includes both the Rayleigh region and wavevectors region: aq ∼ 1, where it is expected the mixing of polarizations to get in. Furthermore, the validity of the Generalized Born Approximation is not restricted to q in a neighbor of ω c i , where c i is the phase velocity of the unperturbed medium for the i-th polarization, thus permitting to describe also features of the average Green's dyadic occurring at frequencies smaller or higher than c i q, as it is the case for the mixing of polarizations. We finally verified that the proposed approximation permits to describe this phenomenon together with the Rayleigh anomalies.
Acoustic modes with mixed polarization have been observed by both IXS and INS as well as by MD simulations in several disordered systems. The phenomenon has never been related, however, to Rayleigh anomalies and quantitatively described as a phenomenon also originating from the disordered nature of the medium. The proposed approximation can permit to reach this goal and to trace the way toward a coherent and experimentally verifiable mathematical description of all the phenomena arising from the elastic heterogeneous structure of an amorphous solid.

AUTHOR CONTRIBUTIONS
MI led research, produced the results, and wrote the paper. GR and SC discussed the results and revisited the paper.