Passive Retrieval of Rayleigh Waves in Disordered Elastic Media

When averaged over sources or disorder, cross-correlation of diffuse fields yield the Green's function between two passive sensors. This technique is applied to elastic ultrasonic waves in an open scattering slab mimicking seismic waves in the Earth's crust. It appears that the Rayleigh wave reconstruction depends on the scattering properties of the elastic slab. Special attention is paid to the specific role of bulk to Rayleigh wave coupling, which may result in unexpected phenomena like a persistent time-asymmetry in the diffuse regime.


I. INTRODUCTION
Whatever the type of waves involved, knowing the Green's function of an heterogeneous medium is the key to many essential applications like imaging, communication or detection. In the last twenty years or so, mesoscopic physics has intensely studied wave phenomena in strongly disordered media: weak and strong localization, radiative transfer and diffusion approximation etc. [1]. But the exact Green's function of a complex medium is not easily tractable, and usually theoreticians study ensemble-averaged quantities like statistical correlations of intensities or wave fields. Moreover, from an experimental point of view, it is not always possible to measure the Green's functions of a complex medium because it requires controllable arrays of sources and receivers that do not perturb the medium. In the context of laboratory ultrasound, this is (nearly) routine; but in other fields of wave physics like seismology, distant sources are not controllable. In that respect, a good deal of publications followed recent works by Weaver. He proposed to cross-correlate the diffuse wave fields obtained at two passive sensors and showed it yields the elastic Green's function between the two receivers, as if one of the receivers was a source [2,3]. That correlations performed on passive sensors should yield the wave travel times is not that new. This principle was applied to helioseismology in the 90's where it provided tomographic images of the Sun's interior [4,5]. Beyond travel time reconstruction, Weaver's experimental retrieval of exact Green's functions was a real breakthrough.
Weaver's work has been followed by different contributions. Some theoretical works are based on the ergodic approximation where time and ensemble average coincide. This approach is very useful for estimating the role of scattering in the correlation asymmetry [6,7]. * eric.larose@ujf-grenoble.fr When several sources are available, another possibility is to average correlations over sources without moving the receivers. This was done in underwater acoustics [8]. Sabra et al [9] showed the possibility of recovering the entire Green's function of the sea waveguide (especially the late contributions of multiples) and proposed a model for estimating the signal-to-noise ratio of the correlations. Derode et al [10] proposed to interpret the Green's function reconstruction in terms of a Time-Reversal analogy, and showed that all the benefits of time-reversal devices in multiple scattering media could be fruitfully applied to "passive imaging" from correlations [11,12]. This idea is based on the mathematical principle of the representation theorem [13], equivalently referred to as the Helmholtz-Kirchhoff theorem [14]: if a perfect series of receivers has been sensing the wave field for ages, then one can mathematically have access to the wave field anytime and anywhere in the area enclosed by the sensors. Very recently, this principle was fully developed and applied to elastic waves in open media by Wapenaar [15].
From the very start of Weaver's work, the correlation principle was successfully applied to seismic waves [17,18,19]. The most energetic part of the Green's function between two seismic stations (Rayleigh and Love wave trains) was retrieved from passive correlations of either coda waves or records of seismic noise. As we mentioned earlier, obtaining impulse responses without a controllable source is of high interest in seismology since it gives the possibility of simulating very energetic and punctual earthquakes everywhere around the Earth, and do imaging without a controllable source ("passive imaging"). Results are quite encouraging although several theoretical problems remain unsolved. They concern issues very specific to the Earth. In particular elastic sources (earthquakes) are naturally never arranged in such a way that they would perfectly surround a couple of seismic stations. In addition, for distant earthquakes, the propagation directions of incident bulk waves are mostly vertical, and in a vertically layered medium they would not couple with Rayleigh waves.
FIG. 1: 10 6 → 1 scaled analogy between the Earth crust and our ultrasonic experiment. In the crust, the scattering mean free path was estimated in Mexico [16] at 1 Hz: ℓ * = 30 km. In the aluminum waveguide we numerically found ℓ * = 5.5 mm at 1 MHz.
So, why do we observe a Rayleigh wave train in the correlation of waves generated by bulk sources? The Green's Function retrieval is linked to equipartition, which is due to wave scattering. So, what is the influence of scattering within the Earth crust in that process?
Here we present experimental results obtained in the lab with laser-induced and laser-detected ultrasonic waves propagating in a heterogeneous slab mimicking the Earth crust. Lots of previous experimental articles applied the "passive imaging" technique to acoustic waves. Here we propose to investigate the emergence of the Green's function in the correlations of elastic waves propagating in a solid heterogeneous medium with a free surface. Because field experiments are tedious and natural environment are mostly unknown (especially the scattering properties), we propose to build an Earth crust model at the scale 1/10 −6 (presented in fig. 1). Waves will be 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 [7], 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. The next section describes the experimental setup and the propagation medium. Section III presents a short theoretical study of the scattering properties (scatteringcross sections σ of a single scatterer and the transport mean-free paths ℓ * P and ℓ * S of the heterogeneous sample). This study is supported by experimental measurements. In section IV, the passive imaging technique is applied to records acquired at the free surface. Time and frequency analysis are proposed, and a brief discussion on the timesymmetry of the correlations concludes the article. The experimental setup is depicted in fig. 2. It was designed to mimic the propagation and scattering of elastic waves through the Earth crust, at ultrasonic frequencies (0.8-3.2 MHz). A duraluminum slab whose dimensions are roughly 10 −6 those of the crust was used. 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 . Ideally, the slab should have had infinite dimensions along x and y. To approach this condition, we stuck a thick layer of dense plasticine on the lateral sides and edges of the duraluminum sample. It was aimed at 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. With the plasticine it decreased to 120 µs (see next section for a detailed discussion on the absorption time). We therefore simulated an open slab in the x and y direction with a free surface at the top. Rigorously, the earth crust is a waveguide that partially leaks energy through the Moho 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 set-up is to mimic an elastic and scattering medium, infinite in the horizontal directions with a free surface on the top.
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). Two regimes of elastic wave generation are possible with a laser source [20]. When the surfacic intensity of the pulse is weak (< 15 MW/cm 2 ), the surface is locally heated up and its dilation creates Rayleigh waves (thermoelastic regime). At higher intensity, the laser evaporates part of the metal. In this ablation regime, both Rayleigh and bulk waves are generated. Theoretical radiation patterns are displayed on fig. 3 for compressional and shear waves. 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.
As to the detection of the free surface motion, it was achieved with a contactless and quasi-punctual device: a heterodyne optical interferometer developed by Royer et al [21] which has the advantage of a very broadband response (20 kHz-45 MHz) and a sensitivity of 10 −4Å / √ Hz. It was mounted perpendicularly to the slab and then provided us with the absolute vertical component of the free surface displacement (top side), with a fine spatial resolution; the size of the laser spot was ∼ 100 µm whereas the typical elastic wavelengths here are ranging 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 only measured the vertical movements of the free surface.
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 the x-z plane (see fig. 1) 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. Scattering at the grain edge is presumably also negligible compared to scattering by the void cylinders. The overall translational symmetry along y 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 in our article 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 absorbed by the plasticine and lost. 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. A brief description of the calculation is given in the Appendix. For a detailed derivation we refer to [22,23,24]. The differential scattering cross-section gives the angular distribution of the scattered surfacic intensity, normalized by the incident surfacic intensity. The total elastic cross-section is σ = ∂σ ∂θ dθ. In 2D it has the dimensions of a 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 noted σ P P , σ P S and σ P = σ P P + σ P S , 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 .
The frequencies ranging from 0.1 MHz to 200 MHz. In average in the frequency band of interest (0.8-3.2 MHz), we obtained σ P = 9.2 mm. This value is comparable to measurements by White [25]. The same calculations were conducted for an incident shear wave, we found an average of σ S = 12 mm in the 0.8-3.2 MHz frequency band.

B. Transport properties
When the elastic wave propagates through the aluminum slab drilled with holes, it undergoes multiple scattering. Let ϕ(t) (resp. ϕ 0 (t)) be the vertical displacements sensed at the free surface through the scattering (resp. 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: cylinder's positions). We underline the difference between the coherent and the ballistic wave (i.e. the first arrival). For a detailed discussion about scattering effect on coherent and ballistic waves, see [26]. 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 the scatterers density n and their elastic crosssection σ: 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. 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 [27]: 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 [28,29,30,31] (also referred to as weak localization) does give a direct estimation of the transport mean free path but 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) [32], under the diffusion approximation. In an infinite slab of thickness H with perfect reflections on both sides, the averaged transmitted intensity reads: 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 ( vp vs ) 2 ℓ * S ) ≈ 40 mm 2 /µs. This assumption is valid after a couple of mean free times, when the equipartition regime [33] is set. Equipartition means that the density of compressional and trans-verse modes equilibrates. Considering the specific velocities of each mode [34], we infer that 80% of the energy is transported by S waves, and only 19% by P waves (and an additional 1% for surface waves). Hence our best fit ( fig. 7) of the intensity decay in the coda gives τ abs = 120 µs±10% and ℓ * = 5 − 20 mm.

IV. TWO-POINT CORRELATION OF DIFFUSE FIELDS
In this section we focus on the experimental reconstruction of the direct Green's function 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 at 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 for 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 was lasting ≈ 800 µs before reaching the noise level (see fig. 4). We underline 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 Green's function. This modal approach is unrelevant to our experiment. The records were correlated and averaged over the 35 available sources. For the scattering slab, this reads: 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:  To enhance the signal-to-noise ratio, each correlation is time symmetrized (= C(+τ ) + C(−τ ) ) and normalized by its maximum. Results are displayed in fig. 8 and 9. 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 enclosed box on 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 [11] 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 is 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 ( fig. 10 and 11). 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 2D scattering medium. The time series were filtered in two frequency bands: 0.8-1.6 MHz and 1.6-3.2 MHz.
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 fig. 10 and fig. 11). 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 fig. 10  and fig. 11). 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 X 0 's side. 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 [6,36] (see right part of fig. 10): 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. 11). 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 [35]. 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 not perturbed (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. 10). 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 thanks 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 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. 8), where lots of cavities are encountered in the first half wavelength, than at higher frequencies where the Rayleigh wavefront can propagate freely ( fig. 10).
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.

V. CONCLUSION
In this paper were presented laboratory experiments of elastic wave propagation in heterogeneous media at ultrasonic frequencies. An aluminum slab was made quasi-infinite by the use of absorbing boundaries to mimic the Earth's crust, in which scattering was obtained drilling cylindrical cavities. A relatively simple theoretical model for wave scattering properties was proposed. Wave field generation and detection was achieved using contactless and quasi-punctual laser devices. 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 Green's function passive reconstruction.
Analysis for different time-windows and frequency bands were conducted. Previous works in acoustics [37] and seismology [17] observed that the Rayleigh wave reconstruction was harder with increasing frequency. Here we found that the Rayleigh wave reconstruction was more efficient in the high frequency band. 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. The trade-off between the two effects should be further investigated. Though our experiment was designed for seismological applications, results should be applicable to other fields of wave physics where both surface and bulk waves are present. 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: 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 [22,23,24]. The scattering crosssections are σ P = σ P →P + σ P →S The corresponding differential scattering cross-sections are given by: