Time reversal of ultrasound in granular media

Time reversal (TR) focusing of ultrasound in granular packings is experimentally investigated. Pulsed elastic waves transmitted from a compressional or shear transducer source are measured by a TR mirror, reversed in time and back-propagated. We find that TR of ballistic coherent waves onto the source position is very robust regardless driving amplitude but provides poor spatial resolution. By contrast, the multiply scattered coda waves offer a finer TR focusing at small amplitude by a lens effect. However, at large amplitude, these TR focusing signals decrease significantly due to the vibration-induced rearrangement of the contact networks, leading to the breakdown of TR invariance. Our observations reveal that granular acoustics is in between particle motion and wave propagation in terms of sensitivity to perturbations. These laboratory experiments are supported by numerical simulations of elastic wave propagation in disordered 2D percolation networks of masses and springs, and should be helpful for source location problems in natural processes.


Introduction
In a non-dissipative medium, the wave equation is symmetric in time. Therefore, for every wave diverging from a pulsed source, there exists in theory a wave, the time-reversed wave, that precisely retraces all its original paths in a reverse order and converges in synchrony at the original source as if time were going backwards. This time-symmetry exists even in a strongly heterogeneous medium where waves are strongly reflected, refracted, or scattered. In the early nineties, an original method for generating such a time-reversed wave was proposed in acoustics [1]: a pulsed wave is sent from a source, propagates in an unknown media and is captured at a transducer array termed a "Time Reversal Mirror (TRM)". Then the waveforms received at each transducer are reversed in time and sent back, resulting in a wave converging at the original source regardless of the complexity of the propagation medium. TRMs have now been implemented in a variety of physical scenarios from hundreds of Hz in ocean acoustics [2] and MHz Ultrasonics [3] to GHz Microwaves [4]. Common to this broad a e-mail: arnaud.tourin@espci.fr b e-mail: xiaoping.jia@espci.fr

1488
The European Physical Journal Special Topics range of scales is a remarkable robustness exemplified by observations that the more scattering the medium, the sharper the focus [2][3][4][5][6][7].
For the last decade the time reversal focusing concept has also been a very active research in seismology, especially for seismic source imaging and source location of seismic events that exhibit no compressional (P-) and shear (S-) wave arrivals, such as tremor, glacial earthquakes and Earth hum [8,9]. In that case the real Earth, i.e., the medium where the wave field is generated and propagates, and the virtual Earth, i.e., the velocity model in which the time-reversed wave is numerically back-propagated, are however different. As a model system for athermal amorphous media or seismic fault gouges, the granular medium constitutes a particular case among strongly scattering systems [10][11][12]. Dry granular media are collections of macroscopic grains that interact through repulsive and frictional contact forces. For given values of macroscopic control parameters, such as packing density and confining pressure, granular media exhibit many microstates characterized by highly heterogeneous contact force networks that can rearrange under driving. These media whose features range from the microscopic scale (grain) to the mesoscopic scale (force-chain) and the macroscopic scale (bulk), may be modelled either as particulate or continuum materials [13].
Elastic waves that propagate from grain to grain provide a unique probe of the contact force networks. Generally speaking, one distinguishes between the longwavelength coherent (P-and S-) waves and the short-wavelength scattered waves scattered by the heterogeneous force chains [11], often referred to as coda waves. The study of the TR focusing of elastic waves in a granular medium raises two challenging issues. First, no wave equation is available at the scale of the force chains. This issue is related to a fundamental question: at what scale is the continuum elasticity applicable in a contact network [14][15][16][17]? Secondly, one may wonder whether time-reversal invariance still holds in a fragile granular medium, beyond a certain wave amplitude where the wave itself not only acts as a probe but also as a pump, leading to the acoustic fluidization of the jammed media via the rearrangement of the contact network [18][19][20]. This situation is fundamentally different from those previously reported where a perturbation was introduced in the continuous medium between the forward and backward propagation steps [6,7].
In this work, we address the above issues by experimentally investigating timereversal focusing of ultrasonic waves in glass bead packings under external load. The robustness of TR invariance is tested with a specifically developed TRM as a function of the source amplitude. A particular attention is paid to the spatial extent of the rearrangement caused by the large-amplitude driving.

Experiments
A sketch of the experimental setup is shown in Figure 1a. Our granular materials consists of dry monodisperse glass beads of diameter d = 1.5 or 3 mm, confined in a cylindrical container of diameter D = 150 mm with rigid walls (i.e., an oedometer cell) which is filled to a height of H ≈ 55 mm with a packing density of about φ ≈ 0.62. A static uniaxial stress P ≈ 85 kPa is applied to the granular packing. To perform the time-reversal experiment, a compressional or shear transducer is placed in contact with the granular packing at the top of the cell and used as a source excited by a 3-cycle tone burst centered at f = 100 kHz. We have developed a specific time-reversal mirror (TRM) with sixteen identical transducers, compressional or shear. Six other transducers surrounding the source (with a pitch of 20 mm between neighbouring transducers) are used to measure the extension of the time-reversed focal spot. The diameters of these transducers are a = 12 mm, which are sufficiently larger than the  Recompression signals recorded at the source location S1 and at the closest neighbouring detectors S2 and S3 after time-reversal and back-propagation of (a) the low-frequency coherent wave and (b) the high-frequency multiply scattered waves in bead packings (d = 1.5 mm). bead size to ensure an efficient detection of transmitted elastic waves. To investigate nonlinear effects, we vary the source amplitude from 5 to 300 V, corresponding to a vibration displacement u 0 ≈ 1-60 nm [19,20].
TR in the linear regime. Figure 1b depicts a typical waveform transmitted through a packing of 1.5 mm-diam glass beads, excited by a longitudinal transducer at small amplitude (u 0 < 5 nm) and measured by one of four transducers located at the centre of the TRM. We clearly identify the early arrival of low-frequency coherent ballistic P wave (f LF ∼ 15 kHz), followed by high-frequency waves (f HF ∼ 80 kHz) resulting from the scattering by the heterogeneous force chains, i.e., coda waves. The compressional wave velocity can be measured by the time-of-flight of P-wave pulse as V P ≈ 500 m/s, which gives a wavelength λ of about 33 mm and 6 mm (∼4d) for the long-wavelength coherent and short-wavelength coda waves, respectively. The observation of a dominant low-frequency P-wave ( Fig. 1b) stems from the large number of beads which are in contact with the ultrasonic detector (a/d ∼ 8); the spatial averaging thus enhances the coherent wave detection and partly cancels scattered waves, i.e., the acoustic speckle [11].
In a first TR experiment, only the coherent wave received by the TRM (16 transducers) is time-reversed and back-propagated. As shown in Figure 2a, the time-reversed signal is focused at the source location around the focal time t ≈ 0 (arbitrary) but with a spot much larger than the source size. In a second series of experiments, only the coda waves are selected, time-reversed and back-propagated. The TR focal spot is found to be finer with the scattered coda waves (Fig. 2b) than with the coherent wave ( Fig. 2a), providing a higher resolution of the source location. In a homogeneous medium, diffraction theory predicts a focal spot of size around λH/D where D is the aperture of the TRM and H the focusing distance. To ascertain that the fine TR focusing provided by coda waves is not simply associated with its higher frequency content, we have performed the same TR focusing experiment with a single-transducer TRM. In that case no focusing is expected for the low-frequency coherent wave but we found that time-reversed coda waves are still focused. Such observation is in agreement with previous TR experiments performed in strongly multiple scattering media [3,7] or chaotic cavity [5]: multiple scattering allows for redirection of the source angular spectrum towards the TRM, which amounts to creating a virtual aperture larger than the actual TRM aperture -a kind of lens effect [3,5]. Notice that TR of the long-wavelength coherent wave is very robust to external perturbations, e.g. by a gentle tapping between the forward and the backward propagation steps, whereas the short-wavelength coda waves are not. Relatively small amplitude of high-frequency components observed in the TR recompression signals (Fig. 2) compared to the injected pulse (inset of Fig. 1b) is due to the spatial filtering by TRM elements mentioned above.
TR in the nonlinear regime. In the following, we concentrate on the source localization by the TR focusing with multiply scattered waves. To this end, we conduct TR experiments in packings of 3 mm-diam glass beads where the elastic wave transmission is dominated by scattered coda waves with λ/d ∼ 2 (Fig. 3a), thanks to a small number of beads in contact with the detector (a/d ∼ 4). The transport mean free path is about l ∼ λ [12]. To measure the fidelity of the wave reconstruction, we follow the pulse recompression at the source location as a function of the driving amplitude either a compressional or shear source. The quality of this reconstruction can be evaluated through the signal-to-noise ratio (SNR), defined as the peak amplitude of the TR recompression signal at the focal time (centered at ≈ 100 kHz as the injected pulse) divided by the standard deviation (RMS) of the symmetrically surrounding side-lobes calculated in a time-window of arbitrary length (rectangular boxes in Fig. 3b). Here it is important to point out that these side-lobe signals are not due to the electronic noises but to the imperfections of the TR focusing [3,5], and therefore the SNR is a clear indicator of the wave fidelity.
We have verified that both the TR signals at the focal time and their side-lobes increase with the source amplitude v 0 in the linear regime when v 0 < 50 mV (or u 0 < 5 nm) (data not shown), giving rise to a constant SNR (taken as 1 in Fig. 3c). However, in the nonlinear regime (u 0 > 6 nm), the peak amplitude starts to increase in a slower way (and tends to saturation) than the side-lobes amplitude, resulting in a decrease of SNR with increasing the source amplitude. Such effect is even more pronounced with a shear source transducer as seen in Figure 3c. The loss in the fidelity of the TR focusing process is likely associated with the breakdown of the TR invariance between the forward and backward propagation steps, caused by the rearrangement of the contact network (without visible grain motion) induced by the large source amplitude [19,20]. This scenario is also consistent with the difference observed between transmitted coda signals excited at small and large source amplitudes, respectively, during the forward propagation step (Fig. 3a).
In order to evaluate the spatial extent of the network rearrangement caused by the strong source vibration, we examine the local change of the contact networks by using correlation of the configuration-sensitive multiply scattered waves, i.e., acoustic speckles [11,21]. More precisely, we measure the normalized correlation function Γ that quantifies the degree of similarity of two successive, small-amplitude multiply scattered waves (used as nondestructive acoustic probe) recorded before and after the large-amplitude driving [19,20]; these waves are transmitted through the different zones of the granular packing surrounding the source S1 (Fig. 4a). We observe that the multiply scattered waves crossing the source location zone, from transducers S2 to S3 or from S1 to D1, exhibit an important decorrelation Γ < 0.4 (Fig. 4b). On the other hand, the scattered waves primarily through the zones away from the source, e.g. from S2 to S6 on the source side or from D2 to D1 on the TRM side, remain highly correlated Γ > 0.9 (Fig. 4c). These results clearly indicate that the vibrationinduced rearrangement of the contact networks in the nonlinear regime takes place primarily in the vicinity of the source and the fidelity loss of the TR focusing is mainly due to the structural change of granular packings between the forward and backward propagation steps.

Numerical simulations
The range of frequency used in this study lies well below the first resonances of a 3 mmdiam glass bead (the shear-like spheroidal mode arises at f res = V S glass /d ∼ 1 MHz). Thus the granular network can be modelled as an effective random network of point masses (beads) and springs [22], which exhibits spatial fluctuations of both density and elastic modulus. However, in the regime of multiple wave scattering, the elastic wave equations deduced from first principles calculation [23] are not available for such amorphous-like granular media. Various numerical simulations using molecular dynamics (MD) or discrete element method (DEM) based on the frictional Hertzian interaction have thus been used out to model the wave propagation through granular packings [24][25][26].
Description of the model. Here, for simplicity, we investigate the TR focusing based on a toy model: a 2D percolation network of point masses connected by springs in which the structural disorder is obtained by randomly placing the masses m on the simple square lattice sites with a fraction p of the sites occupied. Using a uniform spring constant k, such a percolation network in 3D was previously used to simulate the heat diffusion (phonon transport) in amorphous solids where the amount of disorder can be controlled through p (as the number density) [27]. To account for the heterogeneous network of the contact force (and stiffness) in a granular packing [13][14][15][16][17], we add disorder in the distribution of spring constants which are randomly chosen from a uniform distribution between 0 and 2k 0 . Figure 5a and 5b depict sketches of a 2D percolation network where a given mass interacts in general with eight neighbours via different k; if we fill the lattice sites with monodisperse beads (disks) of diameter d (equal to the lattice constant a 0 ) using a p ≈ 0.91, the surface fraction of disks is about φ = 0.72 (φ rcp = 0.82 for the random close packing in 2D).
Unlike previous simulations [22,27], we consider here the full elastic wave propagation (longitudinal and transverse modes) through a percolating network of L × L (L = 70 in Fig. 5b) with the displacement field r i (t) of a given bead at the site i in the plane (x, y). To ensure the linear response to a transverse or shear where r i (t) = u i (t) + r 0 i with u i (t) the dynamic displacement (||u i || < 0.5||r 0 i ||), r ii is the distance vector between the bead i and the nearest neighbour i , a 0 the spring length at rest (lattice constant), K ii the fluctuating spring constants and β a damping constant. Figure 5c shows a typical wave transmission field excited by a source consisting of three beads near the left wall (x = 4d), oscillating along the x-axis u x S (r 0 , t) with one-cycle of sine at frequency f = 0.2f 0 with f 0 = k 0 /m (one takes m = 1 and k 0 = 1). This snapshot of the dynamic displacement u x (r, t) at t = 46t 0 (t 0 = 1/f 0 ) (β = 5 · 10 −3 m/t 0 ) reveals a wave field composed of a weak coherent longitudinal wave (near x = 60 d) followed by an irregular interference pattern likely due to multiply scattered waves; this observation is consistent with the temporal response u x i (t) detected at one site i (Fig. 5d). For the frequency used here, the longitudinal wavelength is λ = V P /f ≈ 6d with attenuation length l e ≈ 10d dominated by scattering (data not shown), which indicates a weak multiple scattering regime (l e < 0.15 L). Recompressed signal by TR at the source r 0. Nonlinear regime: (c) snapshot of a TR wave field u x (r, t) at t ≈ T in the presence of rearrangements. TR focusing around the source location is less efficient than in the case without rearrangement shown in (a); (d) fidelity R decreases as the driving amplitude increases. Inset: comparison of recompressed signals by TR focusing with (in red) and without (in blue) rearrangements.
TR focusing in the linear regime. To test the TR focusing process, we take seventy beads near the right wall (at x = 66 d) as receivers (i.e., TRM). We then time-reverse the displacement signals u D (r, t) detected by each TRM element and back-propagate u D (r, T − t) with T the signal duration. Figure 6a shows a snapshot of the u x field near the focal time t = T (see the whole movie in the supplementary material). We observe that waves add up coherently near the source location but incoherently elsewhere, indicating that time-reversed waves do converge back to the source. We also check in Figure 6b the recompressed displacement signal Ψ(t) = u x S (r 0 , t) at the source location; the sharp peak at t ≈ T confirms the TR focusing of elastic waves at the source.
TR focusing in the nonlinear regime. We now seek to model the irreversible soundmatter interaction effect, i.e., rearrangements of the contact network observed in TR experiments. This phenomenon is presumably related to the contact sliding between grains by the large source driving which modify partly the initial contact network without visible motion of grains [19,20]. To simulate the modified contacts during the forward propagation step, we compare the maximum dynamic displacement ||u max i || at each bead to its static stretching ||r 0 i || and apply a local yield criterion ||u max i || ≥ 0.02||r 0 i ||. The larger the source amplitude, the larger the number of modified contact and the rearranged zone (data not shown). We ascribe to each modified contact a new value k randomly chosen from the uniform distribution of stiffness and build up a new configuration of the network. Hence, we time-reverse the forward propagating signals through the initial network for a given source amplitude and back-propagate them in the accordingly rearranged network (note that in experiments the initial network is modified during the forward propagation step). Figure 6c shows a snapshot of the displacement field u x after the TR process at the focal time in the presence of rearrangements (see the whole movie in the supplementary material). Compared to the case without rearrangement at small source driving (Fig. 6a), we observe that the focusing spot around the source is less intense and that the wave field manifests some leakage of energy from the source location to elsewhere (see also films related to the TR focusing process in the supplementary material). Moreover, an amplitude decrease is found in the recompressed signal Ψ * (t) on the source at t ≈ T (inset of Fig. 6d). Figure 6d shows the ratio R of Ψ * max to Ψ max obtained without rearrangement, calculated as a function of the source amplitude. The decrease of R which may quantify the loss of fidelity in the TR process [6,28] supports well the experimental observations when the source driving is increased ( Fig. 3c and Fig. 4).

Discussions and conclusion
As mentioned in the introduction, the acoustic TR invariance holds in a multiple scattering medium on the shortest wavelength [1][2][3][4][5][6][7]. Consider a scalar displacement u(r, t) that satisfies the wave equation (dissipationless), with ρ(r) the density and K(r) the elastic modulus of the heterogeneous medium.
A source located at r 0 transmits a short pulse u S (r 0 , t) (= δ(t), dirac-like) into the medium, the multiply scattered signals u D (r j , t) = h j (t) are detected by receivers no. j (TRMs) at r j where h j (t) would denote the impulse response in linear acoustics. They are time-reversed u D (r j , T − t) with T the signal duration and retransmitted into the medium. Taking into account the reciprocity principle, interchanging the source and the receiver does not alter the resulting wave field. The signal recreated at the source location can thus be written as s(t) = j h j (T − t) ⊗ h j (t), which is maximum at time t = T indicating to a TR focusing in time at r 0 . In this study, we have shown by experiments (and simulations) in granular packings that the TR process using multiply scattered waves (λ ≈ 2 − 4d) applies correctly in the linear regime. These observations hence suggest that the elastic modulus K and elastic wave velocity V = K/ρ defined in equation 2 within the continuum elasticity still holds on the local scale of a few grain size. This appears consistent with recent wave velocity measurements in stressed granular layers with thickness ∼ 5d where the fluctuation of V P remains less than 10% [29]. However, this TR invariance breaks down in the nonlinear regime with large source driving where the sound-matter interaction becomes irreversible (Fig. 3) due to the contact slipping associated with the nonaffine deformation in granular packings [14-17, 19,20] particularly under shear (see Fig. 3c). Therefore, the propagation medium is different between the forward and backward steps, leading to a modified TR recompression signal s (t) = j h j (T − t) ⊗ g j (t) at the source location where we denote by g j (t) the impulse response of the rearranged granular network (the backward propagation is indeed performed in the linear regime). The decrease of TR focusing, or loss of fidelity, shown in Figure 3c is related to the nonlinear sound-induced structure change, being very different from those observed in other TR experiments where the medium is stationary (as in an ordinary elastic material) and before reaching the shock formation: if the whole harmonics generated during forward propagation are recorded, time-reversed and retransmitted, the initial wave can be reconstructed [30].
Time-reversal has also been used as a diagnostic tool to test and compare the sensitivity of particle motion and wave propagation to perturbations in the initial conditions or in the propagation medium [6,28]. In a multiple scattering system, particles and waves exhibit fundamentally different behaviours: particle motion is chaotic, incapable of returning to the source, while wave propagation is much more stable, despite TR invariance in both Newton's law and the wave equation (Eq. (2)). The physical reason for this can be explained through the concept of ray splitting as follows. A particle follows a well defined trajectory whereas waves travel along all ray directions (a huge number of trajectories) visiting scatterers in all positions within the scale of the wavelength. While a small perturbation can make the particle miss one obstacle or scatterer and completely change its future trajectory, the wave amplitude is much more stable due to coarse graining on the wavelength scale. Granular acoustics described here appears in between particle motion and wave propagation in terms of the sensitivity to perturbations, closely related to the discrete nature of the contact network with a finite coordination number Z ≈ 6 in 3D bead packings (via Eq. 1) in which waves propagate along preferred paths. Further studies are needed to investigate the loss of fidelity as a function of the arrival time of multiply scattered waves [6].
Finally, we think that our study may be useful to those who are using timereversal to seismic source locations or defect detections in fractured materials. In the context of seismic imaging, the application of TR process requires that the difference between the real Earth (forward propagation) and the Earth model (backward propagation) can be neglected at the used wavelength [8,9]. In laboratory experiments, an unknown source may be localized without the use of a numerical model when the reversed wave field is accessible by measurements. This is precisely the case for detecting fissures near the surface of a solid by TR process where harmonics generated by the nonlinear scatters, i.e., defects are selected and time-reversed to the source location [31]. However, our TR experiments in granular materials have shown that an important rearrangement of the medium by a large source driving may decrease the accuracy of the TR focusing and the source location ( Fig. 6) in fractured earth materials.
In conclusion, we believe that our time-reversal investigation in granular media may help to get a better understanding of the instability of wave scattering in nonlinear disordered media [32] and of the imprint of classical chaos on wave (quantum) systems [6,28]. It should also be useful for studying the source locations of seismic events and rupture in heterogeneous materials. We will improve TR focusing process inside granular media, especially to learn to focus high-amplitude ultrasound at a particular position to trigger rearrangements.