# Coherent backscattering of acoustic waves in the near field

L. Margerin, M. Campillo, B. A. van Tiggelen

► To cite this version:

L. Margerin, M. Campillo, B. A. van Tiggelen. Coherent backscattering of acoustic waves in the near field. Geophysical Journal International, 2001, 145, pp.593-603. 10.1046/j.0956-540X.2001.01400.x. insu-03606662

HAL Id: insu-03606662

https://insu.hal.science/insu-03606662

Submitted on 12 Mar 2022

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
# Coherent backscattering of acoustic waves in the near field

L. Margerin,<fnref n="1" /> M. Campillo<fnref n="2" /> and B. A. van Tiggelen<fnref n="3" />

Accepted 2000 November 27. Received 2000 November 10; in original form 2000 April 10

## SUMMARY

Recent studies have shown that interference plays an important role in various phenomena observed for waves propagating through random media. Among these phenomena, the so-called cone of coherent backscattering has received much attention in optics. In this paper, we study analytically and numerically the coherent backscattering of acoustic waves in a seismological context. In particular, we focus on the near-field detection of scattered waves and the effect of transient sources. We show that interference results in an increase of the coda intensity as compared to the prediction of radiative transfer theory. After a transient regime, a spot of backscattering enhancement stabilizes in a sphere of radius half a wavelength centred at the source of seismic waves. Several effects such as absorption, boundary conditions and scattering anisotropy are investigated. Our study demonstrates the robustness of coherent backscattering and may offer a possible means of discriminating single versus multiple scattering in the observed coda.

**Key words:** multiple scattering, seismic coda.

## 1 INTRODUCTION

Coherent backscattering of waves was first observed in electromagnetic wave propagation more than a decade ago by Kuga & Ishimaru (1984). Since then, several experiments conducted in the laboratory have confirmed this phenomenon and clarified its origin (e.g. van Albada & Lagendijk 1985; Maret & Wolf 1985). We briefly describe these studies and refer to Corey *et al.* (1995) for a clear and detailed introduction to coherent backscattering. In optics, experiments are usually performed in the Fourier domain. The source consists of a monochromatic plane wave with wavevector $k_i$ that illuminates a disordered sample, and the backscattered waves are detected in the far field in direction $k_{sc}$. In this case, the backscattered intensity, $I(\theta)$, is measured as a function of the angle $\theta$ between $-k_i$ and $k_{sc}$. When $\theta$ is large, the backscattered intensity varies smoothly with $\theta$ and is predicted extremely well by diffusion theory. When $\theta$ is less than a few milliradians, one observes a sharp increase of the intensity that culminates in the exact backscattering direction. This observation is now well understood as the constructive interference between two waves that follow reciprocal paths in a random medium. Ideally, the interference results in an enhancement of the intensity by a factor of 2, as compared to the prediction of diffusion theory.

In the conventional physical picture for wave propagation in random media, multiple scattering of light is analogous to a classical random walk. The average wave intensity in a disordered medium is described by the radiative transfer equation, analogous to the Boltzmann equation of the kinetic theory of gases. The observation of coherent backscattering shows that phase and interference are crucial to understanding the propagation of classical waves in random media. Coherent backscattering appears in a variety of physical situations such as light scattering by cold atoms (Labeyrie *et al.* 1999), or propagation of pulsed acoustic waves in 2-D random media (Tourin *et al.* 1997). For an interdisciplinary review on the importance of interference effects in disordered media, we refer to the proceedings of the NATO Advanced Study Institute on diffuse waves in complex media (Fouque 1999).

In this paper, we investigate coherent backscattering in a seismological context. Contrary to the experiments in optics and in acoustics, seismic sources are embedded in the medium, and the detectors sit on the surface of the Earth. Moreover, seismic sources radiate energy during a short lapse time and thus the problem is intrinsically time-dependent. Schutz & Toksöz (1993, 1994) have studied numerically the backscattering of elastic waves from a rough interface and found an enhancement of intensity in the

<footnote n="1">¹Department of Geosciences, Guyot Hall, Princeton University, Princeton, NJ 08544, USA. E-mail: margerin@princeton.edu</footnote>

<footnote n="2">²Laboratoire de Géophysique Interne et Tectonophysique, Université Joseph Fourier, Grenoble, France</footnote>

<footnote n="3">³Laboratoire de Physique et Modélisation de la Matière Condensée, Université Joseph Fourier CNRS, Maison des Magistères, Grenoble, France</footnote>
backwards direction. Their results suggest the possibility of an increase of the coda energy resulting from waves diffracted at the Moho. In this paper, we consider scalar wave propagation in an acoustic medium and focus on the near-field detection around point-like sources. The full elastic problem will be addressed in a separate publication.

## 2 LARGE LAPSE TIME ANALYSIS

We illustrate the principle of coherent backscattering enhancement in Fig. 1, where two realizations of an ensemble of random media containing point scatterers are shown. For simplicity, only four scatterers are represented, denoted by letters A, B, C and D. A wave packet (solid line) is radiated from the source, S, and undergoes several scatterings before it is detected at the receiver, R. Its reciprocal counterpart (dashed line) visits the same scatterers in reverse order. The complex amplitudes of the direct and reciprocal waves are denoted by $A_d$ and $A_r$, respectively. Because these two waves encounter the same scatterers in reverse order, the angles of scattering and the total propagation distances are identical, ensuring that these waves have the same phase, geometrical spreading factor and scattering amplitude, so that $A_d = A_r$, in agreement with the reciprocity theorem. Conventional radiative transfer theory assumes that the mean intensity is equal to $\langle|A_d|^2 + |A_r|^2\rangle$, where $\langle\cdot\rangle$ denotes ensemble averaging. This expression is correct when the phases of the waves are uncorrelated. However, when source and receiver coincide, the two reciprocal wave packets are in phase and interfere constructively, independent of the particular realization of the random medium. Interferences will thus persist after ensemble averaging and the true intensity becomes $\langle|A_d|^2 + |A_r|^2 + 2\mathcal{R}_e(A_d A_r^*)\rangle$, where the last term accounts for interference and $\mathcal{R}_e(z)$ denotes the real part of $z$. Since $A_d = A_r$, the true intensity is exactly twice as large as the intensity predicted by conventional radiative transfer theory. If source and receiver do not coincide, the phases of the two reciprocal waves are not necessarily equal and will depend on the particular realization of the random medium. The coherent backscattering effect is thus expected to decrease with source-receiver distance. Coherent backscattering goes beyond the classical transport theory, which considers the two reciprocal wave paths but not the interference effect between them. We conclude this discussion with two important remarks.

(i) No matter how complex the scattering path of the wave, the phase difference between the reciprocal waves will depend only on the position of the first and last scatterings with respect to the source and receiver. This will be helpful in understanding the theoretical expressions below.

(ii) If only one scattering event is involved, no backscattering enhancement is possible since the direct and reciprocal paths cannot be distinguished. Coherent backscattering can therefore be observed only if multiple scattering is dominant.

Our approach to coherent backscattering is based on the theoretical studies by Akkermans *et al.* (1988) and van der Mark *et al.* (1988). These authors showed that the total intensity in a random medium is the sum of two contributions: an incoherent intensity, $I_{\text{inc}}$, that is simply the solution of the radiative transfer equation, plus a coherent intensity, $I_{\text{coh}}$, that takes interference effects into account. To find the expressions for both terms, we need to introduce some notation. $\mathbf{R}_0$, $\mathbf{R}$, $\mathbf{R}_1$ and $\mathbf{R}_n$ are the position vectors of the source, receiver, first scatterer and last scatterer respectively. $G(\mathbf{R}_2, \mathbf{R}_1)$ is the mean Green function for a source at $\mathbf{R}_1$ and a detector at $\mathbf{R}_2$ at frequency $\omega$, and $G^*(\mathbf{R}_2, \mathbf{R}_1)$ denotes its complex conjugate. $P(\mathbf{R}_2, \mathbf{R}_1, t)$, the propagator of the intensity in the random

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Illustration of the enhanced backscattering effect. The path of a wave packet and its reciprocal counterpart are depicted by solid and dashed lines respectively. Scatterers are denoted by letters A, B, C and D. The source is located at S and the energy is detected at the receiver R. Configurations 1 and 2 correspond to two realizations of an ensemble of random media. When source and receiver coincide, the phases of both waves are equal, independent of the particular realization. This results in complete constructive interference.</figcaption></figure>
medium, is proportional to the Green function of the radiative transfer equation for a source at $\mathbf{R}_1$ and a detector at $\mathbf{R}_2$; $t$ is the time elapsed since the energy release at $\mathbf{R}_1$. For simplicity, we consider an isotropic point source of energy embedded in the medium. For the moment, we neglect the presence of interfaces such as the free surface or the Moho. According to Akkermans et al. (1988), the incoherent and coherent components of the intensity can be expressed as

$$I_{\text{inc}}(\mathbf{R}, \mathbf{R}_0, t) \approx \iint d^3\mathbf{R}_1 d^3\mathbf{R}_n G(\mathbf{R}, \mathbf{R}_n) G^*(\mathbf{R}, \mathbf{R}_n) P(\mathbf{R}_n, \mathbf{R}_1, t) G(\mathbf{R}_1, \mathbf{R}_0) G^*(\mathbf{R}_1, \mathbf{R}_0), \quad (1)$$

$$I_{\text{coh}}(\mathbf{R}, \mathbf{R}_0, t) \approx \iint d^3\mathbf{R}_1 d^3\mathbf{R}_n G(\mathbf{R}, \mathbf{R}_n) G^*(\mathbf{R}, \mathbf{R}_1) P(\mathbf{R}_n, \mathbf{R}_1, t) G(\mathbf{R}_1, \mathbf{R}_0) G^*(\mathbf{R}_n, \mathbf{R}_0). \quad (2)$$

These formulae can be given a physical interpretation by reading the integrands from right to left. Eq. (1) corresponds to the following physical picture. The source emits a wave packet that travels to point $\mathbf{R}_1$, where it is scattered for the first time. Its energy is then transported through the random medium to point $\mathbf{R}_n$ according to radiative transfer theory. At $\mathbf{R}_n$, it is scattered for the last time and emits a wave packet to the detector. Similarly, expression (2) represents the interference term between two waves travelling in opposite directions, the direct and reciprocal paths being represented by the Green function and its complex conjugate respectively. The presence of scatterers causes any wave packet to lose energy while propagating through the random medium. As a consequence, the Green's function $G$ decays spatially as $\exp(-r/l)$, where $l$ is the scattering mean free path and $r$ is the propagation length (see Sato & Fehler 1998 for further details).

Formulae (1) and (2) neglect the propagation times from the source and receiver to the first and last scatterings and are therefore valid only when $t \gg \tau$, where $\tau=l/v$ is the scattering mean free time, and $v$ is the wave velocity. For large lapse times, the Green function of the radiative transfer equation $P(\mathbf{R}, \mathbf{R}', t)$ can be approximated by the solution of a diffusion equation (Lagendijk & van Tiggelen 1996). In that case, the propagators take a simple form and the integrals (1) and (2) can be solved analytically. A detailed analysis is given in Appendix A and the final results are

$$I_{\text{inc}}(\mathbf{R}, \mathbf{R}_0, t \to \infty) \sim t^{-3/2}, \quad (3)$$

$$(I_{\text{coh}} + I_{\text{inc}})(\mathbf{R}, \mathbf{R}_0, t \to \infty) \approx I_{\text{inc}} \left( 1 + \frac{\sin^2(k|\mathbf{R} - \mathbf{R}_0|)}{(k|\mathbf{R} - \mathbf{R}_0|^2)} e^{-|\mathbf{R} - \mathbf{R}_0|/l} \right), \quad (4)$$

where $k$ is the wavenumber. For crustal propagation in the 1–15 Hz frequency band, the mean free path roughly ranges from 20 to 200 km (Sato & Fehler 1998) and is therefore much larger than the wavelength of the probing wave. Our analysis is limited to the case $kl \gg 1$, a priori valid in the crust. Recent theoretical and experimental studies (Sheng 1995; Wiersma et al. 1997) suggest that for $kl \sim 1$, interferences can completely block the transport of energy in the medium, a phenomenon known as wave localization.

In Fig. 2, we show the theoretical ratio $(I_{\text{coh}} + I_{\text{inc}})/I_{\text{inc}}$ as a function of source-station distance, assuming that $l \approx 10\lambda$. From eq. (3), we conclude that the distribution of incoherent energy around the source eventually becomes homogeneous as a consequence of energy diffusion. However, Fig. 2 illustrates that interference effects play an important role within a sphere of radius half a wavelength centred at the source. The terminology ‘cone’, widely used in optics, is rather inappropriate for seismology and we will define a ‘spot of coherent backscattering’ or simply a ‘spot’ instead. Because $l$ is assumed to be much larger than $\lambda$, the dependence of coherent backscattering on the mean free path, as predicted by eq. (4), is rather weak.

The asymptotic analysis can be generalized to the case of a slightly anelastic medium by introducing a phenomenological absorption length $l_a$ and a slightly different wavenumber, $k_a$, due to dispersion. We show in Appendix A that eq. (4) still holds true provided one substitutes $1/l$ with $1/l+1/l_a$, and $k$ with $k_a$. The factor of 2 enhancement at the source still applies because the two reciprocal waves are identically attenuated by anelasticity. We note that our results differ markedly from those obtained in acoustics and optics. For example, Tourin et al. (1997) have shown that when the source and receivers are located far outside the scattering medium, the width of the zone of backscattering enhancement varies as $1/k\sqrt{vt}(t \to \infty)$, in sharp contrast with eq. (4), which predicts a stabilization in time. In the next sections, we address in more detail the convergence time of the spot of coherent backscattering towards its asymptotic shape.

### 3 MONTE CARLO SIMULATION OF COHERENT BACKSCATTERING

In this section, we study numerically coherent backscattering using a Monte Carlo method that can cope with various effects such as boundary reflections, scattering anisotropy and time dependence. Our numerical scheme is basically identical to that presented in other papers (Hoshiba 1995, 1997; Margerin et al. 1998). We briefly recall the ingredients of the simulation and address in more detail the calculation of the coherent intensity in Appendix B.

A wave packet or particle is launched at the source and walks randomly in the medium (see Fig. 3). The step length between two scattering events is determined by an exponential probability law $1/l \exp(-r/l)$. At each scattering event, the particle changes its direction. The new propagation direction is selected from the differential scattering cross-section, $d\sigma(\theta, \phi)/d\Omega$, which represents the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of the shape of coherent backscattering predicted by diffusion theory and numerical experiments. The backscattering enhancement ($I_{\text{coh}}+I_{\text{inc}}$)/$I_{\text{inc}}$ is plotted as a function of the source-station distance in terms of the wavelength. The Monte Carlo simulation results have been averaged in a time window extending from 17 to 20 mean free times.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Schematic explanation of the Monte Carlo simulation. The particle leaves the source $R_0$ and walks randomly in the medium. At the receiver $R$ we must take into account the interference between the direct path (solid line) and the reciprocal path (dashed line). This requires knowledge of the first and last scattering positions $R_1, R_n$, as well as the first and last scattering angles $\theta_1, \theta_n, \theta'_1, \theta'_n$.</figcaption></figure>

amount of energy radiated by the scatterer in the direction $(\theta, \phi)$ per unit solid angle and unit incident flux. The angles $\theta$ and $\phi$ refer to latitude and longitude in a spherical coordinate system whose polar axis points in the direction of propagation before scattering. Once these two angles have been selected in the local coordinate system, a rotation allows us to keep track of the motion of the particle in a global coordinate system.

At each scattering, the energy contribution of the particle at the receiver has to be estimated. Usually this is done by calculating the probability of the particle reaching a small area surrounding the detector. However, this approach cannot be generalized to interference effects because one has to estimate a quantity that is not necessarily positive. To cope with this problem, we adopt a point of view that considers the particles as real wave packets, capable of interfering. Appendix B discusses this approach and provides expressions for the estimate of the total intensity, $E[I_t]$, which can be expressed as a product of the classical estimate of the incoherent intensity, $E[I_i]$, times an interference factor,

$$E[I_t] = E[I_i] \left( 1 + \frac{|\mathbf{R}_1 - \mathbf{R}_0| |\mathbf{R} - \mathbf{R}_n| f(\theta'_1) f(\theta'_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}| + |\mathbf{R}_n - \mathbf{R}_0|)]}{|\mathbf{R}_1 - \mathbf{R}| |\mathbf{R}_n - \mathbf{R}_0| f(\theta_1) f(\theta_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R} - \mathbf{R}_n|)]} \right). \quad (5)$$

The geometry and notations are illustrated in Fig. 3. The interference term is a function of the position of the first and last scatterings (denoted by $\mathbf{R}_1$ and $\mathbf{R}_n$) and the corresponding scattering angles for the direct wave (denoted by $\theta_1, \theta_n$) and for the reciprocal wave (denoted by $\theta'_1, \theta'_n$); $f(\theta)$ is a weighting function for the scattered amplitude as a function of the scattering angle $\theta$. Because we assume that the medium is statistically isotropic, the scattering amplitude does not depend on the longitude $\phi$. The function $f(\theta)$ is easily obtained from the solution of the single-scattering problem and is simply related to the differential scattering cross-section, as explained below. The interference factor gives an enhancement of the intensity by a factor of 2 at the source. However, when the receiver is located a few wavelengths away from the source, the interference factor starts to oscillate with $\mathbf{R}_1$ and $\mathbf{R}_n$, indicating that constructive interference disappears in this configuration. The final step of the simulation consists of averaging the results of many independent random walks.

In the numerical simulation, we consider a point-like, isotropic, and instantaneous source embedded at $\mathbf{R}_0=0$ in a uniform random medium. The receiver is located at $\mathbf{R}$. The anisotropy of scattering is described by the normalized phase function, which is related to the differential scattering cross-section by

$$\Phi(\theta, \phi) = \frac{d\sigma(\theta, \phi)/d\Omega}{\int_{4\pi} \frac{d\sigma(\theta, \phi)}{d\Omega} d\Omega} \quad (6)$$
(van de Hulst 1981). In an acoustic medium with velocity fluctuations described by a Gaussian correlation function, the phase function is independent of $\phi$ and assumes the form

$$ \Phi(\theta) = \frac{\mu}{4\pi(1 - e^{-\mu})} e^{-\mu \sin^2(\theta/2)} \quad (7) $$

(Rytov et al. 1989). The parameter $\mu$ is related to the wavenumber $k$ and the correlation length of fluctuations $a$ through the equation $\mu = 2k^2a^2$. We choose this particular function because it depends on only one parameter, $\mu$, and is convenient to simulate. The weighting function $f$ introduced in eq. (5) is related to the phase function $\Phi$ by $f^2(\theta) = \Phi(\theta)$. For $\mu=0$, scattering is isotropic and for $\mu > 0$, scattering is predominantly forward. The strength of the anisotropy can be determined from the anisotropy factor,

$$ \langle \cos \theta \rangle = \int_{4\pi} \Phi(\theta) \cos(\theta) d\Omega \quad (8) $$

(van de Hulst 1981). The simulations have been performed for $\mu=0$ (isotropic scattering), $\mu=3$ ($\langle \cos \theta \rangle = 0.43$), $\mu=6$ ($\langle \cos \theta \rangle = 0.67$) and $\mu=10$ ($\langle \cos \theta \rangle = 0.80$), which corresponds to increasing scattering anisotropy. Unless explicitly stated, we assume a wavelength of 3 km and a scattering mean free path of 30 km, which are realistic values for seismology. In Fig. 2, we show the shape of the spot of coherent backscattering obtained numerically at large lapse times in the case of anisotropic scatterers ($\mu=3$), together with the predictions of formula (4). To get rid of numerical fluctuations, the results of the simulation have been averaged in a time window running from 17 to 20 mean free times. The difference between the analytical and numerical results is a few per cent only, showing consistency between the two approaches.

## 4 TIME DEPENDENCE OF COHERENT BACKSCATTERING

Fig. 4 illustrates the evolution of coherent backscattering with time for isotropic scatterers. The incoherent background intensity $I_{inc}$ and the total intensity ($I_{coh} + I_{inc}$) are shown at different lapse times in terms of the mean free time as a function of the source-station distance in terms of the wavelength. As is clear from Fig. 4, the intensity predicted by conventional radiative transfer theory is almost constant within a distance of two wavelengths from the source. This is in contrast to the true intensity, which exhibits an interference pattern that grows in time and eventually stabilizes around the source, as expected from eq. (4). Although the intensity decays in time by three orders of magnitude, the enhancement effect persists near the source, which may enable a possible observation of coherent backscattering in the coda of local earthquakes.

In Figs 5 and 6, the backscattering enhancement ($I_{coh} + I_{inc}$)/$I_{inc}$ is plotted as a function of time for five receivers spread within one wavelength around the source. Scattering is isotropic in Fig. 5 and moderately anisotropic ($\mu=3$) in Fig. 6. The wiggles on the curves have no physical origin and are caused by an incomplete averaging in the Monte Carlo simulations. The time dependence of the enhancement effect is not simple. In particular, curves in Fig. 6 exhibit an overshoot. Except at $R=0$, the enhancement level rapidly increases, reaches a maximum and then slowly decreases toward an asymptotic value. The position of the maximum moves towards shorter lapse times as the source-receiver distance increases. Very close to the source, the position of the maximum becomes difficult to evaluate. To define a characteristic time of stabilization of the spot as a function of the scattering properties in the medium, we need to understand the origin of the time dependence shown in Fig. 6.

As explained earlier in this paper, singly scattered waves have no reciprocal counterpart and therefore do not contribute to the coherent intensity. However, they do contribute to the incoherent intensity for lapse times shorter than a few mean free times (see e.g. Hoshiba 1991), and therefore the maximum enhancement at the source can be reached only after the energy of the singly scattered waves has become negligible. To understand the role played by single scattering, it is convenient to subtract its contribution from the incoherent intensity. The results are shown in Fig. 7, where the backscattering enhancement obtained for $\mu=3$ is plotted as a function of time. At the source ($R=0$), we observe that the maximum enhancement factor is reached almost immediately. This shows that the time evolution at $R=0$ is mostly governed by the single-scattering term. This property could be used to measure the backscattering coefficient introduced by Aki & Chouet (1975). When the source and receiver coincide, the single-scattering intensity, $I_1(t)$, can be expressed as

$$ I_1(t) = \frac{g(\pi)}{2\pi vt^2} e^{(-vt/l)}, \quad (9) $$

where $g(\pi)$ is the backscattering coefficient, and a unit energy release is assumed. In our notation, we have

$$ g(\pi) = \frac{4\pi\Phi(\pi)}{l}. \quad (10) $$

Apart from the wave velocity (which is usually known), $g(\pi)$ is the only free parameter governing the single-scattering intensity. Therefore, a measurement of the time dependence of the enhancement factor exactly at the source could in principle provide $g(\pi)$. In optics, the influence of single scattering on coherent backscattering has also been investigated theoretically by Mischenko (1992) and experimentally by Wiersma et al. (1995).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Snapshot of coherent backscattering at different lapse times. The total intensity, which includes interference effects (solid lines), and the incoherent intensity (dash-dotted lines) are shown as a function of the source-receiver distance for different lapse times *T* (in terms of the mean free time τ). The unit distance on the horizontal axis is the wavelength.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Backscattering enhancement ($I_{\text{coh}} + I_{\text{inc}}/I_{\text{inc}}$) as a function of time for isotropic scattering. The different curves correspond to different source-station distances *R* in terms of the central wavelength $\lambda$. The time unit is the mean free time of the waves. The spot of coherent backscattering is seen to stabilize after about 15 mean free times.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Same as Fig. 5 for anisotropic scattering ($\mu=3$).</figcaption></figure>

When source and receiver do not coincide ($R \neq 0$), the time dependence is more complicated. The position of the maximum enhancement shown in Fig. 7 is now shifted towards shorter lapse times as compared to Fig. 6. Beyond this maximum, the enhancement curves exhibit a monotonic decay towards their asymptotic value. The comparison of Figs 7 and 6 also shows that the spot stabilizes after about 15 mean free times, whether the single-scattering term is subtracted or not. This leads to the conclusion that the stabilization time of coherent backscattering has no relation to single scattering, except at exactly $R=0$.

We expect our simulations to match the result of the long lapse time analysis only when the diffusive regime is reached. Therefore, it is important to understand how the diffusion constant *D* influences the stabilization time. *D* is related to the anisotropy factor $\langle\cos\theta\rangle$ and to the wave velocity *v* through $D = vI^*/3$, where $I^* = l/(1 - \langle\cos\theta\rangle)$. The diffusion constant is thus fully determined by the wave velocity, the scattering mean free path and the anisotropy factor. $I^*$ is often termed the transport mean free path and physically represents the length beyond which a 'random walker' has lost memory of its initial direction of propagation. We refer to Sheng (1995) and Margerin *et al.* (1998) for further details.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Effect of single scattering on the stabilization of coherent backscattering. This figure is similar to Fig. 6, except that the single-scattering term has been subtracted from the total intensity. At the source, the asymptotic value is reached immediately. The convergence time is unchanged elsewhere.</figcaption></figure>

In Fig. 8 we show the time dependence of the backscattering enhancement at $R=0$, $\lambda/2$ and $\lambda$ for three different values of the diffusion constant. The three cases correspond to the same value of the scattering mean free path but increasing values of the constant $\mu$ defined in eq. (7). As $\mu$ increases, the scattering becomes more and more strongly peaked in the forward direction and $\langle \cos \theta \rangle$ increases. Accordingly, the diffusion constant $D$ also increases with $\mu$. Recall that at $R=\lambda/2$ and $\lambda$, the asymptotic value of the backscattering enhancement is exactly 1. We observe that (except for $R=0$) the speed of stabilization of coherent backscattering increases with $\mu$. This can be physically understood by the fact that the characteristic time to achieve complete isotropy of the wavefield, and thus complete diffusion, increases with the scattering anisotropy. We can roughly estimate this characteristic time as a function of the transport mean free time $\tau^*=l^*/v$. Our simulations suggest that the spot has become stable after five to 10 transport mean free times.

# 5 EFFECT OF THE CRUST

In this section, we briefly discuss the effect of geometry on coherent backscattering. We consider a simplified crustal model with a wave speed jump at the Moho with usual values for S waves in the crust (3.5 km s⁻¹) and mantle (4.7 km s⁻¹). The mantle is assumed to be perfectly transparent, whereas the crust is assumed to be 20 km thick and very heterogeneous with a scattering mean free path $l=20$ km. This set of physical parameters is suggested by previous studies of the coda (Margerin *et al.* 1999) and gives rise to a strong leakage of energy into the mantle. The numerical method to solve the radiative transfer equation in a waveguide geometry was published previously (Hoshiba 1997; Margerin *et al.* 1998). The calculation of the coherent intensity requires one to take into account all the reciprocal paths corresponding to multiple reflections at the boundaries of the medium. The simulation is therefore more complex but the modifications are straightforward to incorporate. In Fig. 9, we show the shape of the spot of coherent backscattering for a lapse time roughly equal to 20 mean free times and isotropic scatterers. It compares fairly well with the analytical solution for the full-space case. This means that the shape of the spot of coherent backscattering is very robust to changes in the boundary conditions. This reduces the number of free variables considerably and favours a possible observation in seismically active regions.

# 6 CONCLUSIONS

We have studied numerically and theoretically the coherent backscattering of waves in a seismological context. Although our investigation is limited to acoustic waves, some important effects such as near-field detection and point-like sources have been taken into account. Our analytical theory provides an exact asymptotic description of coherent backscattering and predicts a coherent intensity that oscillates and decays rapidly with distance from source to receiver. The enhancement persists in time and should be observable as long as a coda is measurable. Coherent backscattering affects a sphere centred at the source of radius half a wavelength, typically 500 m for 3 Hz waves.

Our numerical study gives access to the time dependence of the coherent backscattering effect. For large lapse times, the spot stabilizes as predicted by diffusion theory. This stabilization is a specific feature of near-field detection. In media with a
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Time evolution of the spot of coherent backscattering for strong (solid lines), moderate (dash-dotted lines) and weak (dotted lines) scattering anisotropy. One observes that except at the source, a strong forward anisotropy delays the stabilization of the spot.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Effect of strong energy leakage on the spot of coherent backscattering. The thickness of the crust and the scattering mean free path both equal 20 km. We assume a usual crust (3.5 km s⁻¹)/mantle (4.7 km s⁻¹) velocity contrast at the Moho. Note the robustness of coherent backscattering to changes in the boundary conditions.</figcaption></figure>

constant-scattering mean free path but with increasing scattering anisotropy, the convergence time tends to increase. The estimated time from the simulation is about five to 10 transport mean free times. In the ideal case where detection takes place at the source, the time dependence is influenced by single scattering only. The coherent backscattering effect is stable against changes in the boundary conditions and is preserved in a waveguide geometry.

We finally comment on the possibility of revealing experimentally the existence of coherent backscattering in seismology. Since this effect is measurable only in the multiple-scattering regime and in the vicinity of the source, experiments should be set up in very heterogeneous regions with very shallow earthquakes. Volcanoes seem to be good candidates to fulfil these requirements. Seismologists only have access to one realization of the random medium. Because our calculations are valid in an ensemble average sense only, the measurement of the spot of coherent backscattering requires the ensemble average to be replaced by a time average. The observation of the spot of coherent backscattering would provide direct experimental proof that the coda of earthquakes is caused by multiple scattering of elastic waves.

## ACKNOWLEDGMENTS

We would like to thank S. H. Hung and J. R. Grasso for useful comments. Valuable reviews by J. Kawahara and S. N. Ward helped a lot to improve our manuscript.
REFERENCES

Aki, K. & Chouet, B., 1975. Origin of coda waves: source, attenuation and scattering effects, *J. geophys. Res.*, **80**, 3322–3342.

Akkermans, E., Wolf, P.E., Maynard, R. & Maret, G., 1988. Theoretical study of the coherent backscattering of light by disordered media, *J. Phys. France*, **49**, 77–98.

Born, M. & Wolf, E., 1970. *Principles of Optics*, 4th edn, Pergamon Press, Oxford.

Corey, R., Kissner, M. & Saulnier, P., 1995. Coherent backscattering of light, *Am. J. Phys.*, **63**, 560–564.

Fouque, J.P., ed., 1999. *Diffuse Waves in Complex Media, NATO Sci. Ser.*, Vol. 531, Kluwer, Dordrecht, the Netherlands.

Goodman, J.W., 1985. *Statistical Optics*, Wiley, New York.

Hoshiba, M., 1991. Simulation of multiple-scattered coda wave excitation based on the energy conservation law, *Phys. Earth Planet. Inter.*, **67**, 123–136.

Hoshiba, M., 1995. Estimation of nonisotropic scattering in western Japan using coda wave envelopes: application of a multiple non-isotropic scattering model, *J. geophys. Res.*, **100**, 645–657.

Hoshiba, M., 1997. Seismic coda wave envelope in depth dependent S-wave velocity structure, *Phys. Earth planet. Inter.*, **104**, 15–22.

Kuga, Y. & Ishimaru, A., 1984. Retroreflection from a dense distribution of spherical C. particles, *J. Opt. Soc. Am.*, **A1**, 831–835.

Labeyrie, G., de Tomasi, F., Bernard, J.-C., Mueller, C.A., Miniatura, C. & Kaiser, R., 1999. Coherent backscattering of light by cold atoms, *Phys. Rev. Lett.*, **83**, 5266–5269.

Lagendijk, A. & van Tiggelen, B.A., 1996. Resonant multiple scattering of light, *Phys. Repts*, **270**, 143–215.

Maret, G. & Wolf, P., 1985. Weak localization and coherent backscattering of photons in disordered media, *Phys. Rev. Lett.*, **55**, 2696–2699.

Margerin, L., Campillo, M. & van Tiggelen, B.A., 1998. Radiative transfer and diffusion of waves in a layered medium: new insight into coda, *Geophys. J. Int.*, **134**, 596–612.

Margerin, L., Campillo, M. & van Tiggelen, B.A., 1999. 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–352.

Mischenko, M.I., 1992. Enhanced backscattering of polarized light from discrete random media: calculations in exactly the backscattering direction, *J. Opt. Soc. Am.*, **A9**, 978–982.

Rytov, S.M., Kravtsov, Y.A. & Tatarskii, V.I., 1989. *Principles of Statistical Radiophysics*, Vol. 4, *Wave Propagation through Random Media*, Springer-Verlag, Berlin.

Sato, H. & Fehler, M.C., 1998. *Seismic Wave Propagation and Scattering in the Heterogeneous Earth*, Springer-Verlag, New York.

Schutz, C.A. & Toksöz, M.N., 1993. Enhanced backscattering of seismic waves from a highly irregular, random interface: SH case, *Geophys. J. Int.*, **114**, 91–102.

Schutz, C.A. & Toksöz, M.N., 1994. Enhanced backscattering of seismic waves from a highly irregular, random interface: P-SV case, *Geophys. J. Int.*, **117**, 783–810.

Sheng, P., 1995. *Introduction to Wave Scattering, Localization and Mesoscopic Phenomena*, Academic Press, San Diego.

Tourin, A., Derode, A., Roux, P., van Tiggelen, B.A. & Fink, M., 1997. Time dependent coherent backscattering of acoustic waves, *Phys. Rev. Lett.*, **79**, 3637–3639.

van Albada, M. & Lagendijk, A., 1985. Observation of weak localization of light in a random medium, *Phys. Rev. Lett.*, **55**, 2692–2695.

van de Hulst, H.C., 1981. *Light Scattering by Small Particles*, Dover, New York.

van der Mark, M.B., van Albada, M.P. & Lagendijk, A., 1988. Light scattering in strongly scattering media: multiple-scattering and weak localization, *Phys. Rev. B*, **37**, 3575–3592.

Wiersma, D.S., van Albada, M.P. & van Tiggelen, B.A., 1995. Experimental evidence of recurrent scattering in disordered media, *Phys. Rev. Lett.*, **74**, 4193–4197.

Wiersma, D.S., Bartolini, P., Lagendijk, A. & Righini, R., 1997. Localization of light in a disordered medium, *Nature*, **390**, 671–673.

APPENDIX A: EVALUATION OF THE COHERENT AND INCOHERENT INTENSITIES

We begin with the simple case of a non-dissipative infinite medium. The Green functions and intensity propagators defined in eqs (1) and (2) read

$$G(\mathbf{R}_2, \mathbf{R}_1) = - \frac{\exp(ik|\mathbf{R}_2 - \mathbf{R}_1| - |\mathbf{R}_2 - \mathbf{R}_1|/2l)}{\sqrt{4\pi|\mathbf{R}_2 - \mathbf{R}_1|}}, \quad (A1)$$

$$P(\mathbf{R}_2, \mathbf{R}_1, t) = \frac{v \exp(-|\mathbf{R}_2 - \mathbf{R}_1|/4Dt^2)}{l^2 (4\pi Dt)^{3/2}}. \quad (A2)$$

The Green function in eq. (A1) represents the mean wavefield emitted by a source of unit energy. D is the diffusion constant of the waves related to the transport mean free path through the relation $D = vl^*/3$ as explained in the text. Each wave packet is assumed to be slowly modulated in phase and amplitude, implying a finite frequency band $\Delta\omega$, with central frequency $\omega \gg \Delta\omega$ and central wavenumber $k$. The corresponding coherence time of the wave packet $1/\Delta\omega$ will be large compared to the oscillation time $1/\omega$ so that the interference term will be well approximated by considering monochromatic waves only. We refer to Born & Wolf (1970) and Goodman (1985) for a discussion of interference properties of finite bandwidth signals. When $Dt \gg l^2$, eq. (1) can be approximated by

$$I_{\text{inc}}(\mathbf{R}, \mathbf{R}_0, t) \approx \frac{v}{l^2(4\pi Dt)^{3/2}} \iint \frac{e^{-|\mathbf{R}-\mathbf{R}_n|/l} e^{-|\mathbf{R}_0-\mathbf{R}_1|/l}}{4\pi|\mathbf{R}_0-\mathbf{R}_1|^2 4\pi|\mathbf{R}-\mathbf{R}_n|^2} d^3\mathbf{R}_n d^3\mathbf{R}_1. \quad (A3)$$

After integration, we find

$$I_{\text{inc}}(\mathbf{R}, \mathbf{R}_0, t) = \frac{v}{(4\pi Dt)^{3/2}}, \quad (A4)$$
as expected for diffuse propagation near the source. Again for $Dt \gg l^2$, the coherent intensity reads

$$I_{\text{coh}}(\mathbf{R}, \mathbf{R}_0, t) \approx \frac{v}{l^2(4\pi Dt)^{3/2}} \iint d^3\mathbf{R}_n d^3\mathbf{R}_1 \frac{\exp(-|\mathbf{R}_1 - \mathbf{R}_0|/2l + ik|\mathbf{R}_1 - \mathbf{R}_0|) \exp(-|\mathbf{R}_0 - \mathbf{R}_n|/2l - ik|\mathbf{R}_0 - \mathbf{R}_n|)}{4\pi|\mathbf{R}_1 - \mathbf{R}_0||\mathbf{R}_0 - \mathbf{R}_n|} \\ \times \frac{\exp(-|\mathbf{R} - \mathbf{R}_n|/2l + ik|\mathbf{R} - \mathbf{R}_n|) \exp(-|\mathbf{R} - \mathbf{R}_1|/2l - ik|\mathbf{R} - \mathbf{R}_1|)}{4\pi|\mathbf{R} - \mathbf{R}_n||\mathbf{R} - \mathbf{R}_1|}, \quad (A5)$$

$$\approx \frac{v}{l^2(4\pi Dt)^{3/2}} |J|^2, \quad (A6)$$

where we have introduced

$$J = \int d^3\mathbf{R}_1 \frac{\exp[-(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R} - \mathbf{R}_1|)/2l + ik(|\mathbf{R}_1 - \mathbf{R}_0| - |\mathbf{R} - \mathbf{R}_1|)]}{4\pi|\mathbf{R} - \mathbf{R}_1||\mathbf{R}_1 - \mathbf{R}_0|}. \quad (A7)$$

This last integral is most easily computed by expressing the two Green functions in terms of their Fourier transforms,

$$\frac{G(\mathbf{R}, \mathbf{R}_0)}{\sqrt{4\pi}} = - \frac{\exp(ik|\mathbf{R} - \mathbf{R}_0| - |\mathbf{R} - \mathbf{R}_0|/2l)}{4\pi|\mathbf{R} - \mathbf{R}_0|} \quad (A8)$$

$$= \frac{1}{8\pi^3} \int d^3\boldsymbol{\kappa} \frac{\exp[i\boldsymbol{\kappa} \cdot (\mathbf{R} - \mathbf{R}_0)]}{\boldsymbol{\kappa}^2 - \boldsymbol{\kappa}'^2 + ik/l}. \quad (A9)$$

Inserting this last expression into eq. (A7), we obtain

$$J = \frac{1}{2^4 \pi^5} \int d^3 \mathbf{R}_1 \iint d^3 \boldsymbol{\kappa} d^3 \boldsymbol{\kappa}' \frac{\exp[i\boldsymbol{\kappa} \cdot (\mathbf{R}_1 - \mathbf{R}_0)] \exp[-i\boldsymbol{\kappa}' \cdot (\mathbf{R} - \mathbf{R}_1)]}{(k^2 - \boldsymbol{\kappa}^2 + ik/l)(k^2 - \boldsymbol{\kappa}'^2 - ik/l)}. \quad (A10)$$

Now assuming that the order of integration can be interchanged, we can easily integrate over $\mathbf{R}_1$, leaving us with

$$J = \frac{1}{2\pi^2} \iint d^3\boldsymbol{\kappa} d^3\boldsymbol{\kappa}' \frac{\delta(\boldsymbol{\kappa} + \boldsymbol{\kappa}') \exp[-i(\boldsymbol{\kappa}' \cdot \mathbf{R} + \boldsymbol{\kappa} \cdot \mathbf{R}_0)]}{(k^2 - \boldsymbol{\kappa}^2 + ik/l)(k^2 - \boldsymbol{\kappa}'^2 - ik/l)} \quad (A11)$$

$$= \frac{1}{2\pi^2} \int d^3\boldsymbol{\kappa} \frac{\exp[i\boldsymbol{\kappa} \cdot (\mathbf{R} - \mathbf{R}_0)]}{(k^2 - \boldsymbol{\kappa}^2)^2 + k^2/l^2}. \quad (A12)$$

Making use of the symmetry of the integrand, we perform the integration over $\boldsymbol{\kappa}$ in spherical coordinates with the polar axis oriented along $\mathbf{R}-\mathbf{R}_0$. After integration over the two polar angles, we obtain

$$J = -\frac{i}{\pi} \int_{-\infty}^{+\infty} d\kappa \frac{\kappa e^{ik|\mathbf{R}-\mathbf{R}_0|}}{(\kappa^2 - k^2)^2 + k^2/l^2}. \quad (A13)$$

To evaluate this last integral, we close the contour of integration in the upper sheet of the complex plane with a semi-circle of infinite radius, and apply the residue theorem. The poles $\kappa_1$ and $\kappa_2$ of the integrand are located at

$$\kappa_1 = k(1 + i/kl)^{1/2}, \quad \kappa_2 = -k(1 - i/kl)^{1/2}. \quad (A14)$$

In this last equation, the branch of the square root is chosen such that $Im(\kappa_1), Im(\kappa_2) > 0$, where $Im(z)$ denotes the imaginary part of $z$. For weak scattering ($kl \gg 1$), the poles can be approximated by

$$\kappa_1 = k + i/2l, \quad \kappa_2 = -k + i/2l. \quad (A15)$$

Since the integral over the semi-circle does not give any contribution, the final result reads

$$J = l \frac{\sin(k|\mathbf{R} - \mathbf{R}_0|)}{k|\mathbf{R} - \mathbf{R}_0|} e^{-|\mathbf{R} - \mathbf{R}_0|/2l}. \quad (A16)$$

Collecting together the results of eqs (A4), (A6) and (A16), we finally obtain the formula quoted in the text,

$$\frac{I_{\text{coh}} + I_{\text{inc}}}{I_{\text{inc}}} \approx 1 + \frac{\sin^2(k|\mathbf{R} - \mathbf{R}_0|)}{(k|\mathbf{R} - \mathbf{R}_0|)^2} e^{-|\mathbf{R} - \mathbf{R}_0|/l} \quad (t \to \infty). \quad (A17)$$

These computations can easily be extended to the case where absorption is present in the medium. We introduce phenomenologically an absorption length $l_a$, such as the mean Green functions decay spatially at the rate $1/l + 1/l_a$. The dispersion due to anelastic absorption also changes $k$ into $k_a$. We then see that the new poles are $\kappa_1 = k_a + i(1/l + 1/l_a)$ and $\kappa_2 = k_a + i(1/l + 1/l_a)$. Therefore, we just need to substitute $k$ with $k_a$ and $1/l$ with $1/l + 1/l_a$ in eq. (4).
APPENDIX B: MONTE CARLO ESTIMATE OF THE COHERENT INTENSITY

Let us consider the interference between two reciprocal wave packets such as those represented in Fig. 3. However, now, instead of one configuration, we consider the ensemble average response. As before we call $A_d$ the direct path and $A_r$ its reciprocal counterpart. The ith scattering occurs at $\mathbf{R}_i$, $f(\theta_i)$ is a weighting function for the scattered amplitude in the case of anisotropic scattering, and $\theta_i$ denotes the scattering angle of the ith scattering. For the reciprocal wave, we introduce $\theta'_i$ as shown in Fig. 3. This is necessary since the first and last scattering angles are not exactly the same for the two waves. The direct and reciprocal amplitudes at the receiver are, respectively,

$$A_d = C \frac{f(\theta_1)f(\theta_2)\cdots f(\theta_{n-1})f(\theta_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R}_2 - \mathbf{R}_1| + \cdots + |\mathbf{R}_n - \mathbf{R}_{n-1}| + |\mathbf{R}_0 - \mathbf{R}_n|)]}{|\mathbf{R}_1 - \mathbf{R}_0||\mathbf{R}_2 - \mathbf{R}_1|\cdots|\mathbf{R}_n - \mathbf{R}_{n-1}||\mathbf{R}_0 - \mathbf{R}_n|}, \quad (B1)$$

$$A_r = C \frac{f(\theta'_1)f(\theta'_2)\cdots f(\theta'_{n-1})f(\theta'_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R}_2 - \mathbf{R}_1| + \cdots + |\mathbf{R}_n - \mathbf{R}_{n-1}| + |\mathbf{R}_0 - \mathbf{R}_n|)]}{|\mathbf{R}_1 - \mathbf{R}_0||\mathbf{R}_2 - \mathbf{R}_1|\cdots|\mathbf{R}_n - \mathbf{R}_{n-1}||\mathbf{R}_0 - \mathbf{R}_n|}, \quad (B2)$$

C is an unimportant normalization constant introduced to account for the product of terms such as $1/4\pi$. Now we express the total intensity as

$$|A_d + A_r|^2 = |A_d|^2 + |A_r|^2 + A_d A_r^* + A_d^* A_r. \quad (B3)$$

The first two terms on the right-hand side of this equation are the incoherent intensities. These terms are accounted for by the radiative transfer theory and can be computed with the standard Monte Carlo method. The last two terms correspond to the interference effect. Note that in the simulation, $|A_d|^2$ and $|A_r|^2$ correspond to two different paths that will be simulated independently. Since the Monte Carlo scheme is assumed to exhaust all possible random walks, each path will eventually find its reciprocal counterpart in the simulation. In other words, to estimate the total intensity (that is, the sum of coherent and incoherent contributions), one should not evaluate the whole expression $|A_d + A_r|^2$ because this would be equivalent to erroneously including the same scattering path twice. Instead one needs just to calculate the incoherent intensity $|A_d|^2$ plus the real part of the interference term $A_d A_r^*$ (or $A_d^* A_r$). The estimate of the total intensity at the receiver $E[I]$ can therefore be expressed as the product of the classical estimate of the incoherent intensity $E[I]$ times an interference term,

$$E[I] = E[I] \left( 1 + \frac{A_d}{A_r} \right), \quad (B4)$$

$$E[I] = E[I] \left( 1 + \frac{|\mathbf{R}_1 - \mathbf{R}_0| |\mathbf{R} - \mathbf{R}_n| f(\theta'_1) f(\theta'_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R}_n - \mathbf{R}_0|)]}{|\mathbf{R}_1 - \mathbf{R}_0| |\mathbf{R}_n - \mathbf{R}_0| f(\theta'_1) f(\theta'_n) \exp[i(k + i/2l)(|\mathbf{R}_1 - \mathbf{R}_0| + |\mathbf{R}_n - \mathbf{R}_0|)]} \right). \quad (B5)$$

We can check that when the source and receiver coincide, the estimated intensity equals exactly twice the classical intensity. When the source and receiver are more than a few wavelengths apart, the interference term strongly oscillates and will finally average 1.