Breakdown of equipartition in diffuse fields caused by energy leakage

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 our 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∗ → 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∗≈ 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.

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 our 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 * → 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 * ≈ 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.

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 reverberations (Weaver, 1982), multiple-scattering (Ryzhik et al., 1996) or by continuous excitation of the waves by a random and statistically uniform distribution of sources (Lobkis and Weaver, 2001). Equipartition implies that the average crosscorrelation of the wavefield is proportional to the imaginary part of the Green's function 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 applications 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 scattering theory. A similar observation was made by Paul et al. ( 2005) using earthquake 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 asymmetry 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 spheroidal higher-modes. Sens-Schönfelder et al. (2015) used data from USArray 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 wavefield 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 stabilization 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), 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

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(r,n, t): (∂ t +cn·∇) I(r,n, t) =− 1 τ + 1 t a I(r,n, t)+ 1 4πτ 4π p(n,n )I(r,n , t)dn +S(r,n, t). (1) The specific intensity quantifies the flux of energy directed around the unit vectorn at point r and time t in a scattering medium. The source of intensity is denoted by S(r,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 τ and the absorption time t a . In the regime of weak fluctuations, perturbation theory provides an explicit relation between the phase function p(n,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 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. Zc) and vm (resp. Zm) 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. popular exponential distribution. In what follows, I further assume the medium to be translationally invariant in the horizontal directionsx andŷ. 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 ± (μ) = I(±μ), where μ ∈ [0, 1] denotes the cosine of the incidence angle measured with respect to the vertical directionẑ. Note that, following the seismological practice, theẑ axis points downwards. As a consequence, positive μ 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): Total diffuse reflection (Lambert reflector, top surface in Fig. 1c,d): Full transmission (bottom surface in Figure 1a,c): The European Physical Journal Special Topics Partial specular reflection (bottom surface in Fig. 1b,d): where R(μ) 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: and v c (resp. v m ) denotes the wavespeed in the scattering layer (resp. half-space). In this work, I have only considered constant τ and t a , but depthdependent 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 → ∞) 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: where p(μ, μ ) denotes the azimuthally averaged phase function, normalized such that 1 −1 p(μ, μ )dμ = 2. Note that the result (3) may also be obtained by integrating equation (1) over the horizontal plane. Equation (3) is of the form: ∂ 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, μ) 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 τ −1 l (τ 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, τ/τ * ), with τ * 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 τ * = τ /(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.

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, μ 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: The European Physical Journal Special Topics where the notational convention of equation (2a)-(2d) has been adopted. In equation (4), λ 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). 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 μ = μ 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: where C i are constants to be determined. The first (resp. last) N components of ψ correspond to the down-going (resp. up-going) fluxes. The eigenvalues α i and 2Ndimensional eigenvectors V i may be obtained upon substitution of the ansatz ψ(z) = V e αz into equation (4) using standard techniques from linear algebra. As is clear from equation (5) the eigenvalues and eigenvectors depend on the parameter λ and are generally complex. Note that complex eigenvalues α i always occur as complex conjugated pairs. Additionally, it may be noted that if α i is an eigenvalue then so is −α i . To determine the smallest eigenvalue of the operator L, I employ a simple shooting method. Starting with an initial guess 0 < λ 0 < τ/τ * , one first calculates the set of eigenvalues α i (λ 0 ) and eigenvectors V i (λ 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 Δ(λ 0 ) vanishes. Using a standard optimization algorithm such as Newton's method, it is possible to find a new candidate eigenvalue λ 1 that improves on the initial guess. By monitoring the real and imaginary parts of Δ, and iterating the process with λ 2,··· the sought after nondimensional eigenvalue τ −1 l may be determined to the desired accuracy. The corresponding eigenfunction ψ(z, μ; τ −1 l ) 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 λ 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 τ * /τ ) 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. 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 Bernouilli process. In the case of a Lambert surface, the probability distribution of the particle direction after reflection is given by μ = √ , with 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 μ, which in turn corresponds to a uniform distribution for the reflected specific intensity. After each time step dt/τ = 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/τ = 7.5 to t/τ = 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 . 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).

1362
The European Physical Journal Special Topics  Figure 1. Note that for the thinnest layer a logarithmic scale is employed on the radial axis. The ratio between the minimum (µ → 1) and maximum (µ → 0) value of the specific intensity is larger than 20 and 100 in case (a) and (c), respectively.

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 adminensional eigenvalue of the transport problem is τ −1 l = 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 nondimensionalized transport problem. This correlation between the leakage rate and the anisotropy of the specific intensity will be further illustrated below.

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 ≈ 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): 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: where ω 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 (ω = 2π) 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) 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 m/Zc ≈ 1.34. Left: flat free surface (case b). Right: rough free surface (case d). See Figure 1 for further details. free surface. Although these figures basically convey the same information, they shed light on different aspects of leakage. Figure 4, which represents τ l as a function of H * , will be particularly useful in the following to establish a link between the nondimensional 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 τ 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 τ l to Q c . Media with the same l * but different anisotropy factors g have different mean free time τ , and therefore different Q c according to formula (7). Figure 4 shows that the non-dimensional leakage time τ l depends monotonically on the non-dimensional crustal thickness H * and the anisotropy factor g. Based on the scaling of τ 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 τ 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 τ l ∝ (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 τ 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 τ l ∝ 1 (i.e., τ 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 plays the role of limiting factor which suggests the scaling τ l ∝ 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 τ 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 ≈ 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 * , τ 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 −1 l increases with g.

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 6 for anisotropic scattering. The correlation of the fluctuations is exponential and the anisotropy factor is g = 0.8.
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 τ 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 τ 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.

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.