# Sound propagation in a constrained lattice of beads: High-frequency behavior and dispersion relation

Christophe Coste, Bruno Gilles

► To cite this version:

Christophe Coste, Bruno Gilles. Sound propagation in a constrained lattice of beads: High-frequency behavior and dispersion relation. Physical Review E, 2008, 77 (2), pp.021302. 10.1103/Phys-RevE.77.021302. hal-01406704

HAL Id: hal-01406704

https://u-paris.hal.science/hal-01406704

Submitted on 1 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Sound propagation in a constrained lattice of beads: High-frequency behavior and dispersion relation

Christophe Coste

INSP, Université Pierre et Marie Curie-Paris 6, Université Denis Diderot -Paris 7, CNRS, UMR 7588, Campus Boucicaut,
140 rue de Lourmel, 75015 Paris, France

Bruno Gilles

Université de Lyon, F-69003 Lyon, France; Inserm, U556, F-69003 Lyon, France; and Université Claude Bernard,
Lyon 1, F-69000 Lyon, France

(Received 8 March 2007; revised manuscript received 31 October 2007; published 7 February 2008)

We report on acoustic wave propagation in a regular array of nominally identical beads under isotropic static stress. The weak polydispersity of the beads makes the contact lattice random. Time-frequency analysis of the acoustic signal is performed and allows measurement of the full lattice dispersion relation. Comparison with the theoretical prediction for a perfect triangular lattice gives an indication of the level of randomness in the contact lattice. The results extend, in a consistent way, a previous study restricted to long wavelength propagation [B. Gilles and C. Coste, Phys. Rev. Lett. 90, 174302 (2003)]. The contact lattice is ordered by increasing the stress, and the smaller the wavelength, the higher the stress required to get regular lattice behavior. Measurements involving ballistic propagation of the coherent wave, whatever its frequency, evidence reversible lattice behavior under compression and/or decompression. Nevertheless, correlations of short wavelength incoherent waves are a sensitive probe of disorder, and allow us to exhibit a small irreversible evolution of the lattice.

DOI: 10.1103/PhysRevE.77.021302

PACS number(s): 45.70.—n, 43.58.+z, 43.60.+d

## I. INTRODUCTION

Static granular media exhibit solidlike behavior since they resist their own weight, may hold external loading, and carry acoustical waves, but in a very peculiar way: Even under symmetrical loads, their internal stresses exhibit very large inhomogeneities, evidenced by the large fluctuations of inter-grain forces [1]. Since the contacts between adjacent grains are pointlike, even a very small amount of disorder in the size distribution of the particles leads to a completely random lattice of active contacts [2–4]. A second intricacy arises from the fact that the elastic contact law between grains is nonlinear, because of the local curvature at the contact. It has been calculated by Hertz for elastic spheres [5], and may be generalized to particles of any regular shape using the local radius of curvature.

Acoustic waves may propagate in granular media, and have been studied for a long time [6–14]. Since sound waves are necessarily carried from one grain to another by actual contacts, the intergrain contact law has a strong influence on the propagation. Indeed, the Hertz law implies that the stiffness of an intergrain contact depends on the normal force between the two grains. As a consequence, the sound velocity $c_s$ should vary as $p^{1/6}$ if $p$ is the isotropic pressure applied on the granular material [7].

Experimentally, at low pressure the sound velocity rather follows a $c_s \propto p^{1/4}$ non-Hertzian power law [15,16]. This behavior is observed even for model systems such as regular lattices of nominally identical beads [7,13], the Hertzian power law being recovered at high pressure only. The Hertzian interaction law might be put into doubt, since actual beads are quite far from ideal elastic spheres. Experiments performed on one-dimensional (1D) chains of identical

beads, a geometry where all contacts are obviously activated, have demonstrated that the Hertz law accurately describes the properties of linear [17] and nonlinear [17,18] wave propagation.

Non-Hertzian behaviors have also been observed numerically, and have been attributed to the disorder of the contact lattice [19–24]. Actual beads always display some scatter in their radius, and for two-dimensional or three-dimensional pilings the contact lattice is inevitably disordered. The deformations induced by pressure increments progressively close initially open contacts. This picture is at the core of a recent theory [25], in which an effective medium approach leads to the calculation of the static compressibility of a two-dimensional (2D) triangular lattice of almost identical spheres. The reference state is the high-pressure limit, when all contacts are closed and the lattice is perfectly ordered. The effective medium approximation should be correct for small values of the relative fluctuations of the forces between the grains. The theory is indeed in very good agreement with simulations of the same system by Roux [20], at moderate or high pressure. In the simulations, the pressure may be decreased down to the rigidity threshold of the lattice, and discrepancy is observed at the lowest available pressures. The calculations are also in very good agreement with experiments [13] on low-frequency acoustic wave propagation in a triangular lattice of beads. The only fitting parameter is the uncertainty on the bead diameter, which is set to $\pm 2$ µm in the calculations, to be compared to the nominal value of $\pm 4$ µm in the experiments, for 8 mm beads. It is also shown in [25] that the Hertzian $p^{1/6}$ power law is valid at high pressures, while no simple functional form can describe the dependency of the compressibility on the pressure at lower pressures.
We thus have a fully consistent picture of the lattice behavior at low frequency. The aim of the present work is to study high-frequency wave propagation, when dispersion occurs. Appendix B is devoted to a presentation of some aspects of the dispersion in a triangular lattice. For a discrete lattice, quite significant dispersive behavior should be expected when the difference between the sound velocity $c_s$ and the group velocity $v_g$ is more than 10%. Using the dispersion relation (B1), and the expressions (B2) and (B3) for $c_s$ and $v_g$, it means that the wavelength below which propagation is considered as clearly dispersive is $\lambda_{disp} \approx 5d$, where $d$ is one bead diameter. In our previous study [13], the wavelength was between 8 and 16 bead diameters, which may consistently be called a low-frequency regime. In what follows, wavelengths less than 3 bead diameters are used (see Table II of Appendix B), for which dispersive effects are expected to occur. The focus will be on using such high-frequency waves to determine the level of contact disorder in the regime identified as Hertzian according to low-frequency measurements.

The paper is organized as follows. We describe the experimental setup in Sec. II A, the experimental protocol in Sec. II B, and the wide-band acoustic emitter, specially designed to send well-defined Gaussian bursts of high-frequency waves, in Sec. II C. In Sec. III, we briefly recall the notion of coherent wave propagation in random media. Section IV deals with cross correlation analysis of the received wave forms. In Sec. IV A, we give an experimentally consistent definition of the coherent part of our signals. Then we report on time-of-flight measurements for a low-frequency coherent wave in Sec. IV B. Cross correlations for a low-frequency incoherent wave are provided in Sec. IV C. Considering several different initial configurations of the lattice, we emphasize the behaviors that are configuration dependent, and those common to every configuration. The last paragraph in Sec. IV D is concerned with cross correlations for a high-frequency incoherent wave. Section V is focused on time-frequency data analysis. The dispersion relation is briefly mentioned in Sec. V A. Then, Sec. V B summarizes the basis of time-frequency analysis. In Sec. V C, we represent the coherent wave in the time-frequency domain. Then we calculate the group delay as a function of the frequency, and compare it to the theoretical prediction in Sec. V D. The reversibility of a granular medium behavior during compression and decompression is a subject of fundamental interest [15,20]. We review and discuss all relevant experimental evidence provided in the present work in Sec. VI. Section VII sums up our conclusions. In Appendix A we recall the definitions of correlations used throughout the paper, and in Appendix B we gather relevant information about the perfect triangular lattice dispersion relation.

## II. EXPERIMENTAL SETUP

### A. Bead lattice

The experimental setup is sketched in Fig. 1 (a more detailed description can be found in [13]). Stainless steel beads of diameter $d=8$ mm $\pm$ 4 $\mu$m are held between two horizontal plates of polytetrafluoroethylene, in order to get low

<figure><img src="image_2.png" /></figure>

FIG. 1. Left-hand side: Sketch of the experimental setup. See Sec. II C and Fig. 2 for details about the emitter. Right-hand side: Details of the coupling between the emitter and the lattice. The gray slab is a piezoceramic sandwiched between two one-half steel beads. It constitutes a sensor which allows in situ measurement of the emitted acoustical signal during the experiments.

friction and poor impedance matching between the beads and the framework [17]. They are disposed so that they form a triangular lattice, contained in a hexagonal cell. The hexagon sides are $31d$ in length. The scatter in bead dimensions is so small that the arrangement looks like a perfect lattice, but nevertheless ensures that most beads are not actually in contact with all their neighbors. The lattice of effective contacts is thus disordered, as shown below.

The three bottom sides, in the figure, of the hexagon are fixed, while the others may move independently along their normal directions. The motion of each side is ensured by a computer controlled stepping motor, and the force exerted on the lattice is measured by a static force sensor. A feedback loop allows us to apply a given force to the lattice, which is held constant during an acoustic measurement. In what follows, each moving side exerts the same force on the lattice (within an experimental uncertainty of $\pm 2\%$), in order to apply the stress isotropically.

An important feature of our experiment is that both the acoustic transmitter and the sensor are in contact with only one bead. This configuration favors the observation of disorder effects on the wave propagation [12]. The emitter (see Sec. II C) is specially designed in this way. The receiver is a small piezoelectric pressure sensor (Pinducer Valpey-Fisher VP-1093) of cylindrical shape, 2.5 mm in diameter. This sensor is placed exactly in front of the emitter, so that their relative distance is $31d\sqrt{3}$, i.e., 42.95 cm.

### B. Measurement protocol

Each experimental run proceeds in the same way: first, we arrange the beads so that they form a triangular lattice. Then we apply the static force on the lattice, increasing it without making any measurement: since the lattice can support the stress, the contact network is established. At the end of this first compression, we get a specific contact lattice that depends on the unknown initial distribution of bead diameters. This peculiar lattice will be further referred to as a configuration. Then, we proceed to the measurements by slowly decreasing the stress, dividing the force by a constant factor of 1.14 at each step. A state of the lattice will be given by its initial configuration, the current static force, and possibly the
history of decompression and compression of the lattice. In a given state, we send several acoustic burst waves of a Gaussian envelope (see Fig. 3), with different center frequencies. The deformation due to the acoustic wave is very small and does not disturb the lattice state, so that, for a given applied force, low- and high-frequency waves propagate in the same state of the same configuration. When all waves are recorded, we change the force and resume the process. Having reached the minimum measurable force, above the rigidity threshold, we may increase it again, keeping the same configuration. Experiments on a new configuration can be done by opening the setup, taking off the beads, mixing them, and then starting up again the whole process.

The lattice behaves as a static random medium, described by characteristic length scales rather than time scales. This motivates the use of constant wavelength waves. In this system, the sound velocity changes with the applied force and keeping the wavelength constant requires frequency adjustment. Since the positions of the emitter and sensor are fixed, the wavelength cannot be measured but may be estimated from frequency measurements. To obtain waves of (roughly) 3 cm (i.e., about 4 bead diameters) constant wavelength, we adjust the emitted-wave center frequency using the velocity measured in low-frequency experiments for each static force. Those frequencies are indicated in the relevant figure captions. In Table I we give the characteristics of the different emitted waves. The shortest wavelength available in the experiments is provided in Table II, where it is shown that it is about 3 bead diameters. This should be compared to the cut-off wavelength which is $\sqrt{3d}$.

### C. Acoustic emitter

The emitter is sketched on the left-hand side of Fig. 2. It consists of two piezoceramic elements, sandwiched between two duralumin parts that act as loads to lower the available frequency range [26,27]. The oscillations of the piezoceramic elements are longitudinal and result in the emission of compressional waves. The duralumin parts are screwed in order to prestress the piezoceramic elements, which ensures a good electrical contact. Moreover, the force exerted on the lattice, being small in comparison with the emitter's prestress, does not modify its acoustical response. The reduction of the diameter of the front duralumin load is intended to concentrate acoustic energy [28] in the pointlike region of contact between the emitter and the bead. This design allows us to send burst waves of limited duration, with a Gaussian envelope. To avoid distortion, we work quite far from the emitter resonances: Its spectral response is shown in Fig. 2, and is indeed quite flat.

We show in Fig. 3 the electric excitation that drives the emitter, and the acoustic response. Each burst has a duration of roughly four oscillations. The coupling between the emitter and the lattice takes place on the very small Hertz contact between the flat extremity of the emitter and a bead. The emitted signal is thus recorded with an ad hoc sensor, made of a piezoelectric slab sandwiched between two one-half steel beads, shown in Fig. 1. The slab resonance frequency is very high (typically 2 MHz), and its frequency response in

<table><caption>TABLE I. Characteristics of the emitted acoustic waves (in parentheses, wavelength in bead diameters is given). The first three lines are for constant frequency experiments, and the last line is for constant wavelength.</caption><thead><tr><th>Force range (daN)</th><th>Velocity range (m/s)</th><th>Frequency range (kHz)</th><th>Wavelength range (cm)</th></tr></thead><tbody><tr><td>3 → 150</td><td>400 → 800</td><td>6.5</td><td>6 → 12 (8d → 16d)</td></tr><tr><td>3 → 150</td><td>400 → 800</td><td>15.2</td><td>3.5 → 7 (4d → 9d)</td></tr><tr><td>3 → 150</td><td>400 → 800</td><td>19.5</td><td>2.5 → 5 (3d → 7d)</td></tr><tr><td>3 → 150</td><td>400 → 800</td><td>13.6 → 26</td><td>3 (4d)</td></tr></tbody></table>

our operational range is basically flat. This sensor is in contact with one bead in the lattice, so that the signal in Fig. 3 is recorded in actual experimental conditions. The acoustic response is very similar to the electrical excitation for all frequencies. The Gaussian envelope ensures a slow rise time of the electrical excitation, avoiding high-frequency resonances of the emitter. A finite duration signal allows time-of-flight measurements, as in [13] and Sec. IV B. The signal bandwidth is rather large, which enables time-frequency analysis, performed in Sec. V.

### III. COHERENT AND INCOHERENT WAVES: A HEURISTIC DEFINITION

Wave propagation occurs in the contact lattice, built upon the beads actually in contact. Although the bead center lattice is perfectly regular, the contact lattice is disordered. The increase of the applied stress induces deformation of the beads, compensating their radius inhomogeneities [13], and progressively ordering the contact lattice. This evolution is what we want to characterize.

The acoustic waves are recorded when the lattice is under a constant pressure. This stationary state will not change unless we vary the applied pressure. We thus must contend with a problem of wave propagation in a static random medium. By static, we mean that all characteristics of the random medium are constant in time, in such a way that a peculiar wave form is perfectly reproducible if we do not change the lattice state. When the state is changed, either by changing the static pressure applied to a given configuration, or by changing the configuration itself, the wave amplitude and phase fluctuate randomly. Those fluctuations are of statistical

<figure><img src="image_3.png" /><figcaption>FIG. 2. Left-hand side: Sketch of the emitter, with all dimensions in mm. Right-hand side: Emitter spectral response (in logarithm scale). We indicate the frequencies used in our experiments.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 3. From top to bottom, electrical (upper curve) and acoustic (lower curve) signal for Gaussian bursts of center frequencies 6.5, 15.2, and 19.5 kHz. The acoustic signal is recorded next to the emitter during the experiments, with the help of an *ad hoc* sensor held between the emitter and the lattice (see Fig. 1). The acoustic emitter preserves burst shape.</figcaption></figure>

nature only, and do not depend on time. The signal may be decomposed into the sum of a coherent wave, which is its ensemble average on many statistically independent realizations of the random medium, and of an incoherent wave which is its fluctuating part [29]. The problem being static, only ensemble averages make sense. Recent simulations of Somfaï *et al.* illustrate well those concepts [21]: Averaging wave forms on 1000 realizations, those authors exhibit most

clearly coherent wave propagation in a 2D random packing of spheres.

Experimentally, such a number of different realizations is out of reach. We show in Fig. 4 two series of four signals recorded after propagation in four different initial configurations, together with their ensemble average. Despite the insufficient statistics, much of the tail is indeed suppressed after averaging, while the very beginning of the signal remains almost unchanged. The first oscillations may thus be identified with the coherent signal and their time of arrival is used to define a time of flight and a velocity which is the group velocity of the coherent wave [12]. Another difficulty comes from the fact that our system is confined, and wave reflection on the boundaries cannot be avoided. In this respect, a part of the tail is not incoherent, but corresponds to coherent propagation after reflection on the boundaries. Thus, the incoherent wave cannot be simply identified with the tail of the signal.

In what follows, rather than averaging the very small number of experimental runs at our disposal, we have chosen to explicitly confront those configurations. We exhibit behaviors that do depend on the configuration, and behaviors that are observed whatever the configuration. It means that we must identify the coherent wave without resorting to ensemble averaging. This is done in the next section.

## IV. CROSS CORRELATION ANALYSIS

### A. Correlations: A better insight into the coherent signal

We may define the coherent signal in an approximate way, using cross correlations that are an efficient and sensitive tool to compare the shape of two signals. Since our signals are deterministic, the correlations defined in Appendix A do not imply time averaging [30]. Strictly speaking, they just

<figure><img src="image_5.png" /><figcaption>FIG. 4. Amplitude (arbitrary unit) of the acoustic signal, as a function of time (in µs), for a static force of (a) 154 N, (b) 1030 N. Center frequency is 6.5 kHz. Upper plots correspond to different initial lattice configurations, and the bottom curve is their average.</figcaption></figure>
<figure><img src="image_6.png" /></figure>

FIG. 5. Wave forms recorded in two statistically independent configurations (left-hand ordinate, arbitrary unit), for a static force of 200 N and a center frequency of 6.5 kHz, and their local cross correlation (right-hand ordinate) as a function of time (in µs). It is almost 1 for the first four oscillations, and significantly less after eight oscillations.

quantify the likeness of two deterministic signals.

In Fig. 5, we show two signals from different initial con-
figurations, under the same static force. We calculate local
cross correlations, restricting the signal duration to one
(pseudo) period around each maximum. This gives the dots
indexed by the position of the relevant maximum in Fig. 5.
This correlation is almost one for the first four maxima, then
decreases quickly. Those four oscillations thus constitute the
purely coherent part of the signal, corresponding to the direct
ballistic propagation of the emitted burst. We restrict the cal-
culation of the coherent correlation $C_{\text{coh}}$ (see Appendix A) to
this purely coherent part even if it does not include the entire
coherent signal. Whatever the experimental situation, this co-
herent correlation is always very large (see Fig. 6, Fig. 7, and
Ref. [13]), which confirms the consistency of this choice.

B. Time-of-flight measurements with low-frequency
waves

Cross correlation calculations provide a precise measure-
ment of the time of flight (TOF) (see Appendix A). This
method is valid for nondispersive waves, and we consider
here low-frequency waves only (6.5 kHz, thus a wavelength
of 10 bead diameters at least). We show in Fig. 6 the TOF as
a function of the static force, for three different initial con-

figurations. In [13], the reader may find data for another ini-
tial configuration [31]. According to the Hertz law, the sound
velocity in a perfectly regular lattice of elastic spheres under
a static force *F* scales as *F*¹/₆ [Ref. [7]], and Eqs. (B2) and
(B5) of Appendix B]. It has become customary to plot sound
velocity as a function of the confining force in logarithm-
logarithm scale, in order to test this Hertzian 1/6 power law.
At low force, experimental [6–10,12,13] or numerical
[19,20] data may rather be described by a non-Hertzian 1/4
power law (see [16]).

From Fig. 6 and [13], we see that the main features of the
lattice behavior at large and small static forces does not de-
pend on the initial configuration. The Hertzian behavior is
always recovered at large force, whereas the non-Hertzian
exponent describes the small force behavior quite well. High
force behavior confirms the relevance of the Hertz interac-
tion law for our beads, as previously shown in a 1D system
[17]. Low force behavior evidences that the randomness in
diameter and sphericity of our beads induces disorder of the
contact lattice. The tolerance level of the beads (±4 µm) is
consistent with their deformation on our force range: For a
force *F* on each moving side, a contact is submitted to
√3*F*/31 and the deformation of our steel beads is <1 µm,
≈3 µm, and ≈4.5 µm at respective forces of 100, 500, and
1000 N.

The force beyond which the behavior can be considered
as Hertzian is quite sensitive to the initial configuration. It
may be estimated to occur, from left to right of Fig. 6, at 200
N, 200 N, and 300 N. In the configuration of [13], it occurred
at 450 N. Those values are directly linked to structural modi-
fications of the contact lattice and thus depend on the initial
configuration.

C. Cross correlation for a unique configuration
(low-frequency waves)

To describe the lattice evolution with the force, we calcu-
late the cross correlation between the waves recorded at
force *F* and *F*+Δ*F*, respectively, for a small increment Δ*F*
[see Eq. (A2)]. This is done for each force, following the
evolution of a single initial configuration. As seen in Fig. 6,

<figure><img src="image_7.png" /></figure>

FIG. 6. Time of flight τ₀ (filled disks, left-hand ordinate, logarithm scale; in µs) and cross correlations C~ΔF~(τ₀) (right-hand ordinate, linear scale; open squares for coherent correlations; filled triangles for total correlations) versus the static force *F* on each boundary (abscissa, logarithm scale; in daN), for three different initial configurations and low-frequency (6.5 kHz) waves. The solid line is a *F*⁻¹/₆ fit, the dotted line is a *F*⁻¹/₄ fit. All data are taken during decompression of the lattice, except for the open disks on the right-hand plot that are taken during compression. The arrows indicate the points analyzed in Fig. 11.
<figure><img src="image_8.png" /><figcaption>FIG. 7. Correlations versus the force (in N) for compression and/or decompression of the same lattice $C^c(\bullet)$, and for two different initial configurations $C^\ddagger(\circ)$. Center frequency is 19.5 kHz. Upper (respectively, lower) curves for coherent (respectively, incoherent) waves, with a common scale. For the sake of comparison with low-frequency waves, see Fig. 5 of [13].</figcaption></figure>

the data are recorded at forces equally spaced on a logarithmic scale, hence in geometric proportion. The force increment is such that $\Delta F/F=1.14$.

We calculate either the coherent correlation $C_{\Delta F}^{\text{coh}}$, restricted to the coherent part of the wave (see Sec. III), or the total correlation including the signal tail. Those data are displayed in Fig. 6. The correlation $C_{\Delta F}^{\text{coh}}$ is always very high, whatever the initial configuration. As mentioned in Sec. IV A, this confirms the coherent nature of the first signal periods: From its very definition, the coherent wave should not depend on the lattice initial configuration, as observed.

On the contrary, the total correlation is configuration dependent. The only common feature observed in the three experiments of Fig. 6 is high total correlations (roughly 0.8) with small fluctuations at the highest available forces. But there are strong differences for the three data sets. The total correlation may still exhibit huge fluctuations, showing that strong modifications in the contact network still occur, while TOF measurements are consistent with an ordered lattice behavior (see center and right-hand plots of Fig. 6, and [13]). In theoretical studies [25] and numerical simulations [20,21,32], it is indeed shown that, at forces for which the velocity seems to exhibit the Hertzian power law, the contact lattice is actually far from being completely ordered, if one looks, e.g., for the mean number of effective contacts per grain. Cross correlations including the tail of the signal, and hence its incoherent part, are thus much more sensitive to the disorder than the velocity measurements.

### D. Cross correlation between independent configurations (high-frequency waves)

In this section, we study cross correlations between signals received at the same force for two different initial configurations (see Appendix A). Such measurements have already been discussed for low-frequency waves in [13]. Extending this study to high-frequency waves, we show in Fig. 7 the correlations as a function of the static force, for a

19.5 kHz wave. The correlation restricted to the coherent part of the wave (the first four periods at both center frequencies) is always very high, even at low force. This extends the consistency of our coherent wave definition to the high-frequency and/or short wavelength case.

The correlation restricted to the incoherent signal is very low, even at high force. It is interesting to compare those high-frequency results with previous low-frequency ones, Fig. 5 of [13]. At low frequency, the system exhibits strong evidence of ordering at high force, as seen in Secs. IV B and IV C. Above 500 N, the low-frequency correlation $C_{\text{incoh}}^c$ ranges between 0.6 and 0.75, whereas it is always less than 0.3 at high frequency. At low frequency, such a low correlation is observed only at very low force, when the behavior of the lattice is clearly non-Hertzian (see [13] and Sec. IV B). The high-frequency and/or short wavelength waves emphasize the effects of disorder and show that even in the Hertzian regime, the contact lattice is not perfectly ordered, as mentioned in Sec. IV C.

The randomness of the contact lattice is due to the inhomogeneities of the force distribution carried by the contacts, and obviously decreases when the overall applied force increases. The length scale of those force fluctuations can be estimated from photoelasticity observations [33,34] to a few grain diameters. The scattering cross section of such inhomogeneities should thus be small for long wavelength, but increases for high-frequency waves. The incoherent wave thus corresponds to a greater fraction of the wave energy at high frequency, which is evidenced here by the smallness of the correlation between the tails of high-frequency signals that propagated in different lattices.

### V. TIME-FREQUENCY ANALYSIS

#### A. Dispersion relation

In previous 1D experiments, we studied the propagation of linear and nonlinear waves in a chain of identical steel beads [17,18]. For the same range of static force as in the present study, we showed that the chain may be modeled by point masses connected by Hertzian springs. The present system is different because the 2D contact lattice may be random, whereas all contacts are obviously active in a chain. A way to estimate the randomness of the lattice is a comparison between actual wave propagation and theoretical predictions for the regular lattice.

The latter are recalled in Appendix B. We just mention here the fact that the dispersion relation (B1), the group velocity (B3) as well as the group delay (B4) depend on the force $F$ isotropically applied on the beads through the cutoff frequency $v_c(F)$ only,

$$v_c \approx 11.04F^{1/6} \text{ kHz (with } F \text{ in N).} \quad (1)$$

This formula gives the theoretical predictions for the dispersion relation and the group delay of a perfect lattice that we will use afterwards. The numerical coefficient is provided by former independent experiments [17], with no free parameter left. Group delay measurements require time-frequency analysis of the experimental data, since we must know the
<figure><img src="image_9.png" /><figcaption>FIG. 8. (Color online) Upper left-hand side: Electric input (left-hand side) and received signal (right-hand side) versus time in ms. Time-frequency analysis is done on the thick line signal. Lower left-hand side: Input (thin line) and selected received signal (thick line) power spectra versus frequency in kHz. Right-hand side: Group delay (in ms) versus frequency (in kHz). Solid line is Eq. (B4). The data are taken under a static force of 938 N, at a center frequency of 22.4 kHz. The corresponding lattice configuration is that of the left-hand plot of Fig. 6.</figcaption></figure>

time of arrival of each frequency component. The next section describes the method used for this analysis.

## B. Wigner-Ville transform

A thorough review of time-frequency methods may be found in [35], and in Ref. [36] Yen has considered the application of such methods to acoustical problems. Several experimental studies have demonstrated their relevance for problems involving dispersive waves in homogeneous [38,39] or random media [40].

We will expose here very briefly the relevant basic ideas. Assuming a time-dependent signal described by a complex function $s(t)$, with amplitude and phase information, its Wigner-Ville transform (WVT) is defined as

$$ 
\begin{aligned}
W_s(t, \omega) &= \int_{-\infty}^{\infty} e^{i\omega\tau} s\left(t + \frac{\tau}{2}\right) s^*\left(t - \frac{\tau}{2}\right) d\tau \\
&= \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{i\xi t} S\left(\omega + \frac{\xi}{2}\right) S^*\left(\omega - \frac{\xi}{2}\right) d\xi,
\end{aligned}
\quad (2) $$

where the superscript asterisk means complex conjugation and $S(\omega)$ is the Fourier transform of $s(t)$. The WVT may be thought of as an instantaneous power spectrum [36]. Confirming this picture is the property of time and frequency covariance,

$$ 
\begin{aligned}
s(t) \rightarrow s(t + t_0) \Rightarrow W_s(t, \omega) \rightarrow W_s(t + t_0, \omega), \\
s(t) \rightarrow e^{i\omega_0 t} s(t) \Rightarrow W_s(t, \omega) \rightarrow W_s(t, \omega + \omega_0),
\end{aligned}
\quad (3) $$

showing that a translation in time or frequency of the signal induces the same translation for its WVT. When applied to the received signal, this transform allows us to calculate the time of arrival of each frequency component, that is, its group delay (B4).

Being quadratic, the transform (2) has the very unpleasant feature that it may produce spurious interferences between different frequency components. A convenient way to get rid of those interferences is windowing the signal, defining a pseudo-Wigner-Ville transform (PWVT),

$$ PW_s(t, \omega) = \int_{-\infty}^{\infty} H(\tau - t, \xi - \omega) W_s(\tau, \xi) d\tau d\xi. \quad (4) $$

For a Gaussian window function

$$ H(t, \omega) = \frac{1}{2\pi\sqrt{\Delta T \Delta \Omega}} \exp\left[-\left(\frac{t}{\Delta T}\right)^2 - \left(\frac{\omega}{\Delta \Omega}\right)^2\right], \quad (5) $$

it may be shown that spurious interference terms are suppressed in the PWVT if $\Delta T \Delta \Omega > 1$, recovering in the time-frequency context the classical Heisenberg's uncertainty relation [36].

In what follows, the signal analysis is undertaken with the MATLAB toolbox provided by Flandrin and co-workers [37].

## C. Time-frequency analysis: Experimental results

The analysis is restricted to the purely coherent part of the signal, as defined above (see Sec. IV A). This is illustrated in Fig. 8. Because of dispersion, the burst wave broadens, and its power spectrum shrinks. We also see that the maximum of the power spectrum is not at the same frequency for the electrical excitation and the received signal. The lattice behaves like a low-pass mechanical filter, so that the received signal is centered on a somewhat lower frequency than the emitted one.

We show in Fig. 9 the analysis of signals received after propagation in the same lattice, at four static forces, during its decompression. We plot the amplitude of the normalized PWVT in the time-frequency plane. The time origin is taken at the beginning of the emitted (acoustic) signal. Thus the color of a pixel (increasing from light yellow for null energy, toward dark blue) at a position $(t_0, \nu_0)$ in the time-frequency plane describes the amount of acoustic energy at a frequency $\nu_0$ that arrived on the sensor at time $t_0$. It is seen in Fig. 9 that the signal is correctly localized in the time-frequency plane, as it should, and that windowing has suppressed spurious interference terms. The analysis is done on high-frequency waves. The initial lattice configuration was that of the center plot in Fig. 6.
<figure><img src="image_10.png" /><figcaption>FIG. 9. Normalized amplitude of the Wigner-Ville transform (2) (gray scale, decreasing from black to white) in the time-frequency plane (time in ms, frequency in kHz), for a single configuration (that corresponds to the center plot of Fig. 6). The solid lines are the theoretical prediction (B4) and the cutoff frequency. The relevant forces are indicated in the upper left-hand corners. The center frequencies of the Gaussian burst waves are, from left to right and top to bottom, 13.6, 17.8, 21.5, and 24.3 kHz.</figcaption></figure>

At the highest forces, the density plot in Fig. 9 indicates that for short wavelength waves the lattice is well ordered: Each frequency component arrives at the expected theoretical group delay deduced from (B3), (B4), and (1). We remind the reader that there is no adjustable parameter. Under a force of 154 N, for which time-of-flight measurements exhibit slightly non-Hertzian behavior, there is a clear discrepancy, mostly noticeable for the highest frequency components. At the lowest force, this discrepancy is also clearly observable for low-frequency components. Although such representations are quite illustrative, we will rather discuss group delay measurements that allow a more quantitative analysis.

## D. Group delay measurements

### 1. Group delay as a function of the applied force

From its very definition (B4), the group delay calculations require time-frequency analysis of both the received and emitted signals. Experimentally, the emission and reception times of a frequency component $v$ are defined as the position of the maximum of the relevant PWVT on the time axis, at

frequency $v$. To improve accuracy, we restrict the analysis to frequency components with high enough amplitude, dropping the data corresponding to less than 25% of the frequency spectrum maximum. We display in Fig. 10 the results for waves of center frequencies 6.5 and 19.5 kHz, and waves of constant wavelength $\lambda=3$ cm. The lattice configuration is the same as in Fig. 9.

At low frequency, there is no dispersion, the group velocity is equal to the sound velocity, and the group delay should be equal to the TOF. This is indeed observed for the three forces corresponding to the Hertzian regime. As for the non-Hertzian behavior seen in Fig. 6, we observe at 46 N a low-frequency group delay in the [990, 1150] $\mu$s range, in fair agreement with the TOF measured at the same force (980 $\mu$s). It is not a surprise, but it confirms the results of Sec. IV B and our previous work [13] with completely different data analysis.

For high-frequency data, there is an excellent agreement between the experimental high-frequency group delay, and theoretical predictions at 817 N. At 412 N some discrepancy is evidenced at high frequency, which was hardly seen in Fig. 9: The lattice cannot be considered as ordered at length scales as low as 2.5d (see Table II). This may be understood
<figure><img src="image_11.png" /><figcaption>FIG. 10. (Color online) Experimental group delay, calculated from time-frequency analysis, and theoretical prediction for a perfectly ordered lattice, Eq. (B4). Same configuration as in Fig. 9, under the same static forces, for Gaussian burst waves of the same center frequencies (crosses). Also shown are group delay data for Gaussian burst waves of center frequencies 19.2 kHz (open squares; no available data at 46 N) and 6.5 kHz (open circles).</figcaption></figure>

from the analysis of Thorpe et al. [41–43]. In a model of identical point masses and linear springs, initially arranged in a perfect triangular lattice, they study the evolution of the vibrational density of state when springs are progressively cut. In our experiments, when the confining force is decreased, the active contacts rather exhibit an increasing randomness in their elastic properties than a net decrease of their number. However, both systems are quite similar and may adopt the same qualitative behavior. Randomness is quantified by the fraction *p* of unbroken springs [41]. For a triangular lattice, the rigidity threshold is about *p*=0.65 [44,45]. With *p*=0.85, the peaks in the histograms of a phonon density of states of a regular lattice almost disappeared (see Fig. 4 of [41]). Those peaks are called van Hove singularities, and trace back to the zeros of the group velocity. Their disappearing indicates that the variations of group velocity, if not suppressed, are greatly reduced by even a small amount of randomness. This is what we observe in Fig. 10 at 412 N, and in Fig. 11 at somewhat smaller forces, where it is shown that the experimental group delay variation is less than the theoretical prediction.

At 154 N, the high-frequency waves are strongly filtered; the received signal spectra extends up to 20 kHz, whereas the emitted signal contains frequencies up to 28 kHz. This effect is reinforced at 46 N, where the highest frequency waves are almost completely filtered.

## 2. Group delay versus correlations

In [13] and in Sec. IV C, we argue that the correlations $C_{\Delta F}(\tau_0)$ are more sensitive to disorder than velocity measurements. It is clearly seen in the Hertzian regime, in which there is still an evolution of the correlation (see Fig. 6, center and right-hand plots). From its very definition, Eq. (A2),

<figure><img src="image_12.png" /><figcaption>FIG. 11. (Color online) Same as Fig. 10 (open squares and open circles), for configurations showing very high (left-hand side) and very low (right-hand side) cross correlation at low frequency (indicated by the arrows in the center plot of Fig. 6). The relevant forces are indicated in the upper left-hand corners. The crosses correspond to Gaussian burst waves of respective center frequencies 19.4 kHz (left-hand side) and 21.0 kHz (right-hand side).</figcaption></figure>

$C_{\Delta F}(\tau_0)$ measures the likeness of the contact lattices at two forces that differ by a small increment $\Delta F$. It is thus not a measurement of the absolute level of randomness of the contact lattice. Group delay measurements at a frequency $\nu$, on the contrary, provide an absolute test of the level of randomness at the scale of wavelength $\lambda(\nu)$.

To demonstrate that, we compare two different initial configurations at (roughly) the same static force. We first consider the configuration of the center plot of Fig. 6, for two states that differ greatly in the low-frequency correlation values (indicated by arrows). At a static force of 231 ± 1 N the correlation is 0.8, and at 310 ± 1 N it is less than 0.2. There has been a great structure modification of the contact lattice between those two measurements. The group delay plots for those two forces shown in Fig. 11 exhibit similar discrepancy at high frequency, which means that none of them correspond to an ordered state of the contact lattice. Then we consider the configuration of the left-hand plot of Fig. 6, at typically the same forces (231 ± 1 and 328 ± 1 N). The correlation is almost constant and very high (about 0.8), which means that the contact lattice does not evolve much. We may thus expect that, for this peculiar configuration of the lattice, the contact lattice is ordered at much lower forces than in the two other cases of Fig. 6. This is indeed the case if we look at the group delay measurements for those data, shown in Fig. 12. The agreement with the theoretical prediction is very good at high frequency. We can notice that, in this configuration, the agreement is even better than in the configuration of Fig. 10, at a force of 412 N.

The only way to estimate the level of randomness from correlation $C_{\Delta F}(\tau_0)$ is to consider their variations rather than their absolute value. Since they exhibit large fluctuations, this is very difficult. Hence group delay measurements are a much more reliable tool to this end. However, from the perspective of nondestructive testing, correlations can be an interesting and very precise tool to detect slight modifications of the structure of a granular medium under high load. On a fundamental point of view, disorder may be evidenced by the behavior of either incoherent long waves, which are involved in our correlation measurements, or coherent short waves analyzed with time-frequency methods. Physically, one must compare the travel distance of the wave to its wavelength. Incoherent waves are multiply scattered, in such a way that
<figure><img src="image_13.png" /></figure>

FIG. 12. (Color online) Same as Fig. 10 (open squares and open circles), for the lattice configuration corresponding to the left-hand plot of Fig. 6, which shows almost constant cross correlation at low frequency. The relevant forces are indicated in the upper left-hand corners. The crosses correspond to Gaussian burst waves of respective center frequencies 19.6 kHz (left-hand side) and 20.5 kHz (right-hand side). This figure is to be compared to Fig. 11: The applied force is (roughly) the same, but the agreement with the regular lattice dispersion curve is much better.

their total travel length may become large in comparison to the wavelength. The emitter-receiver distance is sufficient for shorter coherent waves to probe the lattice disorder on their ballistic flight distance only.

## VI. REVERSIBILITY OF LATTICE EVOLUTION

Up to now, we have consistent evidence of the progressive lattice ordering when the applied force increases. A question of interest is the physical process leading to this ordering. A first possible mechanism is the failure of the force network at some static forces and its rearranging into more ramified structures; an example of such a process is the particle chain buckling proposed by Goddard [15]. Another type of mechanism, proposed by Roux [20,32], is progressive closing of opened contacts through elastic deformation of the beads. The latter, at the difference of the other one, is reversible. Hence the question of the reversibility of the lattice behavior is of fundamental importance.

The time-of-flight measurements characterize coherent wave propagation, hence ballistic flight between the emitter and the receiver, of low-frequency waves (wavelength about 10 bead diameters $d$). They are thus our less sensitive probe. The right-hand plot of Fig. 6 shows data obtained during both compression and decompression of the same lattice configuration. The same experiment, done on another initial configuration of the lattice, is shown in [13], and reveals similar features: At a given force, there is good agreement between time-of-flight measurements done during either compression or decompression, and no irreversibility is observed on coherent propagation at length scales of about 10 bead diameters.

Time-frequency analysis gives information on the group delay for high-frequency coherent waves. We plot in Fig. 13 group delay data for the same initial configuration as in the preceding section (right-hand plot of Fig. 6). We consider three forces, high enough for low-frequency time-of-flight measurements to follow Hertzian behavior. The group delay is measured during either compression or decompression of the lattice, for waves of central frequency 19.5 kHz. It is

<figure><img src="image_14.png" /></figure>

FIG. 13. (Color online) Same as Fig. 10, for compression (open squares) or decompression (open circles) of the same lattice configuration (corresponding to the right-hand plot of Fig. 6). Static forces in the upper left-hand corner. The center frequencies of the Gaussian burst waves are 19.5 kHz in each case.

shown in Table II that the wavelength goes down to 2.5 bead diameters. When the static force increases, the agreement between experimental group delay and theoretical expectations for a regular lattice is better at high frequency, as already observed in Fig. 10. The important point here is that even for the lowest force and highest frequency, there is very good matching between data recorded during compression and decompression. Thus ballistic propagation of coherent waves does not show any irreversible behavior of the contact lattice even for wavelength down to $2.5d$.

Those results seem to rule out any mechanism linked to irreversible rearranging of the network such as particle chain buckling. This result is in good agreement with recent simulations of frictional grain packings [24], in which the system response is quite sensitive to the symmetry of the loading: For isotropic compression as the one we use, the simulations are consistent with reversibility of the piling behavior, at forces which are of the same order of magnitude as in our experiments. This is not true when deviatoric stress exists.

However, correlation measurements based on the incoherent part of high-frequency waves reveal some irreversible behavior of the lattice: For low-frequency waves, we showed in Fig. 5 of [13] that the correlation based on the incoherent part of the signals was greater than 0.6 between waves recorded at the same static force, during either compression or decompression of the lattice. The same measurement done with high-frequency waves (typical wavelength about three bead diameters) never gives correlations higher than 0.3, as shown in Fig. 7. This brings to the fore a small irreversibility in the evolution of the contact lattice, which can only be seen using an incoherent part of high-frequency waves, which is much more sensitive to small changes in the lattice, even at very short length scales. Unavoidable changes in mobilization states of grain-grain contact friction between compression and decompression of the same lattice, which takes place at short length scales, could explain low correlations between incoherent parts of short wavelength signals, and would have only a weak influence on group delay measurements.

## VII. CONCLUSIONS

In this paper, we report on what may be the simplest nontrivial example of a granular medium: An array of nominally identical spheres, placed at the nodes of a regular tri-
angular lattice, under isotropic stress. Because of unavoidable scattering in the bead dimensions, the lattice of actually active contacts displays some stress-dependent randomness. This is evidenced by the observation of both a low stress non-Hertzian regime, where the velocity roughly scales as a power 1/4 of the stress, and a high stress Hertzian regime, with exponent 1/6. One of the main objectives of this paper was to go further than this classical picture, and to characterize the level of randomness that persists in the Hertzian regime.

In addition to our previous study [13], we provide results for several different initial configurations, and exhibit typical behaviors, such as Hertzian scaling at high stress, and non-Hertzian scaling at lower ones. On the contrary, the crossover between the two regimes takes place at configuration-dependent forces, ranging from 200 to 450 N. The correlation $C_{\Delta F}(\tau_0)$ [see (A2)] measures the likeness of the contact lattice after a ≈10% force increment $\Delta F$. When it involves the incoherent signal, it evidences lattice ordering [$C_{\Delta F}(\tau_0) > 0.85$] at the highest forces only (more than 800 N). Even in the Hertzian regime, it fluctuates in the [0.2,0.9] range at lower forces, in a manner that greatly depends on the initial lattice configuration. Cross correlation of high-frequency waves ($\lambda \approx 2.5d$, to be compared to the cutoff wavelength of $\sqrt{3}d$ in units of bead diameter $d$) for two different initial lattice configurations is very low, which shows significant randomness at short length scales.

Using time-frequency analysis, we access the lattice dispersion relation. At low frequency, it is fully consistent with previous measurements. For high-frequency waves, at the highest forces (greater than 800 N), there is excellent agreement with the theoretical dispersion relation without any free parameter. This evidences lattice ordering at the wavelength scale (2.5d). Decreasing the force, while staying in the Hertzian regime for long waves, we exhibit more and more discrepancy at high frequency. The actual lattice is less dispersive than a regular one. Such behavior was already observed in numerical simulations of a closely related problem, that of rigidity percolation in a triangular lattice [41].

Measurements for high-frequency coherent waves are consistent with correlation calculations of low-frequency incoherent waves. High and constant (respectively, low and fluctuating) correlation corresponds to good (respectively, poor) agreement with the theoretical dispersion relation. The contact lattice becomes more and more ordered when the applied force is increased, and the less the wavelength, the higher the force required to recover the regular lattice behavior. The force range that allows observation of both randomness at low force, and ordering at higher force, is fully consistent with the width of the diameter distribution of our beads.

Measurements were also done during compression and decompression of the same lattice. For the dispersion relation analysis, we see a good agreement between both data sets. On the contrary, correlations between high-frequency incoherent waves recorded during compression and decompression of the same lattice are very small, which evidences some irreversibility in the lattice evolution under stress modifications. The typical length scale involved here (a few beads) leads us to attribute it to frictional effects, which are

presumably moderate since the applied stress is isotropic.

ACKNOWLEDGMENTS

The authors acknowledge useful talks with C. Caroli, P. Flandrin, and J.-N. Roux.

APPENDIX A: CORRELATIONS

For two signals $s_i(t)$, $i=1,2$, recorded with the same sampling period at instants $t_n, n \in [1, N]$, we define the cross correlation as

$$C(\tau) = \frac{\sum_{n=1}^{N} s_1(t_n)s_2(t_n + \tau)}{\sqrt{\sum_{n=1}^{N} s_1^2(t_n)} \sqrt{\sum_{n=1}^{N} s_2^2(t_n)}}. \quad (A1)$$

In this calculation, time averaging is not performed. Such a process is well suited for deterministic signals [30] that do not include a stochastic component. This is the case of the signals that are recorded during our experiments, since the randomness of the lattice is static, and does not evolve during sound propagation.

When $\tau$ is the delay between signals $s_1$ and $s_2$, $C(\tau)$ is a measure of the likeness of the received wave forms. This tool is used to compare either two signals received at the same static force for different realizations of the experiment, or two signals obtained before and after a small increment of the static force during the same realization of the experiment (see below).

In the first case, the delay is almost zero, and $C(0)$ is taken to compare wave forms. We may define $C_{\text{coh}}(0)$ and $C_{\text{incoh}}(0)$, respectively, restricted to the coherent and incoherent parts of the received signal in the sense of Secs. III and IV A. We also define $C^\ddagger$, if both signals correspond to two different initial configurations, and $C^\ddagger$ if they correspond to the same lattice under compression and decompression. Those definitions make sense if the two signals are recorded at the same static force.

Another use we make of cross correlations is to follow the evolution of a given lattice when the static force is modified. If $\Delta F$ is the increment between two successive measurements at forces $F$ and $F+\Delta F$, we set

$$C_{\Delta F}(\tau_0) = \frac{\sum_{n=1}^{N} s_F(t_n) s_{F+\Delta F}(t_n + \tau_0)}{\sqrt{\sum_{n=1}^{N} s_F^2(t_n)} \sqrt{\sum_{n=1}^{N} s_{F+\Delta F}^2(t_n)}}. \quad (A2)$$

Here $\tau_0$ is the time difference for the arrival of the two signals. Here again, $C_{\Delta F}(\tau_0)$ is a quantitative measure of the likeness of wave forms received at forces $F$ and $F+\Delta F$. Moreover, this provides a very precise way to measure $\tau_0$: We first calculate the correlation $C_{\Delta F}^{\text{coh}}(\tau_0)$ restricted to the coherent wave, taking $\tau_0$ as a fitting parameter. $\tau_0$ is then chosen to maximize the coherent correlation. The time-of-
<figure><img src="image_15.png" /><figcaption>FIG. 14. Sketch of the triangular lattice, showing axes orientation.</figcaption></figure>

flight data of Sec. IV B and [13] are obtained in this way. This method relies on the assumption of no dispersion, and since the incident wave is not monochromatic, it may be applied to low-frequency waves only. For high-frequency waves, we use time-frequency analysis that gives a direct estimate of the group delay (see Sec. V D).

## APPENDIX B: TRIANGULAR LATTICE

In a perfect triangular lattice under isotropic stress, all contacts between the beads are effective and bear the same force $F_0$. The distance $d$ between adjacent bead centers becomes under compression $d-\delta d$, where according to Hertz [5] the distance of approach $\delta d$ depends as a power law on the force: $\delta d \propto F_0^{2/3}$. The sound wave sent in the lattice induces a very small perturbation of the confining pressure, so that the contacts respond like linear springs. Their stiffness $\alpha$ depends on the confining pressure, $\alpha=(\partial\delta d/\partial F)^{-1} \propto F_0^{1/3}$ [see (B6)]. In this simple picture, the lattice may be modeled by point masses $m=\pi\rho d^3/6$, where $\rho$ is the volume mass, linked by springs of stiffness $\alpha$ and length $d$.

A sketch of the triangular lattice in Fig. 14 defines the axes orientation. We are concerned with compressional waves, with polarization along the y axis. According to Fig. 1, the direction of propagation from the emitter to the receiver is also along the y axis. In that case, the dispersion relation can be obtained from the one provided by Garboczi and Thorpe [41] for a perfectly ordered central-force elastic network, and takes the simple form

$$ \omega^2 = \frac{6\alpha}{m} \sin^2 \frac{\sqrt{3}}{4} kd, \qquad (B1) $$

where $\omega$ is the pulsation and $k$ is the wave number norm. It is convenient to introduce the cutoff frequency $v_c = (6\alpha/m)^{1/2}/(2\pi)$. The compressional sound velocity $c_s$ reads as

$$ c_s \equiv \lim_{k \to 0} \frac{\omega}{k} = \pi(\sqrt{3}/2)dv_c, \qquad (B2) $$

and the group velocity $v_g$,

$$ v_g(\nu) = c_s \cos[\arcsin(\nu/v_c)]. \qquad (B3) $$

We will also use the group delay $T_d(\nu)$, which is the time of arrival of a frequency component $\nu$ after propagation along a

<table><caption>TABLE II. Estimates of the shortest wavelength available in the experiments reported in Figs. 8 and 10–13.</caption><thead><tr><th>Force (N)</th><th>Figure</th><th>v<sub>max</sub> (kHz)</th><th>λ(v<sub>max</sub>)/d</th></tr></thead><tbody><tr><td>938</td><td>8</td><td>29.4</td><td>2.7</td></tr><tr><td>46</td><td>10</td><td>15.7</td><td>3.2</td></tr><tr><td>154</td><td>10</td><td>19.5</td><td>3.1</td></tr><tr><td>412</td><td>10</td><td>25.6</td><td>2.7</td></tr><tr><td>817</td><td>10</td><td>27.3</td><td>2.9</td></tr><tr><td>231</td><td>11</td><td>23.3</td><td>2.7</td></tr><tr><td>310</td><td>11</td><td>23.6</td><td>2.8</td></tr><tr><td>231</td><td>12</td><td>21.5</td><td>3</td></tr><tr><td>328</td><td>12</td><td>23.4</td><td>2.9</td></tr><tr><td>329</td><td>13</td><td>25.6</td><td>2.5</td></tr><tr><td>588</td><td>13</td><td>25.6</td><td>2.9</td></tr><tr><td>1050</td><td>13</td><td>25.6</td><td>3.3</td></tr></tbody></table>

distance $L$. It is thus the difference between the reception time $T_r(\nu)$ of frequency component $\nu$, and its emission time $T_e(\nu)$, and is given by

$$ T_d(\nu) \equiv T_r(\nu) - T_e(\nu) = L/v_g(\nu). \qquad (B4) $$

All dependency on the confining pressure is contained in the cutoff frequency, which depends on the force $F_0$ between adjacent beads through the spring constant $\alpha$. The latter may be derived from the Hertz law [5], and the cutoff frequency reads as

$$ v_c \equiv \frac{1}{2\pi} \sqrt{\frac{6\alpha}{m}} = \frac{1}{4\pi} \left(\frac{27}{\pi\rho}\right)^{1/2} \left(\frac{16\sqrt{2E}}{3(1-\sigma^2)d^4}\right)^{1/3} F_0^{1/6}. \qquad (B5) $$

Here $E$ is the Young modulus, and $\sigma$ the Poisson ratio of the bead material (here stainless steel). Let $F$ be the force applied on each moving piston. Then $F_0$ is given by

$$ F_0 = 2 \frac{F}{N} \cos \frac{\pi}{6}, \qquad (B6) $$

where $N=31$ is the number of beads along one side of the hexagonal cell (see Fig. 1). We obtain from a 1D experiment a precise measurement of the sound velocity in a chain of beads [17], which gives the numerical value (1).

From (B1), we obtain an estimate of the shortest measurable wavelength in a given experiment. Let $v_{max}$ be the highest frequency sent in the lattice. In units of the bead diameter $d$, the corresponding wavelength $\lambda(v_{max})$ is

$$ \frac{\lambda(v_{\text{max}})}{d} = \frac{\pi\sqrt{3}}{2} \left[ \arcsin\left(\frac{v_{\text{max}}}{v_c}\right) \right]^{-1}, \qquad (B7) $$

showing that the cutoff wavelength is $\sqrt{3}d$. The relevant numerical values are provided for all experiments in Table II. We insist on the fact that they are just rough estimates, since the calculation assumes perfect agreement with the dispersion relation (B1), which is not correct except at the highest confining pressures.
[1] D. W. Howell, R. P. Behringer, and C. T. Veje, *Chaos* **9**, 559 (1999).

[2] D. M. Mueth, H. M. Jaeger, and S. R. Nagel, *Phys. Rev. E* **57**, 3164 (1998).

[3] F. Radjaï, S. Roux, and J. J. Moreau, *Chaos* **9**, 544 (1999).

[4] J.-N. Roux, *Phys. Rev. E* **61**, 6802 (2000).

[5] L. D. Landau and E. M. Lifshitz, *Theory of Elasticity* (Pergamon, Oxford, 1986).

[6] K. Iida, Bull. Earthquake Res. Inst., Univ. Tokyo **16**, 131 (1938).

[7] J. Duffy and R. D. Mindlin, *J. Appl. Mech.* **24**, 585 (1957).

[8] C. W. Thurston and H. Deresiewicz, *J. Appl. Mech.* **26**, 251 (1959).

[9] P. Talwani, A. Nur, and R. L. Kovach, *J. Geophys. Res.* **78**, 6899 (1973).

[10] N. Warren and O. L. Anderson, *J. Geophys. Res.* **78**, 6911 (1973).

[11] C. H. Liu and S. R. Nagel, *Phys. Rev. Lett.* **68**, 2301 (1992).

[12] X. Jia, C. Caroli, and B. Velicky, *Phys. Rev. Lett.* **82**, 1863 (1999).

[13] B. Gilles and C. Coste, *Phys. Rev. Lett.* **90**, 174302 (2003).

[14] I. Agnolin, J.-N. Roux, P. Massaad, X. Jia, and P. Mills, in *Powders and Grains 05*, edited by R. Garcia-Rojo, H. J. Herrmann, and S. McNamara (Balkema, Leiden, 2005).

[15] J. D. Goddard, *Proc. R. Soc. London, Ser. A* **430**, 105 (1990).

[16] A word of caution is useful here. As shown in the theoretical analysis [25] the functional form of the effective compressibility has no reason to be a power law. This representation is only a convenient way to emphasize non-Hertzian behaviors, and the 1/4 "exponent" should not be taken too seriously.

[17] C. Coste and B. Gilles, *Eur. Phys. J. B* **7**, 155 (1999).

[18] C. Coste, E. Falcon, and S. Fauve, *Phys. Rev. E* **56**, 6104 (1997).

[19] H. J. Herrmann, D. Stauffer, and S. Roux, *Europhys. Lett.* **3**, 265 (1987).

[20] J.-N. Roux, in *Powders and Grains 97*, edited by R. P. Behringer and J. T. Jenkins (Balkema, Rotterdam, 1997).

[21] E. Somfaï, J.-N. Roux, J. H. Snoeijer, M. van Hecke, and W. van Saarloos, *Phys. Rev. E* **72**, 021301 (2005).

[22] J.-N. Roux, e-print arXiv:cond-mat/0405358.

[23] I. Agnolin and J.-N. Roux, in *Powders and Grains 05*, edited by R. Garcia-Rojo, H. J. Herrmann, and S. McNamara (Balkema, Leiden, 2005).

[24] I. Agnolin and J.-N. Roux, *Phys. Rev. E* **76**, 061304 (2007).

[25] B. Velicky and C. Caroli, *Phys. Rev. E* **65**, 021307 (2002).

[26] P. Langevin, *Oeuvres Scientifiques de Paul Langevin* (CNRS, Paris, 1950).

[27] S. R. Rich, U.S. Patent No. 2,947,889 (1960).

[28] L. G. Merkulov, *Sov. Phys. Acoust.* **3**, 246 (1957).

[29] A. Ishimaru, *Wave Propagation and Scattering in Random Media* (IEEE, New York, 1997).

[30] T. J. Cavicchi, *Digital Signal Processing* (Wiley, New York, 2000).

[31] To be consistent with Sec. IV D, we plot the time of flight rather than the sound velocity. Since the distance (42.95 cm) between emitter and receiver is held constant, the conversion is easily done.

[32] S. Ouaguenouni and J.-N. Roux, *Europhys. Lett.* **39**, 117 (1997).

[33] A. Drescher and G. de Josselin de Jong, *J. Mech. Phys. Solids* **20**, 337 (1972).

[34] T. Travers, D. Bideau, A. Gervois, J.-P. Troadec, and J.-C. Messager, *J. Phys. A* **19**, L1033 (1986).

[35] P. Flandrin, *Time-Frequency/Time-Scale Analysis* (Academic, New York, 1998).

[36] N. Yen, *J. Acoust. Soc. Am.* **81**, 1841 (1987).

[37] F. Auger, P. Flandrin, P. Gonçalves, and O. Lemoine, *Time-frequency toolbox*, 1995. The toolbox and its tutorial may be downloaded at http://tftb.nongnu.org/. We refer the reader to this tutorial for the details on the signal processing.

[38] K. Kishimoto, H. Inoue, M. Hamad, and T. Shibuya, *J. Appl. Mech.* **62**, 841 (1995).

[39] W. H. Prosser, M. D. Seale, and B. D. Smith, *J. Acoust. Soc. Am.* **105**, 2669 (1999).

[40] H. Jeong and Y.-S. Jang, *Compos. Struct.* **49**, 443 (2000).

[41] E. J. Garboczi and M. F. Thorpe, *Phys. Rev. B* **32**, 4513 (1985).

[42] Shechao Feng, M. F. Thorpe, and E. J. Garboczi, *Phys. Rev. B* **31**, 276 (1985).

[43] L. M. Schwartz, Shechao Feng, M. F. Thorpe, and P. N. Sen, *Phys. Rev. B* **32**, 4607 (1985).

[44] J. G. Zabolitzky, D. J. Bergman, and D. Stauffer, *J. Stat. Phys.* **44**, 211 (1986).

[45] D. Stauffer and A. Aharony, *Introduction to Percolation Theory* (Taylor and Francis, London, 1994).