Scattering of elastic waves in heterogeneous media has become one of the most important problems in the field of wave propagation due to its broad applications in seismology, natural resource exploration, ultrasonic nondestructive evaluation and biomedical ultrasound. Nevertheless, it is one of the most challenging problems because of the complicated medium inhomogeneity and the complexity of the elastodynamic equations. A widely accepted model for the propagation and scattering of elastic waves, which properly incorporates the multiple scattering phenomenon and the statistical information of the inhomogeneities is still missing. 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 under the first-order smoothing approximation. The model establishes an accurate quantitative relation between the microstructural properties and the coherent wave propagation parameters and can be used for characterization or inversion of microstructures. Starting from the elastodynamic differential equations, a system of integral equation for the Green functions of the heterogeneous medium was developed by using Green’s functions of a homogeneous reference medium. After properly eliminating the singularity of the Green tensor and introducing a new set of renormalized field variables, the original integral equation is reformulated into a system of renormalized integral equations. Dyson’s equation and its first-order smoothing approximation, describing the ensemble averaged response of the heterogeneous system, are then derived with the aid of Feynman’s diagram technique. The dispersion equations for the longitudinal and transverse coherent waves are then obtained by applying Fourier transform to the Dyson equation. The exact solution to the dispersion equations are obtained numerically. To validate the new model, the results for weak-property-fluctuation materials are compared to the predictions given by an improved weakfluctuation multiple scattering theory. It is shown that the new model is capable of giving a more robust and accurate prediction of the dispersion behavior of weak-property-fluctuation materials. Numerical results further show that the new model is still able to provide accurate results for strong-property-fluctuation materials while the weak-fluctuation model is completely failed. As applications of the new model, dispersion and attenuation curves for coherent waves in the Earth’s 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 Mohorovičić 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.

Keywords:Multiple scattering, Mohorovičić Discontinuity, Elastic Continua.

Scattering of elastic waves in heterogeneous media has become one of the most important problems in the field of wave propagation due to its broad applications in seismology, natural resource exploration, ultrasonic nondestructive evaluation and biomedical ultrasound. Nevertheless, it is one of the most challenging problems because of the complicated medium inhomogeneity and the complexity of the elastodynamic equations. A widely accepted model for the propagation and scattering of elastic waves, which properly incorporates the multiple scattering phenomenon and the statistical information of the inhomogeneities is still missing. 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 under the first-order smoothing approximation. The model establishes an accurate quantitative relation between the microstructural properties and the coherent wave propagation parameters and can be used for characterization or inversion of microstructures. Starting from the elastodynamic differential equations, a system of integral equation for the Green functions of the heterogeneous medium was developed by using Green’s functions of a homogeneous reference medium. After properly eliminating the singularity of the Green tensor and introducing a new set of renormalized field variables, the original integral equation is reformulated into a system of renormalized integral equations. Dyson’s equation and its first-order smoothing approximation, describing the ensemble averaged response of the heterogeneous system, are then derived with the aid of Feynman’s diagram technique. The dispersion equations for the longitudinal and transverse coherent waves are then obtained by applying Fourier transform to the Dyson equation. The exact solution to the dispersion equations are obtained numerically. To validate the new model, the results for weak-property-fluctuation materials are compared to the predictions given by an improved weakfluctuation multiple scattering theory. It is shown that the new model is capable of giving a more robust and accurate prediction of the dispersion behavior of weak-property-fluctuation materials. Numerical results further show that the new model is still able to provide accurate results for strong-property-fluctuation materials while the weak-fluctuation model is completely failed. As applications of the new model, dispersion and attenuation curves for coherent waves in the Earth’s 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 Mohorovičić 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.

Keywords:Multiple scattering, Mohorovičić Discontinuity, Elastic Continua.

Scattering of elastic waves is the central topic in the theory of elastodynamics and its applications cover a broad spectrum of physical and engineering disciplines, ranging from seismology, natural resource exploration, ultrasonic nondestructive evaluation to biomedical ultrasound. Multiple scattering is a universal phenomenon which occurs in any materials that exhibit spatial fluctuation of elastic moduli or density. The planet Earth is a typical heterogeneous media. Geophysical studies show the rocks constituting the lithosphere, typically with a characteristic size from several hundred meters to tens of kilometers, exhibit strong fluctuations in both the elastic stiffness and the density. A seismic wave gets scattered by large amounts of such inhomogeneities during its travel from the source, normally tens of kilometers beneath the Earth’s surface, to the observation station located up to thousands of kilometers away from the epicenter. As a consequence, the seismic signal experiences multiple scattering and appears as a wavetrain consists of several dispersive and attenuating direct arrivals followed by a long duration of coda waves. Multiple scattering theory is a necessary tool for the explanation of measured signals and extracting statistical information of the heterogeneous lithosphere. A series of classics in seismology [1-4] are devoted to seismic wave scattering and attenuation and numerous research papers have also been published [5-33]. Geophysicists also make use of man-made earthquakes and borehole explosions to explore oil and gas reservoir [34-36]. Porous rocks saturated by hydrocarbon-enriched liquid mixtures are the major constituent of oil and gas reservoir, of which the mass density and elastic modulus are significantly different from that of tight rocks. By analyzing the unique dispersion and attenuation of the artificial seismic waves, geophysicists can identify the reservoirs, infer its mineral constituents and evaluate reserves of the resources. Ultrasonic waves with a central frequency of several megahertz are extensively used in industry for nondestructive characterization of microstructures in heterogeneous materials, such as hightemperature alloys, composite materials or ceramics [37- 68]. The typical size of the microstructure in these materials, such as pores, grains, crystallites and microtextured regions, lies between a few micrometers to several millimeters, thus the wavelengths are comparable to the characteristic dimension of the microstructures. In this frequency regime, the wave-scatterer interaction is relatively strong and the multiple scattering theory is extremely important for the analysis of measured signals. In analogy to its industrial applications, ultrasound is also used in medicine to diagnose and monitor certain diseases, like cancer and osteoporosis [69-77]. Bone quantitative ultrasound is a typical example of such applications. Cortical bone is a typical two-phase elastic material composed of a porous solid frame saturated by marrow. Clinical studies reveal that the porosity and average pore size are the key microstructural parameters that control the risk of bone fracture. Consequently, there is an increasing interest in quantitatively characterizing the microstructures in cortical bone [70]. X-ray microtomography (XMT) and nuclear magnetic resonance (NMR) analysis show the diameters of the majority of pores lie between 20 μm and 300 μm, while the volume fraction of pores varies between 5% and 30% [78]. These two factors combine to determine that the propagation of ultrasound in cortical bone is dominated by multiple scattering.

From the above discussion, we see an accurate and quantitative description of multiple scattering is of vital importance to the success of ultrasound applications. Nevertheless, theoretical analysis of multiple scattering in heterogeneous media is exceedingly challenging. Complexities of the scattering problem arise from two aspects. First, the mass density and/or elastic stiffness vary randomly with spatial coordinates, thus the governing equations are stochastic differential equations, conventional mathematical physics approaches developed for deterministic equations become invalid for these materials. Second, there exist large amounts of interfaces among different phases, where the boundary conditions, i.e., the continuity of displacement and traction must be enforced. As a consequence, obtaining an accurate analytical solution to the complete set of boundary value problems quickly becomes intractable. Solving for the numerical solution to this problem is equally difficult. The material properties (density and/or elastic moduli) exhibit a sharp jump across the boundaries, and accurate description of the discontinuity is impossible for certain numerical method, like the finite-difference time domain approach. In addition, large quantities of irregular boundaries exist in the medium, for which an accurate depiction in numerical models is also challenging. Moreover, large amount of meshes or elements makes the numerical simulation of the wave scattering extremely time-consuming, especially for the case of strong-property fluctuation materials and high frequency signals. Finally, numerical dispersion and artificial excitation or consumption frequently contaminate the numerical algorithms and result in distorted signals. Therefore, significant efforts have been directed to developing approximate analytical models for wave propagation in heterogeneous media.

During the last century, pioneering scientists have made tremendous efforts in studying the multiple scattering of waves and proposed numerous theoretical models to explain the rich phenomena. Based on the fundamental assumptions and methodologies, theoretical models can be classified into two categories, one is called the macroscopic phenomenological models, the other is known as the scattering models. Phenomenological models are also known as homogenized effective medium models, which focus on searching for different homogenization schemes to obtain the effective constitutive properties. The most representative example of phenomenological models is developed by Biot in his pioneering research on fluid-saturated porous materials [79-81]. In Biot’s theory, two sets of field variables (displacement, strain and stress) are introduced to describe the motion of the solid phase and the fluid phase, respectively. The equations of motion of a representative volume element for the two-phase material were derived from the Lagrange equations of a specific kinetic energy density, in which a set of mass coefficients were introduced to account for the coupling interactions between the two phases. After taking divergence and curl of the governing equations, Biot obtained two independent equations for the dilatational waves and one for the transverse waves, by which he further predicted the existence of two longitudinal propagation modes, known as the Biot waves of the first and the second type, and a transverse mode. In 1980, Plona [82] conducted a series of ultrasonic immersion tests on a slab of porous elastic material. After analyzing the transmitted signals with different incident angles, he identified a slow longitudinal mode in addition to the fast longitudinal wave and the transverse wave. He concluded that these waves correspond to those predicted by Biot. Since then, Biot’s model becomes the standard method for dealing with wave propagation in fluid-saturated porous elastic materials and is widely used in geomechanics, hydrogeology, petroleum exploration and more recently in biomechanics [70,83-88]. However, it is generally acknowledged that Biot’s model have a number of drawbacks. One shortcoming is that the fluid viscosity is regarded as the sole mechanism for attenuation, while the portion contributed by scattering is completely neglected. It is reported that Biot’s model always underestimate the attenuation [86-88]. Another noteworthy defect is that a number of phenomenological parameters, such as the frequency correction factor of the viscosity coefficient and the structural factor are involved, which either have obscured physical meaning or the accurate values are difficult to measure experimentally. Furthermore, Biot’s model cannot predict the anomalous negative dispersion as observed in cortical and trabecular bone [70,73,89], and fluid-saturated sediment sand [90]. Additionally, numerical results show it cannot predict the velocity and attenuation accurately and simultaneously [70]. In order to remedy these deficiencies, researchers have proposed different versions of modified Biot’s model. Dvorkin and coworkers incorporate the effects of local squirt flow into Biot’s model and developed the so-called BISQ model [19-21]. The BISQ model achieves certain degree of success in regards to improving the estimation of seismic attenuation. Muller and Gurevich [26-27] considered the effects of both the waveinduced flow and scattering on the seismic attenuation and derived a statistical smoothing approximation theory based on the Biot equations. Johnson and coworkers introduced the concepts of dynamic permeability and dynamic tortuosity to describe the fluid-saturated pores in the solid frame [91]. It has been applied in analysis of ultrasound dispersion and attenuation of trabecular bones and achieved intermediate degree of success [70]. Another type of macroscopic phenomenological method is known as the dynamic selfconsistent theory. In analogous to the static self-consistent theory [92-95], in which the stress polarization as introduced, Willis and Sabino [59-64] adopted the concept of the momentum polarization for dynamic problems and derived a system of coupled integral equations for the ensemble-averaged displacement and strain. A set of selfconsistent conditions are derived to fix the properties of a homogeneous reference medium. Formulas for waves propagating in general two-phase composite materials are then specialized for the cases of aligned spheroidal inclusions and randomly oriented spheroids. Numerical results for velocity and attenuation shows that their model can successfully capture the negative dispersion, which is a common feature of two-phase materials. An obvious drawback of the model is that at high frequencies, the attenuation predicted by the model approaches zero, which is physically impractical. In light of Willis and Sabino’s work, Zhuck and Lakhtakia [96-97] incorporated the second-order statistics of the random medium and developed a system of renormalized integral equations for the estimation of the effective constitutive properties. In the proposed model, they first considered the singularity of the elasodynamic Green’s tensor, and thus their results are applicable for elastic materials with strong property fluctuation. As pointed out in their work, a major drawback of the homogenization approach is that it is only valid for heterogeneous materials in which the characteristic dimension of the heterogeneities is less than the wavelength of the excitation signals, i.e., at relatively low frequencies. Through the above discussion, it is seen different versions of effective medium models have several shortcomings in common. The most noteworthy is that they are valid either at relatively low frequencies or in a relatively narrow frequency band, none of them can predict the effective properties in the whole frequency range. This puts server limits on practical applications, for instance, analysis of high frequency signals. Another disadvantage is the introduction of phenomenological parameters which require to be measured experimentally or need to be known in advance. These limitations pose challenges for practical applications and consequently, they are rarely used as inversion models for microstructure characterization.

Contrary to phenomenological models, whose attention is concentrated on searching for an equivalent homogeneous medium, scattering models are focused on investigating the interaction of wave fields with randomly distributed small-scale scatterers. Parameters of practical importance, such as backscattering coefficients, scattering cross section, diffusivity and mean free path, in addition to the velocity and attenuation, are extracted through the analysis of the scattering process. At present, a vast variety of scattering theories and approximations have been proposed. All these scattering theories can be classified into several categories following different classification schemes. For instance, according to the frequency range in which they are applicable, scattering theories can be classified into low frequency, stochastic and high frequency theories. Based on the fractional fluctuation of the material properties of the constituent phases, scattering theories are also classified into weak scattering theories [2,16,48,98] and strong scattering theories [99-109]. In view of stochastic medium model, scattering models can also be classified into discrete scatterer models [110-119] and continuum models [2,16,48,98]. In order to obtain a comprehensive understanding of the background of this area, here we give a brief overview of the most popular models and approximations, and compare their advantages and disadvantages. The Born approximation is probably the most famous approximation used in the scattering community. When the fluctuations of the materials properties are weak, and the wavelength is large compared to the typical dimension of the heterogeneities, the field in the inclusions can be approximated by the unperturbed field. This approximation is called the Born approximation. The Born approximation is an ideal tool for dealing with weak scattering problems, and it has been used in nearly all involved disciplines, including seismology [1-3], ultrasonic NDE [37-41] and electromagnetic remote sensing [120-122]. Explicit expressions for measurable quantities like the scattering section and the backscattering coefficient of polycrystalline materials are derived in [37- 41]. Born approximation fails in the high frequency regime and for strong fluctuation materials. Complementary to the Born approximation, the Rytov approximation, the geometric approximation and the parabolic approximation [3,123-124] are developed to analyze scattering behavior in the high frequency regime, although they are still limited to weak scattering problems. When the concentration of the inclusions is low, or in other words, the solid is a dilute solution of the dispersed inclusions, a straightforward approximation is to omit the interactions among different inclusions, and assume that each scatterer interacts with the incident wave independently, this approximation is called the independent scattering approximation (ISA) [18, 33,41]. ISA is applicable to inclusions with weak property perturbation or predict the dispersion and attenuation accurately when the volumetric concentration of the pointscatterers is very small.

Multiple scattering theories account for the sophisticated scattered wavefield faithfully and aim at giving an accurate and consistent description of the propagation behaviors. Rossum and Nieuwenhuizen [125] and Barabanenkov et. al. [126] gave two excellent reviews for the historical development and the current state of general multiple scattering theories. The first multiple scattering theory for elastic waves was developed by Karal and Keller [98]. This model relies on the small-perturbation expansion of differential operators and their inverse operators. After successive iteration of the series expansion of the wave operator, the differential equations are transferred into a system of integral equations, from which the Christoffel equation for the coherent plane waves is then derived. Based on Keller’s approximation, they obtained analytical formulas for the velocity and attenuation which are valid in the whole frequency range. Stanke and Kino [46-47] developed a unified model for polycrystalline materials in parallel to the Karal-Keller approach. This model has been applied in nondestructive characterization of microstructures in polycrystalline alloys. Exact dispersion and attenuation of alloys with cubic crystallites are obtained in the whole frequency range. Weaver [48] developed the first multiple scattering theory based on Dyson’s equation and the firstorder-smoothing-approximation (FOSA) [127-128]. The Dyson equation of the ensemble-averaged Green’s function for polycrystalline medium was obtained by introducing a multiple scattering series expansion of the stochastic differential equations. By invoking the First-OrderSmoothing approximation (FOSA), the explicit expressions of the mass operators were derived for longitudinal and transverse coherent waves. Closed forms of longitudinal and transverse Green’s functions and dispersion equations were then obtained based on the FOSA Dyson’s equation. Because of the extreme complexity of the resulting expressions, the Born approximation was adopted in the last step and the closed-form analytical expressions for longitudinal and transverse wave attenuation were obtained. As pointed out in [48], the final expressions are only valid for weak fluctuation materials and low frequency waves. Since then, Weaver’s model has become one of the most widely used model in the nondestructive evaluation community. Weaver’s original work considered untextured polycrystals with spherical cubic grains only. Turner [57] extended Weaver’s model to incorporate the effects of textures. He studied equal-axed cubic polycrystal with the crystallographic axis aligned along a preferred direction. Green’s functions for transversely isotropic reference medium are incorporated to account for the macroscopically transverse isotropy of the textured polycrystals. Turner and Anugonda further considered two-phase materials with microscopically isotropic components [58]. It is noted that both [57] and [58] focus on the attenuation under the Born approximation, while the exact solution to the exact dispersion equation is still missing. Moreover, it is seen that only one root for attenuation is found for all propagation directions. Calvet and Margerin calculated both the velocity and the Q-factors of polycrystalline materials with cubic and hexagonal crystallites using Weaver’s model, with an ultimate goal to analyze the grain anisotropy of the Earth’s uppermost inner core [13]. All the previous calculations in the framework of Weaver’s model adopt the Born approximation and valid in the low frequency regime only. In order to obtain the propagation characteristics at high frequencies, Calvet and Margerin [12] proposed a spectral function approach, in which the imaginary part of the ensemble averaged Green’s function derived in Weaver’s model is identified as the spectral function. By introducing a simplified spectral function of the coherent waves and invoking the method of least square fitting, the dispersion and attenuation of the coherent wave are reconstructed from the exact spectral function. They first discovered that at high frequencies there exist two sets of longitudinal and transverse modes, although at the end of the paper they questioned the discovery by providing a subtle example for which there should be one propagation mode but their calculations gave two. Calvet and Margerin further extended Weaver’s model to incorporate the effects of grain shape [14]. An interesting phenomenon is that the attenuation and dispersion exhibit obvious anisotropy as a result of the geometric anisotropy. A noteworthy feature of the spectral function approach is that it is only an approximation to the 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-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 T-matrix 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 table 1.

^aWeaver’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.

^bStanke-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 table 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 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-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, including Green’s function of the heterogeneous medium and its integral representation, singularity of the Green tensor, derivation of the renormalized Dyson’s equation using Feynman’s diagram and the first-ordersmoothing-approximation, solution of the Dyson equation and the dispersion equations. In Section III, the exact solutions for the dispersion and attenuation of longitudinal and transverse waves will be discussed. The advantages of the new model will be demonstrated through comparison of the numerical results with that calculated from a multiple scattering theory for weak-property-fluctuation materials. Section IV presents a series of practical applications of the new model in seismology, ultrasonic NDE and bone quantitative ultrasound. To show the power of the model, we calculate the exact dispersion and Q-factor curves of seismic waves propagating in the Earth’s lithosphere, the dispersion and attenuation of ultrasonic waves in porous aluminum and Cu-Al alloy, and the dispersion and attenuation curves of ultrasound in human cortical bone. Through comprehensive parametric study, we show that the new model can give a concise and consistent explanation for a broad spectrum of observed scattering phenomena, such as the Q-factors of seismic waves, negative dispersion in bone, spectral splitting and occurrence of two modes. In section V, we give a brief discussion on how to perform successful numerical and experimental verifications of the model. Future work on extension of the current model to incorporate more complicated microstructures are also presented. Finally, the unique features of the new model are summarized and the major conclusions are highlighted in Section VI.

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 frequency-wavenumber 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.

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:

${[c_{ijkl}(x)u_{k,l}]}_{,j}+\rho (x){\omega }^2{u}_i=0,$ (1)

where ω is the circular frequency, u_i denotes the displacement components, ρ(x) and c_ijkl(x) are the mass density and elastic stiffness. The elastic stiffness has the following symmetries:

$c_{ijkl}(x)=c_{jikl}(x)=c_{ijlk}(x)=c_{klij}(x).$ (2)

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 figure 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 e_a, where e_a. 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:

$c_{ijkl}^{(0)}{G}_{k{a}^{\prime },lj}+{\rho }^{(0)}{\omega }^2{G}_{i{a}^{\prime }}+F{a}_{i{a}^{\prime }}\delta (x-x^{\prime })=0,$ in V_o (3)

$c_{ijkl}^{(1)}{G}_{k{a}^{\prime },lj}+{\rho }^{(1)}{\omega }^2{G}_{i{a}^{\prime }}+F{a}_{i{a}^{\prime }}\delta (x-x^{\prime })=0,$ in V₁ (4)

where V₀ and V₁ represent the volumes occupied by the substrate and the inclusion, respectively, and a_iaʹ denotes the directional cosine of e_aʹ, a_ia=cos(e_a,e_i).

Suppose the inclusion and the substrate are perfectly bonded, so the stress and displacement are continuous at the boundary of the inclusion, thus we have on S.

$c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}(x\to S^{+})N_j=c_{ijkl}^{(1)}{G}_{k{a}^{\prime },l}(x\to S^{-})N_j,$ $G_{k{a}^{\prime }}(x\to S^{+})=G_{k{a}^{\prime }}(x\to S^{-}),$ on S. (5)

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 figure 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:

$c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},lj}^{(0)}+F{a}_{i{\alpha }^{″}}\delta (x-x^{″})+{\rho }^{(0)}{\omega }^2{G}_{i{\alpha }^{″}}^{(0)}=0.$ (6)

Multiplying Eq. (6) and Eq. (3) by G_ia, and G_iα⁽⁰⁾, respectively, and then subtracting the resulting equations, we obtain:

$c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},lj}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(0)}{G}_{k{a}^{\prime },lj}{G}_{i{\alpha }^{″}}^{(0)}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0.$ (7)

Multiplying Eq. (6) and Eq. (4) by G_ia, and G_iα⁽⁰⁾, respectively, and then subtracting the resulting equations, we obtain:

$c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},lj}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(1)}{G}_{k{a}^{\prime },lj}{G}_{i{\alpha }^{″}}^{(0)}+{\rho }^{(0)}{\omega }^2{G}_{i{\alpha }^{″}}^{(0)}{G}_{i{a}^{\prime }}-{\rho }^{(1)}{\omega }^2{G}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0.$ (8)

Appling product rule of partial derivatives to Eqs. (7) and (8), we get:

$$\begin{array}{l}{c}_{ijkl}^{(0)}{[}_{G_{k{\alpha }^{″},l}^{(0)}}-c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime },j}-{[}_{c_{ijkl}^{(0)}}\\ +c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0\end{array}$$ ,(9)

$\begin{array}{l}{c}_{ijkl}^{(0)}{[}_{G_{k{\alpha }^{″},l}^{(0)}}-c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime },j}-{[}_{c_{ijkl}^{(1)}}+c_{ijkl}^{(1)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}\\ +{\rho }^{(0)}{\omega }^2{G}_{i{\alpha }^{″}}^{(0)}{G}_{i{a}^{\prime }}-{\rho }^{(1)}{\omega }^2{G}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0.\end{array}$ (10)

Considering the symmetry of the elastic stiffness, Eqs. (9) and (10) can be simplified as:

$c_{ijkl}^{(0)}{[G_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime }}]}_{,j}-{[c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″}}^{(0)}]}_{,j}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0,$ (11)

$\begin{array}{l}{c}_{ijkl}^{(0)}{[}_{G_{k{\alpha }^{″},l}^{(0)}}-{[}_{c_{ijkl}^{(1)}}+(c_{ijkl}^{(1)}-c_{ijkl}^{(0)})G_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}\\ -({\rho }^{(1)}-{\rho }^{(0)}){\omega }^2{G}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}+F{a}_{i{\alpha }^{″}}{G}_{i{a}^{\prime }}\delta (x-x^{″})-F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })=0.\end{array}$ (12)

Integrating Eq. (11) over the domain V₀, and then applying Gauss’s divergence theorem, we obtain:

$\begin{array}{l}{\displaystyle {\iiint }_{V_0}F{a}_{i{\alpha }^{″}}\delta (x-x^{\prime })G_{i{a}^{\prime }}dV}-{\displaystyle {\iiint }_{V_0}F{a}_{i{a}^{\prime }}\delta (x-x^{\prime })G_{i{\alpha }^{″}}dV}\\ +{\displaystyle {∯}_{S^{\infty }}[c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″}}^{(0)}]}{N}_{j}dS+{\displaystyle {∯}_{S_1}[c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″}}^{(0)}]}(-N_j)dS=0.\end{array}$ (13)

A negative sign is introduced in front of N_j in the integral over S₁ because the out normal of volume V₀ corresponds to the negative of N, which is shown in figure 1. S^∞ represents an imaginary surface that lies infinitely far from the source and the inclusion.

Integrating Eq. (12) over the volume V₁ and applying Gauss’s divergence theorem, we obtain:

$\begin{array}{l}{\displaystyle {\iiint }_{V_1}F{a}_{i{\alpha }^{″}}\delta (x-x^{\prime })G_{i{a}^{\prime }}dV}-{\displaystyle {\iiint }_{V_1}F{a}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}\delta (x-x^{\prime })dV}+{\displaystyle {∯}_{S_1}[c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(1)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″}}^{(0)}]}{N}_{j}dS\\ +{\displaystyle {\iiint }_{V_1}[(c_{ijkl}^{(1)}-c_{ijkl}^{(0)})G_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}-({\rho }^{(1)}-{\rho }^{(0)}){\omega }^2{G}_{i{a}^{\prime }}{G}_{i{\alpha }^{″}}^{(0)}]dV}=0.\end{array}$ (14)

Adding Eqs. (13) and (14) and then invoking the boundary conditions Eq. (5), we have:

$\begin{array}{l}{\displaystyle {\iiint }_{V_0+V_1}F{a}_{i{\alpha }^{″}}\delta (x-x^{″})G_{i{a}^{\prime }}dV}-{\displaystyle {\iiint }_{V_0+V_1}F{a}_{i{a}^{\prime }}\delta (x-x^{\prime })G_{i{\alpha }^{″}}^{(0)}dV}+{\displaystyle {∯}_{S^{\infty }}[c_{ijkl}^{(0)}{G}_{k{\alpha }^{″},l}^{(0)}{G}_{i{a}^{\prime }}-c_{ijkl}^{(0)}{G}_{k{a}^{\prime },l}{G}_{i{\alpha }^{″}}^{(0)}]}{N}_{j}dS\\ +{\displaystyle {\iiint }_{V_1}[(c_{ijkl}^{(1)}-c_{ijkl}^{(0)})G_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}-({\rho }^{(1)}-{\rho }^{(0)}){\omega }^2{G}_{i{a}^{\prime }}{G}_{i\alpha }^{(0)}]dV}=0.\end{array}$ (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 O"x₁"x₂"x₃" and O"x₁"x₂"x₃" coincide with that of the coordinate Ox₁x₂x₃, thus a_ia′ = δ_ia′., Eq. (15) can be simplified as:

$F{G}_{{\alpha }^{″}{a}^{\prime }}-F{G}_{a^{\prime }{\alpha }^{″}}^{(0)}+{\displaystyle {\iiint }_{V_1}[(c_{ijkl}^{(1)}-c_{ijkl}^{(0)})G_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}-({\rho }^{(1)}-{\rho }^{(0)}){\omega }^2{G}_{i{a}^{\prime }}{G}_{i\alpha }^{(0)}]dV}=0.$ (16)

Finally, we obtain the integral representation of the inhomogeneous Green function expressed in terms of the reference homogeneous Green’s function,

$G_{{\alpha }^{″}{a}^{\prime }}=G_{a^{\prime }{\alpha }^{″}}^{(0)}+\frac{1}{F}{\displaystyle {\iiint }_{V_1}[-(c_{ijkl}^{(1)}-c_{ijkl}^{(0)})G_{k{a}^{\prime },l}{G}_{i{\alpha }^{″},j}^{(0)}+({\rho }^{(1)}-{\rho }^{(0)}){\omega }^2{G}_{i{a}^{\prime }}{G}_{i\alpha }^{(0)}]dV}.$ (17)

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,

$c_{ijkl}^{(1)}=c_{ijkl}^{(0)}+\delta c_{ijkl}^{(1)},$ ${\rho }^{(1)}={\rho }^{(0)}+\delta {\rho }^{(1)}.$ (18)

Substitution of (18) into (17) yields:

$$G_{{\alpha }^{″}{a}^{\prime }}(x^{″},x^{\prime })=G_{a^{\prime }{\alpha }^{″}}^{(0)}(x^{\prime },x^{″})+\frac{1}{F}{\displaystyle {\iiint }_{V_1}\begin{array}{l}[-\delta c_{ijkl}^{(1)}{G}_{k{a}^{\prime },l}(x,x^{\prime })G_{i{\alpha }^{″},j}^{(0)}(x,x^{″})\\ +\delta {\rho }^{(1)}{\omega }^2{G}_{i{a}^{\prime }}(x,x^{\prime })G_{i{\alpha }^{″}}^{(0)}(x,x^{″})]dV\end{array}},$$ (19)

where Green’s function of the homogeneous medium is given by:

$G_{i{a}^{\prime }}^0(x,x^{\prime },\omega )=\frac{F}{4\pi }\frac{{\delta }_{i{\alpha }^{\prime }}}{\mu }\frac{e^{i{k}_T\left|x-x^{\prime }\right|}}{\left|x-x^{\prime }\right|}+\frac{F}{4\pi \rho {\omega }^2}\frac{{\partial }^2}{\partial x_i\partial {x^{\prime }}_a}\left[\frac{e^{i{k}_L\left|x-x^{\prime }\right|}}{\left|x-x^{\prime }\right|}-\frac{e^{i{k}_T\left|x-x^{\prime }\right|}}{\left|x-x^{\prime }\right|}\right].$ (20)

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:

$G_{i{a}^{\prime }}^0(x,x^{\prime })=G_{a^{\prime }i}^0(x,x^{\prime })=G_{a^{\prime }i}^0(x^{\prime },x)=G_{i{a}^{\prime }}^0(x^{\prime },x).$ (21)

Applying these relations to Eq. (19), we obtain:

$G_{{\alpha }^{″}{a}^{\prime }}(x^{″},x^{\prime })=G_{{\alpha }^{″}{a}^{\prime }}^{(0)}(x^{″},x^{\prime })+\frac{1}{F}{\displaystyle {\iiint }_{V_1}\begin{array}{l}[-\delta c_{ijkl}^{(1)}{G}_{k{a}^{\prime },l}(x,x^{\prime })G_{{\alpha }^{″}i,j}^{(0)}(x^{″},x)\\ +\delta {\rho }^{(1)}{\omega }^2{G}_{i{a}^{\prime }}(x,x^{\prime })G_{{\alpha }^{″}i}^{(0)}(x^{″},x^{\prime })]dV\end{array}}.$ (22)

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 figure 2.

Obviously, the above derivation can be straightforwardly extended to inhomogeneous medium with multiple inclusions, and the resulting equation is:

$G_{{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })=G_{{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })+\frac{1}{F}{\displaystyle \sum _{n=1}^N{\displaystyle {\int }_{V_n}\begin{array}{l}[\delta {\rho }^{(n)}{\omega }^2{G}_{i{a}^{\prime }}(x,x^{\prime })G_{{\beta }^{″}i}^0(x^{″},x)\\ -\delta c_{ijkl}^{(n)}{G}_{k{a}^{\prime },l}(x,x^{\prime })G_{{\beta }^{″}i,j}^0(x^{″},x)]dV\end{array}}},$ (23)

where δρ⁽ⁿ⁾ and $\delta c_{ijkl}^{(n)}$ represent the material property fluctuation of the n-th inclusion.

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-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:

$G_{{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })=G_{{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })+\frac{1}{F}{\displaystyle {\int }_V[\delta \rho (x){\omega }^2{G}_{i{a}^{\prime }}(x,x^{\prime })G_{{\beta }^{″}i}^0(x^{″},x)-\delta c_{ijkl}(x)G_{k{a}^{\prime },l}(x,x^{\prime })G_{{\beta }^{″}i,j}^0(x^{″},x)]dV}.$ (24)

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.

Integral representation of the perturbed field:Equation (24) is the governing equation for the displacement Green’s function. To complete the description of the elastic field, we also need the strain Green’s function. The strain tensor of a general elastic medium is defined by:

${\epsilon }_{{\alpha }^{″}{\beta }^{″}}=\frac{1}{2}(u_{{\alpha }^{″},{\beta }^{″}}+u_{{\beta }^{″},{\alpha }^{″}}).$ (25)

Consequently, the Green function for strain is given by:

$\begin{array}{l}\frac{1}{2}[G_{{\beta }^{″}{a}^{\prime },{\alpha }^{″}}(x^{″},x^{\prime })+G_{{\alpha }^{″}{a}^{\prime },{\beta }^{″}}(x^{″},x^{\prime })]=\frac{1}{2}[G_{{\beta }^{″}{a}^{\prime },{\alpha }^{″}}^0(x^{″},x^{\prime })+G_{{\alpha }^{″}{a}^{\prime },{\beta }^{″}}^0(x^{″},x^{\prime })]\\ +\frac{1}{F}{\displaystyle {\int }_V\left[\delta \rho (x){\omega }^2\frac{1}{2}[G_{{\beta }^{″}i,{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,{\beta }^{″}}^0(x^{″},x)]G_{i{a}^{\prime }}(x,x^{\prime })-\frac{1}{2}[G_{{\beta }^{″}i,j{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,j{\beta }^{″}}^0(x^{″},x)]\delta c_{ijkl}(x)G_{k{a}^{\prime },l}(x,x^{\prime })\right]dV}.\end{array}$ (26)

Considering the elastic stiffness variation tensor δc_ijkl (x) has the same symmetry as c_ijkl (x), we can write (24) and (26) as

$\begin{array}{l}{G}_{{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })=G_{{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })+\frac{1}{F}{\displaystyle {\int }_V\{\delta \rho (x){\omega }^2{G}_{{\beta }^{″}i}^0(x^{″},x)G_{i{a}^{\prime }}(x,x^{\prime })}\\ -\frac{1}{2}[G_{{\beta }^{″}i,j}^0(x^{″},x)+G_{{\beta }^{″}j,i}^0(x^{″},x)]\delta c_{ijkl}(x)\frac{1}{2}[G_{k{a}^{\prime },l}(x,x^{\prime })+G_{l{a}^{\prime },k}(x,x^{\prime })]\}dV,\end{array}$ (27)

$\begin{array}{l}\frac{1}{2}[G_{{\beta }^{″}{a}^{\prime },{\alpha }^{″}}(x^{″},x^{\prime })+G_{{\alpha }^{″}{a}^{\prime },{\beta }^{″}}(x^{″},x^{\prime })]\\ =\frac{1}{2}[G_{{\beta }^{″}{a}^{\prime },{\alpha }^{″}}^0(x^{″},x^{\prime })+G_{{\alpha }^{″}{a}^{\prime },{\beta }^{″}}^0(x^{″},x^{\prime })]+\frac{1}{F}{\displaystyle {\int }_V\{\delta \rho (x){\omega }^2\frac{1}{2}[G_{{\beta }^{″}i,{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,{\beta }^{″}}^0(x^{″},x)]G_{i{a}^{\prime }}(x,x^{\prime })}\\ -\frac{1}{4}[G_{{\beta }^{″}i,j{\alpha }^{″}}^0(x^{″},x)+G_{{\beta }^{″}j,i{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,j{\beta }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}j,i{\beta }^{″}}^0(x^{″},x)]\delta c_{ijkl}(x)\frac{1}{2}[G_{k{a}^{\prime },l}(x,x^{\prime })+G_{l{a}^{\prime },k}(x,x^{\prime })]\}dV.\end{array}$ (28)

Equations (27) and (28) can be expressed in a matrix form:

$\left[\begin{array}{l}{G}_{{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })\end{array}\right]=\left[\begin{array}{l}{G}_{{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })\end{array}\right]+\frac{1}{F}{\displaystyle {\int }_V\left[\begin{array}{l}{G}_{{\beta }^{″}i}^0(x^{″},x)-{\epsilon }_{{\beta }^{″}ij}^0(x^{″},x)\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}i}^0(x^{″},x)-E_{{\alpha }^{″}{\beta }^{″}ij}^0(x^{″},x)\end{array}\right]\left[\begin{array}{l}\delta \rho (x){\omega }^2{\delta }_{ij}0\\ 0\delta c_{ijkl}(x)\end{array}\right]\left[\begin{array}{l}{G}_{j{a}^{\prime }}(x,x^{\prime })\\ {\epsilon }_{kl{a}^{\prime }}(x,x^{\prime })\end{array}\right]dV},$ (29)

where

${\epsilon }_{{\beta }^{″}ij}^0(x^{″},x)=\frac{1}{2}[G_{{\beta }^{″}i,j}^0(x^{″},x)+G_{{\beta }^{″}j,i}^0(x^{″},x)]$, ${\epsilon }_{{\alpha }^{″}{\beta }^{″}i}^0(x^{″},x)=\frac{1}{2}[G_{{\beta }^{″}i,{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,{\beta }^{″}}^0(x^{″},x)]$ (30)

$E_{{\alpha }^{″}{\beta }^{″}ij}^0(x^{″},x)=\frac{1}{4}[G_{{\beta }^{″}i,j{\alpha }^{″}}^0(x^{″},x)+G_{{\beta }^{″}j,i{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}i,j{\beta }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}j,i{\beta }^{″}}^0(x^{″},x)]$ (31)

${\epsilon }_{{\alpha }^{″}{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })=\frac{1}{2}[G_{{\beta }^{″}{a}^{\prime },{\alpha }^{″}}(x^{″},x^{\prime })+G_{{\alpha }^{″}{a}^{\prime },{\beta }^{″}}(x^{″},x^{\prime })]$, ${\epsilon }_{kl{a}^{\prime }}(x,x^{\prime })=\frac{1}{2}[G_{k{a}^{\prime },l}(x,x^{\prime })+G_{l{a}^{\prime },k}(x,x^{\prime })]$. (32)

Note the homogeneous Green’s function G^o_β"i(xⁿ,x) given by Eq. (20) can be rewritten as G^o_β"i(xⁿ,x), where:

$G_{{\beta }^{″}{i}^{″}}^0(x^{″},x,\omega )=\frac{F}{4\pi }\frac{{\delta }_{\beta i}}{\mu }\frac{e^{i{k}_T\left|x^{″}-x\right|}}{\left|x^{″}-x\right|}-\frac{F}{4\pi \rho {\omega }^2}\frac{{\partial }^2}{\partial {x^{″}}_i\partial {x^{″}}_{\beta }}\left[\frac{e^{i{k}_L\left|x^{″}-x\right|}}{\left|x^{″}-x\right|}-\frac{e^{i{k}_T\left|x^{″}-x\right|}}{\left|x^{″}-x\right|}\right].$ (33)

From Eq. (33) one can easily find:

$\frac{\partial G_{{\beta }^{″}{i}^{″}}^0}{\partial x_j}=-\frac{\partial G_{{\beta }^{″}{i}^{″}}^0}{\partial {x^{″}}_j}.$ (34)

Consequently, we have the following identities:

${\epsilon }_{{\beta }^{″}ij}^0(x^{″},x)=-{\epsilon }_{{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x),$ ${\epsilon }_{{\alpha }^{″}{\beta }^{″}i}^0(x^{″},x)={\epsilon }_{{\alpha }^{″}{\beta }^{″}{i}^{″}}^0(x^{″},x),$ $E_{{\alpha }^{″}{\beta }^{″}ij}^0(x^{″},x)=-E_{{\alpha }^{″}{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x),$ (35)

where

${\epsilon }_{{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x)=\frac{1}{2}[G_{{\beta }^{″}{i}^{″},j^{″}}^0(x^{″},x)+G_{{\beta }^{″}{j}^{″},i^{″}}^0(x^{″},x)],$ ${\epsilon }_{{\alpha }^{″}{\beta }^{″}{i}^{″}}^0(x^{″},x)=\frac{1}{2}[G_{{\beta }^{″}{i}^{″},{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}{i}^{″},{\beta }^{″}}^0(x^{″},x)],$ (36)

$E_{{\alpha }^{″}{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x)=\frac{1}{4}[G_{{\beta }^{″}{i}^{″},j^{″}{\alpha }^{″}}^0(x^{″},x)+G_{{\beta }^{″}{j}^{″},i^{″}{\alpha }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}{i}^{″},j^{″}{\beta }^{″}}^0(x^{″},x)+G_{{\alpha }^{″}{j}^{″},i^{″}{\beta }^{″}}^0(x^{″},x)].$ (37)

Substituting F=1 newton into Eq. (29) and considering the relations given by Eq. (35), we can simplify Eq. (29) as:

$\begin{array}{l}\left[\begin{array}{l}{G}_{{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}{a}^{\prime }}(x^{″},x^{\prime })\end{array}\right]=\left[\begin{array}{l}{G}_{{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}{a}^{\prime }}^0(x^{″},x^{\prime })\end{array}\right]\\ +{\displaystyle {\int }_V\left[\begin{array}{l}{G}_{{\beta }^{″}{i}^{″}}^0(x^{″},x){\epsilon }_{{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x)\\ {\epsilon }_{{\alpha }^{″}{\beta }^{″}{i}^{″}}^0(x^{″},x)E_{{\alpha }^{″}{\beta }^{″}{i}^{″}{j}^{″}}^0(x^{″},x)\end{array}\right]\left[\begin{array}{l}\delta \rho (x){\omega }^2{\delta }_{ij}0\\ 0\delta c_{ijkl}(x)\end{array}\right]\left[\begin{array}{l}{G}_{j{a}^{\prime }}(x,x^{\prime })\\ {\epsilon }_{kl{a}^{\prime }}(x,x^{\prime })\end{array}\right]dV}.\end{array}$ (38)

For the convenience of subsequent discussion, we rewrite Eq. (38) in a more compact form:

Ψ(x^" - x^')=Ψ°(x^" - x^') +∫∫∫_V(x)Γ(x^" - x)Π(x)Ψ(x- x^')d³x, (39)

where

$\Psi (x^{″}-x^{\prime })=\left[\begin{array}{ccc}{G}_{11}& G_{12}& G_{13}\\ G_{21}& G_{22}& G_{23}\\ G_{31}& G_{32}& G_{33}\\ G_{11,1}& G_{12,1}& G_{13,1}\\ G_{21,2}& G_{22,2}& G_{23,2}\\ G_{31,3}& G_{32,3}& G_{33,3}\\ G_{21,3}+G_{31,2}& G_{22,3}+G_{32,2}& G_{23,3}+G_{33,2}\\ G_{11,3}+G_{31,1}& G_{12,3}+G_{32,1}& G_{13,3}+G_{33,1}\\ G_{21,1}+G_{11,2}& G_{22,1}+G_{12,2}& G_{23,1}+G_{13,2}\end{array}\right],{\Psi }^0(x^{″}-x^{\prime })=\left[\begin{array}{ccc}{G}_{11}^0& G_{12}^0& G_{13}^0\\ G_{21}^0& G_{22}^0& G_{23}^0\\ G_{31}^0& G_{32}^0& G_{33}^0\\ G_{11,1}^0& G_{12,1}^0& G_{13,1}^0\\ G_{21,2}^0& G_{22,2}^0& G_{23,2}^0\\ G_{31,3}^0& G_{32,3}^0& G_{33,3}^0\\ G_{21,3}^0+G_{31,2}^0& G_{22,3}^0+G_{32,2}^0& G_{23,3}^0+G_{33,2}^0\\ G_{11,3}^0+G_{31,1}^0& G_{12,3}^0+G_{32,1}^0& G_{13,3}^0+G_{33,1}^0\\ G_{21,1}^0+G_{11,2}^0& G_{22,1}^0+G_{12,2}^0& G_{23,1}^0+G_{13,2}^0\end{array}\right],$ (40)

$\Pi (x)=\left[\begin{array}{ccccccccc}\delta \rho {\omega }^2& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& \delta \rho {\omega }^2& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& \delta \rho {\omega }^2& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& \delta c_{11}& \delta c_{12}& \delta c_{13}& \delta c_{14}& \delta c_{15}& \delta c_{16}\\ 0& 0& 0& \delta c_{12}& \delta c_{22}& \delta c_{23}& \delta c_{24}& \delta c_{25}& \delta c_{26}\\ 0& 0& 0& \delta c_{13}& \delta c_{23}& \delta c_{33}& \delta c_{34}& \delta c_{35}& \delta c_{36}\\ 0& 0& 0& \delta c_{14}& \delta c_{24}& \delta c_{34}& \delta c_{44}& \delta c_{45}& \delta c_{46}\\ 0& 0& 0& \delta c_{15}& \delta c_{25}& \delta c_{35}& \delta c_{45}& \delta c_{55}& \delta c_{56}\\ 0& 0& 0& \delta c_{16}& \delta c_{26}& \delta c_{36}& \delta c_{46}& \delta c_{56}& \delta c_{66}\end{array}\right],$ (41)

$\Gamma (x^{″}-x)=\left[\begin{array}{cc}A& B\\ B^T& D\end{array}\right],$ (42)

and

$A=\left[\begin{array}{ccc}{G}_{11}^0& G_{12}^0& G_{13}^0\\ G_{21}^0& G_{22}^0& G_{23}^0\\ G_{31}^0& G_{32}^0& G_{33}^0\end{array}\right],B=\left[\begin{array}{cccccc}{G}_{11,1}^0& G_{12,2}^0& G_{13,3}^0& G_{12,3}^0+G_{13,2}^0& G_{11,3}^0+G_{13,1}^0& G_{12,1}^0+G_{11,2}^0\\ G_{21,1}^0& G_{22,2}^0& G_{23,3}^0& G_{22,3}^0+G_{23,2}^0& G_{21,3}^0+G_{23,1}^0& G_{22,1}^0+G_{21,2}^0\\ G_{31,1}^0& G_{32,2}^0& G_{33,3}^0& G_{32,3}^0+G_{33,2}^0& G_{31,3}^0+G_{33,1}^0& G_{32,1}^0+G_{31,2}^0\end{array}\right],$ (43)

$D=\left[\begin{array}{l}{G}_{11,11}^0{G}_{12,21}^0{G}_{13,31}^0{G}_{12,31}^0+G_{13,21}^0{G}_{11,31}^0+G_{13,11}^0{G}_{12,11}^0+G_{11,21}^0\\ G_{21,12}^0{G}_{22,22}^0{G}_{23,32}^0{G}_{22,32}^0+G_{23,22}^0{G}_{21,32}^0+G_{23,12}^0{G}_{22,12}^0+G_{21,22}^0\\ G_{31,13}^0{G}_{32,23}^0{G}_{33,33}^0{G}_{32,33}^0+G_{33,23}^0{G}_{31,33}^0+G_{33,13}^0{G}_{32,13}^0+G_{31,23}^0\\ G_{21,13}^0+G_{31,12}^0{G}_{22,23}^0+G_{32,22}^0{G}_{23,33}^0+G_{33,32}^0{G}_{22,33}^0+G_{33,22}^0+2G_{23,23}^0{G}_{21,33}^0+G_{23,13}^0+G_{31,32}^0+G_{33,12}^0{G}_{21,23}^0+G_{22,13}^0+G_{31,22}^0+G_{32,12}^0\\ G_{11,31}^0+G_{31,11}^0{G}_{12,32}^0+G_{32,12}^0{G}_{13,33}^0+G_{33,13}^0{G}_{12,33}^0+G_{32,31}^0+G_{13,23}^0+G_{33,21}^0{G}_{11,33}^0+G_{33,11}^0+2G_{13,13}^0{G}_{12,13}^0+G_{11,23}^0+G_{32,11}^0+G_{31,21}^0\\ G_{11,12}^0+G_{21,11}^0{G}_{22,21}^0+G_{12,22}^0{G}_{23,31}^0+G_{13,32}^0{G}_{12,32}^0+G_{22,31}^0+G_{13,22}^0+G_{23,21}^0{G}_{21,31}^0+G_{11,32}^0+G_{23,11}^0+G_{13,12}^0{G}_{11,22}^0+G_{22,11}^0+2G_{12,12}^0\end{array}\right].$ (44)

One can easily verify that Γ is a symmetric matrix, i.e., Γ^T=Γ.

From Eq. (39) it is seen that the following integral is involved:

${\displaystyle {\iiint }_{V(x)}\Gamma (x^{″}-x)\Pi (x)\Psi (x-x^{\prime })d^{3}x}$,

where Γ −(xⁿ-x) 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 G_β"i"⁰, 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.