# Breakdown of equipartition in diffuse fields caused by energy leakage

L. Margerin

► To cite this version:

L. Margerin. Breakdown of equipartition in diffuse fields caused by energy leakage. The European Physical Journal. Special Topics, 2017, 226, 10.1140/epjst/e2016-60165-6. insu-03677009

HAL Id: insu-03677009

https://insu.hal.science/insu-03677009v1

Submitted on 17 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Regular Article

# Breakdown of equipartition in diffuse fields caused by energy leakage

L. Margerinª

Institut de Recherche en Astrophysique et Planétologie, CNRS – Université Paul Sabatier,
14 Avenue Edouard Belin, 31400 Toulouse, France

Received 19 May 2016 / Received in final form 21 September 2016

Published online 25 May 2017

**Abstract.** Equipartition is a central concept in the analysis of random wavefields which stipulates that in an infinite scattering medium all modes and propagation directions become equally probable at long lapse time in the coda. The objective of this work is to examine quantitatively how this conclusion is affected in an open waveguide geometry, with a particular emphasis on seismological applications. To carry out this task, the problem is recast as a spectral analysis of the radiative transfer equation. Using a discrete ordinate approach, the smallest eigenvalue and associated eigenfunction of the transfer equation, which control the asymptotic intensity distribution in the waveguide, are determined numerically with the aid of a shooting algorithm. The inverse of this eigenvalue may be interpreted as the leakage time of the diffuse waves out of the waveguide. The associated eigenfunction provides the depth and angular distribution of the specific intensity. The effect of boundary conditions and scattering anisotropy is investigated in a series of numerical experiments. Two propagation regimes are identified, depending on the ratio $H^*$ between the thickness of the waveguide and the transport mean path in the layer. The thick layer regime $H^* > 1$ has been thoroughly studied in the literature in the framework of diffusion theory and is briefly considered. In the thin layer regime $H^* < 1$, we find that both boundary conditions and scattering anisotropy leave a strong imprint on the leakage effect. A parametric study reveals that in the presence of a flat free surface, the leakage time is essentially controlled by the mean free time of the waves in the layer in the limit $H^* \to 0$. By contrast, when the free surface is rough, the travel time of ballistic waves propagating through the crust becomes the limiting factor. For fixed $H^*$, the efficacy of leakage, as quantified by the inverse coda quality factor, increases with scattering anisotropy. For sufficiently thin layers $H^* \approx 1/5$, the energy flux is predominantly directed parallel to the surface and equipartition breaks down. Qualitatively, the anisotropy of the intensity field is found to increase with the inverse non-dimensional leakage time, with the scattering mean free time as time scale. Because it enhances leakage, a rough free surface may result in stronger anisotropy of the intensity field than a flat surface, for the same bulk scattering properties. Our work identifies leakage as a potential explanation for the large deviation from isotropy observed in the coda of body waves.

<footnote>ª e-mail: ludovic.margerin@irap.omp.eu</footnote>
# 1 Introduction

Equipartition is a central concept in a number of recent advances in the field of multiple scattering of waves and ambient noise. Following Weaver (1982); Ryzhik et al. (1996) and many others, we say that a field is at equipartition when, in a given frequency band, all the modes of the system are statistically equiprobable. By mode we mean a stationary solution – an eigenfunction – of the underlying wave equation independent of the fact that the system is closed or open. In particular, in infinite translationally invariant media, the modes are plane waves that form a complete continuous basis. In this context, equipartition implies in particular that the angular distribution of energy flux tends to isotropy at long lapse-time. The main contribution of this paper is to show that this conclusion may completely break down when propagation of the waves takes place in an open waveguide. Such a configuration is widespread in seismology but its implications have been largely overlooked so far.

Previous studies have shown that equipartition may be achieved thanks to re-
verberations (Weaver, 1982), multiple-scattering (Ryzhik et al., 1996) or by con-
tinuous excitation of the waves by a random and statistically uniform distribution
of sources (Lobkis and Weaver, 2001). Equipartition implies that the average cross-
correlation of the wavefield is proportional to the imaginary part of the Green's func-
tion in the frequency domain (Barabanenkov and Ozrin, 1991; Van Tiggelen, 2003),
a property which has been widely exploited in recent years for both imaging and
monitoring purposes (Shapiro et al., 2005; Sens-Schönfelder and Wegler, 2006;
Brenguier et al., 2008). In practice, strict sense equipartition is rarely met for a
number of reasons. In ambient noise cross-correlation, the non-homogeneity of the
sources of noise is most often invoked to explain deviations from equipartition. Other
phenomena such as spatial variation of absorption or preferential absorption of some
modes may also hamper equipartition, even if the distribution of random sources or
scatterers is statistically homogeneous (Margerin et al., 2001; Snieder, 2007). In this
paper, I will focus on the equipartition properties of a diffuse field excited by a single
source with a particular emphasis on earthquake coda waves. The word coda refers
to the long time-tail of the seismograms following the ballistic arrivals and which is
composed of waves scattered by medium heterogeneities (Aki, 1969; Aki and Chouet,
1975). Because absorption hampers the detection of coda waves at large epicentral
distance for small to moderate earthquakes, equipartition is often thought to apply
locally, i.e., in a region of finite extent surrounding the source and station. For appli-
cations of the equipartition concept to the noise wavefield, I refer to the review papers
by Sánchez-Sesma et al. (2008), Weaver (2010) and Campillo and Roux (2015).

Sato and Matsumura (1980) were probably the first to address the problem of equipartition in seismology by studying the partitioning of the kinetic energy on the three components of a seismometer installed in a deep borehole in the Kanto plain. With the deployment of seismological networks, it is now possible to analyze the cross-correlation of coda waves across arrays, which provides stringent tests of equipartition. Van Tiggelen (2003) shows that imperfect equipartition manifests itself as a lack of temporal symmetry in the correlation function. The asymmetry is caused by the persistence of an energy flux coming from the source which vanishes only algebraically in time in the diffusion regime. This result clearly puts forward the asymptotic nature of equipartition when a single source is considered. Malcolm et al. (2004) demonstrate experimentally the convergence of the diffuse coda wavefield towards equipartition in a laboratory experiment of ultrasound propagation in a rock sample. They study the temporal asymmetry of the cross-correlation function and show that it depends strongly on the lapse-time in the coda, i.e., on the “age” of the signal used to perform the cross-correlation. The strong asymmetry observed in the “early” coda progressively vanishes as the waves undergo more and more
scattering and surface reflections in the sample, in agreement with multiple scat-
tering theory. A similar observation was made by Paul et al. (2005) using earth-
quake coda waves. These authors show that the level of scattering in the crust in
Alaska is too weak to reach equipartition before the signal enters the noise level.
Emoto et al. (2015) performed an extensive study of the asymmetry of the correlation
function of noise and earthquake coda waves using a large dataset from Japan. Their
analysis demonstrates very conclusively the relation between the temporal asymme-
try of the cross-correlation function and the energy flux radiated from the source.
They further show that an almost isotropic distribution of Rayleigh waves composes
the late coda in the 0.1–0.2 Hz band in Japan.

Some studies have also discussed – either directly or indirectly – the relevance of equipartition for the interpretation of the long-period seismic coda. Using (f, k) analysis, Maeda et al. (2006) have shown that in the 90–180 s period band, the composition of the coda wavefield evolves with lapse time. Array observations of the early coda suggest that it is mostly composed of a sum of scattered fundamental modes Rayleigh waves, whereas the late coda appears to be dominated by high-Q spherical higher-modes. Sens-Schönfelder et al. (2015) used data from USAarray to analyze the angular distribution of energy in the coda of the 2013 Okhotsk earthquake. They showed that multiple reverberations without significant effect of scattering persist at long lapse-time in the coda in certain period bands. Furthermore, they find that the typical time to randomize the propagation direction of deeply steeping body waves approximately equals 10 hours at 40 s periods, which is of the order of the typical duration of the coda before the signal reaches the noise level. This observation underlines the asymptotic character of equipartition.

The primary tool to model quantitatively the convergence of the diffuse wave-
field towards equipartition is transport theory. In particular, the radiative transfer
equation allows one to calculate the modal, spatio-temporal and angular distribution
of energy in the coda as a function of the scattering and absorption properties of
the medium. Margerin et al. (2001) have studied the impact of absorption on the
modal equipartition of multiply-scattered *P* and *S* waves (i.e., without consideration
of the angular or spatial distribution of the energy). They find that equipartition
is stable against perturbations of the anelastic quality factor. In other words, only
strong absorption of one mode with respect to the other can significantly shift the
shear to compressional energy ratio in the coda. Using Monte Carlo simulations of the
transport process, Paul et al. (2005) take an additional step and show that the sta-
bilization of energy ratios in the coda occurs long before the energy flux has become
almost isotropic, i.e., before “ground-truth” equipartition is reached. These numerical
results provide an explanation for the success of equipartition theory in predicting
energy ratios observed in short period seismic data although the condition of perfect
isotropic illumination is still far from being satisfied.

In many of the works discussed above, the key role played in the transport process
by reflection and transmission at the boundaries of the open crustal waveguide is often
disregarded. However, several numerical and theoretical studies have underlined the
impact of stratification and boundary conditions (B.C.) on the multiple scattering
of waves in the heterogeneous Earth (Hoshiba, 1997; Margerin et al., 1998; Bal and
Ryzhik, 2002; Wegler, 2004). In particular, the leakage of energy from the crust to the
mantle has been put forward as a potentially important contribution to the decay of
the coda observed around 1Hz (Korn, 1993; Margerin et al., 1999). Bal and Ryzhik
(2002) consider the impact of surface disorder on guided Love waves and concoct
models of random roughness that hamper leakage. In this work, I revisit the question
of energy leakage and investigate its impact on the asymptotic angular distribution of
intensity in the coda. The critical question to be addressed is the following: To what
extent can coda waves be considered as a superposition of random plane waves coming
from all possible space directions at long lapse time? To answer this question, I use radiative transfer theory for scalar waves to model the multiple scattering process in a waveguide geometry. The role of B.C. is carefully examined and shown to have a profound impact on the depth and angular distribution of the intensity detected in the coda.

It is worth noting that the radiative transfer model does not take into account coherent multiple-scattering effects (Akkermans and Maynard, 1985). This neglect should not be critical in the weak scattering regime considered in this work. Other limitations of the model should be mentioned. Because we use a scalar transfer equation, the coupling between shear and compressional waves is not taken into account in our approach. Weaver (1990), Turner and Weaver (1994, 1995) and Ryzhik et al. (1996) derived transport equations for elastic waves which may be employed in the future to model the energy exchange between different waves types taking polarization of shear waves into account. Because coda waves are thought to be dominated by shear waves, a scalar model may still provide insight into the leakage effect. Another limitation of the present work comes from the fact that boundary conditions are formulated at the level of energy flux conservation laws, thereby neglecting the generation of interface waves. In a seismological context, Rayleigh surface waves are important since they are efficiently excited by shallow sources. Expanding the wavefield into a sum of Love and Lamb waves, Trégourès and van Tiggelen (2002) have established a multi-mode quasi 2-D transport equation in plate geometry, including the Rayleigh surface wave. A notable advantage of this approach is the explicit treatment of boundary conditions at the level of the elastic wave equation. However, because the bounding surfaces are perfectly reflecting, the leakage effect is so far not taken into account in their model. The present work examines the case where the frequency is so high that the contribution of Rayleigh waves to the energy transport in the layer may be considered negligible. In the next section, I introduce the key equations and briefly describe their numerical solution.

# 2 Asymptotic analysis of the transport equation

## 2.1 Basic equations and boundary conditions

To model the transport of energy in a scattering and absorbing medium, I introduce the radiative transfer equation satisfied by the specific intensity $I(\mathbf{r}, \hat{\mathbf{n}}, t)$:

$$ (\partial_t + c\hat{\mathbf{n}} \cdot \nabla) I(\mathbf{r}, \hat{\mathbf{n}}, t) = -\left(\frac{1}{\tau} + \frac{1}{t^a}\right) I(\mathbf{r}, \hat{\mathbf{n}}, t) + \frac{1}{4\pi\tau} \int_{4\pi} p(\hat{\mathbf{n}}, \hat{\mathbf{n}}') I(\mathbf{r}, \hat{\mathbf{n}}', t) d\hat{\mathbf{n}}' + S(\mathbf{r}, \hat{\mathbf{n}}, t). \quad (1) $$

The specific intensity quantifies the flux of energy directed around the unit vector $\hat{\mathbf{n}}$ at point $\mathbf{r}$ and time $t$ in a scattering medium. The source of intensity is denoted by $S(\mathbf{r}, \hat{\mathbf{n}}, t)$. Because we focus on the long time behavior of the intensity, this term will be of no concern to us and will be simply ignored. The physical parameters which enter into equation (1) are: the wave velocity $c$, the scattering mean free time $\tau$ and the absorption time $t_a$. In the regime of weak fluctuations, perturbation theory provides an explicit relation between the phase function $p(\hat{\mathbf{n}}, \hat{\mathbf{n}}')$ – which describes the anisotropy of the scattering – and the power spectrum – the Fourier transform of the spatial correlation function – of the fluctuations in elastic parameters. Analyses of well-log as well as seismic data support the idea that correlation functions of the Von-Karman type represent adequately the randomness of the Earth’s crust (see Sato et al., 2012, for a review). For the numerical applications, I specialize to the
<figure><img src="image_2.png" /><figcaption>Fig. 1. The four scattering models considered in this paper, classified after the combination of B.C. imposed at the bottom and top of the layer of thickness H. (a) Flat free surface with total specular reflection at the top, transparent boundary at the bottom (no impedance contrast). (b) Flat free surface with total specular reflection at the top, partial reflection and transmission at the bottom. $v_c$ (resp. $Z_c$) and $v_m$ (resp. $Z_m$) denotes the wavespeed (resp. acoustic impedance) of the scattering layer and of the underlying half-space, respectively. (c) Rough free surface with total diffuse reflection defined according to Lambert's law, total transmission at the bottom. (d) Rough free surface with total diffuse reflection defined according to Lambert's law, partial reflection and transmission at the bottom.</figcaption></figure>

popular exponential distribution. In what follows, I further assume the medium to be translationally invariant in the horizontal directions $\hat{x}$ and $\hat{y}$. Equation (1) must be supplemented with appropriate B.C. at the surface of the Earth and at depth. We shall consider the four situations depicted in Figure 1. To express the B.C.s, it is convenient to decompose the intensity into up-going (+) and down-going (−) components such that $I^{\pm}(\mu) = I(\pm\mu)$, where $\mu \in [0, 1]$ denotes the cosine of the incidence angle measured with respect to the vertical direction $\hat{z}$. Note that, following the seismological practice, the $\hat{z}$ axis points downwards. As a consequence, positive $\mu$ corresponds to increasing depths in the model. Four different kinds of B.C.s are considered in this work (see Fig. 1). Total specular reflection (top surface on Fig. 1a,b):

$$I^{+}(\mu) = I^{-}(\mu). \qquad (2a)$$

Total diffuse reflection (Lambert reflector, top surface in Fig. 1c,d):

$$I^{+}(\mu) = 2 \int_{0}^{1} I^{-}(\mu) \mu d\mu. \qquad (2b)$$

Full transmission (bottom surface in Figure 1a,c):

$$I^{-}(\mu) = 0. \qquad (2c)$$
Partial specular reflection (bottom surface in Fig. 1b,d):

$$I^{-}(\mu) = R(\mu)I^{+}(\mu) \qquad (2d)$$

where $R(\mu)$ denotes the energy flux reflection coefficient. In the scalar approximation this coefficient depends on the acoustic impedance contrast $Z_m/Z_c$ between the heterogeneous layer and the underlying homogeneous half-space as follows:

$$R = \begin{cases} \left( \frac{1 - \frac{Z_m \mu'}{Z_c \mu}}{1 + \frac{Z_m \mu'}{Z_c \mu}} \right)^2 & \text{for } \mu > \mu_c, \\ 1 & \text{for } \mu \le \mu_c \end{cases} \qquad (2e)$$

where $\mu_c = \sqrt{1 - v_c^2/v_m^2}$ is the cosine of the critical incidence angle, $\mu' = \frac{v_m}{v_c} \sqrt{\mu^2 - \mu_c^2}$, and $v_c$ (resp. $v_m$) denotes the wavespeed in the scattering layer (resp. half-space). In this work, I have only considered constant $\tau$ and $t_a$, but depth-dependent properties may be easily incorporated in the model. The Lambert reflector model is very common in optics but has not been very much in use in seismology. It assumes that the reflected field is perfectly diffuse so that no coherent component survives. It is of course an idealization. In the case of light, ground glass and matte paper have been found to be good Lambertian reflectors (see e.g., Thomas and Stamnes, 2002, chapter 5). The theoretical study by Rogatkin (2004) indicates that a coarse perfectly reflecting surface with a correlation length four times larger than its RMS topography fluctuations is an excellent approximation of a Lambert reflector in the short wavelength limit. In this work, the extreme case of perfectly diffuse reflection is not meant to represent accurately the reflection of waves at the Earth's free surface. However, this extreme model does offer insight into the leakage process as shown below.

Obtaining the full solution of equation (1) is notoriously difficult and has been the topic of a very large number of papers. The goal of this study is much more modest: characterize the asymptotic ($t \to \infty$) behavior of the specific intensity only. As discussed in the introduction, the multiple-scattering process tends to homogenize the energy distribution in phase space. This implies, in particular, that the specific intensity must become independent of both the azimuthal angle of propagation and the position in the horizontal plane at sufficiently long lapse-time in the coda. However, because of the stratification of the scattering properties in the medium, we expect the gradient of intensity in the vertical direction to remain finite, even in the late coda. Clearly, this gradient is caused by the leakage of energy at the bottom of the scattering layer. Taking these properties into account and rescaling the space and time variables by the scattering mean free path and the scattering mean free time respectively, the equation of transfer simplifies to:

$$\frac{\partial I(z, \mu)}{\partial t} + \mu \frac{dI(z, \mu)}{dz} = -I(z, \mu) + \frac{1}{2} \int_{-1}^{1} p(\mu, \mu') I(z, \mu') d\mu', \qquad (3)$$

where $p(\mu, \mu')$ denotes the azimuthally averaged phase function, normalized such that $\int_{-1}^{1} p(\mu, \mu') d\mu' = 2$. Note that the result (3) may also be obtained by integrating equation (1) over the horizontal plane. Equation (3) is of the form: $\partial_t I = -LI$, where $L$ is a linear operator. When the medium is closed, 0 is an eigenvalue of $L$ corresponding to the conservation of energy in the medium. In the case of an open system, characterizing the long-term behavior of the intensity suggests to look for the eigenvalue of $L$ that lies closest to zero. This eigenvalue is interpreted as the inverse residence time of the waves in the layer, using the terminology of Margerin et al. (1999). Note that the operator $L$ is not self-adjoint so that there is in fact no guarantee that the
spectrum of L is real. Nevertheless, numerical simulations of the transfer equation as well as the study of L in the limit of small mean free times (the so-called diffusion approximation) support the idea that the intensity $I(z, \mu)$ asymptotically decays exponentially in a variety of situations which are relevant to our purposes (Margerin et al., 1998). Bowden and Williams (1964) have studied analytically the eigenvalue problem for the operator L in the case of a slab with fully transmitting boundaries on each side. Their results confirm the asymptotic exponential decay of the specific intensity in this geometry. We will therefore take for granted that equation (3) has a well-defined eigenfunction with real eigenvalue $\tau_l^{-1}$ ($\tau_l$ is the leakage time), which encapsulates the long-term behavior of the intensity in the coda. Since absorption is supposed to be uniform, it only shifts the spectrum and will therefore not be considered. Note again that depth-dependent absorption and scattering properties can easily be incorporated in the model. Although this is not necessary from a mathematical point of view, we will restrict our search of the nondimensional eigenvalue to the interval $(0, \tau/\tau^*)$, with $\tau^*$ the transport mean free time of the waves in the layer. From the physical point of view, this ensures that the wavefield can indeed be deemed “diffuse”. We remind the reader that $\tau^* = \tau/(1-g)$ with g the mean cosine of the scattering angle, which will be also referred to as the anisotropy factor.

Equation (3) represents a tremendous simplification compared to the original transport problem (1). As a consequence, the full spatial and dynamical behavior of the specific intensity is lost. In the seismologically relevant case of a localized point source, horizontal gradients will also contribute to the observed angular dependence of the specific intensity and represent an additional source of non-isotropy which is not taken into account in our approach. Diffusion theory (e.g. Malcolm et al., 2004) suggests that these horizontal gradients vanish only algebraically in time, so that deviations from isotropy may be much larger in the observations than in our model due to the effect of the source. Numerical simulations of the transport process for a point source with the aid of Monte Carlo simulations are therefore highly necessary to attest of the relevance of our study to seismic data. They will be presented in Section 3. In what follows, I outline the numerical technique employed to solve approximately the spectral problem.

## 2.2 Numerical solution by the discrete ordinate approach

Following Chandrasekhar (1960), I discretize the specific intensity into a set of directions based on a numerical quadrature scheme to be specified below. As a result, the integral term becomes a finite sum and equation (3) reduces to a system of coupled first-order linear differential equations for the unknown $I(z, \mu_i)$. To ensure that the B.C. are accurately specified numerically, I further decompose the intensity into down-going and up-going components to obtain:

$$ \mu_i \frac{dI^+(z, \mu_i)}{dz} = (\lambda - 1)I^+(z, \mu_i) + \frac{1}{2} \sum_{j=1}^{j=N} w_j p(\mu_i, \mu_j) I^+(z, \mu_i) \\ + \frac{1}{2} \sum_{j=1}^{j=N} w_{-j} p(\mu_i, -\mu_j) I^-(z, \mu_j), \quad (4) $$

$$ -\mu_i \frac{dI^-(z, \mu_i)}{dz} = (\lambda - 1)I^-(z, \mu_i) + \frac{1}{2} \sum_{j=1}^{j=N} w_j p(-\mu_i, +\mu_j) I^+(z, \mu_j) \\ + \frac{1}{2} \sum_{j=1}^{j=N} w_{-j} p(-\mu_i, -\mu_j) I^-(z, \mu_j) $$
where the notational convention of equation (2a)-(2d) has been adopted. In equation (4), $\lambda$ denotes the eigenvalue parameter and the weights of the quadrature scheme $w_j$ are even. In this work, I adopt a variant of the well-known “double Gauss” quadrature formula by splitting the integration domain into 4 intervals: $(-1, -\mu_c)$, $(-\mu_c, 0)$, $(0, \mu_c)$, $(\mu_c, 1)$. This procedure guarantees that the points of integration are clustered where the intensity varies most rapidly and also allows one to take into account the discontinuity of the first derivative of the reflection coefficient at $\mu = \mu_c$. The quadrature points and weights have been obtained from the usual Gauss-Legendre formula by mapping each sub-interval to $(-1, 1)$ through a change of variable. To achieve sufficient accuracy, I use 64 points of integration per interval. As mentioned in introduction, the phase function chosen in this work is based on the assumption that the fluctuations are exponentially correlated. Following the approach of Stamnes et al. (1988), I further approximate the phase function through a Legendre series expansion, which guarantees that no spurious dissipation is introduced by the numerical scheme (Wiscombe, 1977). Analytical expressions for the coefficients of the Legendre series are taken from Le Bihan and Margerin (2009).

As is well-known, the solutions to the homogeneous first-order differential system (4) span a $2N$-dimensional space and are of the form:

$$ \psi(z, \mu_i; \lambda) = \sum_{i=1}^{2N} C_j V_j(\mu_i; \lambda) e^{\alpha_i(\lambda)z} \quad (5) $$

where $C_i$ are constants to be determined. The first (resp. last) $N$ components of $\psi$ correspond to the down-going (resp. up-going) fluxes. The eigenvalues $\alpha_i$ and $2N$-dimensional eigenvectors $V_i$ may be obtained upon substitution of the ansatz $\psi(z) = Ve^{\alpha z}$ into equation (4) using standard techniques from linear algebra. As is clear from equation (5) the eigenvalues and eigenvectors depend on the parameter $\lambda$ and are generally complex. Note that complex eigenvalues $\alpha_i$ always occur as complex conjugated pairs. Additionally, it may be noted that if $\alpha_i$ is an eigenvalue then so is $-\alpha_i$. To determine the smallest eigenvalue of the operator $L$, I employ a simple shooting method. Starting with an initial guess $0 < \lambda_0 < \tau/\tau^*$, one first calculates the set of eigenvalues $\alpha_i(\lambda_0)$ and eigenvectors $V_i(\lambda_0)$. The general linear combination (5) satisfies the system of differential equations (4) but not the B.C. at the surface and at depth. Imposing these B.C. results in a homogeneous linear system of $2N$ equations in the $2N$ unknown $C_i$ which only has the trivial solution unless its determinant $\Delta(\lambda_0)$ vanishes. Using a standard optimization algorithm such as Newton's method, it is possible to find a new candidate eigenvalue $\lambda_1$ that improves on the initial guess. By monitoring the real and imaginary parts of $\Delta$, and iterating the process with $\lambda_2, \dots$ the sought after nondimensional eigenvalue $\tau_l^{-1}$ may be determined to the desired accuracy. The corresponding eigenfunction $\psi(z, \mu; \tau_l^{-1})$ is then easily found by normalizing one of the coefficients, say $C_{2N} = 1$ and solving for the remaining coefficients to match the B.C.s.

For the above procedure to be applicable, the knowledge of an initial guess $\lambda_0$ is necessary. In general, this is a very difficult problem. In this study, however, we are primarily interested in finding the smallest eigenvalue of $L$ for a series of layer thickness $H$. Therefore, we may start with a sufficiently thick layer (say $H/l \gg \tau^*/\tau$) for which the diffusion approximation is applicable. The associated spectral problem is easy to solve numerically and provides an excellent approximation to the eigenvalue of the exact transport equation (see e.g. Margerin et al., 1998, for extensive numerical tests). We may then decrease iteratively the layer thickness by a small increment and make use of the exact eigenvalue obtained at the preceding step as the initial guess for the current step. In this way, we obtain a sequence of eigenvalues of the spectral problem corresponding to decreasing layer thicknesses, as required.
<figure><img src="image_3.png" /><figcaption>Fig. 2. Asymptotic energy density distribution in a diffuse layer with isotropic scattering predicted by Monte Carlo simulations (symbols) and by numerical spectral analysis (solid lines). In both cases, the results have been normalized by their RMS. On the horizontal axis, the unit length is the mean free path. Triangles: case (a). Circles: case (b). Diamonds: case (c).</figcaption></figure>

# 3 Impact of leakage on coda decay and angular pattern of intensity at long lapse-time

## 3.1 Numerical tests of the method

I begin with a numerical verification of the spectral technique described in the previous section using comparisons with Monte Carlo simulations. In Figure 2, I consider the asymptotic behavior of the intensity in a layer of thickness one men free path. In the Monte Carlo approach, the random walk of a large number of particles launched from an isotropic point source positioned at a depth $z = H/2$ is monitored as a function of time. For a flat surface, the B.C.s are imposed through the application of the usual laws of geometrical optics to the particle. The possibility of reflection and transmission is treated as a Bernoulli process. In the case of a Lambert surface, the probability distribution of the particle direction after reflection is given by $\mu = \sqrt{\epsilon}$, with $\epsilon$ a uniformly distributed random number in (0,1). This ensures that the flux reflected in a particular direction by the surface is proportional to the direction cosine $\mu$, which in turn corresponds to a uniform distribution for the reflected specific intensity. After each time step $dt/\tau = 0.025$, the particle density is averaged horizontally in thin layers of thickness $dz/l = 1/80$. The energy density profiles are then averaged temporally from lapse time $t/\tau = 7.5$ to $t/\tau = 15$. Note that the results are independent of the position or radiation pattern of the source, as they should. For the three combinations of B.C.s considered in Figure 1a–c, excellent agreement is found between the discrete-ordinate and Monte Carlo methods. Note that the results have been normalized by their RMS to facilitate the comparison. The nondimensional eigenvalues corresponding to the eigenfunctions shown in Figure 2 are given by $\tau_l^{-1} = 0.297$ (a), $\tau_l^{-1} = 0.135$ (b) and $\tau_l^{-1} = 0.311$ (c) (in units of inverse mean free time). The Monte Carlo and discrete ordinate methods agree to within three digits accuracy. Qualitatively, the gradient of intensity is seen to increase with the rate of leakage. In the case of a Lambert surface the intensity has a maximum inside the layer. Such a phenomenon is not predicted by the diffusion approximation, which only considers flux conservation at interfaces and therefore does not distinguish between a flat and a rough boundary. In the framework of the diffusion approximation, the fact that the intensity has no maximum inside the layer may also be understood as a consequence of the maximum principle (see e.g. Stakgold and Holst, 2011, p. 507).
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Polar diagram of the asymptotic distribution of specific intensity at the surface of a scattering layer with thickness $H/l = 1$ (top) or $H/l = 1/4$ (bottom) for cases (a) (left) and (c) (right) as summarized in Figure 1. Note that for the thinnest layer a logarithmic scale is employed on the radial axis. The ratio between the minimum ($\mu \to 1$) and maximum ($\mu \to 0$) value of the specific intensity is larger than 20 and 100 in case (a) and (c), respectively.</figcaption></figure>

## 3.2 Complete breakdown of equipartition: an example

As a clear illustration of the result put forward in the title of this article, I consider in Figure 3 the case of a layer of isotropic scatterers with either a flat (left) or rough (right) totally reflecting surface at the top and transparent B.C.s at the bottom. I first examine the top plots corresponding to the adimensional layer thickness $H/l = 1$ for which we have already calculated the intensity profile (see Fig. 2). In the case of a flat surface (left), the angular intensity distribution is symmetric with respect to an inversion of the propagation direction as required by equation (2a). In the case of a rough surface, the dowgoing intensity verifies the isotropy conditions in accordance with equation (2b). It is therefore clear that the eigenfunctions shown in Figure 3 verify the B.C.s at the surface. Although not shown here, this property also holds at the lower boundary, which gives us confidence in the accuracy of the numerical results. The calculations demonstrate that, even asymptotically, the angular distribution of energy flux departs significantly from isotropy. The specific intensity is typically larger along the surface than perpendicular to it. This feature will crop up repeatedly in other examples and contrasts remarkably with the prediction of diffusion theory, that the maximum intensity should align with the current vector in the positive z direction.

In the case of a thin layer ($H/l = 1/4$, bottom), the lowest adimensional eigenvalue of the transport problem is $\tau_l^{-1} = 0.905, 0.976$, for case (a) and (c) respectively. In other words the leakage time is only slightly larger than the transport mean free
time in the layer. For such fast leakage rates, the basic Monte Carlo method described above becomes inapplicable. This is admittedly an extreme scenario but it does provide insight into the impact of B.C.s on the angular distribution of intensity in the coda at long lapse-time. In particular, one observes a very large peak of intensity directed along the surface. Typically, the ratio between the specific intensity along and perpendicular to the surface can be as large as 100. Said in loose terms, the waves that remain in the layer at long lapse-time tend to minimize their interaction with the boundaries of the medium, a conclusion which appears intuitively quite reasonable. The calculations shown in Figure 3 also illustrate qualitatively the fact that the anisotropy of the specific intensity increases with the eigenvalue of the non-dimensionalized transport problem. This correlation between the leakage rate and the anisotropy of the specific intensity will be further illustrated below.

### 3.3 Analysis of leakage in open waveguide geometry

As the main application of the theory developed in this work, I revisit the problem of energy leakage from the Earth's crust and its interplay with the persistence of an anisotropic flux of energy at long lapse-time in the coda. This corresponds to models (b) and (d) in Figure 1 with an impedance contrast $Z_m/Z_c \approx 1.34$ between crust and mantle. Leakage is one of the two key processes at the origin of the decay of crustal coda waves, the second being absorption. The coda decay is most often quantified with the aid of a quality factor $Q_c$ (Aki and Chouet, 1975). $Q_c$ gained high popularity in the seismological community because this parameter is easy to measure and conveys information on attenuation which is complementary to the velocity structure. At long lapse-time in the coda, the coda quality factor may be expressed as (Margerin et al., 1998):

$$Q_c^{-1} = Q_i^{-1} + Q_l^{-1} \qquad (6)$$

where $Q_i$ and $Q_l$ denote the absorption and leakage quality factors, respectively. Equation (6) shows that it is essential to remove the leakage contribution from the observed coda quality factor to get a correct estimate of the absorption properties of rocks at depth. It is worth noting that the guiding of the waves between the Moho and the free surface of the Earth is usually neglected in popular analysis methods such as the Multiple Lapse-Time Window Analysis (Fehler et al., 1992). Based on Monte Carlo simulations, Del Pezzo and Bianco (2010) have proposed an empirical a posteriori correction of the outputs of MLTWA for the leakage effect. Using spectral analysis of the transport equation, it is in fact possible to calculate the leakage quality factor of the waves as a function of the transport mean free path in the crust. This offers an alternative approach to the a posteriori correction of the estimated absorption for the leakage effect.

In seismological applications, the thickness of the crust is usually known with much less uncertainty than its scattering properties. To compare models with different levels of scattering anisotropy, it is therefore natural to express the leakage quality factor as a function of the transport mean free path:

$$Q_l = \omega \tau \tau_l(l^*) = \frac{\omega l^* (1 - g) \tau_l(l^*)}{c} \qquad (7)$$

where $\omega$ is the central frequency of the waves, and $c$ is the wave speed. In the examples that follow, I assume a nominal frequency of 1 Hz ($\omega = 2\pi$) and a wave speed of 3.5 km/s, typical of shear waves in the Earth's crust. The results of the calculation of the leakage time and quality factor in a series of anisotropically scattering media are presented, respectively, in Figure 4 and 5 for either a flat (left) or rough (right)
<figure><img src="image_5.png" /><figcaption>**Fig. 4.** Non-dimensional leakage time as a function of non-dimensional thickness in a heterogeneous waveguide. The correlation function of the fluctuations is exponential and the mean cosine of the scattering angle *g* quantifies the scattering anisotropy. The different curves correspond to different levels of anisotropy (see inset). The transport mean free path *l*\* is used as scaling length on the horizontal axis. The impedance contrast between mantle and crust is *Z*<sub>m</sub>/*Z*<sub>c</sub> ≈ 1.34. Left: flat free surface (case b). Right: rough free surface (case d). See Figure 1 for further details.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Leakage quality factor at 1 Hz in a heterogeneous waveguide as a function of the transport mean free path. In this example, a typical crustal thickness *H* = 35 km is assumed and the correlation function of the fluctuations is exponential. The mean cosine of the scattering angle *g* quantifies the scattering anisotropy and the different curves correspond to different levels of anisotropy (see inset). The impedance contrast between mantle and crust is *Z*<sub>m</sub>/*Z*<sub>c</sub> ≈ 1.34. Left: flat free surface (case b). Right: rough free surface (case d). See Figure 1 for further details.</figcaption></figure>

free surface. Although these figures basically convey the same information, they shed light on different aspects of leakage. Figure 4, which represents $\tau_l$ as a function of $H^*$, will be particularly useful in the following to establish a link between the non-dimensional leakage rate and the anisotropy of the energy flux in the coda at long lapse time. Figure 5, which represents $Q_l$ as a function of $l^*$ (for $H = 35$ km), is particularly relevant to seismological applications, where the impact of leakage on the coda decay needs to be assessed. Before discussing in greater details Figures 4 and 5, it is worth noting that models with the same $\tau_l$ and $l^*$ do not necessarily have the same leakage quality factor. This is simply explained by the fact that the scattering mean free time is the key conversion factor from $\tau_l$ to $Q_c$. Media with the same $l^*$ but different anisotropy factors $g$ have different mean free time $\tau$, and therefore different $Q_c$ according to formula (7).

Figure 4 shows that the non-dimensional leakage time $\tau_l$ depends monotonically on the non-dimensional crustal thickness $H^*$ and the anisotropy factor $g$. Based on
the scaling of $\tau_l$ with $H^*$, Figure 4 allows one to define a thick layer ($H^* > 1$) and a thin layer ($H^* < 1$) propagation regime in the waveguide. These two regimes may also be clearly identified in Figure 5. In the thick layer case, the leakage quality factor is a decreasing function of $l^*$ and is independent of the details of the propagation. By contrast, in the thin layer case, the leakage quality factor increases with $l^*$, and both B.C.s and scattering anisotropy play a prominent role. The difference in behavior between $\tau_l$ and $Q_l$ may be traced back to the extra $l^*$ factor which appears in equation (7). The case $g=0$ of Figure 5 (left) was previously examined by Margerin et al. (1998) using the diffusion approximation for $H^* > 1$ and numerical solutions of the transfer equation by the Monte Carlo method for $H^* < 1$. Excellent agreement has been found between their results and the spectral analysis presented in this paper. For the study of equipartition, the spectral technique does offer a number of advantages. It is fast and accurate, even more so when scattering is anisotropic, and allows the computation of both the leakage rate and angular pattern of specific intensity. In what follows, the dependence of the leakage time and quality factor on the scattering properties of the crust will be analyzed separately for the two regimes.

The thick layer regime ($H^* > 1$) has been well studied in the literature and will be described briefly. In this regime, the leakage quality factor is independent of the scattering anisotropy and B.C. at the free surface, and the diffusion approximation does an excellent job at approximating the eigenvalues of the transport problem. The leakage time turns out to be essentially equal to the Thouless time $H^2/D$ which quantifies the propagation time of diffuse waves through the crust. This implies the scaling relation $\tau_l \propto (H^*)^2$ which is indeed approximately verified in Figure 4, independent of the B.C. at the surface. It is worth noting that if one uses the transport mean free time as the time scale in Figure 4, the curves corresponding to different anisotropy factor $g$ collapse on a common master curve for $H^* > 1$ (but not for $H^* < 1$). This confirms that the transport mean free time is the relevant time-scale in the thick layer regime and explains the independence of the leakage quality factor on the scattering anisotropy in Figure 5. The scaling relation for $\tau_l$ indicates that $Q_l$ varies approximately like $1/l^*$ in the thick layer regime.

The thin layer case $H^* < 1$ requires more detailed discussion. In this regime, multiple reflections at the free surface and Moho generate partially guided waves whose properties depend on the angular dependence of the reflection coefficient. These complexities cannot be correctly modeled in the diffusion approximation and call for the use of the radiative transfer machinery. Figure 5 clearly demonstrates that the efficacy of leakage depends to a large extent on the type of bounding free surface (flat vs rough) as well as on the details of the scattering pattern. Further examination of Figure 5 reveals that surface roughness enhances leakage and strongly reduces its sensitivity to the anisotropy of the scattering pattern. The dependence of leakage on B.C.s imposed at the surface may be understood as follows. In the flat surface case and in the absence of scattering, the waves that propagate post-critically remain indefinitely trapped in the crust and dominate the signal at long lapse-time. In the presence of scattering, a randomization of the propagation direction occurs with a time scale typically equal to the transport mean free time. Because the layer is thin, other details of the escape process such as the time taken from the last scattering event to the Moho are negligible. This reasoning indicates that the limiting time of the leakage process is the mean free time and suggests the scaling $\tau_l \propto 1$ (i.e., $\tau_l$ is independent of $H^*$). This, in turn, explains why the curves tend to flatten out at small $H^*$ on Figure 4 (Left). By contrast, in the case of a Lambert reflector at the surface, the wave trajectories get completely randomized after each interaction with the roughness. In other words, even in the absence of scattering in the bulk of the medium, leakage is guaranteed by the presence of the rough surface. In this case, the travel-time of the waves through the crust (itself proportional to $H$) now
<figure><img src="image_7.png" /><figcaption>Fig. 6. Polar diagram of the asymptotic distribution of specific intensity at the surface of a crustal waveguide with thickness $H/l^* = 1$ (top) or $H/l^* = 1/5$ (bottom) for cases (b) (left) and (d) (right) as summarized in Figure 1. The impedance contrast is $Z_m/Z_c \approx 1.34$ and isotropic scattering is assumed.</figcaption></figure>

plays the role of limiting factor which suggests the scaling $\tau_l \propto H^*$. This last relation is approximately verified in Figure 4 (right) and explains the flattening of the $Q_l$ curves for $l^* > H$ in Figure 5. Actually, the fact that $Q_l$ slightly increases as $H^*$ decreases in the thin layer regime indicates that $\tau_l$ increases slightly less rapidly than $H$. In all cases examined, the efficacy of energy leakage, as quantified by the quality factor $Q_l$, is maximum around $H \approx l^*$, i.e., for a transport mean free path in the range 30–40 km. This is approximately the value of the mean free path estimated in the most tectonically active regions in the crust. The results shown in Figure 5 are therefore important to interpret the coda quality factor $Q_c$ measured in these regions. Finally, comparison of Figures 4 and 5 reveals that for fixed $l^*$, $\tau_l$ and $Q_c$ vary oppositely as a function of the anisotropy factor $g$. This difference in behavior is simply explained by the extra $(1-g)$ factor in formula (7). Figure 5 demonstrates that in the thin layer regime the efficacy of leakage is enhanced by scattering anisotropy, i.e., $Q_l^{-1}$ increases with $g$.

## 3.4 Impact of leakage on equipartition in the Earth's crust

In Figure 6 the impact of energy leakage on the angular distribution of intensity in the seismic coda is examined in the case of isotropic scattering for two different layer thicknesses and two kinds of B.C.s. For a moderately thin crust ($H^* = 1$), the specific intensity is only weakly anisotropic. The anisotropy of the energy flow, as quantified
<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Figure 6 for anisotropic scattering. The correlation of the fluctuations is exponential and the anisotropy factor is $g = 0.8$.</figcaption></figure>

by the ratio between the maximum and minimum value of the specific intensity, is of the order of 3/2 only. Like in Figure 3, we remark that the intensity is larger parallel to the surface than perpendicular to it. In the case of a rough free surface, the direction of maximum intensity is not exactly aligned with the horizontal direction but makes a relatively small angle to it. This feature is also noticeable in Figure 3. For a relatively thin crust ($H^* = 1/5$), the imprint of the reflection coefficient at the Moho is manifest in the lack of energy coming from a cone of directions with an aperture of about 45 degrees. The direction of maximum intensity tends to align along the horizontal axis with again a slight deviation in the case of a Lambert reflector at the surface. The anisotropy of the flow is typically larger than 5 and is more pronounced in the case of a rough free surface. This paradoxical result is consistent with the fact that leakage is more efficient in this case (compare the right and left plots in Figure 5). More precisely, our parametric study reveals that the anisotropy of the intensity field qualitatively increases with the inverse of the non-dimensional leakage time $\tau_l$.

In Figure 7 the role of scattering anisotropy is illustrated in the case of a medium with an exponential correlation function and an anisotropy factor $g = 0.8$. For comparison with Figure 6, the choice has been made to use the same non-dimensionalized layer thickness. Qualitatively, Figure 7 presents the same characteristics as Figure 6. The main difference lies in the fact that the pattern of intensity is generally less anisotropic than in the case of an isotropically scattering medium. In the case of a relatively thin layer ($H^* = 1/5$), the anisotropy of the specific intensity can still be rather large, approximately equal to 5. The comparison of Figures 6 and 7 reveals that the transport mean free path alone does not fully determine the intensity pattern
at long lapse time in the coda in a waveguide geometry. Perhaps surprisingly, for the same value of the transport mean free path, perfect isotropic illumination is more severely broken in the case of isotropic scattering in the crustal waveguide geometry. This result is again consistent with the fact that the anisotropy of the field is controlled by the non-dimensional leakage time $\tau_l$. Let me emphasize that the knowledge of $Q_l$ alone does not suffice to characterize quantitatively the anisotropy of the intensity in the coda. This result is rather obvious from Figure 5, which shows that models sharing the same $Q_l$ correspond to the drastically different propagation regimes $H^* < 1$ and $H^* > 1$. The later case corresponds to a very weakly anisotropic intensity field, well described by the diffusion approximation, while the former corresponds to a regime of strongly anisotropic intensity.

## 4 Conclusion

Many of the recent developments in the field of seismic scattering and ambient noise rely on the concept of equipartition, which is often formulated as the perfect isotropic illumination of the target medium. Although it is well understood that this condition is rarely met in practice, the mechanisms responsible for the observed lack of equipartition have not been fully elucidated so far. In this study, I have used a radiative transfer model to quantify the impact of B.C.s on the angular distribution of intensity in the coda. Two propagation regimes have been identified, in which the leakage time and quality factor show different dependence on the crustal thickness, scattering anisotropy and B.C.s at the surface. In the usual propagation regime where the transport mean free path is larger than the crustal thickness, the main conclusion of this work is that the angular intensity pattern can depart very significantly from isotropy, even at long lapse time in the coda. Furthermore, the anisotropy of the energy flux is found to increase with the inverse non-dimensionalized leakage time (i.e., the leakage rate). Because it enhances leakage, surface roughness may result in a more anisotropic intensity pattern than a flat surface, a result which may seem counter-intuitive at first. The potential lack of equipartition in body waves implied by this study may partially explain why these waves have proven more difficult to reconstruct than surface waves in applications of ambient noise cross-correlations at local to regional distance.

I would like to thank the referees for their constructive criticism which helped to improve the presentation of the results. This work was initiated almost 20 years ago during my Ph.D. work at the University of Grenoble, and benefited from an exceptional, truly interdisciplinary environment promoted, among others, by Roger Maynard. Roger is also gratefully acknowledged for introducing me to multiple scattering theory during a summer school in Cargèse in 1995.

## References

K. Aki, Analysis of the seismic coda of local earthquakes as scattered waves, J. Geophys. Res. **74**, 615 (1969)

K. Aki, B. Chouet, Origin of coda waves: source, attenuation and scattering effects, J. Geophys. Res. **80**, 3322 (1975)

E. Akkermans, R. Maynard, Weak localization of waves, J. de Phys. Lett. **46**, 1045 (1985)

G. Bal, L. Ryzhik, Wave transport for a scalar model of the Love waves, Wave Motion **36**, 49 (2002)
Y.N. Barabanenkov, V. Ozrin, Asymptotic solution of the Bethe-Salpeter equation and the Green-Kubo formula for the diffusion constant for wave propagation in random media, Phys. Lett. A. **154**, 38 (1991)

R.L. Bowden, C.D. Williams, Solution of the initial-value transport problem for monoenergetic neutrons in slab geometry, J. Math. Phys. **5**, 1527 (1964)

F. Brenguier, N.M. Shapiro, M. Campillo, V. Ferrazzini, Z. Duputel, O. Coutant, A. Nercessian, Towards forecasting volcanic eruptions using seismic noise, Nat. Geosci. **1**, 126 (2008)

M. Campillo, P. Roux, Seismic imaging and monitoring with ambient noise correlations. Edited by B Romanowicz, A. Dziewonski, *Treatise on Geophysics* (Elsevier, 2015) Vol 1, pp. 256-271

S. Chandrasekhar, *Radiative Transfer* (Dover publications, 1960)

E. Del Pezzo, F. Bianco, Two-layer earth model corrections to the MLTWA estimates of intrinsic and scattering attenuation obtained in a uniform half-space, Geophys. J. Int. **182**, 949 (2010)

K. Emoto, M. Campillo, F. Brenguier, X. Briand, T. Takeda, Asymmetry of coda cross-correlation function: dependence of the epicentre location, Geophys. J. Int. **201**, 1313 (2015)

M. Fehler, M. Hoshiba, H. Sato, K. Obara, Separation of scattering and intrinsic attenuation for the Kanto-Tokai region, Japan, using measurements of S-wave energy versus hypocentral distance, Geophys. J. Int. **108**, 787 (1992)

M. Hoshiba, Seismic coda wave envelope in depth-dependent s wave velocity structure, Phys. Earth Planet. Interiors **104**, 15 (1997)

M. Korn, Determination of site-dependent scattering Q from P-wave coda analysis with an energy-flux model, Geophys. J. Int. **113**, 54 (1993)

N. Le Bihan, L. Margerin, Nonparametric estimation of the heterogeneity of a random medium using compound poisson process modeling of wave multiple scattering, Phys. Rev. E **80**, 016601 (2009)

O.I. Lobkis, R.L. Weaver, On the emergence of the Green's function in the correlations of a diffuse field, J. Acoustical Soc. Am. **110**, 3011 (2001)

T. Maeda, H. Sato, M. Ohtake, Constituents of vertical-component coda waves at long periods, Pure Appl. Geophys. **163**, 549 (2006)

A. Malcolm, J. Scales, B. van Tiggelen, Extracting the green function from diffuse, equipartitioned waves, Phys. Rev. E **70**, 015601 (2004)

L. Margerin, M. Campillo, N. Shapiro, B. Van Tiggelen, Residence time of diffuse waves in the crust as a physical interpretation of coda Q: application to seismograms recorded in Mexico, Geophys. J. Int. **138**, 343 (1999)

L. Margerin, M. Campillo, B. Tiggelen, Radiative transfer and diffusion of waves in a layered medium: new insight into coda, Q. Geophys. J. Int. **134**, 596 (1998)

L. Margerin, B. van Tiggelen, M. Campillo, Effect of absorption on energy partition of elastic waves in the seismic coda, Bull. Seismological Soc. Am. **91**, 624 (2001)

A. Paul, M. Campillo, L. Margerin, E. Larose, A. Derode, Empirical synthesis of time-asymmetrical Green functions from the correlation of coda waves, J. Geophys. Res. **110**, B08302 (2005)

D. Rogatkin, Scattering of electromagnetic waves by a randomly rough surface as a boundary problem of laser radiation interaction with light-scattering materials and media, Opt. Spectrosc. **97**, 455 (2004)

L. Ryzhik, G. Papanicolaou, J.B. Keller, Transport equations for elastic and other waves in random media, Wave Motion **24**, 327 (1996)

F.J. Sánchez-Sesma, J.A. Pérez-Ruiz, F. Luzon, M. Campillo, A. Rodríguez-Castellanos, Diffuse fields in dynamic elasticity, Wave Motion **45**, 641 (2008)

H. Sato, M.C. Fehler, T. Maeda, *Seismic Wave Propagation and Scattering in the Heterogeneous Earth* (Springer, 2012)

H. Sato, S. Matsumura, Three-dimensional analysis of scattered P waves on the basis of the PP single isotropic scattering model, J. Phys. Earth **28**, 521 (1980)
C. Sens-Schönfelder, R. Snieder, S.C. Stähler, The lack of equipartitioning in global body wave coda, Geophys. Res. Lett. **42**, 7483 (2015)

C. Sens-Schönfelder, U. Wegler, Passive image interferometry and seasonal variations of seismic velocities at Merapi Volcano, Indonesia, Geophys. Res. Lett. **33**, 21302 (2006)

N.M. Shapiro, M. Campillo, L. Stehly, M.H. Ritzwoller, High-resolution surface-wave tomography from ambient seismic noise, Science **307**, 1615 (2005)

R. Snieder, Extracting the Green's function of attenuating heterogeneous acoustic media from uncorrelated waves, J. Acoustical Soc. Am. **121**, 2637 (2007)

I. Stakgold, M.J. Holst, Green's functions and boundary value problems (John Wiley & Sons, 2011)

K. Stamnes, S.-C. Tsay, W. Wiscombe, K. Jayaweera, Numerically stable algorithm for discrete-ordinate-method radiative transfer in multiple scattering and emitting layered media, Appl. Opt. **27**, 2502 (1988)

G.E. Thomas, K. Stamnes, Radiative transfer in the atmosphere and ocean (Cambridge University Press, 2002)

N.P. Trégourès, B.A. van Tiggelen, Quasi-two-dimensional transfer of elastic waves, Phys. Rev. E **66**, 036601 (2002)

J.A. Turner, R.L. Weaver, Radiative transfer of ultrasound, J. Acoustical Soc. Am. **96**, 3654 (1994)

J.A. Turner, R.L. Weaver, Ultrasonic radiative transfer in polycrystalline media: Effects of a fluid-solid interface, J. Acoustical Soc. Am. **98**, 2801 (1995)

B. van Tiggelen, Green function retrieval and time reversal in a disordered world, Phys. Rev. Lett. **91**, 243904 (2003)

R.L. Weaver, On diffuse waves in solid media, J. Acoustical Soc. Am. **71**, 1608 (1982)

R.L. Weaver, Diffusivity of ultrasound in polycrystals, J. Mech. Phys. Solids **38**, 55 (1990)

R.L. Weaver, Equipartition and retrieval of green's function, Earthquake Sci. **23**, 397 (2010)

U. Wegler, Diffusion of seismic waves in a thick layer: Theory and application to Vesuvius volcano, J. Geophys. Res. **109**, 7303 (2004)

W. Wiscombe, The delta-M method: rapid yet accurate radiative flux calculations for strongly asymmetric phase functions, J. Atmospheric Sci. **34**, 1408 (1977)

**Open Access** This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.