Micro-slips in an experimental granular shear band replicate the spatiotemporal characteristics of natural earthquakes

Memory effects in seismology—such as the occurrence of aftershock sequences—are implicitly assumed to be governed by the time since the main event. However, experiments are yet to identify if memory effects are structural or time-dependent mechanisms. Here, we use laser interferometry to examine the fluctuations of deformation which naturally emerge along an experimental shear fault within a compressed frictional granular medium. We find that deformation occurs as a succession of localized micro-slips distributed along the fault. The associated distributions of released seismic moments, as well as the memory effects in strain fluctuations and the time correlations between successive events, follow exactly the empirical laws of natural earthquakes. We use a methodology initially developed in seismology to reveal at the laboratory scale the underlying causal structure of this behavior and identify the triggering kernel. We propose that strain, not time, controls the memory effects in our fault analog. Deformation experiments using glass beads generate nano-earthquakes that mimic the spatio-temporal pattern of strain in natural seismicity and suggest memory effects associated with strain fluctuations are strain-dependent rather than time-dependent.

E arthquakes are natural phenomena displaying scale-free statistics 1 . Empirical power laws are observed for the distribution of their moments (Gutenberg-Richter's law), rupture lengths and durations, rupture slips 2 , temporal 3 and spatial correlations between earthquakes 4 , which also express through a decaying rate of aftershocks (Omori's law) 5 , characterized by a scale-free (sub)diffusion 6,7 . Understanding the origin of those laws as well as reproducing them at the laboratory scale remains nowadays major issue. From a fundamental point of view, these scaling laws are reminiscent of nonequilibrium dynamics and critical phenomena 8,9 and raise the question of the existence of a universality class to which earthquakes would belong. Among mechanical systems, possible candidates for such a universality class are those for which deformation occurs through avalanches. The mechanical response of those systems is characterized by intermittent dynamics alternating slow elastic loading and rapid sliding and relaxation, leading to jerky dynamics and/or stress drops.
Power-law distributions of slip sizes or relaxed energies have been evidenced experimentally in various systems, such as the stationary propagation of a fractured front in a heterogeneous material 10 , compression of heterogeneous materials 11,12 , or systems involving frictional sliding on or within a granular media 13,14 . The common basic ingredients underlying the dynamics of those different systems are the existence of the material disorder and the decomposition of the dynamics in elementary events localized both in space and time, coupled together by elasticity. Progressive evolution of avalanche size and duration statistics has been reported for different heterogeneous materials 12 or granular media 15 upon increasing the loading up to a macroscopic yield or failure stress at which scale-free statistics are observed, arguing for a "stress-tuned" critical behavior fundamentally different from a self-organized critical dynamics characterized by steady-state statistics 16 . It has been proposed on the basis of a mean-field model of plasticity that those different systems, as well as deformed microcrystals and earthquakes, could belong to the same class of universality 16,17 . This, however, was mainly addressed from an analysis of the size distribution of avalanches as well as their average shape. In the case of earthquakes, the stress-tuned critical hypothesis was argued 16 on the basis of dependence of magnitude distributions on the slip direction on the fault plane (the rake angle), which gives indirect information about the differential stress acting on the fault 18 . However, if seismic moment distributions appear to be indeed exponentially tapered at very large scales, the associated upper corner magnitude was found to be independent of the region or the depth interval considered, or of the plate velocity, i.e., to be rather "universal" 19 . On the other hand, an important and ubiquitous feature of brittle deformation in the crust is the existence of aftershocks, which occurrences are also governed by scale-free laws. Much fewer attempts have been done to model those spatio-temporal correlations, through the introduction of a memory mechanism such as a slow healing of frictional properties 20 , or viscoelastic relaxation 21 . Similarly, lab experiments reproducing the clustering of events in time and space remain scarce.
As a matter of fact, systems displaying avalanches can have also fundamental differences that limit the pertinence of a universal picture. As already noted, systems in a stationary regime must be distinguished from those for which the spatio-temporal dynamics are evolving (stress-tuned). In the case of frictional granular media, besides a progressive increase of the maximum avalanche size 15 , the spatial distribution of plastic events evolves during loading: initially homogeneously distributed in the bulk of the material, plasticity progressively localizes to form shear bands at macroscopic yield [22][23][24] , in which all the shear rate concentrates afterward while the rest of the system becomes an elastic "solid". Identifying clustering in time is only possible within those shear bands when the spatio-temporal organization becomes stationary. For stationary systems, the dimensionality of the active zone is expected to play a role, at least on the value of the critical exponents. One must thus distinguish tri-dimensional systems (e.g., plasticity distributed in the bulk of an amorphous material), from those where the plasticity is confined to a quasi-2D zone (a fault in the case of earthquakes, a shear band in the case of amorphous granular media), and finally quasi uni-dimensional active zones (e.g., a propagating crack front). Practically, in experimental works pertaining to the plasticity of amorphous media, it is not always clear whether the plasticity is broadly distributed in the bulk of the system or if it is localized along a shear band. Even when the geometry of the active zone is identified, most experimental set-ups are unable to fully resolve the spatio-temporal organization of the avalanches which are solely identified and studied through indirect measurements of their sizes such as acoustic emissions (e.g., ref. 14 ) or stress drops on a loading curve (e.g., ref. 15 ).
To address the challenging issue of reproducing an analog of a fault gouge at the lab scale, a straightforward approach consists of imposing the bi-dimensional geometry in stationary conditions by confining a granular material between elastic plates 14,[25][26][27] . In the vast majority of those experimental studies, quasi-periodical stick-slip events with a typical size are observed, indicating that finite-size effects dominate the dynamics 28 . The slip events then involve the whole length of the shear band and the dynamics lose their universal features. In addition, those macro-slip events are characterized by a reverse asymmetry of the activity compared to earthquakes, with foreshocks of increasing size as approaching the macro-instability, but an absence of aftershocks 14 , likely resulting from the finite-size effects mentioned above. A recent experiment with a quasi-2D shear cell in a stationary regime displayed intermittent dynamics sharing several features of earthquake dynamics, such as the G-R law and a power-law decay of the rate of aftershocks 27 . While giving promising results in terms of the analogy between the dynamics of stationary sheared granular materials and that of earthquakes, it did not give a direct characterization of the localization and the spatial extension of the detected events. Consequently, the aftershock characterization amounts to a time correlation analysis of discrete events, without quantifying the underlying causal triggering. We can thus ask: is it possible to build a laboratory analog of a fault gouge where well-identified events would share all the properties of earthquakes, and more particularly their spatio-temporal, scale-free clustering properties arising from stress transfers and the resulting cascades of triggering [29][30][31][32] .
Here we present experimental results obtained in a 3D granular system in a post-(macro)yield regime displaying a stationary shear band and in which finite-size effects do not dominate the dynamics. Using an interferometric method of measurement of micro-deformations we provide a direct spatially-resolved measurement of the micro-slip events that govern the frictional motion along the shear band. We are able to measure the localization, the spatial extension, and the magnitude of those events, providing the first direct experimental measurement at the laboratory scale of frictional micro-slips along a fault. We show that the statistics of those events display scale-free behavior in close agreement with earthquake phenomenology. Using a methodology developed for earthquake analysis 30 , we go another step further compared to previous experimental studies by quantifying the causal triggering between events. We show that this underlying triggering process can explain the observed spacetime correlations in the dynamics, much like it does for earthquakes. We argue on this basis that a frictional shear band in a granular material represents a formal analog of tectonic faults, with intermittent dynamics probably belonging to the same universality class.
Strain fluctuations inside a shear band Stationary shear band and strain imaging. We use for this study an experimental set-up that consists of a biaxial cell filled with a granular material composed of an assembly of glass beads confined into a rectangular box (see Fig. 1a). The two lateral faces are deformable latex membranes which allow us to impose a confining stress σ 1 to the material. This stress is kept constant during the whole experiment. The material is slowly compressed by moving at a fixed velocity the top face with respect to the bottom face, and the axial macroscopic deformation ε M , and the applied axial stress σ 3 are measured (see "Methods").
The strain fluctuations are imaged using an interferometric technique based on diffusing wave spectroscopy (DWS). For this, the material is illuminated with an extended laser beam, and the speckle images are regularly recorded. We note δε Ã M the macroscopic strain increment between two successive images, and its value is δε Ã M ¼ 5:10 À7 if not otherwise specified. The speckle images are then divided into square zones, and for two successive images, the normalized autocorrelation function of the scattered intensities is calculated for each zone. Associating a color to the value of the correlation at a position, we obtain maps of correlation g I (ε M , r), where r is the position on the observation plane, as shown in Fig. 1d. High correlation g I ≈ 1 (white pixels) indicates that beads are uniformly translated without relative motions, whereas low correlation g I < < 1 (dark pixels) is the signature of bead relative motions. In addition to this interferometric correlation technique, we use a conventional digital image correlation method on the speckle pattern: the displacement of zones of the speckle pattern between different images are measured, giving access to the displacement field. This measure is used to determine the mean relative velocity of blocks when shear bands are formed.
Starting from an initial condition of a material submitted to an isotropic confining pressure, the material is slowly compressed. The beginning of the compression is associated with a plastic flow spatially distributed into the sample and to an increase of the stress difference σ 3 − σ 1 (see Fig. 1b). We analyzed this plastic flow in previous works 22 (see also Supplementary Movie 1) and we do not look further to this initial stage here. When the deformation ε M of the material exceeds a yield strain ε Y ≈ 4.5%, the stress difference σ 3 − σ 1 is roughly constant (see Fig. 1b), and the deformation localizes into the material, with the formation of one or two linear shear bands 33 . The shear is not localized close to a moving mechanical boundary as it is the case for a Couette cell (where the shear band appears at rotor), or in a gouge confined between two rigid blocks. Here, the shear band emerges spontaneously in the system. Its orientation is linked to intrinsic properties of the material (it is linked to the deviatoric stress at failure through the Mohr-Coulomb relationship) and not to geometrical constraints. Its width is the result of the selforganization of the system: the flowing material forming the band and the solid material surrounding it is the same and this separation of phase emerges spontaneously in the system.
Average vs. instantaneous strain. The map of the correlation of the scattered intensity can be linked to the shear motion of the sliding blocks at each side of a band. This may easily be seen qualitatively: for this we consider a correlation map obtained in the stationary regime (see Fig. 1d), i.e., when the stress difference is in a plateau phase, and ε M > 5%. The correlation is close to 1 into the four triangular zones partitioned by decorrelated boundaries. This indicates that the material is split into four rigid blocks separated by deformed zones.
To obtain quantitative information about the shear field inside the band, we assume (this hypothesis will be discussed just below) that the motion of the beads around a point r is mainly a shear γ m (r, ε M ) = ∂u x /∂z, where u is the local displacement of a block with respect to another one, and (x, z) are local coordinates associated to a band (see Fig. 1d for axis definition). By making this assumption, we neglect other components of the strain tensor and uncorrelated motion of the beads. If the beads move accordingly to this shear field, the decorrelation can then be related to the local shear as (see Supplementary Methods 2 with γ 0 = 2.6 × 10 −4 a constant given by the optical properties of the material. The time-averaged local deformation is defined as: The hypothesis of local shear can be quantitatively tested. For this, we integrate (2) along a direction perpendicular to the shear band, and we obtain u x ðzÞ À u Figure 1c shows the displacement field across the shear band, Fig. 1 Imaging shear band fluctuations. a Schematic of the experimental set-up. The material is submitted to a biaxial stress test. The front face of the sample is imaged on a camera. As illumination is done using coherent light, those images display speckles. b Normalized deviatoric stress as a function of macroscopic deformation. ε Y is the yield strain, and the gray zone is the post-yielding zone analyzed in this study. c Relative displacement u x of two blocks separated by a shear band as a function of the direction perpendicular to the shear band. Symbols are experimental data, the plain line is u x ðzÞ ¼ Δu x ½1 þ tanhðÀ2z=wÞ=2, with Δu x = 32 nm and w = 22d. d Map of the correlation between two successive speckle images. The color of the pixel is related to the value of the correlation. Blue arrows represent the directions of displacement of the blocks with respect to the shear bands. The dashed area is the Region of Interest (ROI) for performing speckle analysis. e Schematic of the shear band separating two sliding blocks as composed of discrete shear events of fixed-width w and of size L × L.
demonstrating that the strain is concentrated into a narrow zone of width w = 22d. This value is in agreement with the reported values of shear band thickness in biaxial experiments 34 . The spreading of light due to the diffusive transport of photons inside the granular material is expected to occur on a scale of the order of l *~3 d, which is small compared to our measured value of w. This value of w = 22d is, therefore, representative of the thickness of the actually sheared zone.
For our macroscopic strain increment δε Ã M ¼ 5:10 À7 , the relative sliding of the two blocks Δu This value is close to the imposed value Δu x = 25 nm that we may estimate from the displacement of the top plate and the orientation of the shear band. This agreement confirms the hypothesis that the decorrelation of the speckle pattern is mainly due to the shear motion of the beads. The small difference between the imposed and the measured value of displacements may have different origins. Shearing a granular media usually produces some dilation or compaction, creating a decorrelation of the scattered light. Also, some nonaffine motion of the beads occurring in the sheared zones may produce decorrelation 35 . Finally, our optical model of perfectly spherical beads ignores possible bead rotations.
Memory effect in strain fluctuation. We now analyze the fluctuations of the local shear. Since deformation is located in the shear band, we consider the transverse-averaged shear deforma- Figure 2b shows γ T (x, ε M ) into the (x, ε M ) plane. We can clearly see that the shear is heterogeneous both in space (i.e., along the shear band) and intermittent. In order to analyze those fluctuations, we introduce the normalized spatiotemporal correlation function: where 〈⋅〉 is an average both on deformation and position (see "Methods" for details).

Results
Definition of shear transformation events. As we saw in Fig. 2a, the local shear γ T (x, ε M ) exhibits important fluctuations, alternating activity, and quiescent phases. The concept of shear transformation zone has been introduced to deal with the flow of disordered materials: spatial zones reorganize, creating a local shear. We define such zones from the light scattering data. For this, we apply a threshold to the image γ m (r, ε M ), and we use a particle detection algorithm to obtain individual shear events (see Supplementary Methods 3 and Supplementary Fig. 5 for details about threshold and detection algorithm). Events are numbered, and each event i is associated with a macroscopic deformation ε M ;i at which the event occurs, a position r i (defined as the barycenter of γ m ), a surface Σ i on the image, and a mean microscopic shear γ i . We also define a "seismic moment" as i with μ the shear modulus of the granular assembly, u i the shear displacement, and L 2 i the shear surface (see Fig. 1e). For a shear zone of width w, we have Σ i = wL i and u i = wγ i , and thus The shear modulus of the granular assembly may be estimated from mean-field theory of granular elasticity (equation (14) of ref. 36 ) and is μ = 200 MPa for the pressure of 30 kPa. The shear bandwidth w = 22d is measured from the mean strain (Fig. 1c). The moment magnitude is defined as m w;i ¼ 2 3 log 10 ðM i Þ À 6:07, with M i expressed in 37 Nm. The energy dissipated during an event is where τ is the shear stress along the shear band and is then directly linked to the moment: M i /μ = E i /τ. The value of τ = 33 kPa may be obtained from the principal stresses σ 1 and σ 3 at failure using Mohr-Coulomb construction.
Scaling laws of events. We now look at the statistical characterization of the events. We consider for this the sequence of events occurring on one half shear band for macroscopic deformation 6% ≤ ε M ≤ 10% (gray zone on Fig. 1b). The total number of counted events is N tot ≈ 1.1 × 10 5 . The minimum moment is dependent on the threshold and is M min = 6 × 10 −7 Nm, whereas the largest events have M max = 0.05 Nm. The probability density function of energy dN/dM is plotted on Fig. 3a, and decays as a power of energy dN/dM~M −β , with β ≃ 2.1. Although our moments stand in a range roughly 20 orders of magnitude below that of earthquakes, this behavior is similar to the empirical Gutenberg-Richter's law. Indeed, the number N(M) of earthquakes with a moment magnitude larger than m w is log 10 ðNðm w ÞÞ ¼ a À bm w , leading to dN/dM~M −(1+(2/3)b) . The value of b for faults is usually ≃1.0 1 , leading to a slope β ≃ 1.66. It may be verified (see Supplementary Methods 4 and Supplementary Fig. 6) that the distributions remain the same during the whole post-failure regime, and that they do not depend on the occurrence of large events. In this sense, our model gouge is in a stationary state. The mean deformation of an event of moment M is defined as , for a small dM. As shown on Fig. 3b, this quantity is relatively constant 〈γ〉~M 0.17 . The broad distribution of the values of the moments M i (or equivalently of the relaxed energies E i ) is then mainly due to a broad distribution of sizes L i , but not of deformation γ i . In other words, the stress drop in an event Δτ i = μγ i has always the same order of magnitude. This is consistent with what is generally considered for earthquakes. Indeed, compilations of earthquake data argue for a scaling 38 M~L 3 , while M = μL 2 u = μγL 3 = ΔτL 3 , hence implying a constant stress drop. In our experiment, Δτ i /μ = γ i ≈ 10 −4 , this ratio is relatively close to the one typically obtained for earthquakes where 39 Δτ/μ ≈ 3 × 10 −5 . We should however mention that, in our case, the increase of 〈γ〉 i.e., of Δτ, with the seismic moment, though weak, is significant (Fig. 3b). In the case of earthquakes, a potential similar scaling would be indiscernible, owing to the large uncertainty on the estimation of the average slip and the variety of geophysical contexts. We may also consider the relative fraction of the shear stress which is relaxed during an event: Δτ i /τ = μγ i /τ. For large events, 〈γ〉 ≈ 2.10 −4 , whereas 〈γ〉 ≈ 3.10 −5 for small events, hence showing that the events relax typically~0.1-1 of the mean stress.
We finally look at the ratio of energy dissipated by events of moment greater than M: E ev ðMÞ ¼ ∑ M i ≥ M E i , compared to the total dissipated energy E dis (see Supplementary Methods 5 for the details). Figure 3c shows that typically half of the energy is dissipated in events of moment M > 10 −4 Nm, while all detected events account for a value R s = 90% of E dis . R s is an equivalent of the seismic coupling defined in tectonics 1 . The fact that the value is close to 1 means that only a minor part of the energy is dissipated into a continuous motion, distinct from the detected intermittent shear transformation events. For tectonic faults, the coupling can vary considerably with the geophysical context but is generally strong for interplate continental faults 40 .
Temporal organization of events. The statistical laws governing the succession of shear transformations may be analyzed within the framework of the statistical laws of natural earthquakes.
We first look at the rate of events occurring at the same position as a particular event (so-called mainshock). Figure 4a shows the rate of events occurring after ("aftershocks") of before ("foreshocks") "mainshocks" of magnitude m w ≥ −9 (total number of mainshocks ≈13 × 10 3 ). Only aftershocks or foreshock events occurring at the same position (±15d) are counted. For large delays, the rate of events is constant, corresponding to a background (dN/dε) bg rate of events uncorrelated to the mainshocks, while, at small delays, the rate excesses the background rate. The excess rate of aftershocks (dN/dε) exc = (dN/dε) − (dN/dε) bg decays with the macroscopic deformation as ðdN=dεÞ exc $ δε À1 M (see Fig. 4b). This behavior is reminiscent of Omori's law 5 which states that the rate of seismic events occurring after a mainshock decays with the time t to the mainshock as n(t)~t −1 . Note however that such analysis amounts to blindly characterize time correlations between events. In particular, unlike what is generally done in earthquake analysis, the magnitude of the "mainshock" is not prescribed to be larger than that of its "aftershocks". Such correlation analysis is a signature, but not a formal quantification, of causal triggering, which is explored in details later. An opposite evolution, reminiscent of an inverse Omori's law, characterizes, on average, an increasing rate of foreshocks before mainshocks, though with a smaller rate (Fig. 4a) that expresses an asymmetry of time clustering.
The productivity law describes the number of excess events in response to an event of magnitude m w . For this, we integrate the total number of aftershocks in excess to the background: exc dðδε M Þ. Figure 4c shows the evolution of N exc with the mainshock magnitude m w , and we find that N exc / 10 α:m w , with α ≃ 0.44. This result may be compared with the productivity law for natural earthquakes, where the number of aftershocks n AS / 10 α:m w , with α in the range 0.6-1.2 41 . In striking contrast, the number of foreshocks appears independent of the mainshock magnitude (Fig. 4c). This is in full agreement with a previous analysis of seismic foreshocks showing that such precursory activity before any event is a mere statistical consequence of cascades of triggering 42 . The triggering of deformation events in our system is thoroughly analyzed below. The distribution of events during the loading may be further characterized by considered the first-return deformation probability P(δε M ) which is analogous to the first return time probability for earthquakes. For this, we measure the macroscopic deformation δε M between successive events, occurring at the same position (±15d) along the band. Only events of magnitude greater than m w are considered. Figure 5a shows the distribution of interoccurrence deformations for different moment thresholds. The distributions decay with a power-law, followed by an exponential decay. As shown on Fig. 5b, those distributions may be properly collapsed by considering for every magnitude m w , the normalized deformation of x ¼ δε M =δε M , where δε M is the mean deformation between successive events. As indicated in Fig. 5b, the distribution may be well approximated by a Gamma distribution PðxÞ / x qÀ1 expðÀx=BÞ ð 4Þ with q = 0.6 and B = 1.9. P(x) decays as a power-law with exponent ≈0.4 up to values of x ≈ 1, then exponential decays take place. This behavior is not surprising, as it has been shown to be a mere consequence of triggering dynamics characterized by a GR distribution, the Omori's, and productivity laws 43 . Our results are very similar to observations for tectonic seismicity where q ≃ 0.67 44 , or micro-seismicity q ≃ 0.74 45 . Such scaling laws are also observed in fracture experiments 46 .

Discussion
Shear band viewed as a minimal model of the gouge. Starting from an initially homogeneous assembly of beads, our system organizes spontaneously to reach a stationary regime where all the deformation is concentrated along shear planes. The analysis of the statistical properties of the strain fluctuations along those planes shows a strong quantitative and qualitative analogy with the statistical characteristic of natural earthquakes along a fault: the shear band may be viewed as a simplified gouge. We discuss here why this scale-free organization of deformation along a gouge is not observed in other laboratory systems.  Many mechanical systems other than crustal faults exhibit crackling noise when plastically deformed 10,11,15,27,[46][47][48] . In particular, the statistical properties of deformation or mechanical stress fluctuations follow power laws which reveal the absence of any particular scale in the system at least in an extended inertial (scale) range. Those fluctuations arise from individual deformations which interact elastically to create complex scale-free dynamics. In the case of plasticity of amorphous materials before the macroscopic yielding 15,16,48,49 , the criticality of the system is related to the approach of the yield point. In this case, the plastic events are expected to be initially randomly distributed throughout the bulk of the sample and not localized in structures analog to natural faults. The brittle failure of an amorphous material is another configuration where crackling noise is observed 10 . In this case, the plasticity occurs in a damaged zone close to the propagating crack tip, and stationary plastic deformation cannot be defined. In some experiments such as compression of disordered materials 46 or some granular experiments 13,27 , the spatial extension of events and their localization is unknown.
In order to obtain a zone of intense plasticity in a stationary regime, one can shear an artificial gouge, made of a granular material confined between elastic plates of very different elastic modulus (much stiffer or much softer) 14,25,26 . At first glance, this appears to be a reasonable realization of a natural gouge which consists of highly crushed rocks confined between rigid material. However, in those cases, avalanches of various sizes are not observed, but instead, the dynamics are dominated by macro-slips involving all the sliding interface. De Geus et al. 28 have recently proposed a numerical model of frictional contact consisting of an amorphous layer confined between two elastic blocks, in which scale-free dynamics and large macro-slips events implying the whole interface coexists. Such a competition between an avalanche regime and a periodic stick-slip is reminiscent of the Parkfield segment of the San Andreas Fault, confined between a creeping zone and an unloaded segment, where large and pseudoperiodic earthquakes have been observed 50 . In our experiment, the dynamics of the model gouge are dominated by scale-free avalanches and we do not observe such macro-slips. This difference of behavior may arise from the difference in the confinement of the gouge. When the gouge is confined between elastic plates, there is an important contrast of mechanical properties between the gouge (which is an elasto-plastic granular material), and the plates (which are perfectly elastic). We may then expect that the plates transmit integrally the mechanical stress over all the gouge, leading to macroscopic slip events. In our experiments, materials that compose the fault and the surrounding medium are the same: both consist of the same glass beads. Given the applied pressure in our experiment, we do not expect any bead crushing, and this is in agreement with optical observations. So, the mechanical properties of the material are probably very close inside and outside the shear band, and they are both elasto-plastic. So, the material outside the gouge does not behave as a rigid block transmitting the mechanical stress on all the interfaces. This is probably why we do not observe any large macro-slips but only localized avalanches displaying scale-free dynamics.
In summary, since the stationary shear band emerges from bulk material, we are able to observe a scale-free stationary dynamics occurring in a confined space. The shear band of granular material has the right dimensionality (the 2D shear plane in 3D space) and the right mechanical properties to accurately model a complex gouge at the laboratory scale. As a consequence, we directly observe shear events distributed along the shear plane. The statistical properties of those events are summarized in Table 1, and their size distributions, temporal and spatial organizations, as well as correlations of displacement, are very similar to the ones observed for natural earthquakes. We demonstrate below that the analogy can be pursued one step further through a thorough analysis of triggering.
Triggering of deformation events. Correlations in the deformation field and among deformation discrete events are found both in time and in space and obey power-law regimes that highlight the scale invariance of the system. However, correlation is distinct from causality, which in the present context is equivalent to triggering, i.e., how the occurrence of a deformation event mechanically triggers subsequent deformation events. The underlying causal structure can be inferred from the data using methods that have been developed in seismology 30,31,51 or in social science 52 . We find that triggering obeys a scale-free productivity law, so that the number N of directly triggered events, per mainshock, depends on the magnitude as N $ 10 0:24m w (Fig. 6a), along with an Omori-like kernel, albeit with a relatively steep decay, the density of triggering events decreasing with time t after the trigger as t −p with p in the 1.6-1.8 interval, cf. Fig. 6b. Departure from these power laws is observed for the biggest events, that produce relatively more aftershocks in the early times, but are then followed by a clear activity shutdown, both features being likely due to a finite size effect and exhaustion of the stressed, ready-to-fail patches along with the deformation band after such large events. It is customary in the framework of these models to define a so-called branching ratio, which measures the capacity of a perturbation to sustain itself over potentially an infinite time (if the branching ratio is close to 1) or instead to die off quickly (if it is close to 0). This ratio can be estimated as the number of directly triggered aftershocks per "mainshock" averaged over all the events of the catalog. We here find that the branching ratio is very close to 1, implying that the background (i.e., non-triggered) activity consists of a few % at most so that most of the activity is made of events triggered by preceding events, highlighting the dominance of triggering and therefore clustering in the dynamics of deformation events. This is also fully consistent with earthquakes, for which the branching ratio ranges between 0.8 and 1 53 .
The productivity in 10 0:24m w found here is distinct from the 10 0:40m w scaling observed when stacking the activity past all m w ≥ 9 events as in Fig. 4c; this is due to the fact that, in the latter case, the stacking mixes causally-triggered sequences (e.g., if A triggers B that triggers C, then both B and C will show up in the counting of subsequent activity, while in Fig. 6 B is counted for A while C is counted for B). We thus checked that this mixing does indeed recreate the observations of Fig. 4. To do so, we exploit the fact that the causal structure can be formulated as a linear model, that is simply amenable to simulations 30,31,51 . We thus simulate synthetic datasets of deformation events based on this model and its basic ingredients: (1) seismic moments are independent, identically distributed, and follow the Gutenberg-Richter-like marginal distribution of the experiment; (2) a small proportion (about 5%) of the events occur randomly in space and time, and correspond to "spontaneous" events, i.e., events that are not triggered by previous events; (3) the 95% other events are triggered events from previous "mainshocks", their distribution relative to the time and position of the mainshock following the kernels observed for the experiment dataset (e.g., the temporal kernel of Fig. 6b). Generating such synthetic datasets, we find that the correlations (i.e., stacked rates) seen for the real data are indeed well recovered (they are within the natural fluctuations of simulation outcomes), demonstrating that these correlations effectively emerge from more fundamental triggering kernels, cf. Fig. 6c.
Structural vs. temporal memory effects. Strain correlation functions (Fig. 2), rate excess after main-shocks (Fig. 4a), and the causal triggering kernel (Fig. 6) all indicate the existence of scalefree memory effects in our system. In the context of seismology, memory effects are remarkably revealed by the existence of aftershock sequences which are quantified by Omori's law stating that the number of aftershocks decays as the inverse of the time elapsed since the mainshock 1 . This law assumes implicitly that the time is the physical variable that governs the memory. The origin of such temporal dependence is however unclear. Several mechanisms such as the temporal dependence of microscopic friction law [54][55][56][57] , sub-critical crack growth, the occurrence of afterslip 58 , or poro-elasticity and the evolution with time of pore fluid pressure 59 have been proposed as a possible sources of temporal memory effects controlling earthquake occurrences. However, the direct links between any time-dependent microscopic mechanism and memory effects in seismicity are still debated.
In our experiment, we can test whether time is indeed the right parameter to describe the memory effect that we observe. For this, we performed experiments at different macroscopic deformation rates. Figure 7a shows the normalized correlation functions of the microscopic deformation expressed as a function of the time increment. If every experiment shows a memory effect, the magnitude of memory depends on the strain velocity. At a fixed time delay δt, the correlation function decreases with the velocity. This reveals that time does not seem to be the right parameter to describe memory. This may be evidenced by plotting the correlation functions as a function of the macroscopic strain increment δε M Fig. 7b. In this plot, the curves collapse, demonstrating that the correlation function decays with the strain increment and not with the time increment. This independence on the shear rate may also be evidenced by considering the similarity of the size distributions of events on Fig. 7c.
This laboratory observation is evidently not in contradiction with Omori's law. Indeed, the driving velocity of a given fault is constant on the temporal scale of human observations. So, Fig. 6 Triggering kernel. a Productivity law giving the mean number N of triggered events for a trigger of magnitude m w . The best power-law fit in N $ 10 0:24m w obtained when discarding the last point (biggest events) is shown in magenta. b Triggering kernels in time, conditioned on the moment of the trigger. We consider the same 8 "classes" of seismic moments as in (a). Power law decays in t −p , with p = 1 and p = 2, are shown for visual guidance. c Correlation in time, as in Fig. 4b, for two instances of synthetic datasets and for the real dataset. describing memory effects in terms of time increment or in terms of strain increment is then equivalent for natural faults. The fact that the memory is strain-dependent rather than time-dependent in our system suggests that the memory could be linked to structural/topological rearrangements within the granular medium, inducing a redistribution of local stress and possibly triggering slip events. It is then not surprising that the macroscopic deformation may be the parameter that governs the plasticity around a given position in the material. This also raises the question of the potential role of such geometrical rearrangements in the "time" correlations characterizing natural seismicity. In that case, such mechanisms could combine with truly time-dependent, thermally activated processes such as subcritical crack growth, to explain memory effects in earthquake occurrences. Interestingly, slip-dependent and time-dependent memory effects combine as well in the classical Rate-and-State friction laws 60,61 that remains nowadays a classical framework of interpretation of earthquake physics 55,62 .

Conclusion
By looking at the intermittent strain fluctuations, we showed that a shear band inside an athermal disordered material is an analog of a natural fault: the deformation consists of many micro-slip occurring along a plane, and their collective dynamics is characterized by statistical properties remarkably consistent with the empirical laws of seismology. This analogy with natural faults is obtained when the fluctuations are observed after macroscopic yielding of the granular medium, when a steady-state regime takes place. The laboratory and natural fluctuations observed are then characteristic of a critical behavior after yielding, which is presumably different from the stress-tuned criticality observed for many systems before yielding. Our statistical analysis of microslips also quantifies the causal triggering between events, and reveals that this underlying triggering mechanism is at the root of the space-time correlations in the dynamics, as it has been previously shown for natural earthquakes.
Despite its simplicity, our experimental model recreates the whole set of spatio-temporal characteristics of earthquakes: moment distribution, recurrence time distribution, productivity law, decay of aftershocks rate, and triggering kernel. This simply suggests that those laws may be reproduced by using simple numerical or theoretical systems of frictional particles. Hence, the analysis of such systems should allow understanding the organization of the microscopic stress field close to the shear band. Experimentally, the possibility to control a model fault in the laboratory also opens the road to many studies, such as the effect of mechanical noise on the size distribution of slip events, the influence of the elastic properties of the surrounding material with respect to that of the band, or a study of size effects.

Methods
The experimental setup consists of a biaxial compressive test in-plane strain conditions already described 63 . The granular material is composed of glass beads of diameter d = 70-110 μm and initial volume fraction ≈0.60. The granular material is enclosed between a latex membrane (85 × 55 × 25 mm 3 ) and a glass plate. A pump produces a partial vacuum inside the membrane, creating confining stress −σ 3 = 30 kPa. The sample is placed in the biaxial apparatus where a compression is imposed in the vertical direction (noted 1 on Fig. 8a) while the lateral sides (noted 3 on Fig. 8a) are maintained under a constant stress σ 3 . Two fixed planes ensure that there is no deformation along the direction 2. At the top, a moving plate exerts a compression of the sample at fixed velocity and the bottom plate is fixed. The velocity of the motor is of the order v~μm s −1 (different velocities are used, see main text for precise values), leading to a deformation rate _ ϵ M ¼ v=L $ 10 À5 s −1 where L = 85 mm. The origin of deformation ϵ M = 0 is taken at the beginning of the contact between top plate and sample, but its precise value is without important because analysis considers only deformation increments. The stress applied at the top of the sample is defined as −σ 1 = −σ 3 + F/S, where F is the force measured by a sensor fixed to the upper plate, and S = 55 × 25 mm 2 is the section of the sample that we consider as constant during all the experiment. When the stress plateau is attained, one or two shear bands form into the sample (green and red planes on Fig. 8a).
Deformations are observed through the front glass plate using DWS, a method already described elsewhere 35,64 . A laser beam at 532 nm is expanded to illuminate the entire sample. The light undergoes multiple scattering inside the granular material and we collect the backscattered rays. The latter interfere and form a speckle pattern. The image of the front side of the sample is recorded by a 2352 × 1728 pixels camera. Speckle images are then subdivided into square zones of size 16 × 16 pixels that we call metapixel. The correlation between two images 1 and 2 is then calculated for each metapixel, and a map of correlation of 147 × 108 metapixels is obtained. The size of the metapixel corresponds to 6.0d × 6.0d on the sample. Average