Multiple scattering theory for heterogeneous elastic continua with strong property fluctuation: theoretical fundamentals and applications

In this work, the author developed a multiple scattering model for heterogeneous elastic continua with strong property fluctuation and obtained the exact solution to the dispersion equation derived from the Dyson equation under the first-order smoothing approximation. The model establishes accurate quantitative relation between the microstructural properties and the coherent wave propagation parameters and can be used for characterization or inversion of microstructures. As applications of the new model, dispersion and attenuation curves for coherent waves in the Earth lithosphere, the porous and two-phase alloys, and human cortical bone are calculated. Detailed analysis shows the model can capture the major dispersion and attenuation characteristics, such as the longitudinal and transverse wave Q-factors and their ratios, existence of two propagation modes, anomalous negative dispersion, nonlinear attenuation-frequency relation, and even the disappearance of coherent waves. Additionally, it helps gain new insights into a series of longstanding problems, such as the dominant mechanism of seismic attenuation and the existence of the Mohorovicic discontinuity. This work provides a general and accurate theoretical framework for quantitative characterization of microstructures in a broad spectrum of heterogeneous materials and it is anticipated to have vital applications in seismology, ultrasonic nondestructive evaluation and biomedical ultrasound.


exact solution. Solving for the exact solution to the dispersion equation derived from the FOSA Dyson's equation is still very challenging.
Recently the author conducted comprehensive study on Weaver's model and the spectral function approach and obtained the exact solution to the dispersion equations [129]. By comparing the accurate solution with the results obtained using the spectral function approach, it is recovered that the spectral function method is valid at low and high frequencies only, and in the stochastic transition region it gives incorrect results. Most importantly, the performance of Weaver's model for two-phase materials is quite unstable, in certain cases it gives physically impractical results. Both Weaver's model and Stanke-Kino's model use the Voigt-average material as the homogeneous reference medium. Turner pointed out the Voigt-averaged velocities always overestimate the velocities that measured in experiments [49]. To remedy this discrepancy, he introduced a self-consistent scheme to calculate the properties of the homogeneous reference medium, and incorporated these material properties into Weaver's model [50]. Through comparison of the results obtained by different homogenization schemes, such as the Voigt, Reuss, Hill and self-consistent techniques, he concluded that the self-consistent method significantly improves the predictions of the dispersion and attenuation. It is worth mentioning that all these scattering theories are valid only for weak scattering media because the small-perturbation expansion is adopted for the description of property perturbation.
A number of multiple scattering theories are also developed based on the discrete-scatterer model. One such theory is known as the generalized coherent potential approximation (GCPA) [110][111][112]. Based on the observation that waves propagating in the homogeneous effective medium should have no scattering at low frequencies and have very little scattering in the high frequency range, Sheng [110] proposed that the dispersion equation can be obtained by enforcing the mass operator vanish at low frequency and assuming its minimum in the high frequency range. Further approximations are introduced in order to make the resulting equations solvable, such as the Tmatrix approach. The dispersion curve for liquid suspensions with monodispersed methylmethacrylate spheres are obtained and it shows that in the intermediate to high frequency range, there exist two propagation modes, which was observed experimentally by [113]. Sheng and coworkers developed the theory for liquid suspensions and electromagnetic waves only while the counterpart theory for elastic materials is still missing. Foldy [114] and Twersky [115][116] developed another multiple scattering theory based on the discrete scatterer model. In their model, the complex wavenumber is expressed explicitly using the far-field scattering amplitude of a single inclusion.
The results are applicable when the volume concentration of the scatterers is low. Quasi-crystal approximation [60][61][117][118] is another discrete scattering model applicable for heterogeneous media in which inclusions are located on a regular lattice. The range of validation of different versions of scattering theories are summarized in Tab. 1. WPFM=Weak-Property-Fluctuation-Medium, SPFM=Strong-Property-Fluctuation-Medium a Weaver's model gives reasonable predictions for the velocity and attenuation of polycrystalline alloys in the whole frequency range [129], but it gives unstable predictions for two-phase materials, as shown in this work. b Stanke-Kino's model has been shown to be equivalent to Weaver's model. For polycrystalline alloys, it gives nearly the same results as that given by Weaver's model [50]. However, numerical results for two-phase materials calculated using this model are unavailable.
From Tab. 1 it is seen that the adoption of different types of approximation puts severe limitations on the scope of applications of the existing models. Although each model achieves a certain degree of success, a general multiple scattering theory for heterogeneous 6 materials with strong property fluctuations and valid in the whole frequency range is still missing. As a consequence, the exact dispersion and attenuation behavior of coherent waves have not been obtained yet. Moreover, there are also a number of fundamental problems, such as the dominant mechanism of seismic wave attenuation [9], the existence of the Mohorovičić discontinuity [130][131][132], applicability of the Kramers-Kronig relation to multiple scattered elastic waves [65][66][133][134], and anomalous negative dispersion of cortical bone [70,73] are full of controversy in the scientific community. Intrigued by the limitations of the existing models, in the present work the author aims at developing a most universal theoretical framework that applicable to a large variety of heterogeneous elastic materials, regardless of weak or strong property fluctuation, with low or high volumetric concentration of inclusions, excited by low or high frequency signals, while giving up all the conventional approximations. Based on the exact solutions to the new model, this work will also provide a series of completely new explanations to the longstanding problems. In addition, it will be demonstrated that the new model is a real prediction model that can be used for microstructure inversion and seismic/ultrasonic data interpretation.
The contents of subsequent sections are arranged as follows: Section II presents the rigorous development of the new model,

II. THEORETICAL FUNDAMENTALS
In this section, we present rigorously the theoretical development of the new multiple scattering theory for strong-property-fluctuation elastic continua. It needs to mention that mathematics plays a central role in developing the new multiple scattering model. A series of advanced mathematical techniques, including tensor analysis, system of integral equations, singularity techniques, Feynman's diagram, statistical theory of heterogeneous medium, Fourier transform, theory of complex variables and advanced numerical techniques provide us a precious opportunity to develop an accurate and universal theory while introducing as less approximations as possible.
From the theory of elastodynamics we know, the local wave motion in an inhomogeneous media is governed by a set of partial differential equations (PDEs), known as the Navier-Lamé equations [135]. Boundary conditions are specified by the continuity of stress and displacement across the phase boundaries. Thus, the multiple scattering problem is essentially a sophisticated PDE boundary value problem (BVP). From a mathematical point of view, the integral equation (IE) description and the PDE BVP description comprise two equivalent and complementary representations for the same problem. While the PDE BVP description is frequently utilized to solve scattering problems of a single or a finite number of scatterers with regular geometry [135][136], the IE description has proven to be more convenient for dealing with scattering problems in inhomogeneous media with multiple irregular inclusions or even randomly distributed scatterers [48,98,[137][138]. Transforming the multiple scattering problem from the PDE BVP description into the IE description is the first key step for the development of the new theory. Green's function plays a vital role in this transformation procedure.
As a fundamental solution to the elastodynamic equation, Green's function contains the complete information of the whole system and can be used to analyze more complicated source distributions or time varying excitations. When transferred into the frequencywavenumber domain, Green's function can also be used to extract dispersion behaviors of plane wave components propagating along different directions. Due to its extreme significance for the theoretical development of this work, we first perform comprehensive study on Green's functions of heterogeneous media. 7 Before embarking on the analysis of wave scattering in a general heterogeneous medium, we first consider the most fundamental problem: scattering of elastic waves by a single inclusion in an infinite homogeneous elastic medium.

Green's function of a homogeneous medium with a single inclusion
The time-harmonic wave propagation in a heterogeneous medium is governed by the classic elastodynamics equation [135], given by: (1) where ω is the circular frequency, ui denotes the displacement components, ρ(x) and cijkl(x) are the mass density and elastic stiffness.
The elastic stiffness has the following symmetries: Throughout this work the Cartesian tensor notation is used. A bold-faced letter represents a vector, tensor or matrix, and italic letters with subscript indices represent tensor components or matrix elements. A comma followed by a coordinate index means taking partial derivative with respect to the corresponding spatial coordinate. The Einstein summation convention, i.e., a repeated index implies summation over that index from 1 to 3, is assumed in this work. Consider an infinite homogeneous elastic medium in which an inclusion with different density and elastic stiffness is embedded, as shown in Fig. 1(a). The quantities pertaining to the homogeneous medium and that to the inclusion are discriminated by an attached index "(0)" or "(1)", respectively.
Green's function is defined as the resulting field excited by a time-harmonic unit concentrated force F applied at a generic point Oʹ along eaʹ, where eaʹ represents the coordinate basis of a source-region coordinated system with its origin located at Oʹ. The coordinate basis of the coordinate system Oʹxʹyʹzʹ by no means needs to be the same as that of the coordinate system Oxyz, so at this point we assume they are different from each other. The portions of the field in the substrate and in the inclusion are governed by two different sets of equations, which are given by: (1) (1) 2 , where V0 and V1 represent the volumes occupied by the substrate and the inclusion, respectively, and aiaʹ denotes the directional cosine of eaʹ, i.e., cos( , ) In order to obtain the integral representation of the perturbed Green's function, we now consider the corresponding homogeneous Green's function, which is defined as the field induced in the same medium but in the absence of the inclusion by a unit concentrated force F applied at another point Oʺ along eαʺ, as shown in Fig. 1(b). eαʺ is the coordinate basis of a new coordinate system Oʺxʺyʺzʺ which is used to describe the source distribution in the homogeneous medium. This function satisfies:  (6) Multiplying Eq. (6) and Eq. (3) by ia G  and (0) i G   , respectively, and then subtracting the resulting equations, we obtain:  (7) Multiplying Eq. (6) and Eq. (4) by ia G  and (0) i G   , respectively, and then subtracting the resulting equations, we obtain: (8) Appling product rule of partial derivatives to Eqs. (7) and (8), we get: ijkl k l ia j ijkl k l ia j ijkl ka l i j ijkl ka l i j (10) Considering the symmetry of the elastic stiffness, Eqs. (9) and (10) can be simplified as: Integrating Eq. (11) over the domain V0, and then applying Gauss's divergence theorem, we obtain: A negative sign is introduced in front of Nj in the integral over S1 because the out normal of volume V0 corresponds to the negative of N, which is shown in Fig. 1. S  represents an imaginary surface that lies infinitely far from the source and the inclusion. Integrating Eq. (12) over the volume V1 and applying Gauss's divergence theorem, we obtain: Adding Eqs. (13) and (14) and then invoking the boundary conditions Eq. (5), we have: (15) Throughout this work we assume that all the materials have a small damping, so both the homogeneous and the inhomogeneous Green's functions decay exponentially as the distance between the source point and the field point increases. Consequently, these functions tend to zero at infinity and the integral over S  vanishes. Furthermore, we assume that the coordinate basis of the coordinate systems 1 Finally, we obtain the integral representation of the inhomogeneous Green function expressed in terms of the reference homogeneous Green's function, To further simplify the expression of Eq. (17), we decompose the density and elastic stiffness of the inclusion into two parts, one corresponds to the ensemble averaged value, for this case it equals to the properties of the substrate, and the other corresponds to the property fluctuations, Substitution of (18) into (17) yields: where Green's function of the homogeneous medium is given by: For a detailed derivation of this expression, readers are referred to the author's monograph [129]. It is seen from Eq. (20) that the homogeneous Green's function has the following symmetry: Applying these relations to Eq. (19), we obtain: Eq. (22) is the major result in this section. It establishes the relationship between the perturbed Green's function and Green's function of the homogeneous medium. From a mathematical point of view, Eq. (22) is a system of Fredholm equations of the second kind, which can be solved numerically by introducing certain discretization schemes. At this point it is meaningful to give a physical explanation of Eq. (22). It tells us the perturbed Green's function can be expressed as a summation of the unperturbed Green's function and a perturbation term, and the perturbation term is given by a spatial convolution of the perturbed and unperturbed Green's function weighted by the property fluctuation. When the fluctuation of the material properties is weak, the perturbation of the inhomogeneous Green's function is small compared to the homogeneous one. Based on this observation, M. Born suggested to replace the perturbed Green's function appeared in the kernel by the unperturbed Green's function, and thus an explicit expression for the perturbed function is obtained.
This approximation is known as the Born approximation. Here we do not use this approximation and proceed in an accurate manner.

Green's function for a heterogeneous medium with multiple inclusions
Now we consider a medium with multiple scatterers. A schematic diagram is shown in Fig. 2.

FIG. 2. Schematic diagram of an inhomogeneous medium with multiple inclusions.
Obviously, the above derivation can be straightforwardly extended to inhomogeneous medium with multiple inclusions, and the resulting equation is: One needs to keep in mind that all the boundary conditions are naturally incorporated into this equation. From Eq. (23) it is seen that the perturbed Green's function is expressed as the summation of the unperturbed Green's function and a perturbation term which includes contribution of all the inclusions. Meanwhile, the internal field of each inclusion is dependent on all other inclusions and the unperturbed Green's function. At this point we need to point out that by neglecting the interactions among different scatterers, we obtain the independent scattering approximation (ISA) of the perturbed field [38][39][40][41]. ISA is frequently used in combination with the Born approximation and yields an explicit expression for Gβʺaʹ. If the property fluctuation is weak or the volume concentration of the scatterers is sufficiently low, each scatterer approximately interacts with the incident wave independently and ISA can give reasonable predictions.
However, if these conditions are not satisfied, i.e., material exhibits strong property fluctuation or contains notable volume fraction of inclusions, the interactions among all the scatterers cannot be ignored and the multiple scattering theory is necessary to fully capture the propagation behavior.
For a general heterogeneous material in which the inclusions are densely distributed in the whole space, the material property fluctuations are most conveniently represented by functions of spatial coordinates, and Eq. (23) is rewritten as: Here comes a peculiar problem. Since the inclusions are densely distributed, and may occupy a finite percentage of the total volume of the whole medium, which material should be chosen as the reference medium? The answer to this question will be clear after we find out the singularity of the Green tensor and introduce the renormalization scheme in Section II.C. Contrary to our intuition, the proper choice of the reference material is neither either phase of the constituent materials nor the volumetric average of the component phases.
From Eq. (39) it is seen that the following integral is involved: where ()   Γ xx is composed of homogeneous Green's function and its spatial derivatives up to second order. Meanwhile, we see that the domain of the integration is the whole space of the heterogeneous medium. A question naturally arises here: Since Green's functions and their derivatives are not well defined at   xx , how can we calculate these integrals? As one can see, Green's functions and their derivatives become infinite when x approaches xʺ, this property is called the singularity of Green's functions. Obviously, Green's function and its derivatives of different orders tend to infinity with different speed as x approaches xʺ, this means that they have different degrees of singularity. To properly define and calculate these integrals, we need to introduce the concept of shape-dependent principal value, which will be detailed in the next section.

B. Singularity of Green's tensor and the renormalized integral equation
The singularity of electromagnetic (EM) Green's function was first discovered by Bladel [139]. He found that the field distribution in the source region cannot be calculated using the conventional dyadic Green's function, instead, a term proportional to Dirac-delta function must be subtracted in order to obtain the correct result. Finkel'berg [140] incorporated the singularity of EM Green's functions in the development of an effective medium theory for dielectric mixtures. Kong and Tsang [99][100][101][102][103][104][105][106][107][108][109] further calculated the singularities of EM Green's function of heterogeneous dielectrics with different spatial correlation functions, including the Gaussian, exponential and Von-Karman correlation functions, and derived effective dielectric constants for inhomogeneous media with spherical and ellipsoidal inclusions. They also derived the explicit expressions of backscattering cross section and applied them in satellite remote sensing. Lakhtakia and coworkers extended Kong and Tsang's theory to EM anisotropic media, and developed effective medium theory correspondingly [141][142]. Analogous to the development of EM effective medium theory, Zhunk analyzed the singularity of Green's function of the acoustic waves [139][140] and the elastodynamic Green's tensor [96] and developed an effective medium theory for heterogeneous acoustic and elastic materials. It has been pointed out that any multiple scattering model without considering the singularity of the Green tensor is only applicable to materials with small property fluctuation. In this work, the singularity of the ij E       xx has the δ-singularity which needs to be carefully dealt with. By introducing the concept of shape-dependent principal value, the correct definition of the matrix D is given by: where P.S. D is the shape-dependent principal value of the Green tensor, S is the singularity of the Green tensor, which is also dependent on the shape of the inclusions. For instance, the singularities of heterogeneous media with spherical, spheroidal, ellipsoidal and cylindrical inclusions all have distinct expressions and possess different symmetries. In the current work, we only consider random medium with spherical inclusions, for which the macroscopic properties are isotropic. Correspondingly, P.S. D and S have the following form:   , thus Sijkl is an isotropic tensor. This conclusion is consistent with the assumption that the heterogeneous medium is macroscopically isotropic.
In analogous to the isotropic elastic stiffness tensor, we can rewrite Sijkl as 12 ( ), Similarly, we can give the correct definition of the matrix ()   Γ xx by introducing its shape-dependent principal value . .
Substitution of (51) into (39) yields: Invoking the definition of the Dirac-δ function, we get 03 () Introducing a new quantity, called the renormalized field variable, (55) the integral equation (54) is rewritten as: In this work, we consider heterogeneous media whose component phases are all isotropic materials, for which the matrix Ξ(x) and its elements are given by: One can easily verify the relation: 11 12 44 2.
     Eq. (56) is known as the renormalized integral equation [141][142]. In contrast to the development of effective medium models, in the current work we strive to develop a multiple scattering theory, which should be able to describe the coherent wave propagation in the whole frequency range. To do so, we first perform multiple scattering series expansion of Eq. (56), and then derive the Dyson equation with the help of Feynman's diagram and the first-order smoothing approximation.

C. Feynman's diagram, Dyson's equation and the First-Order-Smoothing-Approximation
Eq. (56) is an integral equation for the unknown quantity ()    Φ xx . We can see that this quantity appears on both sides of the equation and this naturally forms an iteration scheme. By performing iteration successively, i.e., substituting the righthand side of Eq. (56) into the renormalized quantity in the integrand, we get an infinite series Eq. (60). In Eq. (60), the renormalized field is expressed using the unperturbed field variables only. In this work, Eq. (60) is called the multiple scattering series expansion of the perturbed field.
The coherent field, also called the ensemble-averaged field, is defined by averaging the perturbed field over the whole space of realizations of the random medium. Taking ensemble average of Eq. (60) gives: where the angular bracket means the ensemble average of the enclosed quantity. According to the original definition, we need to calculate the random field variables in the whole configuration space, which includes an infinite number of realizations of the random medium.
According to the ergodic hypothesis, the ensemble average of a quantity is equal to the volumetric average of the perturbed quantity, thus the angular bracket can also be understood as taking volumetric average of the enclosed quantity. Since . . ( ) PS   Γ xy is a deterministic quantity, we can take it out from the bracket and get: One need to be informed that in Eq. (62), the subscripts expressed using Latin letters are matrix indices instead of tensor indices.
In Eq. (62), an infinite number of multi-point correlation functions are involved: Multi-point correlation functions, also called higher-order moments, provide key information about the statistical characteristic of heterogeneous media, and they have found broad applications in different subjects, such as characterization of galaxy distribution in 15 cosmology [145], describing random EM media for satellite remote sensing [101,109,122], analyzing small-scale heterogeneities in the Earth's lithosphere [3,29,146], and characterizing microstructures in metal polycrystals [46][47][48]. For a detailed discussion on how to calculate higher-order moments readers are referred to [47,107,[147][148]. No matter what type of statistics the inhomogeneous medium complies with, one can always choose a reference homogeneous medium such that , so we only need to consider the second and higher order moments. Generally speaking, multi-point moments have no simple relations with lower order moments.
However, if the random medium following the Gaussian statistics, i.e., the two-point correlation function can be expressed as a Gaussian function, all the moments of odd order vanish [123,127]: and all the moments of even order can be expressed as products of the second order moments To proceed further, it is convenient to introduce Feynman's diagram technique, which was first developed by Richard Feynman to investigate interactions of elementary particles with quantum many body systems. Rytov adopted this method to study multiple scattering of classic scalar waves [123]. In this work, we extend Rytov's implementation of Feynman's diagram technique to investigate the multiple scattering of elastic waves. The symbolic representations of all the involved quantities are shown in Tab. 2. Table 2. Symbolic representation of field variables used in Feynman's diagram.
Reference field: In addition to Tab. 2, we also adopt the convention that two scattering points lie in the same inclusion, i.e., located in the same bracket, are linked by a dashed line. To proceed further, we have to introduce the concepts of weakly-and strongly-connected diagrams. If a diagram can be divided into two subdiagrams, each of which has two or more scattering points, without breaking any dashed lines, then the diagram is call ed a weakly-connected diagram. Diagrams do not possess this property are called strongly-connected diagrams. For instance, the third term on the righthand side of the equation is a weakly-connected diagram, the second, fourth and fifth terms are strongly-connected diagrams.
For simplicity, we introduce an additional symbol to represent the sum of all the strongly-connected diagrams: This approximation is called the first-order smoothing approximation [123,127], also known as the Bourret approximation [153] and bi-local approximation [122]. Collin [154] compared different selective summation techniques and concluded that the FOSA has the fastest convergence rate. We can write down the explicit equation corresponding to the diagram shown in Fig. 7: Eq. (67) is the so-called Dyson's equation under FOSA, in subsequent discussion, we call it Dyson's equation for simplicity. A primary study on the accuracy of this approximation was carried out in [156] in the context of the homogenized continuum model. It is found the difference between the third-order and second-order (corresponding to the FOSA here) estimates of the homogenized constitutive parameters are very small for all values of characteristic size, even when the constitutive contrast between the component materials is as large as 30. Similarly, it shows that higher order (fourth-order, fifth-order …) approximation does not add significantly to the estimation given by the second-order approximation. A rigorous analysis of the accuracy and error estimation in the multiple scattering scenario is out of the scope of this work. It is noted that different from the homogenization theories, no constraint on the frequency is introduced here, so the conclusions drawn above are assumed to hold in the whole frequency range. This conjecture will be examined and verified by numerical results for a variety of materials, see Sections III and IV. Another point we need to stress is that Eq. (67) is valid for random media with any types of two-point correlation functions, such as the exponential and the Von-Karman correlation functions. This is because all the third-and higher-order moments are completely neglected.
For a statically homogeneous medium, the two-point correlation function is given by [47][48]155]: where () P  yx is called the spatial autocorrelation function (SAF). Heterogeneous materials with different microstructural distributions are described by different SAFs. For example, random media with blurred interfaces are most conveniently described by Gaussian-type SAF, inclusions or scatterers with sharp boundaries follow exponential-type SAF, and Von Karman-type are introduced to describe selfsimilar random media like turbulent media [3]. In this work, we only consider random media with exponential correlation functions since it describes the correlation features of most random media discussed in this work to an accurate or acceptable level. If the inclusions in the medium are of spherical geometry, then the random media is statistically isotropic, and the exponential correlation function is a function of the distance between the two points [3,[47][48]: where a is the correlation length.
The Fourier transform of SAF represents the power spectrum of the medium fluctuation, which is given by: As pointed out before, we can always choose the properties of the reference medium such that the first-order moment vanishes. This choice also ensures the fastest convergence rate of the multiple-scattering series, thus we take: For a two-phase heterogeneous medium, Eq. (72) gives three independent relations [12]: 12 0    : where As can be seen, the mass density of the reference medium is still defined by the volumetric average of the component material s, but the elastic stiffness is obtained by solving a system of complicated nonlinear algebraic equations (74) and (75), which is largely different from the Voigt-averaged values. In Sections III and IV we will see that the Voigt-average estimation always overestimates the quasistatic limit of the velocities of the coherent waves. This conclusion is in agreement with that given by other works [49][50].

D. Solution of Dyson's equation and the dispersion equations
Equation (68) is a system of integral equations of convolution type, which is most conveniently solved by the Fourier transformation technique. In this work, we use the following spatial Fourier transform pair: Applying Fourier transform to the Dyson equation (68), and considering (69), we get: Multiplying both sides of Eq. (78) by , we get: Rearranging this equation, we have: The ensemble averaged response in the frequency-wavenumber domain can be solved as: Simultaneously we obtain the dispersion equation: The solution in the spatial-frequency domain is given by: If the source is a time varying signal F(t), and it is correlated to its spectrum through the following temporal Fourier transform pair: Thus, the complete wavefield in spacetime induced by a general time-varying point source is given by: Eq. (85) gives an explicit expression for the wavefield in the spacetime domain and it has significant implications for seismological applications, for instance, it can be used to synthesize three-component seismogram envelops for realistic earthquake sources [3,176].
However, exploring the solution goes beyond the scope of the current research. In this work, we only focus on the dispersion and attenuation behavior of a plane wave component. In a statistically isotropic medium, the dispersion behavior of a plane wave is independent of its propagation direction. Without loss of generality, we consider a plane wave propagating along the x3 axis, and the The coefficient matrix and the dispersion equations of longitudinal and transverse waves can be expressed explicitly as: 11 The non-vanishing elements of the matrix M are given by:    44  11  1 2  11  11  12  12  44  1 2  11  11  12  12  45 (1)  2  45  12  1 2  11  11  12  12  44  1 2  11  11  12  12  45 (1) The solution to the dispersion equations is a complex propagation constant, Re( ) Im( ) , ω is the circular frequency, V is the coherent wave velocity, and α is the attenuation coefficient of the coherent wave. The propagation characteristics of the coherent wave is fully described by the complex wavenumber. . The correct definition of these integrals is discussed in Sheng [112] and Calvet and Margerin [12], where they gave the explicit expression by splitting the homogeneous Green's functions into the real part and the imaginary part. The real part of Green's functions is defined as the Cauchy principal value at the singular points while the imaginary part is given by the Dirac-delta function. The sign of the imaginary part is determined by the causality of elastic waves, see Sheng [112]. In the numerical implementation, these integrals are broken into two or three parts, and the general forms are and F(s) is a regular function, i.e., has no singularity along the positive real axis.
All the integrands in the above integrals decays proportionally to 1/s 2 as s   . Consequently, all the infinite integrals are convergent. The complex wavenumber is obtained by searching for the roots of the dispersion equations in the complex k-plane. The numerical algorithm is implemented on the platform Compaq Visual Fortran 6.6 for which the powerful IMSL numerical library is integrated. All the numerical calculations are carried out in terms of the dimensionless quantities, i.e., the fractional velocity variation, the adimensional attenuation and the dimensionless frequencies, as given by:

III. COMPARISON WITH AN IMPROVED MULTIPLE SCATTERING MODEL FOR HETEROGENEOUS ELASTIC MEDIA WITH WEAK PROPERTY FLUCTUATION
Before applying the new model to a series of practical problems, it is meaningful to make a comprehensive comparison between the new model and an improved multiple scattering model for heterogeneous elastic media with weak property fluctuation, which was recently developed by the author based on Weaver's model [48]. For the convenience of subsequent discussion, we introduce two acronyms for the two models, where SFMS stands for the Strong-Fluctuation-Multiple-Scattering model developed in this work and WFMS stands for the Weak-Fluctuation-Multiple-Scattering model. We omit the rigorous development of the WFMS model and only list the major results here. Interested readers are referred to [129] for more details.
The dispersion equations for the longitudinal (L) and transverse (T) coherent waves are: The mass operators are complex quantities and the real and imaginary parts are: The Voigt average velocities L V and T V are defined in terms of the Voigt average material properties as: f and 2 f are the Lamé constants, density, and volume fraction of the two component materials. These expressions differ from that given in [12] in that all the cross-terms of mass density and elastic moduli have a minus sign instead of a plus sign. It is exactly due to this "small" change that the results are dramatically different from that given in [12]. It is worth mentioning that the improved model is developed by introducing the concept of exterior product, which has been extensively used in Differential Geometry [157]. It provides an ideal tool for dealing with integrals with multiple independent variables. As pointed out by Calvet and Margerin, the spectral function obtained in [12] gives two peaks in the high frequency domain for Media D and E, but there should be only one longitudinal mode for Media D and E and one transverse mode for Medium E. Contrarily, the spectral functions of Media D and E obtained using the new formula gives spectral curves exactly the same as predicted, i.e., one peak for longitudinal waves in Media D and E and one peak for transverse waves in Medium E in the whole frequency range, which perfectly resolves the raised question. For more details please see [129].

A. Instability of the WFMS model in predicting the dispersion behavior of weak-property-fluctuation elastic media
To keep consistent with [12], we use the same properties and the same numbering system for materials A-E, as shown in Tab     show excellent agreement with that obtained using the SFMS model. The longitudinal velocity starts from the quasi-static limit V0L, and then slightly decreases as the frequency increases. Meanwhile, the dimensionless attenuation αLd increases nonlinearly from zero to 0.1.
As the dimensionless frequency k0Ld approaches 50, a second, faster mode begins to appear, but with much larger attenuation and slight positive dispersion. As the frequency increases further, the velocities of the two modes quickly approach the upper and lower limits VL1 and VL2, respectively, and the attenuation coefficients reach a saturation value near unity. We can see the two models give accurate results in the whole frequency range, from the quasi-static region to the geometric domain. Other scattering models based on different approximations are only valid in a limited frequency range. For example, the Born approximation is only valid in a relatively low frequency domain, corresponding to 0< k0Ld < 70 for this case, while geometric approximation is only valid in high frequency range, corresponding to k0Ld > 150 for this case. In addition to dispersion, the new model can also give the accurate attenuation. All these unique features show that the multiple scattering models have incomparable advantages compared to other scattering models. The WFMS model has been widely used in characterization of polycrystalline alloys, for which only the elastic stiffness is a random variable.
Previous calculations show the WFMS model can give stable prediction for dispersion and attenuation of coherent waves in polycrystals.
Here we will show that the WFMS model gives unstable prediction for the dispersion and attenuation of two-phase materials for which both the elastic stiffness and the density are random variables.
Comparison of velocity (a) and attenuation (c) of medium F predicted by SFMS and WFMS models.
(b) is a semi-logarithmic coordinate plot of (a) and (d) is an enlarged view of (c).
The dispersion and attenuation of longitudinal waves of Medium F calculated using the two models are shown in Fig. 9. It is seen from Figs. 9 (a) and (b) that the velocity predicted by the SFMS model first undergoes negative dispersion and then become positive dispersive, Contrary to Medium A, a second, slower mode appears in addition to the original, faster one at intermediate frequencies.
At high frequencies, they approach the upper and lower limits of the two components, which is similar to Medium A. From these results, we can also observe that the dispersion given by the WFMS model have a similar pattern as Medium A, but it gives negative attenuation at dimensionless frequencies below 5.5, while the SFMS model gives positive attenuation in the whole frequency domain. Negative attenuation is physically impractical since it means the wave amplitude increases with propagation distance while no energy i s input into the system. For velocity dispersion, differences in velocity up to 5% are observed between the results obtained by the two models.
This example shows the performance of the WFMS model is unstable, in certain cases it gives wrong predictions.
Next, we examine the propagation characteristics of transverse waves predicted by the two models. The transverse dispersion and attenuation of Medium A are shown in Fig. 10. It is seen the dispersion curves calculated from the two models exhibit excellent agreement. However, the attenuation predicted by the WFMS model is negative in nearly the whole frequency range and form an image of the curves given by the SFMS model symmetric about the axis αTd=0. Fig. 10(c) shows the detailed variation the attenuation in frequency range (0, 1.0). In this range, the WFMS attenuation first increases with frequency. After reaching its positive maximum, it decreases with frequency and then becomes negative. This oscillation behavior prohibits any attempts to remedy the defect simply by reversing the sign of the attenuation. As an additional example, Fig. 11 gives the dispersion and attenuation curves of transverse waves in Medium G. In this case the velocities predicted by the two models agree well at frequencies lower than 50, and the discrepancy becomes large at high frequencies. Moreover, the velocities calculated from the WFMS model do not converge to its geometric limits. Although it gives positive values, the WFMS model underestimates the attenuation of the slow mode and overestimates that of the fast mode in the whole frequency range. From the results, we can see Medium B has one longitudinal mode at low frequencies and two at high frequencies, while only one transverse mode with negligible attenuation can be found in the whole frequency range. Both longitudinal and transverse waves in Medium C have one propagation mode at low frequencies and two at high frequencies. All these predictions are consistent with that given in [12]. Media D and E provide two subtle examples to examine the reliability of different models. For Medium D, as pointed out in [12], the two component materials have the same longitudinal velocity but different transverse velocities, so there should be one longitudinal mode in the whole frequency range while the transverse waves behave normally, i.e., one mode at low frequencies and two at high frequencies. Contrary to this presupposition, the spectral function in [12] gives two longitudinal modes at high frequencies. As discussed before, this is due to a mistake in the expression of the mass operators. Figure 13 shows results drastically different from those in [12], the SFMS model predicts one longitudinal mode in the whole frequency range, one transverse mode at low frequencies and two at high frequencies, which perfectly confirms the original judgment. Different from that of Medium D, both the longitudinal and transverse waves of the two phases in Medium E have the same velocity. Once again, the predictions given by SFMS model show one longitudinal and transverse mode in the whole frequency range. It is also noted that when there is only one mode, the relative variation of velocities is very small, less than 0.001, which can be regarded as zero. Meanwhile, the dimensionless attenuation is also very small, normally less than 10 -3 . All these predictions are in excellent agreement with the conjectures in [12].

B. Failure of the WFMS model in predicting the dispersion behaviors of strong-property-fluctuation elastic media
In this section, we examine the velocity and attenuation of strong-property-fluctuation materials predicted by the two models. Cortical bone is a typical strong-property-fluctuation material, for which the two component phases: solid bone frame and water/marrow saturated pores, have drastically different material properties, as shown in Tab. 4. As a simplified mechanical model, cortical bone is assumed to be comprised of a homogeneous isotropic solid and water-saturated spherical pores. The pores are uniformly distributed in the solid frame, so the cortical bone is a statistically homogeneous medium. The reference velocities for cortical bones with different porosity are summarized in Tab. 5 and plotted in Fig. 14. Table 5. Material properties of the reference medium of cortical bone with different porosity.
Porosity  27 Figure 15 shows the longitudinal dispersion and attenuation of cortical bone with 10% porosity. According to the prediction given by the SFMS model, the longitudinal velocity starts from the quasi-static limit and then decreases sharply as the frequency increases. This is the famous phenomenon known as the anomalous negative dispersion, which has been observed in experiments and captured in numerical simulations [70,73,89]. After assuming its minimum, the longitudinal velocity increases quickly. At a dimensionless frequency around k0Ld=2, an additional slow mode begins to appear and its velocity slowly increases with frequency. At the dimensionless frequency k0Ld=4.5, the velocity of the fast mode approaches the upper bond, i.e., the velocity of pure bone, and then the two modes quickly separate and the velocity of the fast mode goes beyond the upper limit. This phenomenon could not be explained in the context of scattering, since the fastest propagation speed of a scattered wave is the longitudinal velocity in pure bone. This does not mean that the new model is invalid, instead, it tells us that the coherent wave disappears. The rationale behind this explana tion is as follows: As the dimensionless frequency increases, the wavelength becomes shorter and shorter. When the wavelength is comparable to the size of a single scatterer, the wave is able to discern the structure of individual scatterer and the interaction of the scatterers with the wave gets significantly enhanced. Large property contrast further strengthens the wave-scatterer interaction. As a consequence, the incoherent component occupies a large proportion of the total random wave field. Thus, when the frequency goes beyond a certain limit, the wave losses coherence. In our model, this critical frequency is identified as the frequency at which the coherent wave velocity goes beyond the faster velocity of the two component phases. The disappearance of coherent waves is systematically studied in Sheng [112].
The loss of coherence is also reported in Sato [3, pp. 58-59], where the experimental studies of the coherent wave propagation and dismiss are presented. In seismology, this frequency range is called the saturated scattering regime by Wu [1]. Wu proposed a seismic ray explanation of this phenomenon: In this regime, the wave amplitude fluctuations are saturated and rays split into numerous microrays and interfere with each other, no coherent waves exist. In Fig. 15, no similarities can be observed between the predictions given by the two models. The WFMS model gives wrong initial guess, cannot reveal the negative dispersion, and unable to predict the frequency at which the coherent waves disappear. Moreover, the WFMS model gives negative attenuation at frequencies below k0Ld=5.5. The transverse wave dispersion and attenuation are depicted in Fig. 16. According to the results calculated from the SFMS model, the velocity of transverse waves starts from the quasi-static limit and decreases monotonically as the frequency increases. At the dimensionless frequency around k0Td=2, an additional faster mode appears. This mode has nearly the same velocity as that of the original mode but with much larger attenuation, so it is very difficult to observe in experiments. As the frequency continues to increase, the fast mode undergoes positive dispersion while the velocity of the slow mode still decreases. When the dimensionless frequency k0Td reaches 5, the velocity of the fast wave exceeds the transverse velocity of pure bone, which indicates the disappearance of coherent transverse waves. Once again, the WFMS model fails. It gives positive dispersion at low frequencies and negative attenuation in the whol e frequency range.
Through the above discussion we can draw the following major conclusions: 1) The performance of the WFMS model is unstable, in certain cases it gives results with large errors or physically impractical results, such as negative attenuation, while the SFMS model gives very accurate results and the performance is extremely robust; 28 2) For two-phase media with a unique combination of mass density and Lamé constants such that both phases have the same longitudinal and/or transverse wave velocity, the SFMS model gives only one longitudinal and/or transverse propagation mode in the whole frequency range, the question raised in [12] has been resolved perfectly; 3) The coexistence of two bulk modes at intermediate to high frequency range is a dynamic bifurcation phenomenon, which universally exists in nearly all types of heterogeneous materials, it is not an artefact and the model has no deficiency in this sense; 4) The Voigt average velocity always overestimates the quasi-static velocity for heterogeneous media. This error becomes unacceptable for strong-property-fluctuation materials; 5) The WFMS model completely fails for strong-property-fluctuation materials, while the SFMS model still can give accurate prediction of the dispersion and attenuation behaviors of the coherent waves, it can also indicate the disappearance of the coherent waves by giving a velocity exceeding the upper bond of the two component phases; 6) The spectral function approach developed in [12] is unable to reveal the accurate dispersion and attenuation in the transition frequency range, and cannot distinguish between negative attenuation and positive attenuation, while the new model can accurately capture all these characteristics in the whole frequency range.

IV. PRACTICAL APPLICATIONS OF THE NEW MODEL
In this section, we perform comprehensive studies on the propagation behavior of multiple scattered elastic waves in a series of

A. Applications in seismology
The planet Earth exhibits different degrees of heterogeneity from the inner core [13], to the mantle [15,32,158], and finally to the lithosphere [3,8,[158][159][160]. The apparent attenuation of major arrivals and the existence of coda waves in nearly all the recorded seismograms are direct evidences of the existence of such inhomogeneities. However, due to lack of knowledge of wave propagation in such highly inhomogeneous media, the classical seismology is based on the multilayered model, in which each layer is treated as laterally homogeneous and vertically varying materials as a natural result of the action of gravitation. With the development of modern geophysics and seismology, more and more seismologists realize that the Earth also exhibits considerable lateral heterogeneity and tend to analyze the measured seismograms using scattering theories based on the random medium model [3, 8-9, 11-15, 178-179]. The prominent seismologist Aki first noticed the importance of coda waves and suggested that coda waves contain rich information about the statistical property of the heterogeneous lithosphere [7]. Aki [5][6], Wu [10] and Sato [3] also proposed scattering models to explain the apparent attenuation of the direct longitudinal and transverse arrivals. Up to date, scattering and attenuation of seismic waves have become the central topic of seismology [1][2][3][4]. A brief literature review shows that most previous studies are based on certain approximations, such as the Born approximation [3], the Rytov approximation [3], the single scattering approximation [18], the scalar wave approximation 29 [17,161], or the weak-fluctuation approximation [2,16]. As pointed out in the Introduction, these approximations put severe restrictions on the application of these models. Consequently, a complete understanding of the propagation behavior of seismic waves is still missing.
Many fundamental aspects of seismic waves are still poorly understood, some of which are listed below: 1. Is the attenuation of seismic waves mainly caused by scattering or by other intrinsic mechanisms? 2. Is the single scattering approximation enough to explain the apparent attenuation?
3. Can the Kramers-Kronig relation be used to explain the measured velocity and attenuation? 4. Does the Mohorovičić discontinuity really exist? 5. Is the multi-layered model appropriate to explain the seismic data?
6. Is the classical ray theory applicable to the explanation of seismic data? 7. What are the scattering dispersion and Q-factor in the whole frequency range?
With these questions in mind, we first conduct an in-depth investigation on the fundamental propagation behavior of longitudinal and transverse seismic waves in the Earth's lithosphere based on a realistic random medium model. The numerical results calculated from various combinations of material properties (density and elastic moduli) of the real rocks in the lithosphere will provide us important insights into all these questions. Through comparison of the new results with that observed in local, regional and global earthquakes, we will give novel and consistent answers to a series of longstanding problems.
Evidence from amplitude and phase fluctuation observed by the large scale seismic arrays LASA and NORSAR show that in the lithosphere, there exists a velocity fluctuation about 2% to 10% and the correlation distance of the fluctuation is about 10-20 km [5,8,160]. Moreover, it is also discovered that the velocity fluctuations are remarkably uniform. Based on these observations, the lithosphere is frequently modeled as a statistically homogeneous and isotropic random medium. Two-point correlation function is an important characteristic function to characterize various heterogeneous media. Up to now, random medium model with Gaussian, exponential and Von Karman correlation functions have all been adopted to describe the heterogeneities in the Earth [3,5,[159][160]177]. However, it is reported that the Gaussian correlation function is too smooth to represent the actual inhomogeneities in the Earth. Transverse coherence function (TCF) and angular coherence function (ACF) analysis measured by the LASA seismic array reveal [8] that the exponential or the Von Karman-type correlation functions fit the seismic data better. For the above reasons, in the present work we model the heterogeneous lithosphere as an exponential-type random medium. To model the velocity fluctuation, we consider a twophase continuum composed of a hard phase with relatively large longitudinal and transverse velocities and a soft phase with relatively low velocities. Since the velocity fluctuation in the lithosphere is rather uniform, we assume that the two phases have equal volume fraction, i.e., f1=f2=50%.
Dispersion and attenuation are two prominent features of waves propagating in a heterogeneous medium. However, researches concerning the dispersion of seismic waves are rather limited. This is partially because of the irregular surface of the Earth which significantly disturbs the measurement of dispersion. Another reason is that most seismic signals lie in the high frequency range, or even in the geometric range, in which the seismic waves have very small dispersion. The majority of existing studies are concentrated on the attenuation of direct arrivals and coda waves. It is well known that scattering attenuation and intrinsic attenuation are two major physical mechanisms for the overall attenuation of seismic signals. Intrinsic attenuation transfers elastic energy into heat and thus causes energy loss of the total wavefield. Contrarily, scattering cause magnitude decrease of the coherent waves by exciting coda waves, while the elastic energy of the total wavefield is conserved. The relative contribution of the two mechanisms to the overall seismic attenuation is an open problem in seismology. Nevertheless, many seismologists believe that the scattering attenuation contribute the major part of the total attenuation [3]. It was also pointed out in [9] that multiple scattering theory is necessary to separate the contributions from the two mechanisms. In this work, we assume that the scattering attenuation is the dominant mechanism and all the materials are considered as purely elastic materials. The answer to this question will be drawn after obtaining the exact numerical solution and comparing the numerical results with measured seismic data.
Seismologists usually quantify the attenuation of seismic waves by using the inverse quality factor Q -1 instead of the attenuation coefficient α. For simplicity, we call the inverse Q factor as Q-factor directly. The quality factor is defined as the ratio of the energy dissipated in a cycle of vibration to the total energy supplied. Carcione and Cavallini [162] conducted a rigorous derivation of the Q factor based on this definition and gives the following expression:

30
where k is the wavenumber of the longitudinal or transverse waves.
In textbooks of seismology [2,3], the quality factor is also defined by the asymptotic expression of a spherically outgoing body wave: where Q is the quality factor, f is the frequency, V is the velocity, r is the distance between the source and the field point, and the prefactor 1/r accounts for the geometric spreading.
Considering the definition of attenuation coefficient α: we have the following identity: or equivalently: On the other hand, k can be split into real and imaginary parts as: Substitution of (113) into (108) yields: When the attenuation coefficient is small compared to the real part of k, terms of second and higher order in α can be neglected, thus we Equation (118) gives a rough estimation of the Q-factor ratio for seismic waves in the geometric regime. It tells us that in the geometric regime the Q-factor ratio is frequency-independent and solely determined by the wave velocities. This conclusion will be verified through numerical calculations for a series of material combinations, as shown in the following examples.
Velocity fluctuation can be caused either by density fluctuation or by elastic modulus fluctuation. In the following numerical examples, we study the effects of the two types of properties separately. Seismological and geophysical studies show that the density and velocity of rocks in the Earth's lithosphere have a linear relationship, known as Birch's law [3]. Although the specific coefficients appeared in the Birch law strongly depend on the average atomic weight of the rocks, there is a general conclusion that rocks with larger 31 density also have greater longitudinal and transverse velocities. Geophysicists have proposed different velocity models for the substructure in the lithosphere using various seismic tomography techniques [163][164][165]. It is generally acknowledged that the longitudinal velocity in the lithosphere varies from 5.8 to 8.3 km/s, while the transverse wave velocity varies from 3.5 to 4.7 km/s. In the subsequent calculations, we choose the material properties following these general guidelines. We first study the effects of modulus fluctuation. The material properties of Media I-IV are listed in Tab. 6. It is noted that the material properties of Phase 1 do not change while the elastic moduli of Phase 2 increase gradually, which results in a velocity fluctuation from ±2.5% for Medium I to about ±6% for Medium IV.  additional, faster mode begins to appear. As the frequency continues to increase, the velocity of the fast (slow) mode quickly increases (decreases) and then approaches its geometric limit, i.e., the longitudinal velocity of the fast (slow) phase. At high frequencies, the dispersions of both the fast and slow modes are very small. As a consequence, seismic signals at relatively low frequencies and at high frequencies can propagate for a long distance while maintaining it waveform nearly unchanged. To show the detailed variation at very low frequency and very high frequency, the Q factors are plotted in a double-logarithmic coordinate system, as shown in the second row in Fig.17. At low frequency k0Ld < 1, the attenuation increases from zero following a power law. In the intermediate frequency range, the Q-factor increases with frequency nonlinearly, then it assumes its maximum soon after the faster mode appears. Then the attenuation decreases with frequency following a negative power law. The Q-factor of the fast mode is very large at its emergence and then decreases dramatically. After the slow mode reaching its summit, it decreases following a negative power law similar to that of the slow mode but with a slightly larger value. With the increase of velocity fluctuation, the frequency of the attenuation peak decreases significantly, from k0Ld =41.60 for Medium I to k0Ld =17.12 for Medium IV. The value of the peak increases from about 0.047 for Medium I to 0.10 for Medium IV. The Q-factor in the geometric regime k0Ld > 100 decreases from 0.01 to 0.001.  It is observed that KcrL is always greater than KcrT. At low frequencies, the longitudinal Q-factor is smaller than that of the transverse waves. At frequencies higher than kcrL, longitudinal Q-factors are always greater than transverse Q-factors, and both of them strictly follow a negative power law. The ratios of the longitudinal to transverse Q-factors have a very complicated profile. At low frequencies, the ratio is relatively small, near 0.5. As the frequency increases, the ratio increases gradually and then reaches a peak. As the frequency continues to increase, the ratio gradually decreases back to less than unity and then decreases sharply to its minimum at KcrT. After reaching its minimum, the Q-factors increases monotonically and reaches a stable value between 1.5 and 2.5 at KcrL. Since there are two sets of longitudinal and transverse waves, we have four different Q-factor ratios in the high frequency range. The most prominent feature in the high frequency regime k0Ld > KcrL is that the ratio keeps constant at four different values. The differences among the four values increase notably as the velocity fluctuations increase from Medium I to Medium IV. It is also worth noting that the peak of the Q-factor ratio decreases from 3 to 1 as the velocity fluctuations increase from Medium I to Medium IV. Table 8 gives the Q-factor ratios in the geometric range and the velocity ratios of the two component phases. It is seen the ratios roughly agrees with the approximate relation given by Eq. (118). As the velocity fluctuation increases from Medium I to Medium IV, the relative error increases from 3.24% to 14.6%.    To show their relative changes, the longitudinal and transverse Q-factors and their ratios are plotted in the same frequency scale, as shown in Fig. 22. As can be seen, the fluctuation in density significantly changes the relative magnitude of the longitudinal and transverse Q-factors in the whole frequency range. The ratio starts from about 0.5 and then increases gradually until reaching its maximum, and then decreases rapidly to a minimum about 0.5 at KcrT. The maximum value of the ratio is strongly dependent on the property fluctuations.
As we can see in the figures, it reaches 1 for Medium V and 2.75 for Medium VIII. At frequencies greater than KcrL, the Q-factor ratios of the four combinations approach four different high-frequency limits. To validate the approximate relation given by Eq. (118), the Q-factor ratios and the velocity ratios are listed in Tab. 11. It is seen that the relation is roughly satisfied within an average error about 15%. The We need to point out that the Medium XVI has much larger velocities and will be used as an exemplar model for the next section.  propagation modes are nearly the same as that of the homogeneous reference media. The dispersion of these modes at low frequency k0Ld < 1 or k0Td < 1 is negligible, while in the frequency range 1 < k0Ld, k0Td < 8 the dispersion increases slightly but still very small. In the low frequency range k0Ld < 1 or k0Td < 1, the Q-factors of these media increase from 10 -7 to 0.01 following a power law. In the frequency range 1 < k0Ld, k0Td < 8, the Q-factors increase with frequency following a complicated nonlinear law, and the relative magnitude of the longitudinal and transverse Q-factors strongly depend on the material properties. We can also find that the Q-factors waves is always smaller than that of the transverse waves at low frequencies k0Ld, k0Td < 1, and greater than the transverse Q-factors in the high frequency range k0Ld, k0Td > 30. The critical frequency at which the longitudinal Q-factors reaching their summits is always greater than that of the transverse Q-factors, i.e. KcrL > KcrT. The Q-factor ratios exhibit a rather complicated pattern. At low frequencies k0Ld, k0Td < 1, the Q-factor ratios increase gradually from about 0.5. As the frequency continues to increase, the Q-factor ratios increase rapidly and reaches a summit ranging from 1.0 to 3.0 at k0Td = 2~3. After assuming its summit, the Q-factors decrease rapidly until reaching the minimum at KcrT. The minimum value is also strongly dependent on the material properties, which normally ranges from 0.5 to 1.2. At frequencies higher than KcrT, The Q-factor increases rapidly again. Since in the high frequency range, there are two sets of longitudinal and transverse waves, there are four different combination of Q-factor ratios. As is evident from the figures, the highfrequency limits of the ratios vary from 1.1 to 3.2, for which the range is much larger than that of Media I-VIII. It is worth noting that the high-frequency limits of the four ratios have the following relation: . This conclusion informs us that for high frequency seismic waves and at large epicentral distances, the amplitude of the longitudinal waves is smaller than that of transverse waves. Moreover, the slow transverse wave has the largest amplitude and carries the largest amount of energy, which is consistent with observations in real earthquakes.  Both density and elastic moduli have significant influence on the propagation characteristics of seismic waves. All these theoretical results help us gain important insight into the dispersion and attenuation characteristics of seismic waves measured in local, regional and global earthquakes. In subsequent sections, the theoretical results will be applied to explain seismic data collected worldwide.

Explanation of the observed Q-factors and Q-factor ratios
Attenuation of the direct arrivals of the longitudinal and transverse waves has been the central topic of nemerious seismological studies ever since the beginning of modern seismology, interested readers please see the manograph by Sato and collaborators [3] and references therein. The longitudinal and transverse Q-factors and their ratios of a series of local and regional earthquakes occurred worldwide are summarized in monograph [3]. For the convenience of subsequent diacussion, the results are presented in Fig. 25. measured in local, regional and global earthquakes [3].
It needs to mention that the notations 1 P Q  and 1 S Q  denote the Q-factors of primary and secondary waves, corresponding to the Qfactors of the longitudinal and transverse waves in this work. Before comparing the measured Q-factors with the theoretical predictions, we first give a rough estimation of the dimensionless frequency to which the measured curves correspond. As pointed out before, the average velocities of the longitudinal and transverse waves are about 7 km/s and 4 km/s, and the majority of the measured curves lie between 1 to 100 Hz. The correlation length, i.e., the characteristic dimension of the rocks in the lithosphere lies between several hundred meters to tens of kilometers [5,8], here we use d =10 km for the estimation. Calculation using the above data gives k0Ld = 9-900 and k0Td = 16-1600, thus we can see most seismic signals lie in the intermediate to high frequency ranges. In this regime, both the longitudinal and the transverse Q-factors increase from about 0.001, reaching their peaks lying between 0.05 to 0.2, and then decrease monotonically to about 0.001. In the high frequency regime, the Q-factor curves follows an inverse power law, i.e., decreasing quasi-linearly in the double-logarithmic coordinate system. Comparing the theoretical predictions with the measured curves, we see all these theoretical predictions show excellent agreement with that measured in local and regional events. Both the longitudinal and transverse Q -factors decrease following an inverse power law in the frequency band 1< f < 100 Hz. Meanwhile, the magnitude of the Q-factors lies in the range of 0.1 to 0.001, which is also consistent with the theoretical calculations. In addition, it is observed that the Curve 16 exhibits a sharp peak at around 1 Hz. This curve partially confirms Aki's conjecture that there should exist a peak in the Q-factor curves in the intermediate frequency range [3]. Extending all the measured curves to the intermediate frequencies 0.05 < f < 1Hz, we can infer that the peaks of the longitudinal and transverse Q-factor curves lie between 0.1 and 0.2. Our theoretical predictions, as shown in Figs. 17-24, along with this reasonable inference fully support Aki's conjecture. At this point, we are ready to give an explanation on why the attenuation peaks are rarely captured in practical measurements. From the theoretical calculation, we can see the peak occurs at a critical frequency at which the dispersion curve branches. Correspondingly, the waves with a central frequency near the critical frequency have both very large dispersion and attenuation. As a consequence, the wavepackets in this range quickly broaden and decay as the epicentral distance increases, and finally dismissed in the ambient noises. Therefore, it is an extremely challenging task for seismic recording systems to identify weak signals near the critical frequency. We also note that the curves numbered 1, 2.1, 2.2 and 2.3 lie in the low to intermediate frequency regimes. Due to lack of accurate data, they are marked with thick lines. It was pointed out that attenuation at low frequency is extremely difficult to be measured accurately [3].
The measured Q-factor ratios are shown in Fig. 25(c). At first glance, we may find it is difficult to extract any law by which these curves abide. After careful inspection, we find the curves can be classified into six categories according to their variation tendencies:  17.2, 19. Surprisingly, we can establish a rough correspondence between these curves and the theoretical curves: Curves 1 and 3 correspond to the low frequency range in which the ratio is nearly constant and lies near 0.5; Curves 2.1, 2.2 and 16 correspond to the ascendant stage in which the ratio raises from 0.5 to 1.5, 2 or 3, depending on the specific material property fluctuation; Curves 5, 18 correspond to the descendant stage; Curves 5, 14 correspond to the turning stage in which the ratio decreases to its minimum and then increases; Curves 4, 6, 8, 10, 13, 17.1, 17.2 and 19 correspond to the plateau stage in which the ratio lies between 1 and 3. These curves correspond to the high frequency range, also called the saturated scattering region, in which the ratios are stabilized at four different constants. The shaded region 9 also cover the range of high-frequency limits of the four Qfactor ratios. It need to mention the Curves 17.1 and 17.2 are for the ratios of longitudinal waves to the SV and SH components of transverse waves, respectively, which are less than unity. However, if recalculated for the complete transverse waves, it is highly possible that it will give a value greater than unity. This is because at high frequencies, the attenuation of longitudinal waves is larger than that of transverse waves. The existence of Q-factor ratios with magnitudes less than 1.5 or near unity in the high frequency range strongly suggests that the lithosphere is a heterogeneous medium with strong property fluctuation, as evident in Tab. 14. Through the comparison we see the new model can give an accurate quantitative prediction of the longitudinal and transverse attenuations and their ratios. The rich variety of the tendency of the Q-factor ratios is also perfectly captured in the new model. At this point we need to mention the scattering theory proposed in [10] based on the Born approximation predicts that the Q-factor curve increases monotonically with frequency, which is in contradiction with the measured seismic attenuation. We also note the curves in the high-frequency regime, for which the magnitude of the ratio is a constant near unity, along with the curves in the descendent stage and the turning stage, cannot be explained using the travel-time corrected Born approximation [3].
At last, we try to give a rough estimate of the correlation lengths of the portion of the lithosphere that corresponds to each curve of the Q-factor ratios in Fig. 25(c). Although the information for the velocities is missing, existing seismic data recorded in numerous seismic events reveal that the reference velocities of the transverse waves lie in the range from 3800 to 4000 m/s. For the purpose of a primary evaluation, we use V0T = 3900 m/s. From Figs. 23-24, we can observe that the peak of the Q-factor ratios for most medium models occurs at K0T = 2-3, here we use K0TMAX = 2.5. The minima of the Q-factor ratios occur in the range 8 < K0T < 15, here we choose

Does the Mohorovičić discontinuity really exist?
On October 8, 1909, an earthquake with the macroseismic intensity VIII ˚MCS struck Pokupsko, 40 km southeast of Zagreb, Croatia [166]. Mohorovičić collected the seismograms recorded by numerous European seismological stations and plotted the travel time curves, as shown in Fig. 26(a). He immediately recognized that there are two sets of longitudinal and transverse waves. To clearly indicate these phases, we replot the travel time curves in Fig. 26(b). From the travel-time curve, we can see four different arrivals, denoted by Pn, Pg, Sn and Sg, respectively. The velocities of the four phases can be roughly estimated from the curves, VPn = 8230 m/s, VPg = 5680 m/s, VSn = 4640 m/s, VSg = 3310 m/s. It is also noted that the Pn and Sn phases begin to appear from the epicentral distance d = 300 km and propagate a long distance even longer than d = 2600 km, the Pg and Sg phases exist in the range from the epicenter to an epicentral distance d = 1500 km. and (b) best fitting curves Since then, more and more seismograms recorded in the local, regional and global earthquakes show that two sets of longitudinal and transverse waves can be observed. Figure 27 shows  Fig. 26, all phases propagate a long distance greater than 1500 km. The Sg phase propagates an exceedingly large distance, even longer than 2500 km. It is worth noting that both the Pn and Pg phases are difficult to discriminate in the range from 0 to 500 km, and start to separate from each other at distances greater than 500 km. Similar phenomena can be observed for the Sn and Sg phases.   As the velocity of the material increases with depth, the wave gradually reflected back and gets refracted again at the discontinuity.
After that, the wave propagates in the upper medium again and arrives at the observation station. A schematic diagram for the paths of the seismic rays are given in Fig. 30 (a).
As similar phenomena are observed in the travel-time curves extracted from numerous local, regional and teleseismic earthquakes, as shown in Figs. 27-29, seismologists infer that there must be a global discontinuity existing in the lithosphere. In honor of its discoverer, this discontinuity is named the Mohorovičić discontinuity. Since then, the two-layered model of the lithosphere is accepted worldwide.
The top layer with relatively low velocity is called the crust, and the P-and S-waves propagating in this layer are denoted by Pg and Sg, respectively, while the bottom layer with relatively high velocity is called the mantle, and the waves propagating in this layer are denoted by Pn and Sn [164]. The Mohorovičić discontinuity is identified as the interface between the crust and the mantle. In 1925, Victor Conrad first identified a new set of longitudinal and transverse waves, denoted by P* and S*, when he analyzed the eastern Alps earthquake in 1923 [130,[172][173]. Since the new phases propagate at an intermediate velocity between that of the upper crust and mantle, Conrad postulated that the waves propagate in the lower crust and there must be an additional interface which separates the upper and lower crust. The new phases are also observed worldwide as shown in Fig. 29. Consequently, the interface between the upper and lower crust is named as the Conrad discontinuity. Later on, additional phases are identified in local and regional earthquakes and correspondingly, new discontinuities such as the Riel Discontinuity and the Fortsch discontinuity are identified. Currently, the multilayered model of the lithosphere is widely used in the seismological community for seismic tomography and inversion for the velocity structure [164][165].
Despite its success in explaining the coexistence of two or more sets of body waves, disputes about the multi-layered model along with the existence of the Mohorovičić discontinuity and the Conrad discontinuity never end in the scientific community [130,132,173].
Modern seismology improves the estimation of the velocity distribution and predicts the velocities of primary seismic waves (P-waves) immediately above the Moho are about 6.7-7.2 km/s, and below the Moho they are about 7.6-8.6 km/s [164]. One naturally raises the question: what is the physical mechanism that causes such an abrupt change in velocity, phase transition or change in chemical composition? It is pointed out that phase transition induced by temperature increase results in a transition zone with gradually varying material properties instead of a sharp interface [130]. Thus, geophysicists suggest that the Moho marks a change in chemical composition.
From the knowledge of petrology, we know the longitudinal velocity immediately above the Moho is consistent with that of basalt, while the velocity immediately below the Moho is similar to those of peridotite or dunite. Nevertheless, any believable sorting mechanism gives an irregular surface instead of an even interface [130,132]. Consequently, both the mechanisms encounter difficulties in explaining such an even and sharp discontinuity. Moreover, as discussed above, the amplitudes of all these waves have obvious  Fig. 30(b). If the wave packet emanated from the source contains both low-or intermediate-frequency components and high-frequency components, the original wave packet will split into three branches, as shown in Fig. 29. The phase velocities of the P* and S* phases lie near that of the homogeneous reference medium. In real measurements, the P* and S* phases may be contaminated by the P-coda waves or S-coda waves, so the recorded signal is a superposition of low-or intermedia-frequency phase signal and the high-frequency noises.  Fig. 28(a), Medium XV corresponds to Fig. 28(b), and Medium XVI corresponds to Fig. 29. With these correspondences in mind, we can try to give a consistent explanation for both the dispersion and attenuation characteristics of the measured data. Fig. 25 shows that most seismic signals decrease with frequency following an inverse power law, which indicates that the signals lie in the high-frequency regime. Meanwhile, the majority of the observed travel-time curves exhibit four branches corresponding to the Pn, Pg, Sn and Sg phases, respectively, which 45 also indicates that most seismic signals lie in the high-frequency regime. By adjusting the properties of the component phases, the Qfactor ratios show a variety of different patterns which are sufficient to cover all the values and variations of measured Q-factor ratios, as shown in Fig. 25(c). The travel-time curves in Fig. 28(b) exhibit high-velocity, weak velocity fluctuations and low attenuation, which are related to the existence of the subducting Pacific and Philippine sea slabs. Correspondingly, these distinctive features are also fully captured in the numerical example for Medium XV. From Fig. 24, we can see Medium XV has high velocities, small velocity fluctuations and low attenuation (high-Q). Medium XVI gives a material model corresponding to Fig. 29, from which the P* and S* phases can be explicitly identified. Comparing the velocities in Fig. 29 and that shown in Fig. 24, we can find the velocities of all the phases calculated from Medium XV show excellent agreement with that measured in real seismic events.
In the new explanation, there needs not be any even interface that separates two or more layers with dramatically different material properties. Instead, it only requires that there exists a certain degree of property fluctuation and the fluctuation is statistically uniform, which are more likely to be realized in the heterogeneous Earth.
Through the above discussion we see the new model is capable of providing an accurate and consistent explanation to nearly all the observed seismic data. At this point we are ready to answer several long-standing problems raised in the seismology community, and draw a series of new conclusions: 1) The measured apparent attenuation for short-period seismic waves are mainly caused by scattering in the heterogeneous Earth, the portion contributed by anelasticity mechanisms is less important. This conclusion is also confirmed by the observation that seismic waves can propagates a very long distance (the seismic waves frequently run several times around the earth before their energy is completely dissipated) and last for a long time (about tens of minutes); 2) Multiple scattering theory is a necessary tool for the explanation of the attenuation of the major phases Pn, P*, Pg, Sn, S* and Sg, single scattering theory is insufficient; 3) Elastic waves are tensor waves in their nature, and the causality of the multiple scattered waves is secured by the proper choice of the imaginary part of the homogeneous Green's functions. The Kramers-Kronig relation, which was derived for scalar waves, cannot be applied to the scattered elastic waves [65][66][133][134]. This is because the K-K relation always requires that there exists a single 7) The existence of inhomogeneities that large enough to break the stochastic homogeneous assumption may change the waveform significantly. Thus, the model provides a new imaging mechanism for inversion of the large-scale structures in the lithosphere or beneath the lithosphere [146,174,180].
Finally, we need to stress that the theoretical predictions of the dispersion and attenuation are for the coherent waves. In an idealized case, the coherent waves are obtained by summing up the recorded waves from a number of observation stations that located concentrically about the seismic source. The quantities measured at a single observation station may be distorted by the near-receiver inhomogeneities, as pointed out by Wu [10]. The development of large scale seismic arrays, such as LASA, NORSAR, has significantly improved our ability to obtain such coherent waves by collecting large amounts of single-station signals. These seismic arrays along with modern signal processing techniques provide us a great opportunity to extract the complete set of coherent wave information, including the dispersion and attenuation of each phase. Based on the high-quality data, the new model can perform reliable inversion for the material property fluctuation and characteristic length of inhomogeneities with the aid of advanced optimization algorithms like the least-square method [175]. Furthermore, it is well known that mineral depositions can also influence the material property fluctuation, thus the inversion results can be used for exploring natural resources, such as hydrocarbone reservoirs and coal and iron mines.

B. Applications in the ultrasonic nondestructive evaluation (NDE)
Porous metals have broad applications in aerospace industry [181][182][183][184] and biomedical engineering [185][186][187] due to its light weight and high strength. An additional prominent feature is that its mechanical properties are adjustable and can perfectly match those of biomedical materials like bone, so porous metals like porous Titanium are also used as bone substitution materials [186]. Studying the dispersion behavior of ultrasound in these materials is of vital importance for the nondestructive evaluation and characterization of porous metals [188][189][190][191][192]. Two-phase metals and particle reinforced composites are other categories of heterogeneous materials widely used in industry. These materials have unique physical and chemical properties, such as chemical stability, high strength, tunable mechanical properties [193][194]. The microstructure is the key factor that determines the material properties like fracture toughness, hardness, impact strength, yield strength, and tensile strength. Consequently, there is considerable interest in developing ultrasonic techniques to nondestructively characterize microstructures and finally determine these mechanical properties. It is well known that the ultrasonic attenuation and dispersion are very sensitive to the characteristic parameters of microstructures, like elastic moduli of the component phases, average sizes of the grains, so they provide important parameters for microstructural characterization. In general, the interpretation of frequency dependent velocity and attenuation data requires the use of multiple scattering theory. However, previous research on ultrasonic NDE has been mainly focused on characterization of metal polycrystals such as Titanium alloy, Nicole alloy or iron, which are generally regarded as weak scattering materials [42][43][44][45][46][47][48][49][50]. Very little effort is devoted to characterization of materials with strong property contrast, such as porous metals and two-phase alloys [188][189][190][191][192]. The model developed in this work exactly meet these requirements. As applications, we calculate the velocity dispersion and attenuation of porous aluminum with various porosities and Cu-Al alloys with different phase volume ratios.

Nondestructive evaluation of porous metals
Porous aluminum is an extreme example of two-phase materials. The material properties used in the calculation are shown in Tab. 15,  From the results, we can find the following propagation characteristics of coherent waves in porous metals: 1) The longitudinal velocity first undergoes steep negative dispersion, after reaching the minimum value, it becomes positively dispersive and quickly goes beyond the longitudinal velocity of pure aluminum, which indicates the disappearance of the coherent waves.
Only one longitudinal mode is found in the whole frequency range. This phenomenon tells us that the branch behavior of longitudinal coherent waves is strictly related to the elastic property of the component materials, or in other words, two longitudinal modes can only exist when both of the component phases are elastic materials. The dimensionless attenuation-frequency relation follows a power law at frequencies lower than 1. At dimensionless frequencies higher than 1, the attenuation approaches unity.
2) The coherent transverse waves first undergo negative dispersion with a magnitude smaller than that of the longitudinal waves, then become positively dispersive and approach the transverse wave velocity of pure aluminum. A second mode appears in the intermediate frequency range. The attenuation of the second mode is much larger than the first one at its initial appearance and then decreases quickly.
At moderately high frequencies, the velocity of the fast mode quickly goes beyond the transverse velocity of pure aluminum, which indicates the disappearance of coherent waves. It is also noted that porous Al with low porosity support coherent transverse waves in a relatively wider frequency band.

Nondestructive evaluation of two-phase alloys
Nondestructive characterization of microstructures in multiphase alloys by means of ultrasonic scattering has received extensive studies. The mechanical properties of frequently used metallic phases, including Fe, Al, Cu, Ti, or Nb all exhibit strong property flu ctuation. Our model is naturally suitable for these types of materials. As an example, we consider a two-phase alloy made of Al and Cu [194]. The material properties of Al and Cu are given in Tab. 17. We consider Cu-Al alloys with different volume ratios. Tab. 18 lists the material properties of the reference media.

C. Applications in bone quantitative characterization
Cortical bone is the major supporting structure of human body and most osteoporotic fractures are happened in it. Cortical bone is composed of compact bone and distributed pores saturated with marrow. The pore diameter is normally in the range of 20-300 micrometer and the volume fraction lies between 5% and 15% [70]. However, as a consequence of aging and continuous loss of minerals, the porosity of the middle-aged and elderly people may increase to 30% [78]. Early diagnosis and treatment of osteoporosis are essential for preventing bone fracture. Compared to microXRT or MRI, ultrasound provides an economic and noninvasive approach for the diagnosis of osteoporosis and monitoring the bone status. Different from MRI and XRT, which can only provide information about the porosity and average pore diameter, ultrasonic signals also carry rich information of the mechanical properties of bone, including the density and elastic stiffness. Therefore, quantitative ultrasound (QUS) provides a promising technique for predicting the bone mechanical strength. For this reason, bone quantitative ultrasound has become an active research field [70]. Multiple scattering theory of ultrasonic waves plays a central role in characterizing microarchitecture in bone. However, modeling the multiple scattering in cortical 50 bones progresses very slowly. One reason is that cortical bone is an extremely complicated heterogeneous material with hierarchical microstructures covering a wide length scale. Another reason is that the wave-scatter interactions are very strong since wavelengths of ultrasound excited by most biomedical transducers are comparable to the characteristic dimension of microarchitectures in cortical bone.
Reported experimental and numerical results are also full of controversy. The majority of the experimental research report that wave velocity is negatively dispersive [89], while some other researchers report that velocity have positive dispersion [73]. The interesting and diverse phenomena have attracted extensive attention from both academy and clinical community [70]. Due to the lack of a thorough understanding of the rich phenomena, quantitative ultrasound techniques for cortical bone characterization has been hampered for many years. It is generally realized the single scattering model and Biot's model are insufficient for characterization of bone samples. For instance, the single scattering model can predict dispersion and attenuation when the number density of pores is low, normally less than 5%, and the error increases as the porosity becomes larger. The Biot model cannot predict the anomalous negative dispersion and it involves several phenomenological parameters that difficult to be measured experimentally. In this section, we conduct comprehensive numerical calculation for cortical bones with different porosities and demonstrate the capability of the new model to quantitatively predict the scattering characteristics of cortical bones. In this work, real cortical bones are modeled as a poroelastic medium with spherical pores saturated by water. The dimensionless velocity variation and the attenuation are plotted versus the dimensionless frequency. The properties of the reference medium are listed in Tab. 5.

Dispersion and attenuation of longitudinal coherent waves
The dispersion and attenuation of longitudinal waves are shown in Fig. 35. According to the characteristics of the dispersion curves, the frequency domain can be divided into four regions: Region I (0 < k0Ld < 2): Starting from quasi-static limits, the wave undergoes negative dispersion. The velocity assumes its minimum at a dimensionless frequency near unity. The attenuation increases with frequency following a power law. There is only one propagation mode in this region; Region II: After assuming its minimum, the velocity becomes positively dispersive. The attenuation increases with frequency in a general nonlinear manner. There is still only one propagation mode; Region III: There are two propagation modes in this region, both have positive dispersion. The width of this range is strongly dependent on the volume fraction of the pores, the smaller the fpore, the wider the Region III. At the emergence of the second mode, its velocity is nearly identical to that of the first mode, then the velocity difference becomes larger as the frequency continues to increases. We can also observe that the difference increases with porosity. This result tells us that the two propagation modes are more easily observed in bones with high-volume fractions. The attenuation of the second mode is very large and then decreases rapidly. When the pore volume fraction is smaller than 6%, the fast mode has larger attenuation, when fpore is greater than 7%, the slow mode has larger attenuation; Region IV: The velocity of the fast mode quickly approaches the upper limit, and then goes beyond that of pure bone, which indicates that the coherent wave disappears. The slow mode becomes negatively dispersive and quickly approaches its lower limit.
The above results are presented in the dimensionless parameter space. However, in practical applications we need to analyze the dispersion behavior in the physical space. In doing so, the velocity and attenuation coefficient against the frequency f are replotted in Figs. 36. Figure 36 shows the effects of porosity on the dispersion behavior by considering cortical bones with fixed pore diameter d=100 μm and various porosities, ranging from 1% to 10%. The coexistence of two propagation modes is frequently mentioned in the literature but rarely quantified to the author's best knowledge.
The new model allows us to estimate the condition under which two propagation modes can exist. It is seen from Fig. 35 that when the second mode just emerge, its velocity is extremely close to the first one. Meanwhile, the attenuation of the new mode is very large, so it is extremely difficult to discriminate the two modes. Indeed, the wave observed in experiments is a superposition of the two modes, and the dispersion behavior lies between the two. When the volume fraction of the pores is relatively large, reaching a value between 20% to 30%, the difference between the velocities of the two modes is relatively large. If the average pore size is also not very small, about 200~300 μm in diameter, it is possible to observe the two modes experimentally.

Dispersion and attenuation of transverse coherent waves
The dispersion and attenuation of transverse waves are shown in Fig. 38 [29,71,[195][196][197][198], finite element method (FEM) [67][68][199][200] and spectral element method (SEM) [195]. Compared to experimental approaches, numerical simulation provides an efficient and flexible way to validate different theoretical models. The most attractive feature is its capability of generating finely controlled microstructures. Numerical models for heterogeneous materials with different types of spatial correlation functions [29], with specific material property combination, and with different realizations of random media can be created using advanced preprocessing techniques, such as the Voronoi tessellations [67][68]155]. Moreover, time domain signals with arbitrary waveform and at any frequencies can also be generated, which enables the possibility of extracting the dispersion and attenuation features in a wide frequency band. Another advantage of numerical simulation is it provides a flexible approach to extract signals at any numerical nodes and at any sections at the post-processing stage, which is impractical in experiments. Numerical approaches do not introduce approximations for the scattering process, in other words, it includes all multiple scattering events and thus contains the complete information of the stochastic wavefield. Nevertheless, conventional numerical approaches suffer a number of drawbacks. For instance, discretization with large amount of meshes often introduces numerical dispersion and attenuation, especially for the case of high frequency signals [195]. Additionally, accurate description of the phase boundaries is extremely difficult, approximation to the geometry generates spurious diffraction [198]. Furthermore, the discontinuity of material properties cannot be accurately accounted for. The finite difference method frequently replaces the phase interfaces separating materials with sharp property contrast by an additional mesh with average properties of the adjacent materials [70]. Such homogenization scheme will inevitably change the scattering behavior of the inclusions. It is also noted the derivatives of the displacement is discontinuous at the phase boundaries, which may lead to numerical instability. Consequently, developing new numerical algorithms is necessary to conduct so that each of the signals are completely separated. Additionally, high-quality signal processing techniques are also necessary for the analysis of its spectral content since these signals are highly attenuated.
Researchers have conducted comprehensive experimental studies on the velocity dispersion and attenuation of multiple scattered waves in polycrystalline metals [201][202], particulate composites [203][204] and cortical and trabecular bone [67,[72][73]. These works 56 provide valuable data for evaluating the validity and accuracy of different theoretical models. However, the relevance of the data obtained in the previous experiments to the theoretical predictions in this work is questionable. Accurate experimental verification of the new model is still a very challenging task. A successful experimental verification should meet the following requirements: (1) The microstructure must be finely controlled, including the shape, size and distribution, especially that the spatial autocorrelation functions should match the theoretical hypothesis; (2) No microcracks, interfacial diffusion/debonding or air bubbles should be introduced in the sample; (3) Planar and/or focused transducers are frequently used in practical backscattering and transmission measurements, both of which have finite beam width. Consequently, the diffraction effects must be properly corrected before the dispersion and attenuation of a plane wave component can be extracted [40,[205][206]; (4) Sufficient propagation distance is required to establish a stable multiple scattering signal, so the sample should have a reasonable size. 3D printing is a promising technology for manufacturing such heterogeneous materials, but its capability still needs to verify.
The current work has laid the foundation of the multiple scattering theory for strong fluctuation materials. It can be extended in many aspects to incorporate more complicated microstructures. To the author's best knowledge, the following directions for future developments can be distinguished. First, it is reported that a number of heterogeneous materials have ellipsoidal inclusions. For example, Lotus-Type porous metals have slender cylindrical pores aligned in one direction [181]. High-temperature titanium alloys used in jet engine foils and disks also exhibit ellipsoidal grains due to certain thermomechanical procedures like rolling and casting [209]. To capture these microstructural information, a multiple scattering model for heterogeneous media with aligned ellipsoidal inclusions are in need. Second, the current model only considers materials with isotropic constituent materials. It is well known that polycrystalline alloys are composed of anisotropic crystallites with random or preferred orientations (textures). In order to take the crystallographic information into consideration, a multiple scattering model for alloys with anisotropic grains is necessary. Third, in this work we only considered random media with exponential correlation functions, which is appropriate for the description of inclusions with sharp boundaries. For heterogeneities with blurred boundary, Gaussian correlation function is the most proper choice. The Von-Karman type correlation function is adopted to describe inhomogeneous materials with self-similar microstructures. Thus, the current model can also be extended to study the scattering phenomena occurred in heterogeneous materials with different statistical characteristics. Fourth, this work only considers pure elastic materials, while the anelastic behaviors are completely neglected. Studying the relative contribution of scattering and viscoelasticity to the overall attenuation is important for a more accurate explanation of the measured seismi c data.
According to the theory of viscoelasticity, the effects of viscosity can be accounted for by introducing an imaginary part in the elastic constants. The work provides a theoretical framework naturally suitable for incorporating this type of complexities. Fifth, in the spirit of this work, a system of renormalized Bethe-Salpeter equation can also be developed, which further paves the way for the development of new radiative transfer theory. The new theory is anticipated to give a more accurate description of the energy diffusion in the regime where no coherent waves can exist [3,48,125].

VI. CONCLUSION
This work provides an accurate and universal theoretical framework for investigation of elastic wave scattering in a wide variety of heterogeneous media, such as heterogeneous earth, metallic polycrystals, granular composite materials and biomedical materials like

AUTHOR'S DECLARATION
This work is independently accomplished by the author and he takes full responsibility for the whole work.

APPENDIX. CALCULATION OF THE SINGULARITY OF GREEN'S TENSOR
Green's functions are generally used to calculate the field induced by certain distributed sources, i.e., distributed forces for the case of elastodynamic problems. The resulting field is frequently expressed as integrals of spatial convolution type, as shown in Eq. (39).
Different from the case of scalar wave scattering where only Green's function is introduced, the integral appeared in the elastic wave scattering also contains its first-and second-order derivatives. For the convenience of subsequent discussion, we list the elastodynamic Green's function and its derivatives explicitly: Green's tensor. In the following section, we only consider spherical inclusions and introduce three different but equivalent methods to calculate the singularity of Green's tensor.

Method 1: Calculation of the singularity in the spatial domain
To calculate the principal value of the singular integrals, we assume that a body force of unit magnitude is uniformly distributed in a cylindrical volume with a radius a and height 2b. In accordance with the definition, a small spherical volume centered at the origin is excluded from the cylinder, as shown in Fig. A1. Without loss of generality, we assume that the origin of the coordinate system is located at  x . The principal value of the integral can be simplified as: where ε is the radius of the small exclusion,   R x x is the distance between the source point and the field point.
In the following discussion, we do not distinguish between the longitudinal and transverse waves and denote the wavenumber uniformly by k. When the fundamental results are obtained, we will attach the subscripts "L" or "T" to the wavenumbers of longitudinal and transverse waves again. For the convenience of the following calculation, we denote the spatial variables x1, x2 and x3 as x, y and z. From Eq. (A5) we know that the principal value of the integral involving 33 To further calculate these integrals, the integral domain is subdivided into three parts: VI, VII and VIII. VI is the outer cylindrical shell formed by the inner cylindrical surface r=ε, outer cylindrical surface r=a and the end surfaces z=±b, VII is enclosed by the inner cylindrical surface r=ε, top hemisphere and the top end z=b, VIII is the symmetric counterpart of VII, as shown in Fig. A1(b).  In this section, we will show that the singularity can also be obtained by setting the zero-frequency limit of the convolution between 0 . . ( ) ij P S E  k and the power spectrum function to be zero, i.e.: Substitution of (A33) into (A34) yields: . We assume that the medium has a small damping, so the poles deviate from the real axis, as shown in Fig. A2. Considering the condition that realistic waves must be out-going and with finite amplitude, we need to take the integral path formed by the upper semi-circle and the real axis to complete the integrals, so 2 Res( , ) 2 Res( , )  Comparing the three methods, we see Method 3 is the simplest way to obtain the singularity of the dynamic Green's tensor.