Éric Larose

Mesoscopics of ultrasound and seismic waves: application to passive imaging

published or not. The documents may come

Mesoscopics of ultrasound and seismic waves:

application to passive imaging

Résumé

Physique mésoscopique des ultrasons et des ondes sismiques : application à l'imagerie passive

Cet article de revue rassemble plusieurs aspects fondamentaux et appliqués de la propagation des ondes acoustiques et élastiques dans les milieux hétérogènes, en régime de diffusion simple ou multiple. Après une introduction sur les techniques conventionelles d'imagerie sismique et ultrasonore, nous présentons deux expériences qui mettent en évidence la présence de diffusion multiple dans les enregistrements sismologiques : l'équipartition des ondes, et la rétrodiffusion cohérente (Chap. 2). La diffusion multiple des ondes, qui démontre l'aspect mésoscopique de leur propagation, est une limitation majeure pour les techniques d'imagerie conventionelles (imagerie médicale, sismique réflexion ou réfraction, tomographie...). La deuxième partie du document (Chaps. [3][4][5] est consacrée à une application de cette physique mésoscopique : le principe de l'imagerie passive. En effectuant la corrélation temporelle d'enregistrement de bruit ambiant ou d'ondes diffuses, il est possible de reconstruire la réponse impulsionnelle du milieu entre deux capteurs passifs comme si l'on avait placé une source en lieu et place d'un des capteurs. Cela offre la possibilité de faire de l'acoustique ou de la sismologie sans source. Plusieurs aspects sont présentés dans ce manuscrit : des aspects théoriques et numériques (Chaps. 3,4), ensuite des aspects expérimentaux avec des applications (Chap. 5) à l'échelle des ultrasons (tomographie passive d'un milieu stratifié), et des applications à l'échelle de la sismologie (imagerie du sous-sol de la Californie, et même de la Lune).

Limitations of conventional imaging

During its propagation, an acoustic or seismic wave interacts with the medium: it can be slowed, accelerated, reflected or refracted. From these interactions, one has the possibility to image the medium. In this chapter, we recall the very basics of medical and seismic imaging, and point out several limitations. These limitations will be addressed in the following chapters of the present manuscript, and several solutions will be proposed.

Imaging techniques: the basics 1.Ultrasounds in medical imaging

Piezoelectric transducers are commonly used as both sources and receivers of acoustic waves at ultrasonic frequencies (MHz). Corresponding wavelengths are small enough (mm) to provide images at a resolution compatible with medical diagnosis. Therefore, from the late 1970's, ultrasound has been widely used in medical imaging and sometimes as an alternative to less benign techniques like X-ray imaging. An ultrasound image is based on a three step operation:

• a short pulse (typically one or two oscillations at central frequency of a few MHz) is emitted by one or several sources. The emitted wave might be focused on a given point to be investigated;

• each time the wave encounters a scatterer (an impedance contrast originating from a change in velocity or density), it is partially reflected. This is particularly valid for interfaces between two different tissues, but also within a given tissue (isolated scatterers like kidney stones, disease like artery calcification, etc.);

• backscattered waves (echos) are sensed by the array of transducers. By summing time-delayed records (beam-forming), one can choose to focus the waves on a given target. The intensity of each pixel forming the final image depends on the intensity of the wave backscattered by the focused point (see Fig. 1.1). Ideally, an ultrasonic image is a reflectivity map of the medium. This technique is efficient as long as the backscattered wave field interacts only once with the medium before being sensed by the sounding device. In other words, an ultrasonic image is based on the simple scattering assumption (or Born approximation). This approximation is valid for many soft tissues. When the simple scattering assumption is not valid, interpreting ultrasonic images is very difficult, and sometime impossible. For instance, Derode et al. [1] demonstrated the presence of multiple scattering of ultrasounds in a human trabecular bone. The validity of the simple scattering assumption appears as an essential point that needs to be checked more regularly. This multiple scattering phenomena concerns not only medical imaging but also non destructive testing, for instance in polycrystaline media [2], or in granular materials where several observations have confirmed the diffuse nature of wave propagation [3].

Active seismic experiments

In prospecting seismology where active sources are employed, two kinds of experimental setup can be distinguished. In a reflection seismic setup, like the one at sea displayed in Figure 1.2, the aim is to image the reflectivity of the medium from the backscattered waves. This setup is analogous to the one in ultrasonic medical imaging. The source is typically an explosion (from an airgun at sea, or from chemical explosive inland), and the wavefield is sensed by an array of geophones placed a few hundreds of meters around the source. From the arrival time of reflections, one can distinguish strong reflectors, especially when the soil is stratified. But to spatially localize the reflectors, we also need to know the celerity of the waves in the subsurface. A reflection image is sometime called the "high frequency" image of the subsurface because it enhances the interfaces between different soils, but does not tell anything about the geological layer itself.

In most locations, the Earth is heterogeneous and the wave velocities evolve very quickly from the surface (a few 100 m/s) to deeper layers (several km/s), typically over one order of magnitude. To evaluate the average velocity in a given layer, a refraction setup is needed. It is composed of a controlled source and an array of sensors placed far from the source. From the direct P and S refracted arrival, it is possible to estimate the wave velocity in each layer or geological unit. The obtained image is referred to as the "low frequency" image because it does not reveals the details of the interfaces but only the smoothed vertical dependence of the velocity.

The two images are eventually mixed, or migrated, to localize the reflector and remove artifacts due to reflection hyperbola and wave diffraction. The final image is obtained after an iterative inversion where the reflectivity and velocity parameters of each geological unit are adjusted. The end of the seismic record (on the left, after 5 s) shows a very short spatial coherence that reminds the optical speckle. This seismic speckle might originate from a random distribution of scatterers in the crust in the simple or multiple scattering regime. This part of the seismic record, referred to as the coda, is often neglected and not processed. At higher frequencies (right), the coda is shortened by the absorption. Such a seismic speckle is a clear evidence for the very heterogeneous nature of the Earth's crust.

Ann. Phys. Fr. 31

• N o 3 • 2006
As in medical ultrasonic imaging, the seismic image is obtained assuming simple scattering. In some cases, waves are multiply reflected within a given layer (seismic multiples). It is possible to subtract these multiples from the raw data because they are easily identified. But in most cases, the possibility of observing waves multiply scattered by bulk heterogeneities is neglected. Nevertheless, in very heterogeneous media like volcanoes, complex margin or orogenic zones, multiple scattering is probable, and will degrade the seismic images. In these cases, new imaging techniques are to be developed.

Global seismology

When we try to characterize and image the features of deeper layers in the Earth (from the bottom of the crust to the inner core), active experiments are inappropriate. Propagating waves across the Earth (12 000 km) requires to employ energetic low frequencies sources (from mHz to ≈10 Hz). The energy required to illuminate the deep Earth is incredibly high, roughly comparable to the one delivered by an M w > 5 Earthquake (several kilo-tons of TNT). Only atomic explosions may artificially deliver comparable energy at those frequencies. Therefore, imaging the deep earth relies on the occurrence of natural sources like earthquakes.

The natural vibrations of the Earth are sensed worldwide by different sets of permanent stations like the Global Seismic Network (GSN). The deployment of such broad band seismic network started about 30 years ago. Detection of the main features of the Earth's interior (mantle, core, inner core) and production of tomographic maps of the Earth lithosphere [4] was made possible by computation of arrival times of seismic waves originating from localized earthquakes. In Figure 1.4 is presented the 1D preliminary reference Earth model (PREM) resulting from such global observations. It was proposed in 1981 by Dziewonski et al. [5]. This model must be thought of as an average model of the Earth. Sometime, The aim of ambitious new deployments of seismic arrays, such as the Program for the Array Seismic Studies of the Continental Lithosphere (PASSCAL) and USArray programs, is to improve the resolution of images of Earth's interior by adding more instruments to regional-and continental-scale seismic networks. Nevertheless, traditional observational methods cannot fully exploit emerging array data because they are based on seismic waves emitted from earthquakes, which emanate from select source regions predominantly near plate boundaries (see Fig. 1.5). The imperfect locations of earthquakes is a clear limitation of traditional seismic imaging techniques. Improving the images of the earth will require compensation for the lack of sources in certain aseismic regions.

99% of seismic records are unexploited 2.1. The notion of information grain

Important notions from information theory can be fruitfully applied to acoustic or seismic wave processing. The idea is to evaluate quantitatively the amount of independent information retained in a seismic record of finite duration and frequency band. Very roughly, the elastic impulse response of an heterogeneous medium like the Earth can be approximated by a series of wave packets arriving with random phase and amplitude [6,7]. Each arrival corresponds to a given path, and each impulse represents one information grain. The duration of the initial impulse t is the inverse of the frequency bandwidth, t = 1/Δ f . Therefore, a record of duration T contains N = TΔ f independent and uncorrelated pieces of information. This means that any acoustic or seismic waveform, including diffuse waves and ambient noise, contains deterministic information on the propagation medium. Except for direct or simply reflected waves, this information is hard to extract. Yet, it is there.

Seismic noise

The last twenty years has undergone a wide development of seismic networks. Large on-line databases gather waveforms from hundreds of worldwide seismic stations. Processing the enormous amount of data is a main challenge and, actually, only a small part of the records are eventually used to image the earth or localize the seismicity. The seismic ambient noise that predominantly fills the databases is very often discarded, as the seismic coda (see Fig. 1.6). Only the direct waves are processed, even though they represent at most 1% of the total amount of data.

Ambient noise is made of continuous vibrations initiated by undetectable tremors (or microtremors) and other low amplitude sources occurring at random. Very roughly, high frequency noise (above 1 Hz) is mainly due to human activity like transportation and industry. Below 1 Hz, the sources include wind, storms, and ocean waves. This noise is due to the coupling between the ocean, the atmosphere, and the ground. At very low frequencies where fundamental modes of the Earth are observed, this ambient noise is referred to as "hum". It should be noted that in all frequencies, the exact location of the noise sources and their generation mechanism remain unknown.

In the high frequency band (above 1 Hz), ambient noise is sometimes used to evaluate the subsoil characteristics, which is the first step to seismic hazard assessment studies [8][9][10][START_REF] Bonnefoy | Nature du bruit de fond sismique : implications pour les études des effets de site[END_REF]. At lower frequency, an important effort is made to locate the noise sources and understand their generation mechanism [START_REF] Suda | [END_REF]13]. But most of the time, the noise is not processed. In some cases, it is not even recorded or just deleted from the database.

The diffuse seismic coda

The seismic coda is named after the Latin word for tail. It is composed of long lasting residual vibrations that follow direct waves (cf. Fig. 1.6). In the short period range (0.1-20 s), this coda lasts sometime several tens of minutes, one order of magnitude longer than the direct travel time from the source to the receiver. It contains most of the seismic energy delivered by an earthquake, but has long been neglected and not processed. In his pioneering works, Aki [14,15] showed that the decay of the coda average intensity is well fit in the 1-10 Hz range by

E ω (t) ∼ 1 t γ exp - ωt Q c (ω) , (1.1)
where 1/t γ is due to geometrical spreading (1 < γ < 2) and ω is the pulsation. Q c is the coda quality factor. It is a dimensionless parameter that characterizes the coda duration, and was found to have a regional variability.

To explain the coda decay, two models were proposed [14,15]. In the first one, coda waves originate from simple scattering in the lithosphere where the exponential coda decay is due to absorption. In the second one, waves are multiply scattered and guided in the crust and the coda decay is due to diffusion. In this last regime, the coda Q factor decomposes into

1 Q c = 1 Q a + 1 Q s with Q a = ωτ a and Q s = ωH 2 D (1.2)
where τ a is the absorption time of the energy, D the diffusion (scattering) constant and H the crust thickness. Only recently, numerical works based on the Radiative Transfer Equation (RTE) demonstrated the validity of the second model [16,17].

In Mexico for instance, Margerin et al. [17] found that, around 1 Hz, Q a 1100, Q s 320, which means that scattering is a stronger effect than absorption. Nevertheless, a direct demonstration of the multiple scattering regime is still needed, and the separation of the effect of absorption Q a and scattering Q s is very tricky to achieve. Again, most of the time, coda waves are not processed. Nevertheless, both seismic ambient noise and diffuse coda waves do contain deterministic information on the medium, and should be used to improve the images of the Earth. To that end, we need to develop new imaging and data processing techniques.

Strategy: from ultrasound to seismology

Seismology is based on the observation and processing of natural vibrations. In a passive field experiment where seismic waves originate from earthquakes, scientists are facing two simultaneous problems. They neither know the source location and mechanism nor the medium of propagation. It is therefore very complex to image the source and the medium at the same time. By reproducing in the lab the main features of the seismic propagation and employing controlled sources and sensors, we can focus our efforts on the physics of the wave propagation and develop more comfortably new methods. In small-scale seismology, we control for instance the size of the medium, the scattering properties, and the absorption.

We are then able to adjust one parameter at a time and test the physical models and imaging techniques we develop.

Another reason for carrying out analogous ultrasonic experiments is more tactical: it is related to the order of magnitude of the physical parameters as recalled in Table 1.1. Ultrasonic wavelengths are on the order of a millimeter, meaning that experiments are physically easy to handle. Additionally, the duration of a single ultrasonic experiment is very short (one minute) compared to seismology where we have to wait for earthquakes (year). This characteristic allows us to achieve many more experiments in the lab. and test many parameters over a wide range of magnitude. Several earth sciences laboratories have already decided to develop analogous experiments for methodological developments, as in San Diego-California [18], Colorado [19][20][21], Tsukuba-Japan [22], Paris-France [23,24] and more recently in Grenoble-France [25] to cite only a few. Our article presents some analogous ultrasonic experiments not only devoted to the study of the physics of wave propagation in heterogeneous media, but also to the development of new techniques applicable to seismic waves in geosciences.

Organization of the present book

The first question addressed in the manuscript is the validity of the multiple scattering approach to seismology.

-Are seismic waves in the multiple scattering regime? -How does one quantify the degree of heterogeneity of the medium?

These points are addressed in Chapter 2 with two experiments: equipartition of seismic waves and weak localization.

The following chapters are devoted to a very rapidly spreading idea: the correlation of ambient noise yields the impulse response between two passive receivers as if a source was placed at one sensor. This gives the opportunity of doing ultrasound or seismology without a source. Chapter 3 proposes different theoretical and experimental approaches to show why the Green's function is retained in the time correlation of diffuse fields. The time-reversal interpretation is presented, followed by a statistical work dealing with the speed of convergence When a wave propagates through a scattering medium, it progressively loses its coherence. The energy of the coherent part of the wave is converted to scattered waves that follow longer scattering paths. In seismology, these long-lasting waves constitute the coda [15] that can be observed in the late part of seismograms. Coda waves are not random noise: they are the deterministic signatures of the heterogeneities of the Earth's crust. Particularly, between 1 and 10 Hz, the coda waveforms show clear evidence of strong scattering [26]. Mesoscopic physicists have introduced one parameter to characterize the importance of scattering: the transport mean free path * . Very roughly, * is the distance after which the wave forgets its initial direction. When the distance of propagation is significantly larger than * , the attenuation due to scattering strongly reduces the direct coherent wave while the order of multiple scattering increases, which complicates the Green's function (GF) and makes conventional imaging very difficult. The cancellation of the direct waves due to scattering has long been considered as a dramatic loss of information. This chapter is devoted to the study of the scattering parameters of the Earth. Knowing whether seismic waves in the Earth are in the multiple scattering regime or not is crucial for application like imaging or monitoring. Additionally, quantifying the scattering (or transport) parameters of a natural medium is an original way to evaluate its degree of heterogeneity and, therefore, to characterize the geologic medium. This section starts with a brief theory on wave scattering, including the definition of the scattering cross-section and the scattering meanfree path. Then we will present different observations (field experiments) where multiple scattering of seismic waves is observed, and a measure of * deduced.

Waves in homogeneous media

Elastic wave equation

Let u(r, t) be the displacement vector from which the strain tensor ε is derived: In the present manuscript, we will use the Einstein's notation convention. We also assume that elasticity is linear: strain and stress are linked by a linear law. This is true for most seismic waves over several order of magnitude. We also assume that elasticity is isotropic, which is not always valid in the Earth but yet will simplify all our equations. The stress tensor follows the Hooke's law (see for instance [START_REF] Aki | Quantitative Seismology[END_REF]):

ε ij (r) = 1 2 (∂ i u j + ∂ j u i ). ( 2
σ ij (u) = λε kk (u)δ ij + 2με ij (u) (2.2)
where δ ij is the Kronecker symbol. λ and μ are usualy called the Lamé parameters, they characterize the rigidity of the medium. We also note ρ the mass density. In a homogeneous medium, λ, μ and ρ are constant. The Newton's law applied to the ith component of the displacement u writes:

∂ tt u - λ + μ ρ ∇(div u) - μ ρ Δ u = 0. (2.
3)

It is noteworthy to define two velocities v 2 p = (λ + 2μ)/ρ and v 2 s = μ/ρ. If we note u = ∇φ + Curlψ = u p + u s , then v p and v s are the velocities associated to the compressional u p and shear u s waves, solutions of the wave equation defined by:

∂ 2 φ ∂t 2 -v 2 P ∇ 2 φ = 0; ∂ 2 ψ ∂t 2 -v 2 S ∇ 2 ψ = 0. (2.4)
These two types of waves are crucial in seismology. Compressional P-waves have a polarization (displacement axis) in the direction of propagation, and shear Swaves have a polarization perpendicular to the direction of propagation. In the Earth, we often observe that λ ≈ μ, therefore v p /v s ≈ √ 3. The decomposition of the wavefield into P and S waves is relevant in infinite media. The Earth has a finite size and a free surface, this feature is at the origin of two other types of waves: Rayleigh and Love waves. These waves are also called surface waves, since they propagate ≈one wavelength beneath the free surface. As the P-and S-wave velocities are different, direct P and S wave reach the seismic stations at different times. The difference in travel time is conventionally used to localize the earthquake.

Green's function in homogeneous medium

The Green function (GF) of a given medium between two points A and B is the waveform sensed at B after an impulse force at A. In a homogeneous medium, and in the frequency domain, the acoustic GF obeys the following equation:

-ω 2 v 2 G 0 -ΔG 0 = δ(x).
The form of the solution varies with the dimension:

G 0 (ω, r) = - e ikr 4πr [3-D], G 0 (ω, r) = -i 4 H 0 (kr) [2-D], (2.5)
where H 0 is the Hankel function of first order (and first kind). In the far field, it asymptotically tends to -e ikr+π/4 /2 √ 2πkr. In the time domain, these GF read

G 0 (t, r) = -1 4πr δ t - r v [3-D], G 0 (t, r) = θ(t -r/v) -2π t 2 -r 2 v 2 [2-D], (2.6)
where v is the wave velocity. θ stands for the Heaviside function. These solutions correspond to a cylindrical (resp. spherical) wavefront diverging from the source point.

Waves in disordered media

Classical and quantum waves: one single equation

In a heterogeneous medium, the elastic parameters (λ, μ, ρ) are spatially fluctuating. It is possible to separate mean quantities from fluctuations. In the case of an acoustic wave for instance, the celerity is rewritten as v(r) = v 0 (1 + ξ(r)) and the wave equation becomes

∇ 2 G(r, ω) + k 2 0 G(r, ω) = ω 2 v 2 0 ξ 1 + ξ G(r, ω), (2.7)
where k 0 = ω/v 0 and v 0 = λ 0 /ρ 0 . When fluctuations are weak enough (compared to the average), the left hand side of the equation describes the free propagation of the wave in a homogeneous medium, and the right hand side describes the coupling with heterogeneities (the Born approximation). This last term generates scattering.

For quantum waves, the density probability of the wave Ψ is a solution of the Schrödinger's equation:

i ∂ t Ψ (r, t) = - 2 2m ∇ 2 Ψ (r, t) + V(r)Ψ (r, t). (2.8)
Classical and quantum waves are similar since we can replace 2mω/ by k 2 0 , and 2mV(r)/ 2 (where m is the particle mass and V the scattering potential) by ω 2 /v 2 0 Ann. Phys. Fr. 31

• N o 3 • 2006 (ξ/(1 + ξ)).
Therefore, theoretical and experimental results obtained in quantum physics or optics should fruitfully transpose to acoustic and seismic waves. This is especially true for waves in heterogeneous media, where several concepts have been transfered from one community to the other. Several examples as the radiative transfer equation, equipartition, and weak localization, will be given in the following.

The physical parameters for the wave propagation

Here follows the definition of several parameters. For the sake of simplicity, these parameters are defined for acoustic (scalar) waves. They are generalized to elastic waves later in the manuscript (see Sect. 2.3, Chap. 4).

Elastic and transport scattering cross-sections

When a plane wave impinges a scatterer, i.e. an area where elastic or acoustic parameters differs from the rest of the medium, its trajectory and shape are modified. The scatterer redistributes the wavefield in several directions with different amplitudes and phases. The strength of the scatterer, i.e. its ability to modify the incident wave, is quantified by its elastic scattering cross-section σ. In 2-D, this cross-section is a distance, in 3-D it's a surface. The greater the section, the stronger the scatterer. There exist different definitions for the cross-section, but all compare the intensity of the incident wave per unit of angle to the intensity of the scattered wave per unit of angle.

If the scatterer is localized in the medium, we can define the scattered wave ϕ(ω) as the superposition of the incident wave ϕ 0 (ω) (the wave without the scatterer at pulsation ω) and the (fictive) wave re-emitted by the scatterer ϕ s (ω). The differential scattering cross-section ∂σ/∂θ is the angular intensity of the reemitted wave in the direction θ normalized by the incident intensity. The direction θ = 0 • is usually the incident direction

∂σ ∂θ (ω) = ϕ 2 s (ω, θ) ϕ 2 0 (ω, θ) • (2.9)
It is relatively easy to numerically estimate the scattering cross-section of a scatterer of any shape and size: we simulate a plane wave impinging the scatterer and sense the wave along a circular array (radius R) of receivers placed in the far-field of the scatterer (R λ). The calculation requires two steps:

1. simulation of the reference wave ϕ 0 (t, θ) received on the array (for all angles θ) without the heterogeneity;

2. simulation of the wave field with the scatterer: ϕ(t, θ).

Ann. Phys. Fr. 31

• N o 3 • 2006
The difference between both records is straightforward and yields the exact scattered field ϕ s (t, θ). In 2-D, the differential cross-section reads

∂σ ∂θ (θ) = R ϕ(t, θ) -ϕ 0 (t, θ) 2 dt ϕ 2 0 (t, θ)dt • (2.10)
If the simulation is performed in a finite frequency band, the obtained differential cross-section is the averaged differential cross-section in the frequency band of the incident wave. This simulation has been performed in an acoustic 2-D medium for three cylindrical cavities of radius a (Fig. 2.1). The numerical scheme is from Tanter [START_REF] Tanter | Application du retournement temporel à l'hyperthermie ultrasonore du cerveau[END_REF]. The background velocity is v = 1.5 mm/μs, the central frequency is f 0 = 1 MHz and the bandwidth is 100%. 200 receivers are placed several wavelengths from the scatterer. The incident wave is generated by a finite-aperture linear array of omnidirectional sources (the effects of the array aperture are visible in the figure). In front of the scatterer, the wave is partially attenuated: the incident wave and scattered wave (in opposite phase) cancel each other. The total elastic cross-section characterizes the total strength of the scatterer, it is defined as

σ = ∂σ ∂θ dθ. (2.11)
Another cross-section is of major interest in this manuscript: the transport crosssection defined as

σ * = 2π 0 ∂σ ∂θ (1 -cos (θ)) dθ.
(2.12)

The transport cross-section quantify the ability of the scatterer to diffuse the wave in all directions. It is evaluated in a way comparable to the elastic cross-section, except we now cancel the contributions in the forward direction (θ = 0 • ). The numerical simulation presented in Figure 2.1 has been repeated for three different cylindrical scatterers of various radii: a = 0.025; 0.25; 2.5 mm. Differential cross-sections are presented in Figure 2.2, and resulting total cross-sections are summarized in Table 2.1. These three simulations were performed in order to verify the asymptotic behavior of the cross-sections:

• ka 1: the scatterer is small compared to the wavelength. The differential scattering cross-section is isotropic (Fig. 2.2a), which implies σ = σ * . In general, when the size tends to zero (a → 0), the wave does not "see" the heterogeneity anymore, and the scattering cross-section tends to zero, too. Nevertheless, this is not exactly true for an empty cavity where the normalized scattering cross-section diverges at low frequency;

• ka 1: the scatterer is bigger than a wavelength. Scattering is mainly forward, and the elastic cross-section tends to the double of the geometric cross-section (σ → 4a) and remains much bigger than the transport crosssection (σ σ * ).

Scattering and transport mean-free paths

The scattering mean-free path is the distance that characterizes the intensity decay of the coherent wave. The coherent wave is obtained after ensemble averaging . , which means averaging over the positions of the scatterers for instance. If we note I the intensity transmitted through a scattering medium of thickness L, then I = I 0 e -L/ where I 0 is the transmitted intensity in the homogeneous medium. In optics, this relation is named after Lambert-Beer: where n is the concentration of scatterers, and A 0 is the normalized attenuation of the scattering solution. The elastic mean-free path is, therefore, the distance that the wave must travel to be significantly attenuated. The density of scatterers n, , and σ are related by

I(L) = I 0 e
= 1 nσ • (2.14)
When the heterogeneity of the medium is not a discrete but a continuous distribution it is not possible to locate a single scatterer. The concept of scattering cross-section of a given scatterer is no longer relevant. Yet it is still possible to statistically define the transport parameters of the medium, since these quantities have a statistical meaning. In that case the scattering cross-section is analogous to the medium's ability to scatter the wave per unit of volume. For a scalar wave, the local celerity is written as v(r) = (1 + ξ(r))v 0 . The medium is characterized by the statistics of the fluctuations ξ(r). These fluctuations have two principal characteristics: their correlation length, and their average (rms) intensity. Let Φ(κ) be the spectrum of the heterogeneities (the spatial Fourier transform of ξ(r)):

Φ(κ) = e iκ(r-r') ξ(r)ξ(r') d(r -r'). (2.15)
We can use this quantity to define the scattering mean-free path :

1 2 = π 2 k 2 0 2 2k 0 0 Φ(κ)κdκ = πk 4 0 4 4π Φ(Ω)dΩ (2.16)
where Ω is a solid angle. The shorter the correlation length, the stronger the scattering and the shorter the mean-free path. Respectively, the greater the fluctuation intensity, the shorter the mean-free path. Another quantity can be defined from the transport cross-section: the transport mean-free path

* : * = 1 nσ * , (2.17)
where * is the distance that the wave must travel to forget its incident direction. This length characterizes the ability of the medium to scatter the wave in different directions. As we will see, * is the characteristic length in the diffusion equation. Therefore, * is also named the diffusion mean free path.

Thouless time

Rigorously speaking, the Thouless time τ D is defined as the inverse of the mean distance between two uncorrelated frequencies of an open scattering medium.

In the time-domain, the Thouless time is the average duration of the intensity I(L, t) [7,[START_REF] Akkermans | Physique mésoscopique des électrons et des photons[END_REF]. Because I is solution of the diffusion equation, the duration of the coda is also the time of its maximum:

∂ t I(L, τ D ) = 0. Then it comes τ D = L 2 4D [2-D]; τ D = L 2 6D [3-D] (2.18) Ann. Phys. Fr. 31 • N o 3 • 2006
where D is the diffusion constant. The number of information grains contained in a diffusive waveform of bandwidth Δ f is:

N = τ D Δ f. (2.19)

Three possible equations and approximations

The wave equation: coherent and ballistic waves

The wave equation has the advantage of requiring few assumptions: the linearity of elasticity, reciprocity, and causality. For a given displacement u(r, t), this equation reads:

ρ(r) ∂ tt u i = ∂ i (λ(r) ∂ k u k ) + ∂ j (μ(r)( ∂ i u j + ∂ j u i )).
(2.20)

The ballistic wave propagates directly from the source to the receiver; its time of flight is the crucial information used in tomography, or other imaging applications. The ballistic wave is often considered as the first arrival, it is the only solution of the wave equation in homogeneous media, and is of natural use in weakly heterogeneous media. However, it is more difficult to detect in multiply scattering media since it has a much weaker amplitude than the following scattered arrivals.

The wave equation has the advantage of being rigorous and perfectly deterministic. Nevertheless, solving the wave equation requires perfectly knowing the medium of propagation in addition to the limit conditions (knowing the field everywhere at a given time for instance); this is very hard to achieve in complex media. In addition, numerically solving the wave equation in a heterogeneous medium is computionaly expensive, and formal solutions are tedious. Therefore, another possible approach is to deal with statistical quantities that directly derive from the wave equation: the ensemble averaged wave named the coherent wave for instance. We emphasize the difference between coherent and ballistic waves. The coherent wave has already been introduced in Section 2.2; it suffers an attenuation that depends on the scattering mean free path and on medium thickness:

ϕ(L, t) = ϕ 0 (L, t)e -L 2 +iδϕ , (2.21)
where ϕ 0 is the field obtained through the equivalent homogeneous medium, and δϕ is a possible phase shift coming from the average velocity of the medium. ϕ(L, t) is then the wave that has propagated in an homogeneous effective medium.

The radiative transfer equation (RTE)

The RTE was empirically established by Chandrasekhar in the early 20th century [START_REF] Chandrasekhar | Radiative Transfer[END_REF]. It describes the specific intensity of the wave I(x, k) at position x and propagating in direction k. Margerin [31] proposes a detailed review on the RTE. The propagation of the specific intensity is assumed to be analogous to the propagation of particles thrown in a scattering medium, the phase information is lost. Each particle at x and propagating in direction k is possibly scattered in the direction k . The RTE describes the mean specific intensity after ensemble averaging:

k • ∇I(x, k) = - I(x, k) + 1 4π 4π I(x, k )p(k , k)dk + e(x, k) (2.22)
where k is a unit vector. The link between the mean specific intensity and the stationary wave equation has been explored in references [START_REF] Rytov | Principles of Statistical Radiophysics I: Elements of random process theory[END_REF][START_REF] Tourin | Diffusion multiple et renversement du temps des ondes ultrasonores[END_REF], but the derivation is tricky and fastidious. Historically, the RTE was introduced empirically, and equation (2.22) derives from the conservation of energy:

• the left hand side is the particle advection;

• -I(x, k)/ is the loss due to scattering into another direction. It is related to the probability of the particle propagating in direction k to be scattered in any other direction k . This scattering probability is eventually related to the elastic mean free path ;

• (1/4π) 4π I(x, k )p(k , k)dk is a term of gain. p(k , k)
is the density of probability that a particle k be scattered in direction k;

• e(x, k) is the source term: the probability to create a particle at x propagating in direction k.

In order to establish the dynamic RTE (the response of the medium after an impulse), one need to empirically add the following term:

1 c ∂I(x, k, t) ∂t + k • ∇I(x, k, t) = - I(x, k, t) + 1 4π 4π I(x, k , t)p(k , k)dk + e(x, k, t).
(2.23) By noting 1 c

∂ ∂ t + k • ∇ = 1 c d dt ,
the RTE is simply a dynamic balance between gain and loss of specific intensity in an elementary volume [START_REF] Tourin | Diffusion multiple et renversement du temps des ondes ultrasonores[END_REF][START_REF] Margerin | Diffusion multiple des ondes élastiques dans la lithosphère[END_REF][START_REF] Lacombe | Propagation des ondes élastiques dans la lithosphère hétérogène : modélisations et applications[END_REF].

The diffusion equation

The diffusion equation derives from the RTE. It no longer deals with the specific intensity in a given direction, but on the energy density ρ:

ρ(x, t) = 1 c 4π I(x, k, t)dk. ( 2.24)
On one hand we have lost the information about the preferential direction of propagation of the particles, but on the other hand, as we will see, the equation for ρ is much simpler. Let's denote the local current density as J(x, t) = 4π I(x, k, t)kdk.

Multiple scattering processes tend to uniformize the angular dependence of intensity because each scattering event distributes energy over the whole solid angle. Hence, after a sufficiently large number of scattering events, we may expect the intensity to differ only slightly from isotropy. The local direction of maximum energy is directed along the current vector. The physical idea of the diffusion approximation is to write the intensity as a sum of two terms: first the angular average, and second a term which takes into account the slight deviation from isotropy expressed in terms of the current vector. In terms of J and ρ the intensity is thus expanded as:

I(x, k, t) = v 4π ρ(x, t) + 3 4π k • J + ...
Let's now apply to the RTE the two following operations. The first one is an integration over all direction of propagation 4π (.)kdk. The RTE becomes:

1 c ∂ t J + v∇ρ = -J/ + 1 4π Jp(θ) cos(θ)dk. (2.25)
The transport mean free path * naturally appears from the last two terms: * = 1 4π 4π p(θ)(1cos(θ))dk In a quasi-static regime, the temporal derivative can be neglected and we obtain the Fick's law:

J = - v * 3 ∇ρ.
(2.27)

We now apply to equation (2.23) the operation 4π (.)dk and obtain:

∂ t ρ + ∇ • J = δ(x)δ(t).
(2.28)

By combining equations (2.27) and (2.28), we obtain the well-known diffusion equation:

∂ t ρ -D∇ 2 ρ = δ(x)δ(t); D = v * 3 • (2.29)

Mesoscopics of seismic waves

To deal with waves in heterogeneous media, we have so far introduced four characteristic lengths: the absorption length l a , the size of the medium L, the transport mean free path * , and the wavelength λ. The comparison of these length allow to define several regimes. Two regimes are probably the most interesting.

In the ballistic regime, the wave propagates with little absorption and scattering:

λ < L < l a and * . (2.30)
In this regime, direct waves are well resolved and the simple scattering approximation applies. It is therefore the ideal regime for applications like imaging, monitoring, or communication. Another regime is of high interest: the mesoscopic regime [START_REF] Akkermans | Physique mésoscopique des électrons et des photons[END_REF][START_REF] Sheng | Introduction to wave scattering, localization and mesoscopique phenomena[END_REF][START_REF] Trégourès | Approche Mésocopique des Ondes en Milieu Complexe : des Micro-Ondes aux Ondes Sismiques[END_REF], that obeys the following relations:

λ < * L l a . (2.31)
In this regime, waves may be multiply scattered before reaching the receiver or being absorbed. Wave trajectories are very complex and well described by random walks. Nevertheless, as long as the medium does not move, the wave equation is still fully deterministic and phase effects like wave interferences can occur. In this regime, both wave (microscopic) and Boltzmann (macroscopic) equations are valid. This regime has undergone intensive investigations during the last thirty years, especially in optics and quantum mechanics (electronic conductance). New applications were found in acoustics, and, as we will see in this section, mesoscopics also applies to seismic waves.

Equipartition of seismic waves

Theory

The multiple scattering assumption has an important consequence: at each scattering event, waves are converted. Therefore, after several scattering, a stationary Ann. Phys. Fr. 31

• N o 3 • 2006
regime is set and mode conversion equilibrates: the initial energy spreads with equal probability over all modes. This equipartitioned regime was fully investigated by Weaver in the 1980's [START_REF] Weaver | [END_REF][39][40]: he calculated the energy partitioning of ultrasonic waves in an elastic body with and without a free surface. In an infinite 3-D medium, the energy density of bulk P and S waves is ρ p,s = ω 2 /v 3 p,s (this formula is simply derived from modal density). The equipartition ratio is the ratio between potential energy:

E s E p = 2v 3 p v 3 s • (2.32)
In the Earth crust (where v P ≈ 6 km and v S ≈ 3.5 km) the expected ratio is 10.4 [41,42]. It is important to emphasize that this ratio is reached only in the multiple scattering regime. In the simple scattering regime, this ratio mainly depends on the source mechanism. For instance, an isotropic explosion (divergent excitation) emits mainly P waves, and the ratio tends to 0. On the contrary, just after an earthquake (a double couple force), this ratio tends to ∞.

Observation of equipartition

An important effort was made by Shapiro and co-workers to experimentally measure this ratio. Their aim was to measure the energy of compressional and shear potential energy:

E P = λ 2 + μ (div u) 2 ; E S = μ 2 (Curl u) 2 .
(2.33)

To that end, in 1998, they dispatched in Mexico an array of seismic stations in a geometry that allows the measurement of the spatial derivatives of the seismic wavefield u (see Fig. 2.4). Spatial derivatives are well-evaluated if the distance between the three-components stations d is shorter than the wavelength λ d:

∂u i ∂x = u 3 i -u 1 i d , ∂u i ∂y = u 2 i -u 1 i d , i = x, y, z.
(2.34) The vertical spatial derivative is easily evaluated by means of the traction free condition at the surface:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∂u y ∂z = - ∂u z ∂y ∂u x ∂z = - ∂u z ∂x ∂u z ∂z = - v 2 S v 2 p ∂u x ∂x + ∂u y ∂y • (2.35)
They observed that, in the coda, the energy between P and S wave equilibrates, and the ratio fluctuations are weak, as predicted by theory. Nevertheless, they found an average ratio of 7.30 ± 0.72, which is not in agreement with the 3-D elastic theory. Actually, theory was initially developed for elastic waves without the free surface. The surface waves drastically modify the partition ratio, and their omission was not justified. Therefore, at the free surface, all modes (bulk and surface waves) have to be taken into account in the equipartition ratio. The modified theoretical equipartition ratio at the Earth surface was later found to be 7.19 [26]. Theory and experiment then (almost) matched. The observation of the stabilization of the energy partition ratio, and the absolute value of the equipartition, are two strong indications that coda waves, at least in Mexico (Fig. 2.5), are made of multiply scattered waves.

A new model for the Earth crust

From these equipartition consideration, along with the work of Margerin on the radiative transfer equation, a new model of the Earth crust was proposed [17,26,[START_REF] Margerin | Diffusion multiple des ondes élastiques dans la lithosphère[END_REF]43]. The continental crust, whose thickness varies from 30 to 70 km, is a very heterogeneous slab overlying a rather homogeneous mantle. The impedance contrast at the Moho (the crust-mantle discontinuity) induces waveguiding for Schematic view of the earth crust as a possible model for seismic waves propagation as proposed by Margerin et al. [44], and confirmed by the observation of equipartition.

large angles of incidence and partially traps the seismic energy in the crust. For other incidences, the waves leak into the mantle and are lost. Scattering properties of the crust around 1 Hz vary over several orders of magnitude from one place to another. Nevertheless, in average, one can assume that the mean-free path * is about the crust thickness. This picture is presented in Figure 2.6.

Weak localization of seismic waves

The seismic coda is not always processed because it is believed not to contain any structural information that is easily extractable using standard imaging techniques. During the past two decades, radiative transfer was successfully introduced to model the energy decay of coda waves [START_REF] Sato | Seismic Wave Propagation and Scattering in the Heterogeneous Earth[END_REF]. It describes the transport of the wave energy in space and time, but does not take into account phase information. Is the phase information preserved in the multiple scattering regime? And, moreover, can we demonstrate the presence of multiple scattering in the seismic coda? These points are addressed in this section where we observe the Weak localization (WL) effect.

Weak localization: an interference effect

Weak localization is a manifestation of interference of multiply scattered waves in disordered media and was first discovered 20 years ago in quantum physics. It was recognized to be the origin of novel features in the electronic magnetoresistance at low temperatures [START_REF] Yu | [END_REF][47][START_REF] Kawabata | Anderson Localization[END_REF], and a genuine explosion of mesoscopic physics followed. The discovery of WL constituted the desired counterexample of the intuition that multiple scattering of waves destroys wave phenomena, reducing it conveniently to classical radiative transfer, where waves are treated similarly to hard spheres colliding with obstacles.

Weak localization finds its origin in the constructive interference between long reciprocal paths in wave scattering [START_REF] Al'tshuler | [END_REF][50][START_REF] Van Tiggelen | Wave Propagation in Complex Media[END_REF]. This enhances the probability to return to the source by a factor of exactly 2, which results in the local energy density enhancement by the same factor. In optics [START_REF] Tsang | [END_REF][53][54] and in acoustics [55,56] is better known as coherent backscattering, and it was shown to be an accurate way to measure transport mean-free paths * or diffusion constants. However, the WL effect has never been observed in seismology. Such an observation is presented in the following. The aim of this work is to show the relevance of mesoscopic physics to seismology and its necessity to interpret seismic observations.

In seismic experiments, we expect WL to appear as an enhancement of seismic energy in the vicinity of a source [57,58]. The main features are:

1. the energy enhancement factor around the source is two, 2. the spatial extension of this enhancement is the Rayleigh wavelength, 3. and the characteristic time for its onset is the scattering mean-free time.

Figure 2.7 proposes an empirical interpretation of WL, or coherent backscattering effect. In the simple scattering regime, waves randomly interfere with each other. For two waves of unity amplitude, the intensity of their interference is 2 in average, which results in a constant energy background backscattered around the source. In the multiple scattering regime, this picture is still valid for distant source and receiver. But if source and receiver are close together, we have to take into account coherent interference of reciprocal paths. Each trajectory SR made of unit amplitude wave adds up constructively with its counterpart. The resulting intensity is 4; twice as much energy as for other backscattered waves. This constructive interference occurs over one wavelength. The detailed calculation of the spatial enhancement is given in next section.

Theory for the spatial enhancement

Since seismic waves sensed at the surface are mostly made of 2-D surface waves, we propose here a 2-D scalar formulation for the spatial energy enhancement. In addition, the 2-D scalar calculation is much simpler than the full 3-D elastic case. Of course, we have to keep in mind that this assumption is not rigorous: at least part of the scattering occurs in the bulk.

2-D scattering open medium

Let an acoustic wave of pulsation ω propagate in a 2-D scattering medium of average velocity v. It is emitted by S and recorded at R. Let D 1 denote the first scatterer. Statistically, the wave propagates over a distance (the scattering mean free path) before the first scattering event. Similarly, D n denotes the last scattering event, and is at a distance from the sensor R.

The propagation from S to D 1 is equivalent to the free propagation in the effective medium where the wavefront suffers from an attenuation equivalent to e -|S-D 1 |/2 . The wave re-emitted by the last scatterer undergoes an equivalent attenuation. The corresponding GF reads:

G(S, D 1 ) = -i 4
H (1) 0 (k • (S -D 1 )) e -|S-D 1 |/2 , (2.36)
where k = ω/v is the wavenumber.

The field sensed at R is the sum of all ray paths in the 2-D scattering medium. This sum can be expanded into a double sum D 1 ,D n G(S, D 1 , D n , R) where D 1 and D n are the first and last scattering event. One has to keep in mind to take into account reciprocal paths

D 1 ,D n G(S, D n , D 1 , R).
It is not possible of course to describe in details all the ray trajectories between the first and last scatterer. Fortunately, we do not need to know the exact GF between D 1 and D n . Indeed, a description in terms of mean intensity is enough: the phase information of the exact waveform of that part of propagation is lost when processing the intensity at R. This propagator, noted P(D 1 , D n , t), is the probability density for the wave at D 1 to get at D n after a lapse time t. The square of this propagator is for example the solution of the RTE, or more simply solution of the diffusion equation ∇ 2 P 2 -D∂ t P 2 = δ(t)δ(r): The field at R writes:

P 2 (D 1 , D 2 , t) = 1 4πDt e -|Dn-D 1 | 4Dt 2 . ( 2
G(S, R, t) = G(S, D 1 )P(D 1 , D n , t)G(D n , R) + G(S, D n )P(D n , D 1 , t)G(D 1 , R) dD 1 dD n . (2.38)

Coherent and incoherent intensity

The total intensity I(t) received at R at time t is:

I(t) = G(S, D 1 , D n , R) + G(S, D 1 , D n , R)dD 1 dD n × G * (S, D 1 , D n , R) + G * (S, D 1 , D n , R)dD 1 dD n (2.39)
where * stands for complex conjugate. After an expansion and a simplification of equal terms, this formula splits into two contributions: one is called the coherent intensity I coh , the other is the incoherent intensity I inc :

I inc = G(S, D 1 )G * (S, D 1 )P 2 (D 1 , D n , t)G(D n , R)G * (D n , R)dD 1 dD n , I coh = G(S, D 1 )G * (S, D n )P 2 (D 1 , D n , t)G(D n , R)G * (D 1 , R)dD 1 dD n .
Though I coh and I inc have similar form, they have totally different physical interpretation. In the incoherent part, there is no interference between two different paths, the incoherent part is spatially constant. In the coherent part interference between reciprocal paths occurs: D 1 → R interferes with D n → R. In average, this interference cancels the coherent intensity. But when the receiver is close to the source Δr = |S -R| < λ, then the two couples D 1 R -D n R and D 1 S -D n S exactly match, which results in the doubling of the intensity. The spatial extension of the enhancement around the source S(Δr) is given by:

S(Δr) = 1 + I coh I inc = 1 + |J| 2 (2.40)
where

J = G(D 1 , S)G * (D 1 , R)dD 1 . (2.41)
To evaluate I coh and I inc , we substitute G and P by equations (2.36) and (2.37), then we carry out the double integration over D 1 and D n . Note that the two GF G(D 1 , S) and G(D 1 , R) can be expressed in terms of their Fourier transform:

G(D 1 , R) = 1 4π 2 e iκ•(D 1 -R) k 2 -κ 2 + ik/ d 2 κ, (2.42
) The integration over D 1 leaves:

G(D 1 , S) = 1 4π 2 e -iκ •(D 1 -S) k 2 -κ 2 + ik/ d 2 κ . ( 2
J = δ(κ -κ )e -iκ •(S-R) (k 2 -κ 2 + ik/ )(k 2 -κ 2 + ik/ ) d 2 κd 2 κ . (2.44)
The integration over κ yields:

J = 1 4π 2 J 0 (κΔr)κdκ (k 2 -κ 2 ) + k 2 / 2 • (2.45)
By using the equality J 0 (κΔr) = H

(1) 0 (κΔr)+H

(2) 0 (κΔr), and noticing that H

(2) 0 (κΔr) = H

(1) 0 (-κΔr + iε) with ε → 0, and applying the Cauchy theorem we find:

J = π 2k H (1) 0 kΔr + iΔr 2 -H (1) 0 -kΔr + iΔr 2 • (2.46)
To summarize, the spatial extension for weak localization, in a 2-D scalar case and in the asymptotic ( Δr) regime, takes the form

S(kΔr) = 1 + J 2 0 (kΔr), (2.47)
which confirms that the backscattered energy is doubled around the source with a spatial extension characterized by the wavelength. Note that because some 3-D shear and compressional waves are also in the seismic record, and because their respective wavelength are slightly larger (λ P , λ S > λ R ), we expect the WL width to be slightly larger than the one calculated in the Rayleigh 2-D approximation.

Experimental setup

The observation of WL requires four conditions to be fulfilled:

1. some receivers must be placed less than one wavelength from the source (interference condition); We measured the vertical ground motion using a linear array of 23 geophones separated by 2.5 m. The ground motion is the result of a sledgehammer strike at time t = 0 on a 20 cm × 20 cm aluminum plate which was repeated 50 times for each location.

The individual strikes produce forces that may fluctuate slightly in direction around the ideal vertical direction, and may have slightly different frequency contents. By stacking the records of 50 repeated strikes, we simulate a vertical point force which can be considered as a narrow impulse in the 15-30 Hz frequency range. Because the receivers are placed at the free surface, the detected waves are both bulk waves (with either compressional or transverse polarization) and surface waves (Rayleigh waves with elliptical polarization), each propagating at its own velocity. The wavelengths λ roughly range from 9 m (30 Hz Rayleigh waves) to 40 m (15 Hz compressional waves). A typical record is presented in Figure 2.10.

The first 0.5 s of the 3 s signal is composed of direct and simply reflected waves, which are traditionally used in seismic prospecting. In this work we will process the average energy of the subsequent seismic coda. For a given source, the intensity in the coda shows spatial and temporal fluctuations that are similar to the optical speckle. Because of its random nature, the speckle is suppressed by a configurational average while the deterministic WL effect survives. The only average conceivable in seismology is one over source and receiver positions for a fixed source-receiver distance Δr. To this end, we kept the receiver array fixed and we placed the source next to a receiver every 5 m along the array, which provides a total of 12 different configurations. For a diffuse field, the correlation length is λ/2 [61], which implies that each of these measurements corresponds to an independent source-receiver configuration for wavelengths of 10 m or less. For larger wavelengths, statistical correlations still persist, which may degrade the ensemble-averaging process.

Experimental concerns

Studying actively generated coda waves is an experimental challenge. Because of the very fast coda decay, seismic signals must be recorded over a huge dynamic range without reaching the ambient seismic noise level, or even the electronic noise. For this reason, recent 24 bit seismic acquisition devices are crucial. Yet, it's not enough. To increase the signal-to-noise ratio in the coda, three additional steps are required:

1. the seismic signal must be high in the frequency range of interest (here: 10-100 Hz). We increased the earth-geophone coupling by burying the sensors at a depth of 20 cm in the ground. In addition, we repeated the source 50 times for each position of the source. Of course, an explosive source would have been more convenient, but it does not meet the reciprocity condition for WL (condition number 2);

2. the ambient noise must be weak. We carried out experiments during quiet times: by night and under anticyclonic conditions, and in a quiet place: away from roads, towns, and industrial areas. We noted that even a weak wind in the nearby trees might be a noticeable source of seismic noise. Additionally, burying the geophones reduces the effect of aerial ambient noise;

3. the coda must be long lasting (high Q factor). To that end, we chose a highly heterogeneous geological site: a volcano. We also tested several other places, like limestone sediments or moraines. In the first case, scattering was not very high, in the latter absorption was more important. Nevertheless, we stress that larger experiments followed by meticulous analysis could also reveal WL at these sites.

The identification of WL must be accompanied by a close study of different kinds of noise that contaminate the seismic record. In the following discussion, we separate the ambient noise from that generated by the operator of the hammer, and identify the mesoscopic regime where noise is negligible. The operator noise comes from the person manipulating the hammer who is subject to residual movements just before and after the hammer strike. This noise is difficult to separate unambiguously from the signal, because it is local and non-stationary and could be misinterpreted as WL. Ambient background noise is stationary and random, but the operator noise might not be... Fortunately, biophysical studies have revealed that the reproducibility of human motion is limited to frequencies lower than 10 Hz [62]. This suggests that the noise produced by the operator can be considered as random in our frequency band. In order to verify this property and to test quantitatively how efficient the stacking is in increasing the signal-to-noise ratio, we study the sum of M signals s i (t) produced by repeated strikes at the same location. Each signal results from N = 10 strikes that were automatically stacked in the field. We expect both the ambient and the human noise to add up incoherently (∝ √ M) while the seismic signal deterministically generated by the impacts should add up coherently (∝ M). Figure 2.11 confirms these assumptions. For a more quantitative study, we analyze the time evolution of the signal-to-noise ratio using the cumulative index R:

R(t) = 1 M M i=1 s i (t) 2 M i=1 s 2 i (t) • (2.48)
The brackets . denote an average over one oscillation period T = 40 ms. The ratio R(t) takes its maximum value 1 for a perfectly deterministic signal and equals 1/ √ M for pure random noise. confirms the randomness of the operator noise (t < 0) and the deterministic nature of the seismic signal. Between 0 and 2 s, R(t) always exceeds 90% which enables the processing of the coda with an excellent signal-to-noise ratio.

Results

To evaluate the spatial enhancement of energy S(Δr), we normalize the average energy E C around the source by its measured average value E D sufficiently far away (15 m) from the source where the energy density is independent of the source-receiver distance Δr. The theoretical prediction for S(Δr) at the free surface of an elastic body was obtained in [59] and Section 4.2. The vertical force generates both bulk and Rayleigh waves, which undergo mode conversions resulting in equipartition in phase space [26,[START_REF] Weaver | [END_REF]. While both waves play a part in the dynamics of scattering, the Rayleigh waves dominate the local energy at the free surface once equipartition is established. As a result, the rigorous expression obtained in reference [59] can be approximated by the profile predicted for 2-D random media [63,64]:

S(Δr) ≡ E C E D ≈ 1 + J 2 0 (2πΔr/λ) (2.49) Ann. Phys. Fr. 31 • N o 3 • 2006
where λ is the wavelength of the predominant Rayleigh waves and J 0 is the Bessel function. Note that for the near-field regime the size of the WL spot is independent of elapsed time t, contrary to the far field regime [55,56]. The energy distribution E(t) at each sensor is integrated over one sliding window of one cycle duration. The dynamics are studied by analyzing the signals in non-overlapping time windows of 0.4 s duration. In each window, E(t) is normalized at each time t by the maximum over the array, and then averaged over the 12 configurations with equal Δr. This procedure compensates for the exponential decay of the total energy, and provides an unbiased average over the different strikes. Finally, we integrate the normalized, averaged energy E(t, Δr) over the entire time window. S is then computed from equation (2.49).

In Figure 2.12, we plot the seismic energy around 20 Hz measured in the coda as a function of source-sensor distance, and for three specific 0.4 s windows. Around 0.3 s only simply reflected waves are recorded and no energy enhancement is visible around the source. The remaining fluctuations are ascribed to the incomplete suppression of speckle. As from 0.7 s, WL is observed with a gradually increasing enhancement factor at the source. After 1.7 s, the profile including the enhancement factor 2 has stabilized, as predicted by the theory for WL in the near field. Therefore, we attribute this enhancement to WL. According to equation (2.49), the spot has a spatial extent equal to the wavelength λ. This gives the estimate v = 260 m/s for the phase velocity of the Rayleigh waves around 20 Hz. Since at least two scattering events are necessary to generate the enhancement effect, the rise of the enhancement factor corresponds to the transition from the simple to the multiple scattering regime. It was verified in numerical studies [57] that the characteristic time governing the rise of the enhancement factor is the scattering mean-free time t * . We thus conclude that this important time scale is of the order of 0.7 s around 20 Hz. For a velocity v = 300 m/s, this implies a scattering mean-free path * ≈ 200 m. We emphasize that this parameter is very difficult to measure with traditional techniques based on attenuation studies because absorption is hard to separate from scattering effects. Finally, we have studied the frequency dependence of WL. To this end, the seismograms were filtered in three consecutive frequency bands, and the energy profiles were computed as above, though now averaged over the entire coda that exhibits the stabilization of the enhancement S(Δr) (Fig. 2.13). Three different WL widths are observed. The values for the wavelengths estimated from a fit to equation (2.49) have been indicated. We have separately measured the wavelength of Rayleigh waves from a dispersion analysis of direct arrivals in the original records. Both estimates of the wavelength are consistent and indicate a significant dispersion due to the depth dependence of elastic properties. As a result, the spatial width of WL depends nontrivially upon frequency. Future studies might even reveal the frequency dependence of the scattering mean free path * , which would provide precious information on the nature of the heterogeneity.

Conclusions and implications for imaging

We have here presented two observations that demonstrate the presence of multiple scattering in the seismic coda. One is the seismic equipartition, which result in an equilibration between the shear and compressional potential energy in the coda. The other is weak localization of seismic waves in a shallow volcanic structure. WL was observed both in space and time, and was found to be in good agreement with the near-field theory which predicts a size of one wavelength for the enhancement spot.

The study of the WL dynamics turns out to offer a unique opportunity to measure the scattering mean-free time without the bias of absorption. We found an estimate of t * ≈ 0.7 s for the mean-free time for seismic waves around 20 Hz. This time characterizes the transition from the single scattering regime to the multiple scattering regime, which has an important implication for standard seismic imaging.

• In a reflection seismic survey, t * is the maximum time to be processed. Since waves arriving later than t * contain more multiply than simply scattered waves, their processing will not enhance the quality of the seismic images, and may even degrade them.

• In a refraction seismic survey, we suggest that direct waves are not taken into account if they arrive later then 10t * . Indeed, scattering along the ray path results not only in attenuation, but also in fluctuations in arrival time.

Again, the seismic image could be degraded.

All the geological media are potentially multiple scattering. The question of whether or not they are multiply scattering now appears irrelevant. Their degree of heterogeneity is connected to their ability to scatter seismic waves, which can be quantified with the scattering mean free time t * . Only the comparison of t * to the intrinsic absorption of the medium t a , and to the record length T, allow conclusions to be made on the validity of the simple scattering approximation. Observing multiply scattered seismic waves does not only depend on the heterogeneity of the medium, but also on our ability of recording low-amplitude late arrivals.

The scattering mean free time can vary over several orders of magnitude from one site to another (see Tab. 2.2). This clearly opens up new routes for evaluating and characterizing natural media like the Earth. To conclude more generally, our experiments reveal the mesoscopic nature of seismic waves that have traveled many kilometers in the earth crust. As has been the case in nanophysics and in colloid physics, mesoscopic physics may open up new fields of investigation and application in seismology. Precisely, one potential application of mesoscopics of seismic waves is presented in next three chapters.

Correlation of diffuse waves and noise

The physics of wave propagation in complex media covers various areas of research, ranging from quantum mechanics to classical waves like optics [START_REF] Sheng | Introduction to wave scattering, localization and mesoscopique phenomena[END_REF]53,54,[69][70][71][72], ultrasound [56,73], acoustics and room reverberation, seismology [15,26], astrophysics, or ocean acoustics [74,75]. This diversity gave rise to fruitful transdisciplinary approaches. From the 1980's on, huge improvements have been achieved in understanding and modeling wave propagation in inhomogeneous, random, or reverberant media. Many domains of applications are concerned, e.g., imaging [START_REF] Fink | Imaging of Complex Media with Acoustic and Seismic Waves[END_REF], detection, or communication [START_REF] Derode | [END_REF] in a complex environment. Usually, the first, essential step is to know the Green's function of the medium under investigation. When possible, the Green's function (or impulse response) h AB between two points A and B is determined by a direct transmit/receive measurement.

The purpose of this chapter is to develop another way to recover the GF of a medium. We will take advantage of the phase information contained in the seismic coda or in the seismic noise. We propose a passive technique that thoroughly exploits the seismic records to produce an image of the medium. This will be done by correlating records obtained at two different points. Within the last 2 decades, strong attention was paid to correlation of scattered waves. In optics, short and long range intensity correlations of speckle patterns were thoroughly investigated [START_REF] Sheng | Introduction to wave scattering, localization and mesoscopique phenomena[END_REF]69,70]. Time-varying correlation of scattered fields [73,78] (diffusive wave spectroscopy) have given new insight for monitoring changes in complex media. What we propose here is different. When averaged over time and/or sources, this correlation yields the impulse response between the two receivers, as if a source was placed at one of them. This is the basis of the very rapidly spreading Passive Imaging technique. Different theoretical approaches are presented in the following. Because a real seismic experiment is time consuming and never easy to setup, we will illustrate this principle with numerical simulations and ultrasonic experiments.

The pioneers

In his pioneering works in 1957, Aki [8] noticed that the seismic wavefield originating from different noise sources is spatially coherent over one wavelength.

More precisely, if the seismic surface waves are incoming from all directions, meaning that the noise is diffuse enough, the spatial correlation takes the form of the Bessel function ϕ(x, ω)ϕ(x + r, ω) * = |F(ω)| 2 J 0 (kr). He suggested that the observation of such a spatial correlation could provide an accurate way of measuring the wavelength of surface waves at several frequencies, and, therefore, through an inversion of the dispersion curve, could give the structure of the subsurface at the given site x. This Spectral Auto-Correlation method (SPAC) assumes a perfectly diffuse seismic noise made of 2-D surface waves. It is now widely used to constrain the elastic properties of the first km beneath the Earth's surface, especially for seismic hazard assessment [START_REF] Bonnefoy | Nature du bruit de fond sismique : implications pour les études des effets de site[END_REF].

Later on, Claerbout stressed that this observation could be extended for distances r much larger than the wavelength. He conjectured that the spatial correlation of seismic waves, averaged over time, should yield the actual impulse response of the medium. This property was demonstrated theoretically in a 1-D layered medium [79], but was not at that time generalized to 3-D media and was found to have only feeble experimental confirmations.

During the early 1990's, Duvall and co-workers developed a similar technique to image the sun's interior [74]. From the correlation of surface velocity records obtained by Doppler measurements, they computed the acoustic impulse response of the sun between two points at the surface and, from the arrival time and the ray theory, imaged the sun's interior. More recently, Weaver and Lobkis demonstrated that the time-averaged correlation of diffuse ultrasounds or ultrasonic noise yields the elastic Green's function of the medium. Their approach was quite original: it took advantage of the full waveforms and, more precisely, of the mesoscopic nature of diffuse fields. If the wave is sensed by two passive receivers at A an B, they showed both theoretically and experimentally that the time derivative of the correlation is the exact Green's function of the medium, as if a source was place at A or B:

∂ τ C AB (τ) = ∂ τ ϕ A (t)ϕ B (t + τ)dt = G + (A, B, τ) -G -(A, B, -τ) (3.1)
where G + and G -are the causal and anti-causal Green's function. This result provided the opportunity of doing ultrasonics without a source [80], and now finds considerable application in seismology.

Theoretical approach 2.1. The fluctuation-dissipation theorem

The idea of retrieving some deterministic information on the medium from the correlation of a fluctuating variable is not recent; it is actually the underlying idea of the Fluctuation-Dissipation Theorem (FDT), developed by Einstein for the Brownian motion of particles in a viscous fluid [81]. The FDT was generalized by Kubo [82], who formally linked the fluctuations of any given field ξ(r, t) to the linear response G of the system (dissipation for instance), provided ergodicity.

The system is ergodic as long as time and ensemble averaging are equivalent. This last assumption sometime deserves attention. More recently, Rytov et al. [START_REF] Rytov | Principles of Statistical Radiophysics III: Elements of Random Fields[END_REF] proposed to extent the FDT to elastic and electromagnetic waves under thermal equilibrium. Their demonstration is essentially based on the assumption that sources are continuously distributed in the whole scrutined volume (thermal equipartition). The central part of the FDT is the averaging, which can be realized over time, over space, over disorder, or over sources.

Each averaging might yield a different result. In the following, we will perform source and/or time averaging. This averaging is probably adapted to practical applications like imaging, but it is not always the easiest way to perform analytical calculation. As we will see, we will sometimes need to invoke ensemble averaging, especially to understand the role of scattering and the dynamics of correlations in the next chapter [20,[START_REF] Van Tiggelen | [END_REF].

Correlation in a reverberant elastic body

We propose in the following a theoretical approach based on a modal expansion valid in closed reverberant bodies, which perfectly applies to the aluminum cavities employed by Weaver and Lobkis in their experiments [80,85,86]. We will point out the role and importance of the Heisenberg time (break time) and address the question of the possible time-derivative of the correlations. Note that the dimension of each quantity is presented in square brackets.

First, we recall the link between the elastic wave equation and the local displacement u due to a volumic force f 0 :

-ρ(x)∂ tt u i + ∂ i (λ(x) ∂ k u k ) + ∂ j (μ(x)( ∂ i u j + ∂ j u i )) = f i (x, t) (3.2)
where λ and μ are the Lamé coefficients, G [m s -1 N -1 ] is the Green function of the system: the impulse response of the medium, which allows the evaluation of any displacement u(x, t) [m] resulting from any force f(y, t ) [N m -3 ]:

u i (x, t) = V t -∞ G ij (x, y, t -t ) f j (y, t )d 3 ydt = G ⊗ f (3.3)
where ⊗ represents convolution, and V is the volume of the body. We introduce the spatial derivative operator L x defined as:

L x u i = ∂ i (λ(x) ∂ k u k )+∂ j (μ(x)(∂ i u j +∂ j u i )).
In the frequency domain, the wave equation reads:

L x + ρ(x)ω 2 G ω (x, y) = f 0 δ(x -y) + Traction free condition at the surface. (3.4) Let φ n (x) [kg -1/2
] be the eigenmodes of the cavity, associated with the eigenvalues

-ρ(x)ω 2 n : L x φ n (x) = -ρ(x)ω 2 n φ n (x).
Ann. Phys. Fr. 31

• N o 3 • 2006
The φ n form a complete set of orthogonal normalized functions:

ϑ ρ(r)φ n (r)φ m (r)d 3 r = δ nm .
This relation is derived from the energy conservation. The GF then reads:

G ij+ ω (x, y) = n φ i n (x) φ j n (y) ω 2 -ω 2 n + iε [s 2 kg -1 ] ε → 0 + . (3.5)
Applying the inverse Fourier transform gives the causal (exponent +) GF:

G ij+ (x, y, t) = θ(t) n φ i n (x)φ j n (y) sin(ω n t) ω n [m N -1 s -1 ], (3.6)
where θ(t) is the Heaviside function. The correlation of the wave sensed at A and B, averaged over a time window ΔT (the record length for instance), and for components j and k, reads:

C ΔT AB, jk (τ) = 1 ΔT ΔT 0 G ij (S, A, t)G ik (S, B, t + τ)dt. (3.7)
Using the modal expansion, and after integration over time, the correlation reads:

C ΔT AB, jk (τ) = n,m G ij (S, A, ω n )G ik * (S, B, ω m )2sinc[(ω n -ω m )ΔT]e -iω m τ e i ωn-ωm 2 ΔT . (3.8)
The sinc function is the central point of the correlation. A characteristic time now shows up: the Heisenberg time T H , defined as the inverse of the average distance between two eigenfrequencies:

T H = 2π ω n -ω m • (3.9)
The statistical distribution of eigenfrequencies in a chaotic reverberant body is a Wigner distribution that follows from the random matrix theory. It results in mode repulsion and a characteristic distance between two consecutive eigenfrequencies.

In a large elastic body, the Heisenberg time depends on the compressional and shear wave velocities v P and v S , the size L of the body, and the central pulsation ω: Comparing T H to the integration time ΔT gives the asymptotic behavior of the correlations. If the integration (or record) time ΔT is much greater than the Heisenberg time ΔT T H , then the sinc function resembles a Dirac function, and the correlation simplifies into:

T H = 3π 3 L 3 ω 2 ⎛ ⎜ ⎜ ⎜ ⎜ ⎝ 1 v 3 P + 2 v 3 S ⎞ ⎟ ⎟ ⎟ ⎟ ⎠ • For instance,
C yz, jk (τ) = n φ i n (S)φ i n (S)φ j n (A)φ k n (B) e iω n τ ω 2 n ,
where one should recognize G AB (τ):

∂ τ C AB (τ) = G(A, B, τ) × G(S, S, τ). (3.10)
This equation was named the cavity equation by several authors [23,87,[START_REF] De Rosny | Milieux réverbérants et réversibilité[END_REF]. As G(S, S, τ) is essentially a peak at time τ = 0, the correlation is very close to the GF between the two receivers. If we need to obtain the GF more rigorously and with one single source S, then a tricky deconvolution is required. Another idea is to average the correlations over a few sources to cancel late time fluctuations of G(S, S, τ) and only retain the peak at zero. This later operation is probably the most efficient, and yet elegant. We would like to address now the following question: do we always need to take the time derivatives of the correlations to estimate the GF? The time derivation was emphasized by several authors [START_REF] Roux | [END_REF][90][91], while neglected by others without affecting the quality of the GF reconstruction [23,[92][93][94]. Actually, a dimensional analysis is enlightening: one just has to refer to the dimensions of the modal decomposition presented earlier in this section. If the receivers are recording the displacement, and if we want to reconstruct the response of the medium to an impulsive force, then the time derivation is required. Nevertheless, if we measure other quantities, like the velocity or the acceleration, but still want to reconstruct the force impulse response, then no derivation is required. An integration is even needed in the latter case. To conclude, the issue of the derivation of the correlations depends of the dimensions of the correlated fields and of the expected GF. In some cases, where the phase of the wave field is studied, the derivation step will be needed.

Thermal equipartition

In most experiments, it is very hard to obtain records whose lengths are comparable to or greater than the Heisenberg time (except in isolated low-absorbing bodies). Weaver and Lobkis proposed another idea: average simultaneously over time and noise sources. They showed that the thermal fluctuations of the displacement, though very weak, are recordable by conventional piezoelectric transducers and, when correlated, yield the GF of the medium. The idea is that a thermal phonon has a very long life time (coherence time). It can therefore reverberate several times and hit both receivers before disappearing. In other words, they proposed to take advantage of the mesoscopic nature of thermal ultrasonic phonons.

The responses recorded in A and B are:

u(A, t) = n A n φ n (A) cos(ω n t) (3.11) u(B, t) = m A m φ m (B) cos(ω m t).
(3.12)

The time-averaged correlation then reads:

C AB (τ) = n,m A n φ n (A)A * m φ m (B)e i(ω n -ω m )t-iω m τ thermal .
The thermal excitation is random in phase and amplitude, so that the A n 's [m kg 1 2 ] are random Gaussian variables. In addition, we assume that the thermal noise is ergodic: averaging over time and noise realization are equivalent here:

A n A m thermal = 2kT ω 2 n δ nm [J s 2 ],
where 2kT/ω 2 n is the thermal noise power spectrum. Therefore, the time derivative of the correlations averaged over thermal noise, which mainly consists in a succession of sources random in time and position, does contain the reference GF:

∂ τ C AB (τ) = n 2kT ω n φ n (A)φ n (B)e -iω n τ thermal = 2kTG(A, B, τ). (3.13)

Correlation of acoustic and elastic waves in homogeneous media

We now focus on open media, where modes are now continuously distributed in the frequency domain. Since the Heisenberg time is now infinite, the modal expansion of the wavefield will be of little help. Nevertheless, we propose different derivations that connect to the now well established result that averaged correlations yield the GF of the medium.

Sources in the far-field

Let a series of sources be placed in the farfield of two receivers A, B. Assuming the incoming plane wave is a very broadband, and θ the incident direction (see Fig. 3.1), the wavesfield sensed at

A is G A (t) = δ(t), and at B is G B (t) = δ (t -(D/v) cos(θ)).
The correlation for a given angle theta reads: An averaging over the sources position is equivalent to an averaging over the angle θ. In the 3-D case, this integration yields:

C AB (τ, θ) = δ τ - D v cos(θ) = 1 2π ∞ e iω(τ-D v cos(θ)) dω. ( 3
C AB (τ) = 1 2π ∞ e iω(τ-D v ) -e iω(τ+ D v ) Dω/v dω. (3.15)
After the inverse Fourier transform, we find again that the time-derivative of the averaged correlation gives the causal and the anti-causal GF:

∂ τ C AB (τ) = v D (δ(τ -D/v) -δ(τ + D/v)) . (3.16)
In the case of an elastic medium, one additional condition is needed. Not only the waves have to come from all directions to ensure a perfect angle averaging, but also do they have to be fully equipartitioned [95]. Equipartition has been discussed in Section 2, it results in an equilibrium between the energy of longitudinal and transverse waves (noted E P and E S ): E S /E P = 2v 3 P /v 3 S . The derivation proposed by Sanchez-Sesma and Campillo [95] establishes clearly this result.

Continuous distribution of sources: the stationary phase approximation

Snieder [96] proposed to extend this result to a homogeneous source distribution, without invoking any far-field approximation but the stationary phase theorem. This means that one additional averaging of the correlation is made, so that the sources might be in the near-field of the receivers. Let the sensor's position be x = ±D/2, y = 0, z = 0. The source S position (x, y, z) is random, the sourcereceivers distances are

r A = (x + D/2) 2 + y 2 + z 2 , r B = (x -D/2) 2 + y 2 + z 2 .
The 3-D acoustic GFs are:

G S,A or B (t) = v 4πr A or B δ(r A or B -vt).
The correlation for a given source S in the frequency domain reads:

C S AB (ω) = v 2 4 3 π 4 r A r B e iω r A -r B v .
As previously, we average this correlation over any position of S with a density of probability n: ∞ nC S AB (ω)dxdydz. This integration can be performed invoking the stationary phase approximation, which is recalled in the following. Let f (r) be a C ∞ function; if there exists a unique point r 0 at which f (r 0 ) = 0 and f (r 0 ) 0, then for ω → ∞:

g(r)e iω f (r) dr = g(r 0 )e i f (r 0 ) 2π ω f (r 0 ) •
This theorem is applied twice to the averaged correlation with f (r) = (r Ar B )/v and g(r) = 1/r A r B for r = y and then r = z at y = 0 and z = 0. The average correlation the reads:

C AB (ω) → ∞ nv 2 4 3 π 4 2π v ω 1 |x A -x B | e i ω v (x A -x B ) dx with x A = |x + D/2| and x B = |x -D/2|.
The correlation simplifies into:

C AB (ω) = v 3 32π 3 ω e ikD + e -ikD D (3.17)
which leads again, in the time domain, to

∂ τ C AB (τ) ∝ G + AB (τ) -G - AB (τ)
. Note that Roux et al. [START_REF] Roux | [END_REF] proposed an elegant change of variables that allows a rapid evaluations of the 3-D space integrals and leads to comparable results.

The Helmholtz-Kirchhoff's theorem

In the two previous derivations, we assumed either a perfect azimuthal distribution of distant sources or a homogeneous distribution over the whole space. This means that in a 3-D experiment, correlations need to be average over a 3-D distribution of sources to converge to the GF. As noticed by several authors, for time-reversal purposes [97] (see next section) and later for correlations [98,99], this condition is not absolutely necessary, especially if we take advantage of all the properties of the wave equation. If A and B are included in a domain D with boundary ∂D (a surface in 3-D, a closed line in 2-D), then the GF between A and B can be estimated by the correlation averaged over the surface

∂D ∞ ϕ i (S, A, t)ϕ j (S, B, t + τ)dtd 2 S = G ij (A, B, τ) + G ij (A, B, -τ). (3.18)
The equality is valid provided reciprocity is ensured: G(S, A, τ) = G(A, S, -τ).

We have here defined new a condition for the correlations to yield the GF: they must be averaged over a continuous distribution of sources forming a surface that encloses the medium being imaged and the two passive receivers. This condition is less restrictive than a distribution of sources in the whole volume.

The time-reversal interpretation: numerical simulation

In the previous theoretical derivations, correlations are found to converge to the Green's function under important assumptions. They rely on a perfect source distribution resulting in a perfect diffuse and equipartitioned wavefield. What happens when there is only a finite number of sources in the medium and when the field is not completely equipartitioned? And most of all, what is the physical meaning of correlations? In the following we propose an interpretation of correlations based on an interpretation with Time-Reversal [23]. To illustrate this interpretation, we present 2-D numerical experiments of acoustic scattering on rigid inclusions randomly located either in a closed cavity or in a open medium.

The cavity equation

To begin, let us consider two receiving points A and B and a source S placed amongst a random collection of scatterers, as represented in Figure 3.2. The scatterers are in water; only lossless acoustic waves are considered here. The wave equation is solved by a finite differences simulation (centered scheme); the boundary conditions are implemented following Collino's work [100]. Naturally, a finite-difference scheme shows numerical dispersion. However, the essential point is that the fundamental symmetries of reciprocity and time-reversal still hold in the numerical experiments. At the edges of the grid, the boundary conditions may be either perfectly reflecting (Dirichlet) as in a closed cavity or absorbing (open medium). The signal transmitted by S is a pulse with a center frequency 1 MHz and a Gaussian envelope (σ = 0.7 μs). We will note G IJ (t) as the impulse response between I and J, i.e., the wave field sensed at I when a Dirac δ(t) is sent by J. If e(t) is the excitation function at S, then the wave field ϕ A and ϕ B received at A and B will be respectively e(t) ⊗ G AS (t) and e(t) ⊗ G BS (t), ⊗ representing convolution. The cross-correlation of the fields received at A and B is then

C AB (t) = Φ A (t)Φ B (t + θ)dθ = G AS (t) ⊗ G BS (-t) ⊗ f (t) (3.19)
with f (t) = e(t) ⊗ e(-t). f (t) depends only on the excitation imposed at the source, whereas the information regarding the impulse response between A and B is hidden in G AS (t) ⊗ G BS (-t). Indeed, the impulse responses of a closed cavity satisfy a remarkable property, as shown by [87], which he termed the "cavity equation":

G AS (t) ⊗ G BS (-t) = G AB (t) ⊗ G SS (-t). (3.20)
For this relation to hold, the cavity must be lossless, and its eigenmodes must not be degenerate. Note that, in practice, the correlations cannot be performed over an infinite time interval (the ring time of a cavity is infinite if it is lossless); therefore, the cavity equation can be compared to experimental results if the integration time ΔT is sufficiently large compared to 1/Δω, with Δω the characteristic distance between modes, so that the modes are resolved. 1/Δω is sometimes referred to as the Heisenberg time, or break time. Figure 3.3 illustrates the validity of the cavity equation; here the impulse responses have been recorded during an integration time of 80 ms (2 × 10 6 time steps), and the Heisenberg time is ≈5 ms. From Draeger's cavity equation, the correlation between the fields received at A and B is:

C AB (t) = G AB (t) ⊗ G SS (-t) ⊗ f (t). (3.21)
Therefore, similarly to Weaver's results [80,85,86] the direct Green's function G AB is present in the correlations of the field within a closed cavity and can be recovered from C AB provided that the G SS term can be properly deconvolved, at least in the frequency domain limited by the spectrum of f (t). The estimation of G SS is possible after averaging C AB over all A and B positions.

Comparison between C AB (t) (thick continuous line) and

G AB (t) ⊗ G SS (-t) ⊗ f (t) (dotted line).
The overall correlation coefficient between the two waveforms is 98.7%.

Simulation of perfect time reversal

Is this valid in an open medium? We have conducted the same numerical experiment, with the same distribution of scatterers, but with absorbing instead of reflecting boundary conditions. As a result, the cavity equation is no longer valid, and the correlation of the scattered field C AB shows no resemblance whatsoever with the Green's function (the correlation coefficient function between wave forms represented in Figure 3.4 is less than 0.5%).

However, a physical argument indicates that the Green's function can still be recovered from the correlations C AB , even in an open medium, if several judiciously distributed sources are used instead of a single point S. To that end, we propose to analyze the experiment in terms of time-reversal symmetry. Indeed, there is a strong link between correlations of a diffuse field and time reversal [101].

Because the scatterers do not move and there is no flow within the medium, the propagation is reciprocal, i.e. G IJ (t) = G JI (t). When we cross-correlate the impulse responses received at A and B, the result G AS (t) ⊗ G BS (-t) is equal to G SB (-t)⊗G AS (t). Now, imagine that we do the following time-reversal experiment: B sends a pulse, S records the impulse response G SB (t), time-reverses it and sends it back. The resulting wave field observed at A would then be G SB (-t) ⊗ G AS (t), which, because of reciprocity, is exactly the cross-correlation G AS (t) ⊗ G BS (-t) of the impulse responses received at A and B when S sends a pulse. We would like the direct Green's function G AB to appear in the cross-correlation. But in the most general case, G SB (-t) ⊗ G AS (t) has no reason to be equal to G AB , as was shown in Figure 3.4. Yet we can go further: imagine now that we use several points C to perform the time-reversal operation, and that we place them in such a way that they form a perfect time-reversal-device, with no loss of information. Following from the Helmholtz-Kirchoff theorem, such would be the case if the sources S were continuously distributed on a surface surrounding the scattering medium. Then the time-reversal operation should be perfect. During the "forward" propagation, B sends a pulse that propagates everywhere in the medium (including at A where the field received is G BA (t)). This pulse may be scattered many times and is eventually recorded on every point of the time-reversal device, with no loss. When the field is time-reversed, since nothing of it has been lost, it should exactly go backwards in time (and refocus on B at time t = 0) everywhere in the medium, which implies that the field received at A after the time-reversal is exactly G BA (-t), the time-reversed version of the direct Green's function. Then, once the wave has refocused on B (at time t = 0), it does not stop since there is no "acoustic sink" at B [102]: the wave diverges again from B and gives rise, at times t > 0, to G BA (t) at A. Thus, if we use a collection of sources C arranged in such a way that they form a perfect time-reversal device, we should have

S G AS (t) ⊗ G SB (-t) = G BA (t) + G BA (-t).
(3.22)

A more detailed analysis, taking into account the monopolar or dipolar nature of the source/receivers, is given by [97]. Equation (3.20) implies that the impulse response G BA (t) can still be recovered from the correlation of a diffuse field, even in an open medium, provided that the sources S are distributed judiciously, and all the correlation functions are summed over the source positions. Unlike the case of a closed medium, no additional deconvolution by G SS is needed. From this time-reversal interpretation, we deduce a condition for the Green's function to emerge from cross correlations in open media: the sources S must be placed so that they form a perfect time-reversal device.

We have checked this in the numerical experiments depicted in Figure 3.5. The results are in excellent agreement with equation (3.22), as is shown in Figure 3.6: the degree of correlation between waveforms is 97.4%. Of course, when the sources are not placed as a perfect mirror, as presented in Figure 3.5b, the results are less good (the degree of correlation between waveforms is 81.9%) because one part of the waves are not recorded by the time-reversal device due to the presence of scatterers outside the sources. Yet the main features of the Green's function can still be recognized, even at late times. If the number of sources is decreased, the reconstruction of the Green's function is less satisfactory, as shown in Figure 3.7. With only 50 sources (instead of 250 previously) regularly spaced every ≈ 5λ/3 on the perimeter of a circle as in Figure 3.5a, the correlation coefficients between waveforms is 70%. However if the 50 sources are gathered together in a 72 • angular sector (pitch λ/3), it drops to 53%. Indeed, since the coherence length of a diffuse wave field is λ/2, it is useless to place the sources closer.

Application to noise

So far, we have considered that the origin of the field measured at A and B was an active and coherent source transmitting a short pulse (or a collection of such sources). What if there are no such sources in the medium, but diffuse continuous noise? The physical origin of this noise may be thermal vibrations [80]. In seismology, noise in the seismograms comes from a variety of different sources (traffic, sea waves, weather, human activity, etc.) continuously and (allegedly) randomly pumping energy into the earth and essentially exciting surface waves. In ocean acoustics, noise may originate from boats, surf, wind, animals, etc. By definition, the noise sources cannot be controlled. In the light of the discussion above, in order to recover the Green's function from the cross-correlation of the noise received at A and B, the most favorable situation would be that in which noise can be considered as coming from virtual point sources S randomly distributed everywhere in the medium and continuously generating uncorrelated white noises n S (t). In that case, the cross-correlation between waveforms sensed at A and B would be

S S G AS (-t) ⊗ n S (-t) ⊗ h S B (t) ⊗ n S (t). (3.23)
If the observation time ΔT is long enough compared to the correlation time of the noise, then n S (t) ⊗ n S (-t) will converge to δ(t)δ(SS ). Moreover, if the virtual noise sources S are distributed everywhere in the medium (in other words, if each degree of freedom is excited randomly and independently) then the S-points would necessarily constitute a perfect time-reversal device, so equation (3.21) should be verified again.

We have carried out a numerical experiment based on this idea. Two hundred forty sources where distributed at random inside the open scattering medium shown in Figure 3.2, and 240 uncorrelated white noises, convolved by e(t), were transmitted by these sources during 40 ms. The resulting wave forms are received at A and B. Their cross-correlations C AB (t) is compared to the direct Green's function G AB (t) ⊗ f (-t): the agreement is still very good (61% correlation coefficient) even at late times.

In this section, we have argued that recovering the Green's function was possible in closed and open multiple scattering medium and we have proposed a criterion based on reciprocity, time-reversal symmetry, and the Helmholtz-Kirchoff theorem: if sources are placed as if they were to form a perfect time-reversal device, then the Green's function can be recovered by summing the cross correlations. This has been validated by numerical experiments. The reduction of the number of sources was also discussed, and the possibility of using noise sources was illustrated. There is still much food for thought, particularly regarding the role of scatterers in the reconstruction of the Green's function. The argument we developed here is valid for any medium (homogeneous, high-order multiple scattering, reverberant, etc.) where reciprocity and invariance under time-reversal hold. The field does not need to be thermally diffuse for the Green's function to emerge from the correlations, as long as there are enough well-positioned sources. Another approach is to consider the scatterers as secondary sources which are necessary to truly randomize the wave field emanating from a single original source. Correlating the late part of the coda would then permit us to reconstruct at least the early arrivals of the Green's function. This will be discussed in the next section.

Correlation and the time-reversal interpretation:

ultrasonic experiment

Perfect source distribution

In the following, we report an ultrasonic experiment designed to retrieve the Green's function of an open medium beyond the ballistic front. We aim at retrieving the signature of two localized and well-separated scatterers. A perfect average over sources is realized, following the Helmholtz-Kirchoff theorem and similarly to a perfect time-reversal experiment.

The experiment [129] is carried out in clear water (sound speed c = 1.480 mm/μs, average wavelength λ ≈ 1.5 mm). Two reflectors A and B are placed in the vicinity of a single ultrasonic piezoelectric transducer R (Fig. 3.8). This sensor is nearly omni-directional. Its (reciprocal) response r(t) has a flat spectrum between 0.5 and 1.5 MHz. Its width is 0.5 mm, and its height ∼10 mm. The reflectors are two empty aluminum cylinders of radius 4 mm. The experimental setup is quasi-2D (Fig. 3.8). The source S sends a broad-band pulse (typically two oscillations at 1 MHz). The resulting wave R(t) is recorded at the receiver R and contains the direct wave S → R and scattered contributions (S → A → R, S → B → R, S → A → B → R and so on...). The impulse response between S and R is noted h SR (t) so the recorded signal is R(t) = r(t) ⊗ h SR (t) ⊗ s(t), ⊗ meaning convolution. For a given source S the auto-correlation function reads: In general, this correlation has no reason to equal the Green's function h RR (t).

C S (τ) = R(t)R(t + τ)dt (3.24) = r(t) ⊗ r(-t) ⊗ h SR (t) ⊗ h SR (-t) ⊗ s(t) ⊗ s(-t). ( 3
Yet different approaches showed that if the correlation is averaged over a large number of sources perfectly surrounding the medium, the exact impulse response h RR (t) can be fully retrieved [23,[START_REF] Roux | [END_REF]96,97,99]. Here such a source distribution is obtained by rotating the source with a step-motor (Fig. response F(t) (continuous line). This reference is the convolution of r(t)⊗h RR (t)⊗r(t) (pulse-echo measurement) with s(t) ⊗ s(-t) (the source auto-correlation). The two waveforms show a degree of resemblance of 95.4%.

This is an experimental proof that the averaged correlation C(τ) reveals the detailed structure of this open medium. The wavefront labeled (b) in Figure 3.12 corresponds to the reflection on the first cylinder. It is followed by a weak pulse: the first cylinder's creeping wave. The second strong peak (c) is the reflection on the second cylinder, and the latest identified contribution (d) corresponds to a more complex ray path: R → A → B → R, i.e. a multiply scattered wave. This experiment is a direct illustration of the Helmholtz-Kirchhoff's theorem, that states that the exact GF of the medium between two passive sensors is reconstructed by averaging correlations obtained for a continuous set of sources perfectly surrounding the medium to image and the two receivers. By reciprocity, this experiment is also equivalent to a perfect time-reversal experiment as presented earlier.

Imperfect source distribution

This section is now devoted to a more realistic configuration, where a perfect source distribution is not achieved. We address the issue of causality and show the role of multiple scattering for the reconstruction of the Green's function. Ultrasonic experimental results are presented to illustrate the argument. We particularly address two issues:

1. physically, the GF G AB is causal, but the correlation between the wave fields received at A and B may be non-causal; therefore, should one keep the causal part, the anti-causal part of the correlation, or both to estimate G AB ?

2. in an inhomogeneous medium, what is the role of scattering in the reconstruction of the GF from field-field correlations?

We also present experimental results to support the argument.

As developed in the previous section, correlations are analogous to a two step time-reversal operation. During the time reversal experiment, (1) the wave is emitted at A and sensed in several points S perfectly surrounding the medium. Then (2) the field is time-reversed in the medium and (3) is finally sensed at B. The wavefield sensed at B is the actual GF G(A, B, t) provided that S are indeed perfectly distributed, i.e. placed so that they would form a perfect TR device. This is therefore valid in the correlation experiment where A and B are passive receivers, and S are independent sources: the impulse response G AB (τ) can be retrieved from either the causal (τ > 0) or the anti-causal part (τ < 0) of the sum of field-field correlations C AB (τ).

In real experiments, whatever the type of waves involved, this condition is hard to meet. In seismology, for instance, the displacement field at the earth surface is recorded by seismic stations (A, B) but the sources (S) of the earthquakes are far from being arranged as a perfect TR device: they are mostly aligned along faults. Yet, the elastic GF can be partially retrieved using correlations of the late seismic codas produced by distant earthquakes [92]. Why is this possible when the TR criterion is not fulfilled? A laboratory experiment can help us find an answer and shed light on the role played by multiple scattering (Fig. 3.13). The coherent field is totally extinct in the received waveforms G AB and G AS . Since L * , the correlation length of the field emerging from the slab is ∼ λ/2, and, from the Van Cittert-Zernike theorem the fields sensed at A and B are spatially incoherent since the transverse size of the slab is 25 cm while the distance between B and the slab varies between 15 and 25 cm.

Here, the GF between A and B is a well-defined pulse arriving at time |AB|/v, followed later by lower amplitude reflections from the rods. The experimental results show that the emergence of the GF from C AB (t) highly depends on the position of B and on the number of sources employed. With only one source (Fig. 3.14a) C AB (t) is too noisy to see the emergence of the GF. But at the same point B, with 21 sources (Fig. 3.14b) instead of one, C AB (t) shows a strong peak at time t = 33.5 μs, which is exactly the travel time |AB|/v. Yet for a different position of B (Fig. 3.14c), C AB (t) shows a peak at time t = -33.35 μs, the negative of the expected travel time. And for a third position of B (Fig. 3.14d), even with 21 sources, the GF does not emerge at C AB (t). So it appears that the GF (at least its first arrival) can indeed be recognized in the correlation C AB (t), but only at certain times (causal or anti-causal), and for certain positions of B relative to A. Why is that?

Role of sources and scatterers positions

The TR interpretation gives the answer: for this particular setup, in a fictitious TR experiment where A would be the source and the 21 S-points a finite-size TR device, the time-reversed pulse would hit B at times t < 0 only if B is between A and S (i.e. x B > 0), and at times t > 0 only if it is behind A (x B < 0). Consequently, when one cross-correlates two wavefields in order to reconstruct the GF of an unknown medium, one has to know the location of the receivers relative to the sources and to the scatterers in order to keep only the relevant part of C AB (t).

It also appears that G AB cannot be properly reconstructed for any position of B, as was shown in Figure 3.14a. Figure 3.15 compares C AB (t) to the theoretical travel times for 61 positions of B (x B = -50 mm, y B = -15 to 15 mm), with and without the rods: the curvature of the GF emerges only in the former case. This emphasizes the role of multiple scattering: the region for which the arrival time is well retrieved is much smaller in a homogeneous medium (water) than through the forest of rods. This too can be interpreted via the TR interpretation: a TR experiment works better (meaning that it better reconstructs the "initial scene") through a multiple scattering medium [6] than in a homogeneous medium (here If no scattering, the GF is reconstructed only in the red (B 1 ) area. With the scattering device, the area is enlarged and represented in green (B 2 ). In both cases, because of the finite aperture of the array of sources or scattering device, the GF is not recovered in the other area (B 0 ). The + andsigns denotes the areas where the causal or the anti-causal GF is reconstructed.

the "initial scene" would be the propagation of a spherical pulse emitted by A, Fig. 3.16).

In our experiment, we do not satisfy the equipartition condition nor the "perfect TR device" condition, but we approach it using several sources instead of one. Our experimental results show that with a limited number of sources it is possible to estimate at least the first arrival of the GF, not everywhere at every time but in a limited area whose size is larger in the presence of multiple scattering.

Statistical approach and convergence rate

In the following, we propose a third possible approach for the correlations, which differs from the time-reversal interpretation or the analytical derivations presented previously. For practical application like imaging, one need to know precisely how much data is required to obtain a "good" image. In other words, how long do we need to record the diffuse field or the noise, and how many sensors do we need? We already noticed that:

• C AB converges to G AB with increasing source number (source averaging);

• C AB converges to G AB with increasing time record (time averaging);

• C AB converges to G AB with increasing bandwidth;

• the rate of convergence is greater if sensors are closer together;

• the rate of convergence is greater at low frequency.

What is the physical interpretation of these observations? Can we propose a model for the rate of convergence? We address in these issues in the following section.

Numerical simulations in a 2-D homogeneous medium

Let's start with two receivers A and B placed in a 2-D homogeneous medium (Fig. 3.17), D is the distance between A and B. S is a source randomly placed in the medium, sending a broadband pulse (Δω ≈ ω) whose autocorrelation is f (t). The acoustic wave sensed at A and B are correlated for each source position S. The correlation is numerically estimated using the Fourier transform of the acoustic GF:

C S AB (τ) = 1 2π G 0 (ω, SA)G * 0 (ω, SB)e iωτ dτ ⊗ f (τ), (3.26)
where G 0 (ω, r) = (-i/4)H 0 (kr) is the 2D acoustic GF, and H 0 is the first-order Hankel function of the first kind. k is the wave vector, corresponding to the wavelength λ and the celerity v. This calculation is repeated for different source positions, and correlations are eventually averaged. Note that no other processing, like normalization or 1-bit digitization, is applied. This is very important to recover the exact phase and amplitude of the actual GF. Correlations are compared to a reference time-symmetrized GF, RG AB (τ), defined as [G AB (τ) + G AB (-τ)]⊗ f (τ). Note that no time derivation is needed here. Correlations C AB (τ) N are averaged over different source numbers N. Results are plotted in Figures 3.18 and 3.19 for, respectively, D = 1λ, 3λ, and 10λ. In each figure, correlations are averaged over N = 1, 10, 100, and 5 000 sources. The correlations are compared to the reference GF G AB (τ) (red dotted line). For a few sources only, correlations have no reason to equal the GF. But for a number of sources as large as N = 5 000, correlations were found to have well converged to the GF.

The common way to estimate the quality of a given signal is to evaluate the "Signal to Noise Ratio" (SNR). In the present manuscript, this denomination is slightly hazardous for at least three reasons:

1. seismic or acoustic noise might be used as the actual signal to correlate; 3. practically, when correlations have not converged to the GF yet, they show uniform fluctuations (SNR 1). In that case, how does one evaluate the level of the signal?

In this manuscript, the "signal" is the reference GF RG AB and the "noise" is made of the remnant fluctuations of C AB . We choose to quantify the rate of convergence with the normalized correlation coefficient CC between the correlations C AB and G AB :

CC = ∞ RG AB (τ)C AB (τ)dτ ∞ RG 2 AB (τ)dτ ∞ C 2 AB (τ)dτ • (3.27)
If the convergence estimator CC = 1, the averaged correlation perfectly matches the reference symmetrized GF, and remnant fluctuations are canceled. If fluctuations are not negligible, CC < 1. If correlation shows a purely random waveform, than CC ∼ 0. By increasing the averaging, we expect CC to rise from 0 to a value as close to 1 as possible. Studying the rate of convergence of the correlations will consist of formulating CC as a function of the frequency f , the distance D between A and B, and the number of sources N.

The concept of coherent zones

The concept of coherent zones is based on a geometrical interpretation of correlations. If a source S is near the axis A -B but outside this segment, it will contribute to the reconstruction of the GF (or at least the direct wave) in the correlation. If this first assertion is not valid, the source will not contribute to the GF but noise in the correlation. The ensemble of points contributing to the GF in the correlation is named the "coherent zone" (Fig. 3.20). The idea was suggested by Snieder [96], and derived by Roux et al. [75] and Sabra et al. [91] who named the coherence zones "end-fire lobe".

Here follows a geometrical definition of the coherent zones. To start with, let's consider the GF G AB of an homogeneous medium and the correlation C AB , both restricted to a frequency bandwidth Δω and approximated by:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ G AB (t) = ω+ Δω 2 ω-Δω 2 cos ωD v -ωt dω C S AB (t) = ω+ Δω 2 ω-Δω 2 cos ω v |SA -SB| -ωt dω.
(3.28)

Note that only the phase of the GF and of the correlation is studied, not the amplitude. We stress that the coherent zones are defined by phase interferences in the correlations, and that amplitudes do not play any role in their definition. If the source S is on the AB axis, then |SA -SB| = D and correlations are equivalent to the direct GF. If S is moved slightly from this axis, it will partially contribute to the GF, and partially to noise. To quantitatively evaluate how much a source S might contribute to the GF depending on its location, we introduce the correlation coefficient CC(S) between C AB and G AB :

CC(S) = ω+ Δω 2 ω-Δω 2 cos ω v |SA -SB| cos ωD v dω. (3.29)
This formula expands:

CC(S) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 2Δω cos ω v (|SA -SB| + D) sinc Δω 2v (|SA -SB| + D) [a] +2Δω cos ω v (|SA -SB| -D) sinc Δω 2v (|SA -SB| -D) . [b]
(3.30) The degree of resemblance CC(S) depends on the arguments of the sinc functions (Fig. 3 are more rapid then the fluctuations of the sinc term. The correlation coefficient CC(s) is therefore not negligible until the cos term drops the first time to zero (i.e. from 0 to π/2). In average, the contribution of the cos term is noticeable for 0 < ω(|SA -SB| -D)/v < π/3. To conclude, S contributes to the GF if:

|SA -SB| < D - λ 6 • (3.31)
This last equation defines the coherent zones: the ensemble of points S that meet this inequality. These zones are delimited by two hyperbola defined by SA -SA = D -λ/6 with focal points at A and B (Fig. 3.22). These coherent areas might be thought of as an anti-Fresnel zone, in reference to the seismic Fresnel zones between a source and a receiver. The two coherent hyperbolas, in a 2D medium, are defined as:

x 2 -y 2 D D -λ/6 2 1 = D 2 4 1 - D -λ/6 D 2 , (3.32)
which can be rewritten in polar coordinates:

r(θ) = D 2 -(D -λ/6) 2 sin 2 (θ) - D D -λ/6
θ 0 = ±arctan D D -λ/6 2 -1,
which can be approximated in the case of D λ by θ 0 → ± √ λ/3D. This derivation is similar to the one proposed by Roux et al. [75]. In their far-field model, the source S is defined by the angle of incidence θ. The correlation coefficient then reads:

CC(θ) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 2Δω cos ωD v (1 + cos (θ)) sinc ΔωD 2v (1 + cos (θ)) [a] +2Δω cos ωD v (1 -cos (θ)) sinc ΔωD 2v (1 -cos (θ)) . [b]
(3.34) Again, CC is dominated by the term [b], which after a second order expansion leads to:

CC(θ) = 1 - θ 4 D 2 8v 2 (ω 2 + Δω 2 /12) + o(θ 4 ). (3.35)
The coherent angle is therefore

θ 0 = 2 √ 2v/Dω ≈ √ λ/2.
2 D, which is equivalent to the angle proposed earlier. An important point to notice is the dependence of the coherent zone on the non-dimensional parameter λ/D. A series of examples of coherent zones are displayed in Figure 3.23 for different λ/D ratios. The longer the distance, or the smaller the wavelength, the smaller the coherent zones. This quantitatively confirms the observations at the beginning of this section: it's more difficult to reconstruct the GF through correlations at high frequency, or if the receivers are far apart than at lower frequency, or for closer receivers.

Convergence rate in homogeneous medium

In the following, we propose a statistical model of the convergence rate based on a geometrical approach. We assume that the source has a random position in a 2-D space:

• if the source is in the coherent zone around A, we assume that it will yield a coherent contribution to G A→B (τ) = G + AB ;

• if the source is in the coherent zone around B, we assume that it will yield a coherent contribution to G B→A (τ) = G - AB ;

• if the source is in the incoherent area, it will contribute to the remnant noise in the correlation.

-2 0 2 -1 0 1 D = 1 × λ -2 0 2 -1 0 1 D = 1 × λ -2 0 2 -1 0 1 D = 3 × λ -2 0 2 -1 0 1 D = 5 × λ -2 0 2 -1 0 1 D = 7 × λ -2 0 2 -1 0 1 D = 9 × λ Figure 3.23.
Coherent zones ϑ for different ratios of D/λ. The greater the distance D between the receivers, the smaller the coherent zones, and the slower the convergence of the average in the passive reconstruction to the GF.

Contributions from different sources within the coherent zones will add up linearly, resulting in constructive interferences. Contributions from sources in the incoherent area will add incoherently, resulting in random interferences. For N normalized sources evenly distributed in the whole surface or volume, N c is the number of sources in the coherent zone ϑ around A or B, and N i is the number of sources in the incoherent zone ϑ. Then the expected intensity of the GF in the correlation C N AB is proportional to N 2 c , and the intensity of the remnant noise is proportional to N i . CC(N) is the convergence estimator of the correlations averaged over N sources. Given the coherent/incoherent contribution of sources, the convergence rate reads:

CC(N) = Nα N 2 α 2 + N(1 -α 2 ) , (3.36)
where α is the statistical contribution of one source. This contribution is the ratio of coherent to incoherent surfaces, weighted by a geometrical factor that arises from the geometrical spreading between the source and the two receivers:

α = ϑ 1 SA • SB dS ϑ 1 SA • SB dS • (3.37)
Note in the far field approximation, α is just a the surface ratio that simplifies into α = 2θ 0 /2π with θ 0 = ± √ λ/3D. For D λ, in that case, the probability for a given source to be in the coherent zone is (2/2π) √ λ/3D, and the signal to noise ratio of the correlations can be approximated by This result is close to the one proposed by Weaver et al. [106], and also Sabra et al. [107]. In the case of real ambient noise, the medium is continuously excited by random sources. The number of independent sources in a record of length ΔT and bandwidth Δω is roughly ΔTΔω. The signal-to-noise ratio in the case of random noise is therefore:

RSB = 2 2π Nλ 3D [2-D]. ( 3
RSB = 2 2π ΔTλ 3DΔω • (3.39)

Verification of the convergence rate using numerical simulations

Correlation of the wave field originating from a random source S is evaluated numerically using the same calculation as Section 5.1:

C S AB (τ) = 1 2π G 0 (ω, SA)G * 0 (ω, SB)e iωτ dτ ⊗ f (τ). (3.40)
Then correlations are averaged over N random sources, and they are compared to the reference GF. N varies from 1 to 5 000 sources. The correlation coefficient CC(N) is additionally averaged over the realizations of the source distributions.

For our theoretical model of convergence rate, the coefficient α is also evaluated numerically using 10 6 positions. Curves for the convergence rate resulting from our theoretical model and from simulations are displayed in Figure 3.24 for different distances D between A and B ranging from λ to 10λ. An excellent agreement is found between simulations and our statistical model. We conclude that the convergence rate of correlations in a homogeneous medium reduces to the evaluation of the geometrical factor α that coincides with the normalized contribution of one source. As α is linked to the coherent areas, we found that this model quantitatively describes the dependence of the convergence rate on the frequency f and the distance D: convergence is slower at high frequencies and for long distances between receivers.

Role of scattering and time-symmetry in the correlations

In this chapter, we discuss quantitatively the impact of the anisotropy of energy flux on the time asymmetry of the Green's function (GF) recovered by crosscorrelations. To that end, we will propose a theoretical approach followed by numerical simulations, ultrasonic experiments, and finally an application to real seismic data.

Multiple scattering and isotropy 1.Theory

For sake of simplicity, we consider here the theory for scalar waves. The results are not expected to be different for the energy of elastic waves as it will be discussed later. The angular distribution of energy at position R and time t is described by a specific intensity I R ( p, t) defined as the energy flux in space direction p per unit solid angle. The specific intensity is the solution of a radiative transfer equation, which can be derived from an ensemble average of the wave equation [2,16]. The radiative transfer equation is a local detailed energy balance which describes the transport of energy through a multiple scattering medium. In general, it is an integro-differential equation that can be solved only numerically. However, after many scattering events, the initial distribution of energy in phase space tends to be homogenized, which implies that the angular dependence of the specific intensity departs only slightly from isotropy. When this assumption is valid, I R can be expanded as follows:

I R ( p, t) = 1 4π ρ(R, τ) + 3J(R, τ) • p + • • • (4.1)
where ρ denotes the energy density, and J is the energy current vector which points in the direction of maximum energy flow. The dots denote higher order multipoles that are neglected. Equation (4.1) forms the basis of the diffusion approximation, which should apply at t τ * , where τ * denotes the transport mean free time. We recall that the transport mean free time τ * is the typical time after which the scattered energy of a wave in a particular direction is spread over all directions. We introduce the "ideal" field-to-field correlation function C between two stations located at R + r/2 and Rr/2 as

C R (r, τ) = +∞ -∞ u(R + r/2, t + τ/2)u * (R -r/2, t -τ/2)dt, (4.2)
where R represents the mean distance between source and stations and r the interstation distance. We will assume now that the averaging over time and over the ensemble can be interchanged. The diffusion approximation can be used to derive an asymptotic (t → ∞) relation between C and the ensemble average GF of the medium [START_REF] Van Tiggelen | [END_REF]:

C R (r, τ) = ρ(R, τ) ∂ ∂τ [ G B (r, τ) -G B (r, -τ) ] -3J(R, τ) • ∇ [ G B (r, τ) -G B (r, -τ) ] .
(4.3) In equation (4.3) brackets • denote an ensemble average, and G B is the retarded causal GF filtered in frequency band B. Equation (4.3) applies to ensemble averaged quantities only and is, therefore, not expected to apply strictly in the seismological case. However, it establishes a relatively simple relation between the field-to-field correlation function from a single source and partial derivatives of the (time-symmetric) GF.

It can be inferred that the partial derivatives ∂ τ and ∂ r acting on the brackets [•] yield even and odd functions of time respectively. Thus, as long as the dipolar (J) and isotropic (ρ) terms are of the same order, a time asymmetry is expected to persist.

In the diffusive regime, the energy density is the solution of a simple diffusion equation:

∂ t ρ(R, t) -D∇ 2 ρ(R, t) = δ(t)δ(R) (4.4)
and is related to the energy current by Fourier's law:

J(R, t) = -D∇ρ(R, t). (4.5)
The diffusion constant of the medium is denoted by D and is related to the transport mean free time τ * by

D = v 2 τ * /3. (4.6)
In equation (4.4) the delta functions represent idealized source terms for small earthquakes. Note that equations (4.4) and (4.5) are valid for coupled elastic waves. In that case, ρ must be interpreted as the sum of P and S energy densities and relation (4.6) takes a more complex form [108]. For scalar waves in a simple infinite scattering medium with homogeneous background, the diffusion equation is easily solved, and the ratio Γ between the causal and anti-causal parts of the correlation function can be written as [20] Γ

(R, t) = 1 + 3R/2vt 1 -3R/2vt • (4.7)
This relation is easily obtained by noting that J/ρ = R/2t and ∂ r = v∂ τ for a propagating wave of the general form h(t -R/v). Equation (4.7) shows that for a single source the convergence of the ideal correlation function toward time symmetry is algebraic, of order t -1 only. This result has to be used with some caution, since it relies on the assumption that the angular dependence of the specific intensity can be described by equation (4.1).

Monte-Carlo simulation of radiative transfer

To assess the range of validity of this expansion, we numerically solved the full elastic radiative transfer equation using the Monte Carlo method of [109]. Our analysis is limited to elastic body waves in an infinite space. The code previously developed to evaluate energy densities has been adapted to calculate the angular distribution of flux. Numerical solutions of the transport equations and analytical results of the diffusion approximation are shown in Figure 4.1. The medium is composed of spherical inclusions with slight (5%) deviations of density, P, and S velocities from the homogeneous background. The product of shear wavenumber and sphere radius is k S a ≈ 2. The scattering mean free time of shear waves is roughly 8 s and the detector is located 80 km away from a point-like shear source. Figure 4.1a demonstrates the rapid mixing of P and S modes. After about 50 s (six mean free times), the P to S energy ratio has stabilized. However, Figure 4.1b shows that the energy flux is strongly anisotropic at the same lapse time: the ratio between forward and backward fluxes is still larger than 4. It is therefore important to carefully make the distinction between the stabilization of the P to S energy ratio, the validity of the diffusion approximation, and equipartition.

In the equipartition regime, all directions of propagation and all polarizations are supposed to be equally represented. Under these condition, the theoretical value of the energy ratio v 3 S /2v 3 P (where v P and v S denote the shear and compression velocities) can be obtained from a simple mode counting argument [START_REF] Weaver | [END_REF]. Note that this result is valid only for an open space and can be generalized to the case of the half-space [26,40]. Figure 4.1 shows that the anisotropy of the field remains at very large times and, therefore, that equipartition occurs asymptotically (t → ∞). On the contrary, the P to S energy ratio stabilizes at the theoretical value v 3 S /2v 3 P after a finite time. At this stabilization time, the diffusion approximation wrongly predicts both energy density and flux anisotropy. After about 100 s, radiative transfer and diffusion theory agree extremely well, but the residual anisotropy of intensity is still of the order of 2 and decays algebraically as predicted by equation (4.7). The total energy density is correctly predicted by the diffusion approximation only when the dipolar term describes with sufficient accuracy the anisotropy of the specific intensity. The calculations prove:

1. that the stabilization of the P to S energy ratio is a rapid phenomenon; It is important to notice here that formal equipartition would imply perfect isotropy (all modes, i.e. directions, equally represented). The stabilization of the P to S ratio occurs before the equipartition which is the asymptotic behavior of diffuse waves for large times. Therefore, one may expect large time asymmetries of the field-to-field correlation functions as shown by equation (4.7), provided the sources are located in the same distant region. The asymmetry is expected to disappear both in the average cross-correlations of late coda signals, or with an isotropic distribution of sources around the seismic network. In the latter case, each source produces an asymmetric correlation function, but the antisymmetric terms from many sources are expected to average out. In next section, we illustrate these points with a different numerical simulation scheme. This simulations will solve the real wavefield and will allow field-field correlations.

Numerical simulation of wave propagation and correlation

Another set of numerical simulations is conducted in a 2-D acoustic medium. This configuration is chosen because it is a simple way to describe wave propagation at the surface of the Earth. We solve the wave equation using a finite difference code [START_REF] Tanter | Application du retournement temporel à l'hyperthermie ultrasonore du cerveau[END_REF]110]. The field produced by each of several sources S is computed at each point in the medium. We consider a weakly scattering medium, where the distance of propagation is smaller than the transport mean free path * of the waves. The scattering is caused by a distribution of small rigid scatterers with radius a. The background velocity is 3.3 km/s. The product of the wavenumber k by the radius a, ka = 1. Following the time-reversal interpretation developed by Derode et al. [23] (and in Sect. 3, Chap. 3), we choose to place the sources S all around A (the reference point at the center of the grid marked with a cross in Figure 4.2a in order to form an equivalent of a perfect time reversal mirror. Each source S sends a broadband pulse with 0.1 Hz central frequency. The correlations are computed between the field G SA (t) at the reference point A and the field G SR (t) at any other point R(x, y) of the grid. The correlation is averaged over the entire set of sources S. The wave field reconstructed by correlations is displayed in Figure 4.2 for correlation times -30 s, 0 s and 30 s. Time τ = 0 is the central time of the correlations, when all the energy is focused at A as if A were a source. We observe a converging wavefront at negative times, and a diverging wavefront at positive times. These wavefronts correspond to the causal (positive times) and anticausal (negative times) parts of the GF between A and any point R in the medium. The nearly perfect reconstruction of the GF (including the converging and diverging wavefronts) is due to the quasi-ideal distribution of sources around A, the length of the coda (as long as allowed by the numerical schemes: 200 oscillations), and the absence of absorption. Let's now illustrate the effect of a non-homogeneous distribution of sources and finite durations of recording, in conditions closer to seismology. The conditions of computation are similar to those used to produce Figure 4.2 but we now consider a configuration which mimics a set of earthquakes along a fault. 40 sources S are aligned in the x-axis direction along a 400-km-long segment located 450 km away from the receivers as depicted in the sources S, corresponding to the anti-causal part of the GF. For τ > 0, only the causal part of the GF is reconstructed in the region opposite to the sources.

This numerical experiment confirms that the spatial distribution of the sources controls the time symmetry of the correlations. If the source distribution is asymmetric, the time-symmetry of the correlations can be broken. This is particularly true in weakly scattering (or homogeneous) media. This can be understood as well in terms of the time-reversal interpretation (see Sect. 4, Chap. 3, [94,103]). The uneven distribution of earthquakes in a limited region has a similar effect as a limited aperture of a time-reversal mirror. At long lapse times, the field becomes diffuse, and the argument given above (Eq. (4.7)) holds. We therefore expect that even with an inhomogeneous distribution of sources, scattering restores the broken time symmetry of the correlations. We checked this expectation with the same numerical experiment as before, but we correlated only the waves of the late coda. Late times in the coda are defined as τ > t * , where t * is the mean free time, t * = * /v. After t * , the waves have traveled more than * , and the field is evolving toward isotropy. The result is shown in Figure 4.4a and proves that the time symmetry of the correlations is restored by using large enough lapse times.

As a consequence of the discussion above on the isotropization of the field, we expect stronger scattering to improve the reconstruction of the GF even for an inhomogeneous distribution of sources. This effect was already confirmed by laboratory experiments [103]. We performed numerical simulations in a strongly scattering medium (Fig. 4.4) and processed the complete coda window. Again, the converging wavefront is clearly symmetric despite the uneven distribution of sources.

Besides the question of the source distribution, the distribution of scatterers can be a source of time asymmetry as well. In the case of an inhomogeneous distribution of scatterers, one may expect preferential directions of arrival of scattered waves. As discussed above, and whatever its origin, such an anisotropy of the energy flow results in a temporal asymmetry of the correlations.

Passive reconstruction of Rayleigh waves in an ultrasonic analogous crust

In this section, we present experimental results obtained in the lab with laserinduced and laser-detected ultrasonic waves propagating in a heterogeneous slab mimicking the Earth crust. Our aim is to illustrate the role of scattering in the time-symmetry of correlations, and investigate the role of mode coupling (bulkto-surface conversion) in the Rayleigh wave passive retrieval. In order to achieve small-scale seismology, we propose to investigate the emergence of the GF in the correlations of elastic waves propagating in a solid heterogeneous medium with a free surface. Because field experiments are tedious and natural environments are mostly unknown (especially in the scattering properties), we build an Earth crust model at the scale 1/10 -6 (presented in Fig. 4.5). Waves are sensed at ultrasonic frequencies at the free-surface of an elastic open medium. In our experiment surface waves are not initially excited. This is different from the work of Malcolm et al. [20], where ultrasonic Rayleigh waves were generated on the same surface they were measured. In addition they used a finite cylindrical medium with possibly round trip wave trains whereas our experiment is conducted in a nearly open medium. In our configuration, we would not expect the Rayleigh wave to be reconstructed in the correlation, except if scattering is present and mode conversion between Rayleigh and bulk waves occurs. To verify this assumption and study the role of mode conversion, we used two samples of identical dimensions, one with scatterers the other without.

Laser generation and detection of ultrasound

The laser devices have the advantage of being sharp, broad-band, and noncontacting (remote sources and detectors). In Figure 4.6 are presented the wavefronts propagating in an homogeneous medium in a source-receiver configuration. At time t after the source, the bulk wavefronts are localized along a sphere of radius v P t and v L t. On the free surface (z = 0), the perturbation is localized at x = ±v R t. Generation of ultrasound in solids with a pulsed laser has been subject to a great deal of publications (see for instance [START_REF] Mason | Physical Acoustics[END_REF][START_REF] Scruby | Laser ultrasonics -Techniques and Application[END_REF][START_REF] Royer | Génération et détection optiques d'ondes élastiques[END_REF]). Very schematically, the laser beam heats the surface of the elastic body. Depending on the intensity of the incoming beam, two regimes are possible (see Fig. 4.7). In the thermoelastic regime, the heat induces dilation of the focal spot just beneath the surface. This thermal expansion mainly creates lateral shear waves (essentially Rayleigh waves). If the laser intensity is increased (I > 15 mW/cm 2 for aluminum), the heat is strong enough to vaporize a thin layer of metal. The laser then mimics a vertical force. The laser detector used in the following is a heterodyne interferometer, and has been developed by Royer et al. [114]. It has the advantage of a very broadband response (20 kHz-45 MHz) and a sensitivity of 10 -4 Å/ √ Hz (for a complete description, see also [START_REF] Royer | Génération et détection optiques d'ondes élastiques[END_REF]).

Elastic waves in a homogeneous slab

Our first experiment takes place in a homogeneous medium: an aluminum slab without holes. Our aim is to mimic seismic waves propagating in the earth. We reproduce in this small-scale experiment the main features of the Earth crust. A duraluminum slab (25 mm thick) with dimensions roughly 10 -6 of those of the crust is used (Fig. 4.8). The source we employed to simulate earthquakes was a Q-switched Nd:YAG laser that shot 24 ns pulses at the bottom side of the slab (each pulse energy: 9 mJ). In the ablation regime, both Rayleigh and bulk waves are generated. Rigorously, such a source is not truly reproducible since the laser impact can damage the surface. To make the experiment as reproducible as possible while staying in the ablation regime, the shot intensity was no more than 280 MW/cm 2 . A 1 ms record was acquired 100 times without any observable change. In the following experiments, for a satisfactory signal to noise ratio, each impulse response was averaged over 100 consecutive shots.

The detection of the free surface motion is achieved with a non-contacting and quasi-punctual device: a heterodyne optical interferometer. It is mounted perpendicularly to the slab and then measures the absolute vertical component of the free surface displacement (top side), with a fine spatial resolution; the size of the laser spot is ∼100 μm whereas the typical elastic wavelengths here range between 1 and 10 mm. This is similar to seismology, where sensors are nearly punctual compared to the wavelengths considered (several kilometers at 1 Hz). However seismic sensors usually provide time records of the three components of the displacement field. Here the interferometer measures only the vertical movements of the free surface.

Ideally, the slab should have infinite dimensions along the directions x and y. To approach this condition, we attach a thick layer of dense plasticine on the lateral sides and edges of the duraluminum sample with the aim of creating absorbing boundary conditions and avoiding the generation of Rayleigh waves by mode conversion at the edges. The energy decay time τ abs was initially found to be 23 000 μs (see Fig. 4.9). With the plasticine it decreased to 120 μs. We therefore simulated an open slab in the x and y directions with a free surface at the top. Rigorously, the earth crust is a waveguide that partially leaks energy through the Moho discontinuity to the underlying mantle. To perfectly match this feature we should have placed a infinite medium with a different impedance at the bottom side of the aluminum slab. Yet we think that our results and conclusions do not suffer too much from this omission: the central point of our setup is to mimic an elastic medium, infinite in the horizontal directions with a free surface on the top. In an elastic body, three different kinds of wave polarization are possible. Compressional (or longitudinal) waves are analogous to acoustic waves in fluids (velocity v P = 6.32 mm/μs in duraluminum). Shear (transverse) waves have two possible polarizations (velocity v S = 3.13 mm/μs): one we call SV (vertical) in Ann. Phys. Fr. 31 • N o 3 • 2006 the x-z-plane (see Fig. 4.5) and one SH (horizontal) in the x-y-plane. SH waves have no contribution in the z-direction and, therefore, will not be detected by the interferometer. In addition to bulk waves, surface waves exist but here only Rayleigh waves (v R = 2.9 mm/μs) will be taken into account, since the others cannot be detected (no vertical displacement). The shortest wavelength in the aluminium slab is 0.9 mm, which is much greater than the duraluminum alloy grain size. Since the orientation of the grains is random, we consider the alloy to be isotropic for elastic waves in the frequency band of interest.

Elastic waves in a heterogeneous slab

The experimental setup depicted in Figure 4.10 is designed to mimic the scattering of elastic waves through the Earth crust, at ultrasonic frequencies (0.8-3.2 MHz) [22,115]. Fifty-four cylindrical holes (radius a = 2 mm) were drilled at random along direction y so that the waves propagating through the slab could undergo multiple scattering. The 2-D spatial Fourier transform of the hole positions was calculated and is almost perfectly flat in the range of ultrasonic wavelengths involved in the experiment, which confirms the absence of spatial correlation between the holes. The density of scatterers was n = 0.0105 mm -2 . Experimental setup: the propagation medium is an aluminum block drilled with 54 vertical cylindrical holes (diameter: 4 mm). A Q-switched Nd:YAG laser shoots on the bottom side 100 pulses (24 ns duration and 280 MW/cm 2 intensity each). On the top side a heterodyne interferometer senses the vertical displacement of the aluminum/air interface. This measure is repeated at 7 different locations (X 0 = 0 mm → X 6 = 60 mm) for each source position. Plasticine was stuck on the edges and the four lateral sides to mimic absorbing boundary conditions and avoid the generation of Rayleigh waves by mode conversion at the edges. The overall translational symmetry along the y-axis of both the free surface and the cylindrical scatterers avoid any coupling between SH mode and the other SV and P modes. Therefore SH waves will not be considered here and SV waves will be referred to as S (shear) waves. The wave propagation in our experiment will be treated as 2-dimensional (and quasi 1-D for the surface Rayleigh waves). Waves initially propagating in the y-direction are rapidly lost by absorption by the plasticine.

The laser source and the laser interferometer can be translated independently: 35 sources and 7 sensing positions were used during the experiment, resulting in a set of 35 × 7 impulse responses. A typical waveform is depicted in Figure 4.11. Around 1 MHz it lasts nearly 800 μs, and shows a long diffusive decay comparable to the seismic coda. Due to the strong scattering on the cylindrical cavities, no top-bottom reflection was observed in the data. This confirms the highly diffusive nature of the propagation in the scattering slab. The relevant scattering properties are discussed and evaluated in next section.

Wave scattering and transport properties

Scattering cross section of an empty cylinder

In order to evaluate the amount of scattering and mode conversion, we calculated the differential cross-section ∂σ/∂θ and the total scattering cross-section σ of a cylindrical void in a elastic medium excited by a compressional or shear plane wave. For a detailed derivation we refer to references [START_REF] Pao | Diffraction of Elastic Waves and Dynamic Stress Concnetration[END_REF][START_REF] Faran | [END_REF][118]. The differential scattering cross-section gives the angular distribution of the scattered surfacic intensity, normalized by the incident surfacic intensity. The total elastic crosssection is σ = ∂σ/∂θdθ. In 2-D it has the dimensions of length. It corresponds to the scattering strength of an object at a given frequency. In an elastic medium, mode conversion can occur and different cross-sections must be considered. In the case of an incident compressional wave, they are denoted as σ PP , σ PS and σ P = σ PP + σ PS , respectively for the P to P, P to S and total P elastic cross-sections. We also calculated the elastic cross-sections for an incident shear wave (S), σ SP , σ SS and σ S = σ SP + σ SS .

Here follows a brief calculation of the wave field scattered by a cylindrical cavity insonified by a plane compressional wave. Three displacement potentials are relevant: Φ i is the displacement potential of the incident P wave, Φ s is the scattered P wave potential and Ψ s is the scattered S wave potential. They can be expanded as:

Φ i (r, t) = ∞ m=0 ε m i m J m (k P r) cos(mθ)e -iωt , Φ s (r, t) = ∞ m=0 A m H (1) m (k P r) cos(mθ)e -iωt , Ψ s (r, t) = ∞ m=1 B m H (1) m (k S r) sin(mθ)e -iωt ,
where J m and H (1) m are respectively the Bessel and Hankel functions both of first kind and of order m, and ε m is the Neumann factor. Taking into account the null traction condition at the surface of the cylinder allows the calculation of the A m and B m coefficients. Those coefficients are given in [START_REF] Pao | Diffraction of Elastic Waves and Dynamic Stress Concnetration[END_REF][START_REF] Faran | [END_REF][118]. The scattering cross-sections are

σ P→P = 2 k P ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ 2|A 0 | 2 + ∞ m=1 |A m | 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , σ P→S = 2 k S ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∞ m=1 |B m | 2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , σ P = σ P→P + σ P→S .
The corresponding differential scattering cross-sections are

∂σ P→P ∂θ (θ) = 4 k P ∞ m=0 A m i -m cos(mθ) 2 , ∂σ P→S ∂θ (θ) = 4 k S ∞ m=1 B m i -m sin(mθ) 2 .
The differential cross-sections plotted in Figure 4.12a have been computed at 1.2 MHz and 2.4 MHz frequencies. The elastic scattering sections are plotted in Figure 4.12b for frequencies ranging from 0.1 MHz to 200 MHz. The average in the frequency band of interest (0.8-3.2 MHz) was σ P = 9.2 mm. This value is comparable to measurements by [119]. The same calculations were conducted for an incident shear wave. The average was σ S = 12 mm in the 0.8-3.2 MHz frequency band.

Transport properties

When the elastic wave propagates through the aluminum slab drilled with holes, it undergoes multiple scattering. Let ϕ(t) (respectively ϕ 0 (t)) be the vertical displacements sensed at the free surface through the scattering (respectively homogeneous) medium. Classically, this field is split into two contributions: the coherent and the incoherent part. The coherent wave is the ensemble-averaged field ϕ(t) (averaged over disorder configurations, here: cylinders positions). We emphasize the difference between the coherent and the ballistic wave (i.e. the first arrival). For a detailed discussion about scattering effects on coherent and ballistic waves, see [6]. Away from resonances, the coherent wave can be roughly thought of as an attenuated version of the direct wavefront ϕ 0 (t). When there is no intrinsic dissipation, the energy of the coherent wave decays with the slab thickness H as e -H/ , where is the elastic mean free path. Assuming a dilute set of scatterers, the elastic mean-free path is simply related to density of scatterers n and their elastic cross-section σ by

= 1 nσ •
From the theoretical scattering cross section calculated above, we find P = 10.3 mm. In order to measure the mean-free path experimentally, we used two aluminum slabs of exactly the same dimensions. The first served as a reference and provided measurements of ϕ 0 (t) for different source-sensors positions. The second one was drilled with holes. By translating the source-receiver device along the slab, we achieved something very similar to a configurational averaging and measured the energy of the coherent wave ϕ(t) 2 . Between 0.8 and 3.2 MHz, we obtained P = 9 ± 2.5 mm from these experiments.

Ann. Phys. Fr. 31

• N o 3 • 2006
The intensity of the incoherent part was also studied. The time evolution of the averaged incoherent intensity I(t) = ϕ(t) 2 is governed by another parameter: the transport mean free path * . In an elastic body, transport quantities have been theoretically defined by [108]:

* P = 1 n σ S -σ * SS + σ * PS (σ P -σ * PP )(σ S -σ * SS ) -σ * PS σ * SP (4.8) * S = 1 n σ P -σ * PP + σ * SP (σ P -σ * PP )(σ S -σ * SS ) -σ * PS σ * SP (4.9)
with σ * = (∂σ/∂θ) cos(θ)dθ. It was evaluated numerically: * P ≈ * S = 5.5 mm. In an experiment, this parameter is very hard to measure with a reasonable precision. The coherent backscattering effect [53,54,56,66] (also referred to as weak localization) does give a direct estimation of the transport mean free path. However, our experimental configuration did not allow this special measurement, since we could not place a laser sensor in the vicinity of the laser source. Yet we checked that the experimental intensity decay I(t) gives an order of magnitude for * that is consistent with the theoretical value.

For the sake of simplicity we propose a 2-D scalar wave model for I(t) [59], under the diffusion approximation. In an infinite slab of thickness H with perfect reflections on both sides, the averaged transmitted intensity reads:

I(X, t) = I 0 ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 1 2H √ πDt + ∞ n=1 (-1) n H √ πDt e -n 2 π 2 Dt H 2 ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭ e -X 2 4Dt -t τ abs
with τ abs the absorption time (taking into account the intrinsic absorption in the aluminum and the lateral leaking due to the plasticine) and D the diffusion constant. X is the lateral distance between source and receiver. This formula is obtained using a modal decomposition of the diffusion equation in the z direction. The intensity I(t) is a mix of compressional and shear waves, each mode traveling with its own parameters (velocity, * , diffusion constant) and interchanging their energy through scattering events. In our experiment * S and * P are of the same order. The diffusion constant was approximated by

D = 1 2(1 + 2v 2 P /v 2 S ) v P * P + 2v S v p v s 2 * S ≈ 40 mm 2 /μs.
This assumption is valid after a couple of mean free times, when the equipartition regime [26] has been reached. Equipartition means that the density of compressional and transverse modes equilibrates. Considering the specific velocities of each mode [START_REF] Weaver | [END_REF], we infer that 80% of the energy is transported by S waves, and only 19% by P waves (and an additional 1% by surface waves). Hence our best fit (thick line in Fig. 4.9) of the intensity decay in the coda gives τ abs = 120 μs ± 10% and * = 5-20 mm. The early part of the intensity is significantly sensitive to * , whereas the late part mainly depends on τ abs .

Two-point correlation of diffuse fields

In this section, we focus on the experimental reconstruction of the direct GF from passive correlations. The main idea is to correlate diffuse fields sensed at two different locations on the top side when a source generates bulk waves on the bottom. Since we record the vertical component of the surface displacements, the two-point correlations should simulate a vertical source at the surface, which mainly generates surface waves. Indeed, the experimental correlations we obtained reveal a wave packet that travels at the speed of a Rayleigh wave. We insist that our sources do not generate surface waves at the top side of the slab. Moreover if a surface wave happened to be generated anywhere, it would be completely absorbed by the plasticine. Under these conditions no Rayleigh wave should travel on the top surface, and no Rayleigh wave should be passively retrieved by correlations. Why then should passive imaging give rise to a Rayleigh wave train in our experiment?

We propose first to examine the role of scatterers in the emergence of the direct Rayleigh wavefront in the correlations. To that end we separately correlated coda records obtained through two different aluminum slabs: the first drilled with holes, the second without. Each impulse response lasted ≈ 800 μs before reaching the noise level (see Fig. 4.11). We emphasize that these record lengths are far from the Heisenberg time (break time) at which the modes of the aluminum block would be resolved (here T H ≈ 10 6 μs) and correlations would naturally converge to the GF. This modal approach is irrelevant to our experiment. The records were correlated and averaged over the 35 available sources. For the scattering slab, this reads:

C ij (τ) = 35 S=1 t=600 μs t=0 μs ϕ(S, X i , t)ϕ(S, X j , t + τ)dt,
where X i and X j are the sensors points (running from X 0 = 0 mm to X 6 = 60 mm along the array). And for the homogeneous slab:

C 0 ij (τ) = 35 S=1 t=600 μs t=0 μs ϕ 0 (S, X i , t)ϕ 0 (S, X j , t + τ)dt.
To enhance the signal-to-noise ratio, each correlation is time symmetrized (= C(+τ) + C(-τ) ) and normalized by its maximum. Results are displayed in Figures 4.13 and 4.14. A propagating wavefront (traveling at the Rayleigh wave velocity v R = 2.9 mm/μs) is clearly visible in the presence of scatterers, whereas it does not appear in the homogeneous slab. We also summed the 6 normalized propagating peaks after having delayed each signal according to the Rayleigh wave travel time. The summation is displayed in the inset in each figure. In the scattering slab, its amplitude nearly corresponds to the coherent addition of 6 pulses. In the homogeneous device, the amplitude of the summation is ≈2.5 (incoherent addition of 6 uncorrelated fields). We conclude on the necessity of mode coupling due to scattering for the Rayleigh wave to emerge from the passive correlations of diffuse fields generated by bulk waves sources. This is especially relevant for applications to seismology. Therefore, it appears once again [103] that the role of scattering is crucial in passive imaging. Firstly, because of multiple scattering, at late times the equipartition regime can be attained whatever the sources/receivers positions. Secondly, because of mode conversions due to the scatterers, a Rayleigh wave emerges from the passive correlations even though no Rayleigh wave was generated by the sources. Note that the bandwidth in the upper band record is a little wider than in the lower band. This was done to compensate for the coda shortening (<600 μs) so that the product TΔ f was kept constant.

We can go a little further and catch a glimpse of the time symmetry properties. We performed the correlations into two consecutive time windows: from 0 to 45 μs and from 45 μs to 600 μs, and did not time-symmetrize the correlations (Figs. In the first time-window, from 0 to 45 μs, correlations are asymmetric in time. This means that the causal part (τ > 0) and the acausal part (τ < 0) of the correlations are different (see left part of Figs. 4.15 and 4.16). In the causal part, a Rayleigh wavefront is clearly visible whereas noise is dominating the acausal part. This is due to the preferential direction of Rayleigh wave propagation (waves traveling from X 0 to X 3 in our experiment). There is a net flux of energy from X 0 to X 1 , X 2 and X 3 (distances 10, 20 and 30 mm in Figs. 4.15 and 4.16). This flux is due to the uneven distribution of sources in comparison to the receiver couples X 0 -X 1 , X 0 -X 2 and X 0 -X 3 : most of the sources are on the side of X 0 . At the early times of the coda (from 0 to 45 μs), the diffusion regime is not yet attained.

Later in the coda, from 45 μs to 600 μs, the wave field in the bulk of the aluminum slab is very likely to be equipartitioned. In the low frequency band (0.8-1.6 MHz), the time-symmetry of the correlation is indeed restored [START_REF] Van Tiggelen | [END_REF]120] (see right part of Fig. 4.15): Rayleigh waves travel in all directions. Nevertheless and surprisingly, the asymmetry persists in the high frequency band (from 1.6 MHz to 3.2 MHz, see right part of Fig. 4.16). To interpret this observation, we have carefully studied the location of the scattering sources around the array. In the high frequency band, the Rayleigh wavelength is ∼1.5 mm. The generation of Rayleigh waves by scattering necessarily occurs in the first half-wavelength beneath the free surface [START_REF] Maeda | Modeling of conversion scattering between rayleigh and body waves in a randomly inhomogeneous half space based on the born approximation[END_REF]. In our scattering slab, one hole was nearly showing on the surface (position X < 0), another one was 0.61 mm beneath (position X > 60 mm), the others being located much deeper. Rayleigh wave trains are mainly generated by the hole nearest to the surface, then propagate along the array of receivers (from X 0 to X 6 point). These waves are almost unperturbed (attenuated) until they reach the edges of the slab and the absorbing plasticine. They contribute to a very clear propagating pulse in the positive part of the correlation. The weak coupling due to the deeper hole (z = -0.61 mm) on the X 6 side contributes to a smaller pulse propagating from X 6 to X 0 in the negative part of the correlations.

This interpretation is in agreement with observations in the low frequency band (0.8-1.6 MHz), where the average Rayleigh wavelength is 3 mm (Fig. 4.15). At least 5 holes are present in the first half wavelength and should cause significant scattering of Rayleigh waves and coupling between surface and bulk waves. This time, the holes are evenly distributed along the sensor array. In the late coda, correlations X 0 × X 1 , X 0 × X 2 and X 0 × X 3 are nearly symmetric. The time symmetry is obtained due to scattering by a symmetric distribution of scatterers. Under these conditions, a global equipartition among bulk and surface waves is guaranteed. Furthermore the reconstructed surface wave is strongly attenuated along its path because it senses many scatterrers along the array. In addition, these scatterers contribute to the signals around τ = 0 in the correlations, which degrade the reconstruction. We think this interpretation explains why the symmetric wavefront is much more noisy at low frequencies (Fig. 4.13), where a lot of cavities are encountered in the first half wavelength, than at higher frequencies where the Rayleigh wavefront can propagate freely (Fig. 4.15).

Finally, we comment on the possible misidentification of the waves that are reconstructed in our experiment. Indeed, if many scatterers are present within a wavelength, the overall wave velocity may be different (effective medium). In the heterogeneous plate, a reduced-speed shear wave might propagate with the same wavespeed as a Rayleigh wave in the bulk aluminium plate (i.e., without cavities). In our experiment, the hole interspacing is 10 mm on average, which is larger than the largest ultrasonic wavelength. Thus, it is reasonable to assume that the shear waves propagate in the aluminium plate with the same wavespeed as in the bulk. The measured wave velocity is 2.9 mm/μs at all frequencies (from 0.8 MHz to 3.2 MHz) and indeed corresponds to the velocity of the Rayleigh wave.

To conclude, the cross-correlation of diffuse fields was performed, allowing us to retrieve passively the Rayleigh wave between two sensors only when scattering was present. Without scatterers, and in the case of bulk wave generation and surface detection, no Rayleigh wave was reconstructed. This illustrates the role of scattering and mode conversion in the GF passive reconstruction. Analysis for different time-windows and frequency bands confirms observations in acoustics [106] and seismology [92] that the Rayleigh wave reconstruction was harder with increasing frequency (see also Chap. 3 on the convergence rate). Here we found that the Rayleigh wave reconstruction was more efficient in the high frequency band (Fig. 4.17). In addition we observed that even in the late coda where waves are expected to be equipartitioned, asymmetry in the correlations may remain. Both observations are due to the very specific coupling between bulk and Rayleigh waves, which occurs if scatterers are present in the first wavelength beneath the free surface. We emphasize that equipartition of bulk waves does not always mean equipartition of surface waves. On the one hand, our experiments show the need of scattering to passively retrieve the impulse response between two sensors, on the other hand they show that scattering occurring between the sensors degrade the reconstructed Rayleigh waves.

Correlation of seismic coda in Alaska

In seismology, two kinds of fields are usually considered as random: the seismic noise and the scattered waves of the coda. Seismic noise has the advantage to be easy to record and to exist even in regions without earthquakes. Shapiro et al. [START_REF] Shapiro | [END_REF] have shown that coherent broadband dispersive wave trains emerge from the cross-correlation of ambient seismic noise records between stations many hundred kilometers apart and that dispersion curves can be measured in the absence of earthquake records. However, the question of the origin of seismic noise is still open. Without a clear understanding of the source of the seismic noise, the relevant propagation regime (ballistic, diffusive, etc.) remains unknown. It is thus difficult to better understand the properties of the emergence of the GF from diffuse wave fields using only noise records.

In present section [120], we concentrate on coda waves since they are produced by a known source, with precise location and origin time, and they have been shown to result from multiple scattering in the Earth lithosphere. Campillo and Paul [92] showed that extracting the GF from field-to-field correlation of scattered waves is a valid approach not only in the controlled and favorable conditions of the laboratory, but also with natural signals such as actual seismograms produced by earthquakes. They used coda waves produced by moderate earthquakes in Mexico. The part played by multiple scattering and diffusion in the coda of the seismograms has been a subject of discussion since the pioneering papers of [14,15]. The importance of multiple scattering for the seismic coda is attested to by the success of the use of radiative transfer theory to explain the observation of energy decay [17,44,123,124]. Furthermore, the diffuse character of coda waves has been demonstrated by the observation of the stabilization of different energy ratios [26,41,[START_REF] Campillo | Waves in complex media[END_REF], a phenomenon associated with the principle of equipartition in the phase space of a random field (e.g. [START_REF] Weaver | [END_REF]). Equipartition means that all the possible modes are excited at the same level of energy.

In present context, the word modes refers to all body waves or normal modes of surface waves that are potentially excited in the region probed by the diffuse field and can be coupled by local scattering. Equipartition is a property of diffuse elastic waves. It implies a stabilization of the ratio of S to P wave energy with time, independently of the source. It is an important point in the context of the present paper, since it indicates that the modes which make up the GF are represented in the diffuse field. In the following, we first summarize theoretical and experimental arguments that clarify the conditions for the emergence of the GF from field correlations. We then present an application with data from a temporary broadband network in Alaska. Finally, we discuss the role played by the source distribution and the time evolution of the wave-field structure toward complete randomness.

Data processing

Temporary networks of seismic stations installed in regions with a high level of seismic activity provide useful data sets to study the properties of the emergence Ann. Phys. Fr. 31 • N o 3 • 2006 of GF from the correlations of coda waves. They include numerous stations with identical instrument characteristics that make it possible to compute crosscorrelations between many station couples separated by a large range of distances and with different orientations. Time-distance seismic sections can be constructed from the correlation signals where the propagating waves can be easily identified even with a poor signal-to-noise ratio. Actually, only limited averaging can be performed with coda records, and the deterministic signal is mixed with the remnant random fluctuations of the diffuse field.

We present here an application of the extraction of GF from coda waves to the data of a temporary network in Alaska. During the Broadband Experiment Across the Alaska Range (BEAAR), a network of 36 broadband seismographs was operated during 2.5 years ending in August 2001 [126]. Most stations recorded simultaneously about 100 regional earthquakes with magnitude larger than 3.3. Locations of stations and epicenters are shown in Figure 4.18a. Note that the distribution of earthquakes is not even and that they are mostly concentrated to the southeast of the network. The broadband seismograms were first band-passfiltered. For each couple of stations, the horizontal components of the seismograms were rotated assuming the inter-station great circle path to be the radial direction. We used coda records starting 20 s after the arrival of the S wave and ending when the signal-to-noise ratio was smaller than 4. An example of a record is shown in Figure 4.18b. Because of the exponential decrease of coda amplitude with time, we cannot perform a simple cross-correlation between the coda signals recorded at the two stations without strongly overweighting the earliest part of the coda. To avoid this problem, we followed [127] by disregarding the amplitudes completely and considering only one-bit signals. Campillo and Paul [92] checked that this procedure leads to the same results as compensating for amplitude attenuation with time in successive time windows. For each couple of stations, the cross-correlation of one-bit signals was computed for each earthquake, normalized by the maximum amplitude, and averaged over the entire set of events. This processing was performed for all combinations of components, such as vertical to vertical (Z/Z), vertical to transverse (Z/T), radial to radial (R/R), etc. The results of [92] suggest that these different field-to-field correlations contain the different terms of the elastic Green's tensor. Wapenaar [99] demonstrated the retrieval of the elastic Green's tensor from surface displacement fields produced by a distribution of sources on a closed surface in the medium. Even if the source distribution does not fulfill such a condition in actual experiments, he gives a firm theoretical argument for the reconstruction of the polarized elastic response.

Figure 4.19 shows the results of the processing in the frequency band 0.08-0.3 Hz. The nine polarization permutations correspond to the terms of the Green's tensor, that is for example the R/Z correlation corresponds to the vertical displacement produced by a force in the radial direction. Since we disregard the amplitudes of the coda waves by using one-bit signals, the relative amplitude of the reconstructed signals between the components is lost. Traces with a maximum amplitude at negative time have been time-reversed to concentrate the large-amplitude pulses at positive times and have a better view of the presence, or the absence, of symmetry in time. Note that the location of the maximum amplitude pulse at positive or negative times depends only on the order of the signals in the cross-correlation, which is arbitrary. Clear propagating wave trains can be observed on the Z/Z, Z/R, R/Z, R/R, and T/T components. The Z/T, R/T, T/Z, and T/R components contain only noise, as expected from the symmetries of the Green's tensor. This figure confirms the conclusions of [92] for numerous paths with different azimuths.

When a propagating branch can be seen at positive times, a symmetric one is more or less clearly visible, depending on the component, at negative times. This symmetric wave train is particularly clear on the Z/Z and T/T components at short offsets. We discuss the question of the time symmetry in more detail in next section.

To test the reliability of the extraction of the GF from the data, we computed synthetic GF in a 1-D velocity model derived from the seismic profiles recorded by [128] in the neighbor Tanana terrane. The computations are performed for a vertical point source acting at the free surface. The receivers are also at the free surface in a configuration similar to our station-to-station measurements. The synthetic section in the frequency band 0.08-0.3 Hz is shown in Figure 4.20 in the same distance range as in Figure 4.19. It is dominated by low-frequency surface waves (Rayleigh waves in this case), as the cross-correlation records. As expected for a point force source at the free surface, the body wave contributions to the GFs are negligible. We performed the same computation in a frequency band centered on 1 Hz (Fig. We are able to reconstruct only the surface waves because the GF is expected to be dominated by surface waves. A high level of noise remains after the limited averaging performed. It prevents an unambiguous identification of body waves. The reconstruction of the high frequencies is difficult because absorption limits the durations of the available coda records. These limitations could be overcome with a larger data set.

Time symmetry and isotropy of diffuse wave fields

In Figure 4.19, we arbitrarily chose the direction of time so that the maximum amplitude of the average cross-correlation is at a positive time. Swapping the two stations in the cross-correlations gives a time series which is exactly the timereversed of the original by definition of cross-correlation. Since the identification of the waves relies only upon the variations of the travel time with distance, this choice has no consequence. On the other hand, Figure 4.19 exhibits a strong asymmetry of the cross-correlation time functions. One must remember that theoretically we expect to see both causal (retarded) and anti-causal (advanced) parts of the GF. Nevertheless we observe that for most station pairs, a single direction of time is favored. A similar observation was made by [92] in Mexico. They used earthquakes along the subduction and found that the cross-correlation functions are asymmetrical for stations inland with an orientation perpendicular to the coast and symmetrical for stations located along the coast within the earthquake source region. They suggested that time asymmetry indicates a preferential direction of flow of the energy of coda waves, a property that could appear as paradoxical in the context of random fields. In the following, we will discuss this issue in relation with multiple scattering and the uneven distribution of epicenters which are mostly located to the southeast of the network in Alaska.

Consider the Z/Z correlation profile of Figure 4.19. Using different bandpass filters, we attempted to see whether or not the level of asymmetry varies with frequency. However, the signal-to-noise ratio does not allow any convincing visual comparison. To achieve a more quantitative analysis we performed a slantstack of the cross-correlations. For each trace, the arrival time of the Rayleigh wave was measured in the positive times, where, by construction, the largest amplitude is found. We then shifted the traces of the time of the maximum and stacked all the arrivals in the positive times. We stacked the signals in the negative times in the same way using time shifts opposite to the ones measured in the positive times. This operation enhances the coherent arrivals and takes into account the variations of wave velocity among the different paths. Finally we computed the ratio between the maximum amplitude of the stacked Rayleigh pulses in the positive time to the maximum amplitude of the stack in the negative time. The stacks were performed for all traces corresponding to interstation distances between 20 and 120 km which have the best signal-to-noise ratio. This processing was applied to the sections of the Z/Z correlation profile of Figure 4.19 after bandpass filtering. We obtained an amplitude ratio of 3.8 in the band 0.12-0.3 Hz and 4.0 in the band 0.08-0.15 Hz indicating that the symmetry is slightly stronger for high frequency.

However, we also expect the late coda to behave differently from the early coda. Due to anelastic absorption, the early coda includes more high frequencies. As a consequence, the influences of both lapse time and frequency on the symmetry are mixed in this analysis when we consider cross-correlations computed on the whole coda length.

To clarify this point, we studied the dependence of the causal to anticausal amplitude ratio with the lapse time in the coda window. We formulate as a first order hypothesis that the larger the lapse time, the more isotropic the coda. According to this argument, we expect the correlations computed in later time windows to be more symmetric in time since the corresponding diffuse wave field is more isotropic, and all directions of propagation are closer to being equally represented. To verify this effect, we compared the amplitude ratio of positive and negative times for correlations computed for different lapse times and the same two frequency bands as before. We considered first the correlations computed from the first 300 s of coda (early coda) filtered in the low frequency band (0.08-0.15 Hz). We found a ratio of stacked amplitudes of 5.6. We repeated the measurement for correlations computed from record windows starting 300 s after the beginning of the coda (late coda). The amplitude ratio between the causal and anticausal signals is 3.9. This indicates that the correlations are more symmetrical when measured from the late coda than from the early coda. In other words, the longer time the scattered waves propagate, the more isotropic they are.

We performed the same analysis in the band 0.12-0.3 Hz. For the early coda the amplitude ratio is 4.3 while it is 3.1 in the late coda. This confirms that the pulses emerging from the cross-correlations of the late coda are more symmetric in time than those computed from the early coda. Comparing ratios obtained in the two frequency bands, we find that the values are smaller for the band 0.12-0.3 Hz than for 0.08-0.15 Hz. This indicates that the high frequency waves evolve faster toward isotropy. In the following, we investigate theoretically and from numerical simulations what is the expected evolution of the net flux of energy with time along a seismogram and if this evolution accounts for the observations.

To conclude on our experiments on real earthquake records in Alaska, the emergence of the GF is clearer at low frequency (0.1-0.3 Hz) than at higher frequency. However, the correlations exhibit a clear time asymmetry (Fig. 4.19). This asymmetry proves that the forward flux dominates and that the excitation of the propagation modes that make up the diffuse field is anisotropic. Indeed the earthquakes we used are concentrated in the southeast of the center of the network as shown by Figure 4.18a. Due to absorption and noise, the coda records are not long enough to reach the isotropy of the diffuse field. The resulting average flux of energy is responsible for the observed time asymmetry of the cross-correlation. The relative amplitude of causal and anti-causal waves, and its evolution with time, can be used to measure scattering properties of the medium such as transport mean free path. the wave fields Φ A and Φ B received at A and B will be e(t) ⊗ G AS (t) and e(t) ⊗ G BS (t), ⊗ representing convolution. The cross correlation C AB of the fields received at A and B is then

C AB (τ) = Φ A (t)Φ B (t + τ)dt (5.1) = G AS (-t) ⊗ G BS (t) ⊗ f (t) (5.2)
with f (t) = e(t) ⊗ e(-t). It was argued in previous sections that this two-point correlation C AB of a perfect diffuse field, averaged over time or averaged over a series of sources perfectly surrounding the medium, would yield the GF between the passive receivers [23]. Such a perfect distribution of sources surrounding sensors A and B has been tested in numerical simulations [23] where A and B were deep inside the scattering medium, and with ultrasound (see Sect. 4, Chap. 3 and Ref. [129]).

In practice, whatever the type of waves involved, expecting a perfect distribution of sources is unrealistic for at least two reasons: first the limited number of sources, second their uneven distribution. In seismology, for instance, seismic stations (A, B) record the displacement field at the Earth's surface, but the sources (S, the earthquakes) are mostly aligned along faults.

Experimental setup

In order to test and illustrate the possibility of imaging based on the correlation of coda waves with a limited number of sources, we have designed laboratory experiments (Fig. 5.

1).

The experiments take place in a water tank. We use a 118-element ultrasonic array to simulate 118 earthquakes. For each run of the experiment, one of the elements sends a short pulse (3 MHz center frequency) that traverses a highly scattering medium. This medium consists of a random arrangement of vertical steel rods (average density 29.5 rods/cm 2 , diameter 0.8 mm). The sample's mean free path is * = 3.5 mm (this was estimated by the coherent backscattering effect [56]), whereas its smallest dimension is 15 mm; therefore, the waves undergo many scatterings before they emerge from the scattering slab and reach the receiver.

Behind the slab, we place the medium that we want to image. It consists of four liquid layers (alcohol, oil, water, sugar syrup) with different sound speeds. The receiver is a 0.39 mm piezoelectric transducer. It is translated downwards along the z-axis, and records the scattered signals that are generated each time one of the elements on the array fires a pulse. Those records are equivalent to seismograms, except they are made only of compressional waves. A typical signal is plotted in Figure 5.2a. It lasts more than 300 μs, i.e., 300 times the initial pulse, and shows a high degree of multiple scattering. It should also be noted that the receiver is facing downwards; therefore, it cannot record direct waves coming from the sources, but purely multiply scattered contributions. Moreover, following the Van Cittert-Zernike theorem and given the typical distances involved here, the spatial coherence of the field measured by the receiver is ≈ λ. No coherent wavefronts are propagating between the receiver's positions. The 118 sources fire successively, and 118 "seismograms" are generated and recorded for each position z of the receiver. The scattered waveforms are G Sz (t), with S the index of the earthquake and z the position of the receiver.

Normal correlations

Next, we choose a pair of receiver positions (z 1 , z 2 ), and we cross correlate the field due to the source S C S (z 1 , z 2 , t) = G Sz 1 (t) ⊗ G Sz2 (-t).

(5.3)

A typical result is shown in Figure 5.2b for S = 60, z 1 = 36 and z 2 = 16 mm. Nothing seems to emerge from the correlation. Then, we repeat this for the 118 sources and sum the correlations to get

C(z 1 , z 2 , t) = 118 S=1 C s (z 1 , z 2 , t). ( 5

.4)

A typical result is plotted in Figure 5.2c. This time, a strong peak emerges from the correlation at time t = -13.6 μs, corresponding to the time of flight |z 1z 2 |/c. This is the signature of the direct GF between z 1 and z 2 . In this example, the distance between the receivers is 20 mm, so from the arrival time t = -13.6 μs we get an estimation of the sound velocity between z 1 and z 2 : 1.47 mm/μs. The peak appears only at negative times and not at positive times due to the location of the receiver, i.e. the scatterers and its limited directivity [103]: the receiver sees only what is coming toward its front face.

Here, the sound velocity between z 1 and z 2 could be estimated only because the peak was sufficiently above the surrounding noise. Actually, this term is improper since it is not noise per se, but fluctuations of the correlation between the waveforms outside the peak. In order to see the peak emerge from the signals presented in Figure 5.2, we had to average the cross correlations on a large number of sources (up to 118). The evolution of the signal-to-noise ratio (SNR) in the averaged cross correlation versus the number of sources N is presented in Figure 5.3. The SNR was calculated as follows. The signal level is the value of the envelope of the peak, and the noise level is the mean amplitude of the envelope of C(z 1 , z 2 , t) outside the peak, in a 25 μs-long time window. As the number of sources N is increased, the peak emerges more and more. The SNR seems to grow like N 0.4 , whereas a classical N 0.5 dependence would have been expected. This indicates that the contributions from different sources are not fully decorrelated (indeed, the array pitch is slightly smaller than the wavelength) and/or do not contribute equally to the signal, but the range of N is too small to conclude that the exponent 0.4 is really meaningful.

Even though the sources were not arranged as a perfect time-reversal device, the experimental results show that the main feature of the direct GF can be retrieved from the correlations. This is achieved due to multiple scattering in the random sample. Once again we can refer to the TR interpretation. As we have argued, in a time-reversal experiment, there is a forward step (propagation from A to B, record the field in S) and a backward step (time reverse and send back the field in S, observation of the resulting field at point B). If the time reversal were perfect, then the backward step would be identical to the forward step. Similarly, the key question to the retrieval of the GF is: if this were a TR experiment, would the backward step be identical to the forward step? Although this is almost never the case, except in a thought experiment, past experiments have already shown that time-reversal focusing is more efficient (i.e., the backward step and the forward step are more alike) in the presence of strong multiple scattering or reverberation [6,87,130,131]: the focused peak is narrower in space, indicating that the angular spectrum of the wave field is precisely recovered. Strong multiple scattering or reverberation virtually enhance the size of a TR device, i.e., the number of sources involved.

1-bit correlation

Interestingly, this technique also works with one-bit correlations: instead of recording the entire waveforms G Sz (t), we record and cross correlate only their sign, i.e. a series of +1/ -1 as shown in Figure 5.4. And, one-bit correlations seem to give even better results than "normal" correlations (see Figs. 5.2b and 5.4b for a comparison). Again, it has already shown [127,132] that one-bit time reversal in a multiple scattering or reverberating medium gives a higher SNR than a classical time reversal, since it gives more importance to the longest scattering paths, artificially reinforcing multiple scattering. The benefit of one-bit correlation will be effective only if the recording time is significantly larger than the decay time of the multiply scattered signals. For a multiple scattering slab with thickness L, the typical decay time is the Thouless time L 2 /D, with D the diffusion constant (here, D ≈ 2.6 mm 2 /μs). As long as the recording time is larger than the Thouless time, the effect of one-bit digitization is to reinforce the weight of the longest and most diffracted scattering paths relative to early arrivals. Experimental results show the interest of one-bit correlations: in Figure 5.3, the SNR is systematically higher by ≈4 dB with one-bit correlation; therefore, a smaller number of sources can be employed which is interesting for any practical implementation of the method.

It should also be noted that the retrieval of the GF takes advantage of the large frequency bandwidth available in pulsed ultrasound (here, the spectrum spans from 2 to 4 MHz). Suppose we had only one source, working in a narrow frequency band: then, the retrieval of the GF would completely fail. Once again, the time-reversal interpretation is enlightening. Indeed, in a one-channel time-reversal experiment performed in a multiple scattering medium on a single realization of disorder, focusing cannot be achieved if the frequency band is too narrow [6]: the reemission of the phase-conjugated monochromatic wave just creates a speckle pattern that is not focused back at the source. However, as the frequency bandwidth Δω is progressively enlarged, it has been shown that a TR device manages to refocus the wave through the multiple scattering slab, even with only one source [6,87,130]. The underlying idea is that we take advantage of frequency averaging as soon as the bandwidth Δω is larger than the correlation frequency δω of the scattering medium. In a homogeneous and lossless medium, δω Δω. But, in a multiple scattering slab, the correlation frequency Δω is inversely proportional to the Thouless time (Δω is also often referred to as the Thouless frequency). Since there are roughly Δω/δω decorrelated frequencies available in the spectrum, the SNR can be expected to rise as √ Δω/δω, if the power spectral density is flat. Hence, using a large frequency bandwidth is another way of increasing the SNR. This is illustrated in However, it should be noted that enlarging the frequency band is not a substitute for source averaging. If we want to retrieve all the details of the exact GF, the only solution is to have sources surrounding the medium. But, if we are satisfied with a simple estimate of the first arrival of the GF, then enlarging the frequency band helps because it increases the peak-to-noise ratio, at least as long as dispersion in the medium to be imaged between z 1 and z 2 can be neglected. A different approach consists of extending the bandwidth to infinity (Δω δω) and invoking a self-averaging property to retrieve an ensemble-averaged GF with a single source [START_REF] Van Tiggelen | [END_REF]. Naturally, this holds only if the ensemble-averaged GF itself does not depend on frequency (i.e., the scattered wave field expressed as a function of frequency has to be a stationary, or even ergodic, random process) so that an ensemble average can be replaced by a frequency average.

As a preliminary conclusion, the time-reversal interpretation indicates that in order to retrieve the exact GF from the correlations, the sources of the earthquakes should be placed in such a way that they completely surround the medium and the sensors. But, the time-reversal interpretation (based on previous works on Ann. Phys. Fr. 31 • N o 3 • 2006

Passive imaging of a dynamic medium

The same experimental procedure was also applied to a two-layer medium (oil/sugar syrup). Initially the liquids are at rest, and the image of velocity profile clearly shows the two layers (Fig. 5.7). Then, the medium is scrambled to form an emulsion: the velocity profile we obtain shows an apparently homogeneous medium with a sound speed of 1.57 mm/μs. The experiment is repeated on the same sample while each liquid progressively separates from the other. After 12 hours the separation is complete. Thus, the images obtained from correlation of the scattered fields were able to monitor the evolution of a medium undergoing a slow structural change. This amounts to "passive diffusive wave spectroscopy", in interpretation with diffusive wave spectroscopy (DWS) or coda wave interferometry (CWI) [19,133]. However, in the situation depicted here, the change in the scattered wave fields is not due to a movement of the scatterers, but to a change in the medium we image through the scatterers. Sound-speed profiles obtained as the medium under investigation changes. Initially, there are two wellseparated layers (oil/syrup). Then, they are mixed together to form an apparently homogeneous emulsion. Progressively, the two phases of the emulsion separate again, and the process can be monitored by the sound-speed profiles. Twelve hours later, the separation is complete. To conclude this section, the TR interpretation tells us that when the TR criterion (perfect sources distribution) is only partially fulfilled, multiple scattering helps achieve a correct estimate of the arrival time of the GF. It was also shown that one-bit correlations (i.e. with processing only the sign of the wave field) can give similar estimations with fewer sources, because one-bit processing tends to give more weight to the longest and most diffracted paths compared to normal processing, which enhances the role of multiple scattering. In next section, we will aim at imaging localized reflectors.

Passive imaging of interfaces in open scattering media

The previous experiments showed the possibility of reconstructing the direct (ballistic) wave between two passive sensors from the time-correlations of acoustic We now aim at imaging an aluminum/water interface with only a few sources. The experimental setup (Fig. 5.8) is comparable to the one in Section 1 [94]. We use N = 118 aligned sources S to illuminate a scattering medium made of vertical steel rods randomly distributed. The interface to be imaged is behind this complex structure. A typical example of such a coda waveform is shown in Figure 5.8. Again, the whole setup is 2-D.

As in the previous experiment, a conventional ultrasonic transducer (size ∼ λ) is used. Unfortunately the transducer is not omni-directional. This is compensated for by rotating the sensor at least twice at each position. Then the receiver is translated along the x-axis. The emission/reception sequence is repeated for 12 different locations i. This provides a set of 118 × 12 records for each sensor angle. For each pair (z i , z j ), we compute the time-and sourceaveraged correlation:

C ij (τ) = 1 N N S=1 G Sz i (t)G Sz j (t + τ)dt. (5.6)
The resulting traces are displayed in Figure 5.9 where both the direct and the reflected hyperbolic-shaped wavefronts are visible. Note that one single correlation C S ij (τ) does not contain any clear arrival: a source averaging is needed. The measurement of the arrival times along the sensor array provides a precise measurement of the interface position: it was parallel to the array and 33.2 mm from the x-axis. The estimated error on the reconstructed travel times is of the order of one record sample: 0.05 μs. The signal-to-noise ratio of the reconstructed waveforms roughly corresponds to Weaver's theoretical prediction [106] These results show that it is possible to retrieve the details of the GF of an open medium beyond the first (ballistic) arrival. Note that it differs from [START_REF] Bakulin | Virtual source: new method for imaging and 4d below complex overburden[END_REF] in that the whole coda (high-order multiple scattering contributions) is processed and taken advantage of; the signal-to-noise ratio is noticeably improved. These results have been obtained using distant sources, but can be generalized to the case of purely random noises in any field of wave physics. In seismology, these observations suggest that the application range of Passive Imaging would extend to seismic prospecting or imaging of deep Earth.

Application to seismic noise

Passive tomography of California using ambient seismic noise

To image the Earth interior, the traditional technique is to deploy dense arrays of receivers and record the seismic wavefield emanated from distant earthquakes. Nevertheless this traditional observational method cannot fully exploit emerging array data because they are based on seismic waves coming from select source regions predominantly near plate boundaries and are observed at stations far from the source regions, such as most locations within the United States. An original idea to cross-correlate noise to passively image the medium has been applied to ambient seismic noise [90,[START_REF] Shapiro | [END_REF] at different places in North America. Two teams, one lead by Shapiro [93], the other by Sabra [START_REF] Sabra | [END_REF], took advantage of the dense array of seismic stations (see map, Fig. 5.10a) available in California to reconstruct a series of Rayleigh wave paths and perform a high resolution tomography of the subsurface.

The seismic noise is not as fully diffused as thermal ultrasonic noise [80], but the distribution of the ambient sources randomizes when averaged over long times. Ambient seismic noise is additionally randomized by scattering from heterogeneities within Earth [26]. Surface waves are most easily extracted from the ambient noise because they dominate the GF between receivers located at the surface and also because ambient seismic noise is excited preferentially by superficial sources, such as oceanic microseisms and atmospheric disturbances [13]. The seismic noise field is often not perfectly isotropic and may be dominated by waves arriving from a few principal directions. To reduce the contribution of the most energetic arrivals, the first team disregarded the amplitude by correlating only one-bit signals [92,94] before the computation of the cross-correlation, whereas the other clipped (saturated) the highest amplitudes (to lower their statistical weight).

The number of total stations and the exact processing technique differs from one team to another. Nevertheless, both teams found similar maps of the subsurface, confirming the strong robustness of estimating the Rayleigh wave arrival from correlation of ambient noise. The authors produced detailed images of the crystal structure that delineated sedimentary basins from igneous complexes, and even fault lines that offset different rock types. One example is presented in Figure 5.10b for Rayleigh waves filtered around 7.5 s.

Lunar subsurface investigated from correlation of seismic noise

Passive imaging (Fig. 5.11) has been achieved in seismology using seismic coda [92,120] as well as seismic noise [93,[START_REF] Shapiro | [END_REF][START_REF] Sabra | [END_REF] for frequencies ranging from 0.025 to 0.2 Hz. On Earth, these studies have shown that a statistical treatment of seismic noise yield the GF when averaged over sufficiently long time. Microseismic background noise is mainly excited by surface sources like oceanic and atmospheric perturbations [13]. Our aim now is to take the passive imaging method beyond the Earth and apply it to extraterrestrial planets on which neither an ocean nor an atmosphere exists. In particular, our aim is to investigate the method in the case where the microseismic background noise is made of fully diffusive waves, at frequencies higher than previously investigated [93]. 30-minute passive listening periods following the missions, with the weekly listening modes terminated on December 7, 1973, thereafter only monthly operational checks were performed (ALSEP status report, 1976). The main objective was to infer the velocity structure of the uppermost part of the crust down to a few km depth using traditional seismic refraction techniques [141].

Apollo 17 touched down on the floor of the Taurus-Littrow Valley near the southeastern rim of the Serenitatis Basin, which consists of irregular, heavily cratered regolith developed on lava flows partly filling an embayment between massifs 2 km high [START_REF] Heiken | Lunar Sourcebook -Users Guide to the Moon[END_REF]. The Apollo 17 Lunar Seismic Profiling Experiment (LSPE) was deployed on December 14, 1972, at a distance of about 180 m WNW of the lunar module. Four geophones, with natural resonant frequencies of 7.5 Hz, labeled G1 to G4, were deployed in a triangular array (see Fig. 5.12) and recorded the vertical ground velocity in the frequency range from 3 to 30 Hz during the landing mission in 1972 [START_REF] Kovach | Lunar seismic profiling experiment[END_REF] and then again from August 1976 to April 1977. The geophones were simultaneously connected to the central station where the seismic signals S 1...4 (t) were sampled at 118 Hz and digitized before being telemetered to the Earth. Prior to conversion into a 7-bit digital format (counts 0 to 127), the signals were conditioned using a logarithmic compression, which was done to increase the dynamic range. A typical four-channel record is displayed in Figure 5.12. We processed continuous data from August 15, 1976, to April 24, 1977. They are time-windowed into 2016 record samples, each lasting about T = 3 hours for a given record d. The correlated traces are calculated directly from the raw data without applying any decompression or correction filter. The absence of decompression tends to equalize the noise amplitude. For a given pair of geophones i, j, and the record d, a cross-correlated trace is determined as Why should this pulse emerge from the correlations? Seismic waves propagating in the direction j → i add up coherently in the correlation, and contribute to this well-defined pulse for times τ > 0 (causal part of the correlations). Seismic waves propagating in the opposite direction (i → j) contribute to the well defined pulse for times τ < 0 (acausal part of the correlations). Waves propagating in other directions add up incoherently and contribute to the residual fluctuations in the correlations [75].

C d ij (τ) =
The reconstructed pulses are interpreted as Rayleigh waves between geophones i and j, i.e. the ground velocity response of the subsurface at j (i) to a vertical impulsive force at i ( j). They are found to propagate with an average velocity of ∼50 m/s, and are clearly dispersive, as shown in Figure 5.13a. The group velocity is defined as the envelope arrival time of the Rayleigh pulse. This pulse is filtered around 14 different frequencies, ranging from 3.6 Hz to 11.4 Hz (bandwidth 33%). Velocities were calculated when the signal-to-noise ratio (SNR) of the correlated traces was >2. When possible we averaged velocities over the three central pairs of geophones. We plotted the dispersion curve in Figure 5.13b, which is used to image the subsurface velocity structure by inverting the Rayleigh wave group velocity [START_REF] Herrmann | Computer Programs in Seismology: Surface Waves, receiver functions and crustal structure[END_REF]. In our inversion, a four layer model was assumed with variable thickness and velocities. The uncertainties on observed velocities are estimated to be 5%. The obtained shear-wave velocity profiles are displayed in Figure 5.13c. The indicated jump below 5 m depth is taken to be the base of the lunar regolith, consisting of impact breccias. Cooper et al. [141], using travel time inversion of refracted P -waves, found that the compresionnal wave velocity beneath the Apollo 17 site was around 100 m/s from 0 to 4 m depth and then 327 m/s from 4 to 32 m depth. Our shear wave profile thus complements and is consistent with previous observations [141,145]. Thus, by comparing these profiles, we can conclude that the v P /v S ratio of the Lunar regolith is ∼2. The correlations were not only calculated for the 3 central pairs, but also for pairs 1-2, 2-4 and 4-1. However, the amplitude of the reconstructed Rayleigh wave was found to be very low, which can possibly be explained by the fact that as the distance between those sensors is almost twice that of the central pairs, the Rayleigh pulse consequently suffers a greater amount of geometrical spreading and scattering attenuation. Another explanation is due to Sabra et al. [107] and Weaver et al. [106], who showed that for a given distance between the sensors, the efficiency of the reconstruction also depends on the wavelength, with longer wavelengths resulting in better GF reconstruction. This idea was illustrated by numerical simulations earlier in the manuscript (see Sect. 5, Chap. 3). We have also observed that the Rayleigh wave reconstruction from correlations between the central pairs filtered above 11 Hz is less efficient with increasing frequency.

As a curious feature of our correlations we found that the Rayleigh wave emergence fluctuates with time of the year. To quantify these variations, we estimated the SNR of the correlated traces as a function of record date d [106]. The signal level corresponds to the amplitude of the (causal or acausal) peak. The noise is defined as the standard deviation of the subsequent fluctuations from τ = ±2 s to τ = ±4 s (Fig. 5.14a). Rigorously speaking, this part of the correlated trace, representing apparent group velocity below 30 m/s, should obscure the rest of the impulse response (for instance reflections of the Rayleigh pulse by surrounding heterogeneities). These contributions are often much weaker than the direct pulse, and probably negligible compared with the remaining noise (the seismic signal that does not yield the GF in the correlation and the electronic noise). The SNR displayed in Figure 5.14b is calculated for both the causal (right) and acausal (left) part of the correlation. It is found to fluctuate with a periodicity of exactly 29.5 days which is also the lunation period. The SNR increases during daytime (displayed as thick horizontal bars in Fig. 5.14b), and decreases at night, pointing to an origin of the Rayleigh wave with solar illumination. The following is proposed as an explanation for the coupling between solar radiation and the Rayleigh wave generation. When the sun heats the lunar surface, temperature increases from -170 • C to +110 • C [START_REF] Langseth | Heat flow expermient[END_REF], leading to very high vertical thermal gradients, and resulting in slumping or cracking of the lunar surface material. This interpretation is in line with what was observed during the Apollo era where a large percentage of the seismic events were found to be very small Moon-quakes, termed thermal Moon-quakes, occurring with great regularity [147]. Their activity starts abruptly about 2 days after lunar sunrise and decreases rapidly after sunset.

The SNR was also calculated for correlated traces averaged over an increasing number N of time-windows d. From N = 1 to 100 (corresponding to a total record length of 300 hours), the SNR is found to increase as √ N, meaning that waves not contributing to the Rayleigh pulse add up incoherently. For longer record lengths (higher N), the SNR saturates, which might be due to remnant sources. Indeed, identical waveforms were observed from thermal Moon-quakes occurring at different times. Another reason for the SNR to saturate might be cross-talking between acquisition channels.

The correlated traces displayed in Figures 5.12 and 5.14a are not found to be symmetric in time. Though the pulse arrival times are symmetric, waveform amplitudes may change due to preferential direction of propagation of the incident wave field [START_REF] Van Tiggelen | [END_REF]120,148]. The peak amplitude for each part of the correlations, and for each couple of geophones, is evaluated. Since the reconstructed peak amplitude is proportional to the intensity of the incident waves, we interpret this asymmetry as an anisotropy of the wavefield. Let us assume that the symmetric part of the traces is generated by an isotropic flux Φ 0 (waves coming from all directions onto the array), and furthermore that the asymmetric part, corresponding to the excess of amplitude observed on one part of the traces, is generated by an anisotropic flux Φ (preferential wave incidence). The peak amplitude is therefore proportional to Φ 0 ±(1/2)Φ cos(θ i ) where θ i is the incident angle of the anisotropic flux onto the pair of geophones i -3 (Fig. 5.15). To quantify the asymmetry, we estimate the following ratio:

P + 3i -P - 3i P + 3i + P - 3i = Φ Φ 0 cos(θ i ) (5.8)
where P + 3i is the peak amplitude of the C 3i (τ) positive part. Note this ratio for pair 3-4 can be deduced from Figure 5.14a. It should be noted that this flux is not a remnant wave train (which would appear as a peak at different times τ in the three pairs). It corresponds to the preferential direction of propagation of the diffuse field intensity. Using this array technique, we evaluated quantitatively the incident direction of the seismic anisotropic flux to be +22 • relative to G4 -G3 (Fig. 5.15), pointing in the S-E direction to Steno crater, as either a dense area of scatterers or an area with intense thermal Moon-quake activity.

To conclude, we have shown that even in a quiet seismic environment with feeble sources, it is possible to obtain a Rayleigh wave dispersion curve from correlation of ambient seismic noise. From the inversion of this dispersion curve we estimated the shear wave velocity profile, which provided new information on the velocity structure of the regolith at the Apollo 17 site. Moreover, we have established the Sun as an active generator of lunar seismic noise. In particular the SNR of the reconstructed Rayleigh pulse was found to be strongly dependent on solar illumination. The seismic activity originates from the strong thermal gradients induced during lunar day as well as night. We have applied an array technique for locating noise sources. The weak anisotropy of the seismic diffuse waves reveals a preferential direction of propagation, possibly originating from an area with dense scatterers or increased thermo-seismic activity.

The results presented here establish the method of extracting the Rayleighwave GF by cross-correlating seismic noise, as indeed extendable to extraterrestrial planets that not only differ in size, evolution, and consequently seismic activity, but also in nature of origin of the noise from that of the Earth. This Ann. Phys. Fr. 31 • N o 3 • 2006 provides a novel avenue for future seismic exploration of the planets on which quakes might occur infrequently and are most probably also inhomogeneously distributed, of which Mars might be cited as an example [149]. In addition it holds the potential of increasing the scientific return, as noise between events can be successfully used rather than being discarded as has traditionally been the case. Specifically, to probe deeper into the subsurface on future seismic missions, the following points should be kept in mind. Distances between stations should preferably be between one and a few tens of wavelengths, as scattering attenuation might become important for greater distances. The depth probed being roughly 1/3 of the Rayleigh wavelengths implies going to lower frequencies which necessitates the deployment of broad-band seismometers [150].

Figure 1 . 1 .

 11 Figure 1.1. Left: the simple scattering regime. There is a perfect matching between time of flight and the source-scatterer distance. Right: multiple scattering regime. The wave trajectory is analogous to a random walk. The equivalence between time and distances is no longer valid, rendering impossible any imaging technique.

Figure 1 . 2 .

 12 Figure 1.2. Schematic view of an at sea seismic experiment. The source is emitting a short impulse that propagates in the water before penetrating into the ground. Reflected waves are sensed along the streamer made of several hydrophones. The setup is translated horizontally by the boat to probe different areas along the seismic section.

Figure 1 . 3 .

 13 Figure 1.3. The seismic intensity (logarithmic scale) backscattered by the crust in the Marmara sea plotted versus time (vertical axis) and versus the source-receiver distance.The end of the seismic record (on the left, after 5 s) shows a very short spatial coherence that reminds the optical speckle. This seismic speckle might originate from a random distribution of scatterers in the crust in the simple or multiple scattering regime. This part of the seismic record, referred to as the coda, is often neglected and not processed. At higher frequencies (right), the coda is shortened by the absorption. Such a seismic speckle is a clear evidence for the very heterogeneous nature of the Earth's crust.

Figure 1 . 4 .

 14 Figure 1.4. The preliminary reference Earth model [5] presents the average P and S-wave velocity versus depth, along with the density.

Figure 1 . 5 .

 15 Figure 1.5. (a) Distribution of M > 5 seismic events (circles and triangles) around the globe. Earthquakes are localized along main active faults and at the border of tectonic plates, leaving large aseismic areas. (b) Distribution of seismic stations (black symbols) [maps are from www.iris.edu].

Figure 1 . 6 .

 16 Figure 1.6. Example on a one day seismic record (velocity versus time in second) that includes a large earthquake. About 95% of seismic data are made of seismic noise, 4% of coda waves, and only 1% is used in traditional seismic imaging.

  Ann. Phys. Fr. 31 • N o 3 • 2006 of correlations to the impulse response. Chapter 4 studies the time-symmetry of the correlations and the role of multiple scattering. Applications to passive imaging (imaging with ambient noise) are presented in Chapter 5. Experiments include: -ultrasound passive tomography; -passive imaging of California with seismic noise; -passive imaging of the Moon's subsurface with seismic noise collected during the Appolo 17 mission.

(a) t = 1 2 . 1 .

 121 μs (b) t = 10 μs Figure Snapshot of the wave simulated with a finite difference scheme. (a) Early time: the plane wave is generated. (b) After 10 periods of propagation, the scattered wave is clearly visible.

1 Figure 2 . 2 .

 122 Differential cross-section ∂σ/∂θ of an empty cylinder of radius (a) a = 0.025 mm, (b) a = 0.25 mm et (c) a = 2.5 mm (λ = 1.5 mm). θ = 0 • marks the incident direction.

Table 2 . 1 .

 21 Parameters for the numerical simulations and evaluation of the scattering and transport cross-section of an empty cylinder. λ (mm) Radius a (mm) σ (mm) σ * (mm)

Figure 2 . 3 .

 23 Figure 2.3. Three possible approaches to describe a scalar wave through a heterogeneous medium: (a) the exact wave equation: the field is a sum of several arrivals interfering with each other. (b) Radiative transfer equation: wave packets propagate as energy quanta with an incident velocity. The phase information (the oscillation of the wave) is lost, but the packets retain the memory of their trajectory. (c) The diffusion equation: the intensity isotropically diffuses as a Brownian motion.

Figure 2 . 4 .

 24 Figure 2.4. Experimental setup for measuring the intensity of the shear and compressional part of the seismic wavefield.

Figure 2 . 5 .

 25 Figure 2.5. (a) Example of seismogram recorded in Mexico. (b) Energy partition ratio E S /E P versus time.

Figure 2

 2 Figure 2.6.Schematic view of the earth crust as a possible model for seismic waves propagation as proposed by Margerin et al.[44], and confirmed by the observation of equipartition.

Figure 2 . 7 .

 27 Figure 2.7. Illustration of the enhanced backscattering effect. (a) Single scattering regime. (b) Multiple scattering far from the source. (c) Coherent interference of reciprocal paths at the source.

Figure 2 . 8 .

 28 Figure 2.8. Schematic view of the propagation of waves from the source S to the receiver R. Dashed line is for reciprocal paths. If the source-sensor distance Δr is smaller than the wavelength, reciprocal paths interfere constructively. Crawling arrows are for the intensity diffusion between the first and the last scatterer.

Figure 2 . 9 .

 29 (a) Seismic source in the middle of an array of geophones. Solid and dashed arrows illustrate reciprocal scattered wave paths. (b) Map and picture of the volcano. The line AB marks the position of the array. Map from IGN.

Figure 2 . 10 .

 210 Figure 2.10. (a) Example of vertical ground motion signal s i (t) at the source location. (b) Zoom into the coda. (c) The cumulative ratio R as a function of time, calculated from equation (2.48), indicates whether or not the noise generated around the source by the operator is negligible. Figure reprinted with permission from [É. Larose et al., Phys. Rev. Lett. 93, 48501 (2004)]. Copyright 2007 by the American Physical Society.

Figure 2 .

 2 10c shows an example of R(t), computed for M = 5 signals recorded at the source position. R ≈ 1/ √ 5 indicates that the record is dominated by random noise, whereas R = 1 indicates that the record is strongly dominated by deterministic waves produced by the impact. This plot

Figure 2 . 11 .

 211 Figure 2.11. Top: incoherent summation of human noise. Middle: coherent summation of the seismic response. Bottom: incoherent summation of ambient background noise. Each record is made of N = 10 strikes.

Figure 2 . 12 .Figure 2 . 13 .

 212213 Figure 2.12.Energy ratio S(Δr) around 20 Hz as a function of source-receiver distance Δr for three different lapse times. The WL effect sets in at a time of roughly 0.7 s, and is fully stabilized at 1.7 s. Figure reprinted with permission from [É. Larose et al., Phys. Rev. Lett. 93, 48501 (2004)]. Copyright 2007 by the American Physical Society.

. 14 )

 14 Ann. Phys. Fr. 31 • N o 3 • 2006

Figure 3 . 1 .

 31 Figure 3.1. Plane wave incoming two receivers.

Figure 3

 3 Figure 3.3.Comparison between C AB (t) (thick continuous line) andG AB (t) ⊗ G SS (-t) ⊗ f (t) (dotted line).The overall correlation coefficient between the two waveforms is 98.7%. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America.

  Figure 3.3.Comparison between C AB (t) (thick continuous line) andG AB (t) ⊗ G SS (-t) ⊗ f (t) (dotted line).The overall correlation coefficient between the two waveforms is 98.7%. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America.

Figure 3 . 4 .

 34 Figure 3.4. Comparison between C AB (t) (thick continuous line) and G AB (t) ⊗ G SS (-t) ⊗ f (t) (dotted line) in the open scattering medium for one single source S. The impulse responses were recorded during 800 μs until they became negligible. The overall correlation coefficient between the two waveforms is 0.48%. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America.

Figure 3 . 5 .

 35 A and B are receiving points. Two hundred fifty source points are placed regularly on the perimeter a circle with radius 18.7 mm, 100 scatterers being inside the circle, which (a) completely surrounds the medium, (b) or only partially. The boundary conditions on the edges of the grid are absorbing (open medium), in both cases. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America. Ann. Phys. Fr. 31 • N o 3 • 2006

Figure 3 . 6 .

 36 Figure 3.6. Comparison between C h AS (t) ⊗ h BS (-t) ⊗ f (t) (dotted line) and G AB (t) ⊗ f (t) in the open scattering medium surrounded by 250 sources S as depicted in Figure 3.5a, at early times (a) and in the late coda, 360 μs later (b). The overall correlation coefficient between waveforms is 97.4%. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America.

Figure 3 . 7 .

 37 Figure 3.7. Correlation coefficient between S G AS (t) ⊗ G BS (-t) ⊗ f (t) (dotted line) and G AB (-t) ⊗ f (t) versus the number of sources employed. Figure reprinted with permission from [A. Derode et al., J. Acoust. Soc. Am. 113, 2973 (2003)]. Copyright 2007, Acoustical Society of America.

Figure 3 . 8 .

 38 Figure 3.8. Sketch of the first experiment. A source S sends a broadband pulse s(t) in an open medium. It is sensed by a fixed omni directional receiver R. Two 8 mm-diameter empty cylinders are placed nearby (approximative distances displayed). The source is rotated all over the setup for 2 000 different positions, along a perfect circle of radius 80 mm. Figure reprinted with permission from [É. Larose et al., Appl. Phys. Lett. 88, 104103 (2006)]. Copyright 2007, American Institute of Physics.

. 25 )Figure 3 . 9 .

 2539 Figure 3.9. Photo of a pin-transducer developed by Montaldo in the LOA-Paris. The black part is the active area composed of a reciprocal rectangular shaped piezoelectric. The transducer has a noticeable response in the 0.5-2 MHz frequency range. The transducer lateral directivity is displayed, and show a nearly omnidirectional (actually quadrupolar) pattern.

Figure 3 . 10 .

 310 Figure 3.10. Picture of the experiment. The central transducer is the fixed receiver. The second transducer is a moving source that is rotated over 2 000 different positions. The medium to image is composed of two aluminum cylinders immersed in clear water.

Figure 3 . 11 .Figure 3 . 12 .

 311312 Figure 3.11. Individual contribution of each of the 2 000 sources to the total Correlation C(τ), on a log scale. As noticed by Snieder (2004) based on a stationary phase approach, and Roux et al. [104], sources efficiently contributing to the Green's function are in the receiver-reflector direction. The wavefronts labelled (a-b-c-d) correspond to the direct (a) and secondary (b-c-d) arrivals in the retrieved Green's function (Fig. 3.12). Figure reprinted with permission from [É. Larose et al., Appl. Phys. Lett. 88, 104103 (2006)]. Copyright 2007, American Institute of Physics.

Figure 3 . 13 .

 313 Figure 3.13. (a) Experimental setup. The receiver A is fixed at the origin, the experiment is done for various position of B ranging from x B = -50 mm, y B = -15 mm to x B = 50 mm, y B = 15 mm. 21 sources S are used (size: 0.39 mm, pitch: 0.42 mm, frequency 3.1 MHz). The distance between A and S is 35 cm. (b) Waveform received by A when a 1 μs pulse is sent by one of the sources. Figure reprinted with permission from [A. Derode et al., Appl. Phys. Lett. 83, 3054 (2003)]. Copyright 2007, American Institute of Physics.

Figure 3 . 14 .

 314 Figure 3.14. Cross-correlation C AB (t) for: (a) 1 source, x B = -50 mm, y B = -5 mm; (b) 21 sources, x B = -50 mm, y B = -5 mm; (c) 21 sources, x B = 50 mm, y B = -5 mm; (d) 21 sources, x B = 25 mm, y B = 15 mm. Figure reprinted with permission from [A. Derode et al., Appl. Phys. Lett. 83, 3054 (2003)]. Copyright 2007, American Institute of Physics.

Figure 3 . 15 .

 315 Figure 3.15. Cross-correlation C AB (t) for 21 sources, the position of B is x B = -50 mm, y B ranging from -15 mm to 15 mm, (a) with and (b) without the multiple scattering slab. The thick line shows the theoretical arrival time of the GF. It is properly reconstructed (within 0.05 μs, which is the sampling time) in an angular sector of 28 • with the multiple scattering slab, versus only 6 • in water. Figure reprinted with permission from [A. Derode et al., Appl. Phys. Lett. 83, 3054 (2003)]. Copyright 2007, American Institute of Physics.

Figure 3 . 16 .

 316 Figure3.16. Schematic view of the area where the GF is reconstructed by correlation between B and A. If no scattering, the GF is reconstructed only in the red (B 1 ) area. With the scattering device, the area is enlarged and represented in green (B 2 ). In both cases, because of the finite aperture of the array of sources or scattering device, the GF is not recovered in the other area (B 0 ). The + andsigns denotes the areas where the causal or the anti-causal GF is reconstructed.

  Ann. Phys. Fr. 31 • N o 3 • 2006

Figure 3 . 17 .Figure 3 . 18 .

 317318 Figure 3.17. Sketch of the numerical simulation of 2-D acoustic waves in a homogeneous medium. Two receivers A and B, distant by D, record the wave emitted by S. Correlations are then averaged over N different sources position.

Figure 3 . 19 .

 319 Figure 3.19. Same as Figure 3.18 except D = 10λ.

Figure 3 . 20 .

 320 Figure 3.20. Schematic view of coherent and incoherent zones. If a source S is inside the coherent zones, it will contribute to the GF. If not, it will contribute to noise.

Figure 3 . 21 .

 321 Figure 3.21. The (normalized) correlation coefficient CC(S) is plotted for D = 5λ versus the location of S: [a]+[b] (thick line). Each term [a] and [b] of CC(S) are evaluated separately to compare their importance. The broken line [a] is clearly always close to zero, confirming that CC(S) is well approximated by [b].

Figure 3 . 22 .

 322 Figure 3.22. Coherent hyperbola and their focal points A and B. Dotted lines: asymptotes of these hyperbola, defining the far field coherent angle θ 0 .

2 - 1

 21 . Fr. 31 • N o 3 • 2006 Assuming far-field sources (SA D, SA D and D ≥ λ), coherent zones are defined by the asymptotic branches of the hyperbolas:

Figure 3 . 24 .

 324 Figure 3.24. Scattered dots: convergence of the correlations CC(N) versus the number of sources N for different distances D between A and B. This convergence is estimated through the normalized correlation coefficient between the reference GF and the correlations. Thick lines: theoretical curve from equation (3.36).

2 .Figure 4 . 1 .

 241 that this stabilization does not imply isotropy of the wave field; Ann. Phys. Fr. 31 • N o 3 • 2006 Comparison between numerical (Monte Carlo) solutions of the radiative transfer equation, and analytical solutions of the diffusion equation. (a) Energy density (top) and P to S energy ratio. (b) Angular distribution of elastic energy flux. The dashed and solid lines show the results of the diffusion approximation and radiative transfer equation respectively. The energy flux decreases monotonically from θ = 0 (forward direction) to θ = π (backward direction), where θ denotes the angle between the propagation direction and the source-observer vector. The results for θ = π/4, π/2 and 3π/4 are also plotted. Figures from [120], reprinted with permission from AGU. Ann. Phys. Fr. 31 • N o 3 • 2006 3. that the diffusion approximation may largely underestimate the anisotropy of the energy flux, even in the multiple scattering regime.

Figure 4 . 2 .

 42 Figure 4.2. Numerical simulation of the reconstruction of the causal and anticausal parts of the GF from cross-correlations. (a) Configuration of the numerical experiment. 1000 sources S (×) are surrounding the reference point A (+). The black dots indicate the point scatterers. (b) Snapshot of the cross-correlation between the field in A with the field at location (x, y) after averaging over the sources S for correlation time τ = -30 s. The weakly diffusive medium is characterized by the transport mean free path * = 640 km which is larger than the distance between the points where the correlations are computed. A converging wavefront is well-defined and constitutes the anti-causal part of the GF. (c) Snapshot for correlation time τ = 0 s: the wavefront is focused on A. (d) Snapshot for τ = 30 s: the diverging wavefront corresponds to the causal part of the GF. Figures from [120], reprinted with permission from AGU.

Figure 4 . 3 .Figure 4 . 3 .

 4343 Figure 4.3. Numerical simulation of the asymmetry of the reconstructed GF. (a) 40 sources S are aligned along the x-axis (crosses). The reference point is at the center of the plot, indicated by a "+". (b) Snapshot of the cross-correlation between the field in A with the one at location (x, y) after averaging over the sources S for correlation time -30 s. The converging wavefront is only partially reconstructed in the direction of the sources. (c) Snapshot for correlation time t = 0 s: the wavefront is focused on A. Note the high level of remaining fluctuations. (d) Snapshot for t = 30 s: the diverging wavefront is defined only in the direction opposite to the source region. Figures from [120], reprinted with permission from AGU.

Figure 4 . 4 .

 44 Figure 4.4. (a) Same experiment as in Figure 4.3 for time t = -30 s in the weakly diffusive medium. Only the late coda was processed, corresponding to lapse times larger than the mean free time. This part of the acoustic signal is made of multiply scattered waves propagating in all directions. The wavefront converges on A. Once again diffusion has restored the wavefront isotropy and the time-symmetry of the GF. (b) Same experiment as in Figure 4.3 for time τ = -30 s for a strongly scattering medium with mean free path * = 120 km. The converging wavefront is isotropic, and, therefore, the time symmetry of the GF is restored. Figures from [120], reprinted with permission from AGU.

Figure 4 . 5 .

 45 Figure 4.5. 10 6 → 1 scaled interpretation between the Earth crust and our ultrasonic experiment. In the crust, the scattering mean free path was estimated in Mexico [44] at 1 Hz: * = 30 km. In the aluminum waveguide, * = 5.5 mm at 1 MHz. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 . 6 .

 46 Figure 4.6. Schematic view of optical generation and detection of ultrasonic waves in an elastic body.

Figure 4 . 7 .

 47 Figure 4.7. Two possible mechanisms of elastic wave generation with a laser source. In the thermoelastic regime (low-intensity incident beam), surface waves are mainly generated.If the laser intensity is increased, the source turns to the ablation regime where the top of the metal is partially vaporized. Rayleigh waves are still emitted preferentially, but bulk waves are stronger than in the thermoelastic regime.

Figure 4 . 8 .

 48 (a) Experimental setup in transmission configuration. (b) Example of the signal detected on the top after a source on the bottom (0.8 MHz high-pass filter). Several arrivals are clearly visible: direct P and S wave, followed by multiply reflected/converted paths (P 3 , PPS, P 5 ...).

Figure 4 . 9 .

 49 Decay with absorbing layer (and scattering) (a) Averaged intensity I(t) obtained in transmission through the homogeneous aluminum slab without plasticine. (b) Averaged intensity I(t) with absorbing lateral layer, and theoretical fit (see Sect. 2.4) for scalar wave diffusion in a 2-D semi-infinite slab. Note that at late time, scattering does not play any role in the intensity decay. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 .

 4 Figure 4.10.Experimental setup: the propagation medium is an aluminum block drilled with 54 vertical cylindrical holes (diameter: 4 mm). A Q-switched Nd:YAG laser shoots on the bottom side 100 pulses (24 ns duration and 280 MW/cm 2 intensity each). On the top side a heterodyne interferometer senses the vertical displacement of the aluminum/air interface. This measure is repeated at 7 different locations (X 0 = 0 mm → X 6 = 60 mm) for each source position. Plasticine was stuck on the edges and the four lateral sides to mimic absorbing boundary conditions and avoid the generation of Rayleigh waves by mode conversion at the edges. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

  Figure 4.10.Experimental setup: the propagation medium is an aluminum block drilled with 54 vertical cylindrical holes (diameter: 4 mm). A Q-switched Nd:YAG laser shoots on the bottom side 100 pulses (24 ns duration and 280 MW/cm 2 intensity each). On the top side a heterodyne interferometer senses the vertical displacement of the aluminum/air interface. This measure is repeated at 7 different locations (X 0 = 0 mm → X 6 = 60 mm) for each source position. Plasticine was stuck on the edges and the four lateral sides to mimic absorbing boundary conditions and avoid the generation of Rayleigh waves by mode conversion at the edges. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 . 11 .

 411 Figure 4.11.Top: typical waveform obtained through the scattering slab around 1 MHz. The sensitivity of the heterodyne interferometer is 10 -4 Å/ √ Hz corresponding the optimal level of optical reflection on the sensed surface. After averaging over 100 records we reached the precision of 10 -2 Å. Bottom: intensity averaged over several source/sensor positions. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

  Elastic cross-section Figure 4.12. (a) Differential scattering cross-sections of a cylindrical cavity calculated for a compressional P or shear S incident plane wave. Each pattern is normalized by its maximum. (b) Elastic scattering cross-sections. Between 0.8 and 3.2 MHz, scattering is stronger for shear waves. At high frequencies the elastic cross-sections tend to the limit of twice the geometrical diameter. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 . 13 .Figure 4 . 14 .

 413414 Figure 4.13. GF reconstruction for different pairs of receivers. Correlations are averaged over the 35 sources and filtered in the 0.8-1.6 MHz frequency band. (a) 7 cross-correlations C ij (τ) in the diffusive aluminum plate. (b) 7 cross-correlations C 0 ij (τ) obtained in the homogeneous aluminum block, where the mode conversion possibly occurring at the edges was suppressed by the plasticine. The insets show the summation of the 6 wavefronts from X 0 = 10 mm to X 6 = 60 mm after a time-reduction based on the Rayleigh wave speed (dotted line, v R = 2.9 mm/μs). Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 . 15 .

 415 Figure 4.15. Evolution of the asymmetry in the reconstructed GF in the low frequency regime (0.8-1.6 MHz), for early (a) and late (b) times. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

  4.16 and 4.15). 45 μs is twice the time after which the diffuse energy spreads homogeneously along the array of receivers (length L = 30 mm) : L 2 /4D ≈ 22 μs in an open 2-D scattering medium. The time series were filtered in two frequency bands: 0.8-1.6 MHz and 1.6-3.2 MHz.

Figure 4 . 16 .

 416 Figure 4.16. Evolution of the asymmetry in the reconstructed GF in the high frequency regime (1.6-3.2 MHz), for early (a) and late (b) time windows. Figure reprinted with permission from [É. Larose et al., Phys. Rev. E 72, 046607 (2005)]. Copyright 2007 by the American Physical Society.

Figure 4 . 17 .

 417 Figure 4.17. (a) Homogeneous slab: no mode conversion. (b) At late times + high frequency, mode conversion to Rayleigh waves mainly occurs at one scattering cavity, leading to an anisotropic energy flux, and an asymmetric correlation. (c) At late times + low frequency, mode conversion to Rayleigh waves occurs everywhere, leading to an isotropic energy flux, and a symmetric correlation.

Figure 4 . 18 .

 418 (a) The paths between couples of stations for which average cross-correlations have been computed are shown inset. (b) Example of broadband record from a regional earthquake. Note that the signal remains well above the noise level for several hundred seconds after the S wave arrival. Figures from[120], reprinted with permission from AGU. Ann. Phys. Fr. 31 • N o 3 • 2006

4

 4

Figure 4 . 19 .

 419 Figure 4.19. Average cross-correlations as a function of interstation distance. The correlations have been computed for every permutation of components of motion, vertical (Z), radial (R) and transverse (T).Figure is from [120], reprinted with permission from AGU.

Figure 5 . 1 .

 51 Figure 5.1. setup. Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 2 .

 52 Figure 5.2. Typical experimental results. (a) Waveform G(S = 60, z 1 = 36 mm, t) received by the passive sensor at depth z 1 = 36 mm when source #60 fires a pulse. (b) Cross correlation C S=60 (z 1 = 36 mm, z 2 = 16 mm, t). (c) Cross correlation C(z 1 = 36 mm, z 2 = 16 mm, t) averaged over the 118 sources.Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

  Figure 5.2. Typical experimental results. (a) Waveform G(S = 60, z 1 = 36 mm, t) received by the passive sensor at depth z 1 = 36 mm when source #60 fires a pulse. (b) Cross correlation C S=60 (z 1 = 36 mm, z 2 = 16 mm, t). (c) Cross correlation C(z 1 = 36 mm, z 2 = 16 mm, t) averaged over the 118 sources.Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 3 .

 53 Figure 5.3. Experimental results for the evolution of the signal-to-noise ratio (in dB) versus the number of sources. The solid line represents the one-bit correlations, the dashed line for "normal" correlations. The observation points were at z 1 = 36 mm and z 2 = 16 mm. Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 4 .

 54 Figure 5.4. Same as Figure 5.2, except that only the sign (+1/ -1) of the scattered wave forms has been recorded and cross correlated (one-bit correlations). Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 5 .

 55 Figure 5.5. Experimental results for the evolution of the signalto-noise ratio (in dB) versus the number of sources for different frequency bandwidths. The observation points were at z 1 = 36 mm and z 2 = 16 mm. The scattered waveforms were one-bit digitized. The central frequency is 2.9 MHz. The frequency bandwidth Δ f is indicated in the legend, as well as α the percentage of energy within each frequency band. Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 5 ,

 55 where the SNR has been plotted versus the number of elements and the bandwidth. The proportion α of energy within a frequency band (ω 1 , ω 2 ) is denoted by ω) the power spectral density of the scattered signals. In Figure5.5, the bandwidths are indicated by the values of α.

Figure 5

 5 Figure 5.7.Sound-speed profiles obtained as the medium under investigation changes. Initially, there are two wellseparated layers (oil/syrup). Then, they are mixed together to form an apparently homogeneous emulsion. Progressively, the two phases of the emulsion separate again, and the process can be monitored by the sound-speed profiles. Twelve hours later, the separation is complete. Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

  Figure 5.7.Sound-speed profiles obtained as the medium under investigation changes. Initially, there are two wellseparated layers (oil/syrup). Then, they are mixed together to form an apparently homogeneous emulsion. Progressively, the two phases of the emulsion separate again, and the process can be monitored by the sound-speed profiles. Twelve hours later, the separation is complete. Figure reprinted with permission from [É. Larose et al., J. Appl. Phys. 95, 8393 (2004)]. Copyright 2007, American Institute of Physics.

Figure 5 . 8 .

 58 Figure 5.8. Experimental setup and typical waveforms. The sources S shoots a 1 μs pulse one after another. The receiver's directivity pattern is represented (lobes). It is rotated twice at each position. Figure reprinted with permission from [É. Larose et al., Appl. Phys. Lett. 88, 104103 (2006)]. Copyright 2007, American Institute of Physics.

20 (Figure 5 . 9 .

 2059 Figure 5.9. 118 source-averaged correlations C z i z j (τ) for 12 distances z iz j . The direct and reflected signals are superimposed. Dotted lines represent the arrival times for the direct and reflected paths. The scatterer is a plane interface positioned 33.2 mm away from the sensor array. Figure reprinted with permission from [É. Larose et al., Appl. Phys. Lett. 88, 104103 (2006)]. Copyright 2007, American Institute of Physics.

Figure 5 . 10 .

 510 Figure 5.10. (a) Paths where Rayleigh wave group velocity measurements were obtained from cross-correlations of ambient seismic noise. White triangles show locations of USArray stations. (b) Group-speed maps constructed by cross-correlating 30 days of ambient noise between USArray stations filtered from 7.5-s-period Rayleigh waves. Black solid lines show known active faults. Figure is from [93], all rights reserved.

Figure 5 . 11 .

 511 Figure 5.11. Passive imaging with seismic noise: beyond the Earth... The only solar system body other than the Earth from which we have seismic observations pertinent to its interior properties is the Moon. From 1969 to 1972 the US Apollo program installed one short-lived and four long-lived seismic stations on the Moon. All but the first station were operated until 1977. The data collected by the Apollo seismic network provided the basis for a number of studies of lunar seismicity and internal structure published during the 1970s, early 1980s, and again recently [136-140]. In addition to the above passive seismic experiment, active seismic experiments were also carried out on missions 14, 16 and 17. The active seismic experiments performed at the Apollo 14 and 16 sites were turned on only during the landing missions in active mode and weekly

Figure 5 .

 5 Figure 5.12. Configuration of the Apollo 17 LSPE experiment. It was composed of four geophones (velocimeters): G1, G2, G3 and G4, all connected to a central station for signal sampling, digitization and telemetering. Sample of raw seismic data from each geophone are displayed in the enclosed box (no logarithmic decompression), for August 16, 1976 (UTC time). Typical correlations C 13 (τ), C 23 (τ), C 34 (τ) are plotted between the corresponding pair of sensors. Crosses in the center of each correlation trace mark τ = 0. NASA photo A17.S72_37259. Figure from [24], reprinted with permission from AGU.

  Figure 5.12. Configuration of the Apollo 17 LSPE experiment. It was composed of four geophones (velocimeters): G1, G2, G3 and G4, all connected to a central station for signal sampling, digitization and telemetering. Sample of raw seismic data from each geophone are displayed in the enclosed box (no logarithmic decompression), for August 16, 1976 (UTC time). Typical correlations C 13 (τ), C 23 (τ), C 34 (τ) are plotted between the corresponding pair of sensors. Crosses in the center of each correlation trace mark τ = 0. NASA photo A17.S72_37259. Figure from [24], reprinted with permission from AGU.

  t + τ)S j (t)dt.(5.7) Correlated traces are then filtered in the frequency range of interest (4-12 Hz). The central part of each trace |τ| < 0.5 s (cross talking between geophones contribute to a peak at the correlation time τ = 0) is not displayed. Each trace is normalized by its maximum value. Finally, we computed the average of the correlations over the set of records d. Typical results are displayed in Figures 5.12, 5.13a and 5.14a. A well-defined pulse is observed for both positive and negative correlation times τ.

Figure 5 . 13 .Figure 5 . 14 .

 513514 Figure 5.13. Dispersion analysis of the Rayleigh pulse. (a) The wave packet is filtered in two non-overlapping frequency bands (around 4.5 Hz and 9 Hz). (b) Observed dispersion curve of the Rayleigh wave group velocity. Dots: observations. Crosses: calculated data from the profile shown in (c). (c) Result of the inversion of the dispersion curve showing the upper 10 m of the shear wave velocity profile. Figure from [24], reprinted with permission from AGU.

Figure 5 . 15 .

 515 Figure 5.15. Schematic view of the anisotropic flux.

Table 1 . 1 .

 11 Comparison of the physical parameters between seismology and ultrasound.

		Wavelength Frequency	Size	Duration of an experiment
	Seismology	km	mHz -Hz 6 000 km	month-year
	Ultrasound	mm	MHz	m	min

  the effect Ann. Phys. Fr. 31 • N o 3 • 2006

Table 2 . 2 .

 22 Example of scattering mean free times.

	Site	t *	frequency reference
	Volcano (Merapi)	0.1 s	10 Hz	[65]
	Volcano (Auvergne)	0.7 s	20 Hz	[66]
	Crust (Mexico)	10 s	1 Hz	[17]
	Crust (France)	60 s	2 Hz	[67]
	Mantle	2000 s	1 Hz	[68]

  in an aluminum body of dimension 10 × 10 × 10 mm, and around 1 MHz, T H = 100 μs, representing about one hundred oscillations. T H increases linearly with the volume: for bigger cavities it is therefore rapidly reaching several ms, rendering its observation very hard due to absorption and energy loss.

Ann. Phys. Fr. 31 • N o 3 • 2006

Ann. Phys. Fr. 31 • N o 3 • 2006

Acknowledgments

This manuscript includes contributions from M. Campillo, D. Clorennec, A. Derode, J. Derosny, M. Fink, R. Hennino, A. Khan, O. Lobkis, L. Margerin, G. Montaldo, Y. Nakamura, A. Paul, N. Shapiro, M. Tanter, A. Tourin, N. Trégourès, B. van Tiggelen, and R. Weaver. The author wishes to thank all these collaborators, and especially A. Derode and M. Campillo who supervised the Ph.D. thesis from which this paper is inspired. The author acknowledges an "Allocation couplée" grant from the French Ministry of Research for funding his PhD, and a Lavoisier grant from the French Ministry of Foreign Affairs along with the NSF EAR 0543328 grant for funding his stay at the University of Illinois at Urbana-Champaign as a post-doc scholar.

É. Larose

based on a flat-layered structure must be moderated by the fact that surface waves at these frequencies are strongly diffracted by topography and shallow structures which are neglected in the 1-D model. It is nevertheless a good indication that the assumption that the prominent arrivals in the actual surface/surface GF are surface waves is reasonable for the frequency and distance ranges considered here.

The reconstructed arrivals can be observed at distances as large as 200 km, confirming the existence of strong long range correlations in coda waves. As already observed by [92], and demonstrated by the comparison with synthetics, the dominant signals correspond to the fundamental Love and Rayleigh waves. A simple phase velocity measurement on the Z/Z and T/T components of Figure 4.19 confirms that the velocities of the prominent wave trains are in a good agreement with the expected dispersion curves of Rayleigh and Love waves in the simple model derived from [128]. The signal-to-noise ratio decreases with increasing distance. This is a natural consequence of the spatial decay of the GF. The deterministic signal decreases rapidly while the physical fluctuations of the diffuse fields remaining after partial averaging are independent of the epicentral distance. We assume here that the fluctuations are proportional to the square of the amplitude of the diffuse field. This amplitude is known to vary weakly with distance in the diffuse regime. It has been shown to be almost independent of distance for the late coda (e.g. Fig. 3 in [67]). These observations suggest that even amplitude characteristics could be reconstructed from field correlations. However, a reliable measurement of the amplitude decay with distance would require that the same set of earthquake records and time windows are used for each pair of stations. It is not the case here, since all stations were not operating exactly at the same time.

Application to passive imaging

In this chapter, the passive imaging technique through correlation of diffuse waves or noise is applied:

1. to make the tomographic image of a stratified medium with diffuse ultrasound (Sect. 1);

2. to localize an interface with diffuse ultrasounds (Sect.

Imaging from one-bit correlations of wideband diffuse wave fields

We present an imaging technique based on correlations of a multiply scattered wave field. Usually the Green function (GF) G AB between two points (A, B) is determined by direct transmit/receive measurement. When this is impossible, one can exploit another idea as proposed in Chapters 3 and 4: if A and B are both passive sensors, G AB can be retrieved from the cross correlation of the fields received at A and B. Though the principle is applicable to all kinds of waves, it is illustrated here by experiments performed with ultrasound in the MHz range. These experiments must be thought of as an example of "small-scale seismology": we will build an image of a layered medium through a highly scattering sample via the correlations of the ultrasonic coda waves. We will emphasize the role of broadband multiple scattering and source averaging in the efficiency of the method, and we will show the benefit of performing one-bit correlations. We will also apply this imaging technique to the case of a slowly changing medium, similarly to diffuse wave spectroscopy (DWS). Let's recall the basic notations used in our experiment and processing. A and B are two receiving points, and S is a source. We will note G IJ (t) as the wave field sensed at I when a Dirac δ(t) is sent by J. If e(t) is the excitation function in S, then time reversal) also indicate that when this condition cannot be fulfilled, one can still estimate the main features of the GF (at least an arrival time) using several tricks:

1. take advantage of a multiple scattering medium; 2. use a frequency bandwidth significantly larger than the correlation frequency of the scattering medium;

3. digitize the waveforms over a single bit;

4. average the results over all available sources. This is what was done in Figure 5.4 to retrieve the arrival time of the direct GF between the two observing points located at z 1 = 36 mm and z 2 = 16 mm.

Passive tomography

Now, we can repeat the same procedure for every pair of neighboring observation points (z 1 , z 2 ) and estimate the velocity profile (i.e., build an image) of the layered medium (Fig. 5.6). We have done so along the full length of the vertical axis, with a displacement step of 2 mm. Note that, at room temperature, it is difficult to distinguish between the water and the oil layers, since their estimated velocities are 1478 and 1472 m/s, respectively. Yet interestingly, if we repeat the experiment after heating the sample by 8 • C, we see that the measured velocity of the water layer increases to 1492 m/s, whereas that of the oil layer decreases to 1444 m/s, which is consistent with what is known of these two liquids, and the two layers are better separated on the profile (Fig. 5.6). The velocities estimated by passive imaging were found to coincide with those obtained by direct pulse-echo measurements in the four liquids to within 2%.