Interevent time distributions of avalanche dynamics

TL;DR

This paper investigates the interevent time distributions of avalanches in stick-slip systems, revealing their ability to distinguish different mechanical states and indicating large correlations, unlike size and duration distributions.

Contribution

It demonstrates that interevent time distributions can differentiate system states and reveal correlations, providing new insights beyond traditional avalanche size and duration analyses.

Findings

Interevent time distributions vary with system state.

Large packing ratios show earthquake-like interevent times.

Interevent times reveal correlations not seen in size/duration.

Abstract

Physical systems characterized by stick-slip dynamics often display avalanches. Regardless of the diversity of their microscopic structure, these systems are governed by a power-law distribution of avalanche size and duration. Here we focus on the interevent times between avalanches and show that, unlike their distributions of size and duration, the interevent time distributions are able to distinguish different mechanical states of the system, characterized by different volume fractions or confining pressures. We use experiments on granular systems and numerical simulations of emulsions to show that systems having the same probability distribution for avalanche size and duration can have different interevent time distributions. Remarkably, for large packing ratios, these interevent time distributions look similar to those for earthquakes and are indirect evidence of large space-time correlations in the system. Our results therefore indicate that interevent time statistics are more informative to characterize the dynamics of avalanches.

Equations

$$G(t_i)=\frac{C}{t_m}(\frac{t_i}{t_m}{)}^{(\gamma -1)}\exp (-\frac{t_i}{t_m})$$

$$R\sim \frac{2\delta \sqrt{N}}{L}$$

$$A_{sup}(t)=su{p}_i[A_i(t)]$$

$$E_r=[-\sigma \frac{d\sigma }{dt}|A_{sup}]$$

$$E_r\sim A_{sup}^{1/2}$$

$${\psi }_r\equiv \frac{1}{r^2}{\int }_{B(r)}dxdy{\left[\frac{A_i}{⟨A⟩}\right]}^2$$

$$S\sim t_E^2\sim {\int }_{t_E}dt[-\sigma \frac{d\sigma }{dt}]\sim t_E{E}_r\sim t_E{A}_{sup}^{1/2}\sim A_{sup}$$

Full text

Introduction

Many physical systems subject to small external driving forces exhibit complex burst dynamics in space and in time Fisher ; Sethna ; Zapperi ; uhl_etal:2015 ; maass_etal:2015 . Burst events are the signature of energy release in the system and, in many cases, they have successfully been described in terms of avalanche dynamics. Both theoretically and experimentally, the statistical distribution of the size, $S$, and the time duration, $t_{E}$, of avalanches have been shown to satisfy well-defined scaling laws. In the case of plasticity of soft glasses and fracture dynamics of amorphous solids, the scaling exponents are mostly independent of the details of the microscopic interactions, suggesting some form of universality BenZion ; Goldenfeld ; Wyart ; uhl_etal:2015 ; Denisov17 . Theoretically, the scaling exponents have been explained by a number of different mean- and non mean-field theories, which take into account the basic physical properties of the system and its intrinsic “randomness” Fisher ; BenZion ; Wyart . This apparent universality might also explain the observed scaling properties seen in earthquake dynamics, e.g. the well-known Gutenberg-Richter law gb:1954 ; BenZion ; uhl_etal:2015 . The question naturally arises whether some statistical properties of avalanches are able to discriminate between different states of these systems. One interesting quantity is the interevent time, $t_{i}$, (also referred to as recurrence time, return time or waiting time) between two consecutive avalanches maass_etal:2015 . Most of our information on the interevent time distribution $P(t_{i})$ comes from the analysis of earthquake catalogs. There is a general consensus that for earthquakes, $P(t_{i})$ can be characterized by two different time scales: at short times $P(t_{i})\sim 1/t_{i}$ follows Omori’s law utsu_etal:1995 , while at very long time scales $P(t_{i})\sim exp(-\alpha t_{i})$. The exponential behaviour is explained by the (reasonable) assumption, that for long time scales, earthquake events are likely to be independent. Corral Corral2004 , however, showed that $P(t_{i})$, for long $t_{i}$, is better fitted by a Gamma distribution, leading to a time scale defined by the average interevent time. Interestingly, this Gamma distribution is observed for many different geographic regions, as well as for the whole earth, once $P(t_{i})$ is re-scaled to the regional or global average interevent time Corral2004 . A detailed investigation of seismicity induced by mining and fluid injection revealed the same Gamma distribution for interevent times Davidsen ; Davidsen2 as well as acoustic emissions of various other systems Niccolini ; Ribeiro . If the interevent time distribution is independent of any size threshold or region, the only possible shape for $P(t_{i})$ is exponential, unless complex space and time correlations are present and the system is close to criticality Molchan ; Corral2006 . Another option is that the Gamma distribution emerges as a superposition of independent probability distributions, one with an exponential tail and one with the Omori short-time behaviour Saichev2006 ; Touati . In the case of superposition, $P(t_{i})$ does not reveal any new physics besides the well documented avalanche scaling laws (Gutenberg-Richter and Omori). Finite detection thresholds Santucci2016 have also been suggested to explain the emergence of Gamma distributions. It is thus a fundamental question to know whether $P(t_{i})$ informs us on long-range correlations in the system or not. In the past, studies have tried to settle this question by analyzing the fit of data to various functional forms for $P(t_{i})$. Without the knowledge of the underlying state of the system, best fit arguments are obviously difficult to make. We will instead analyze systems where we control the mechanical state (e.g. rigidity) and investigate whether $P(t_{i})$ informs us on physics beyond that of $P(S)$ and $P(t_{E})$.

We report on a systematic study of the distribution of interevent times, $P(t_{i})$, across different model systems to show that it provides crucial information on the mechanical state of the material, namely the packing ratio or its confining pressure. We complement experimental measurements on granular systems with numerical simulations on emulsion-like models, and compare the resulting interevent time distributions with those for earthquakes. Our results uniquely show that unlike the statistical properties of avalanche sizes and durations, the probability of interevent times $P(t_{i})$ strongly depends on the packing ratio or confining pressure of the material. For relatively low packing ratios, the interevent time distribution is exponential with an Omori behaviour at small time scales. For high packing ratios, $P(t_{i})$ is similar to that observed for earthquakes at long time scales Corral2004 , with an Omori scaling at small time scales. This implies that the functional shape of $P(t_{i})$ depends on the mechanical properties of the system and that, in addition to avalanche size and duration, it provides a crucial measure to distinguish avalanche-like relaxation mechanisms. We show that spatio-temporal correlations are responsible for these different regimes, and therefore interevent time distributions provide insight into the nature of correlated deformations in dense suspensions and earthquakes.

Results

We use granular systems under well-controlled normal forces, apply slow shear strain rates to induce avalanches and monitor them with high temporal resolution (see Methods section). The recorded force signal from our shear cell exhibits strong intermittency: force increases are followed by sudden force drops that demarcate energy release events (Fig. 1 top). We measure the force drops with high temporal resolution to resolve the dynamics of both the large and small avalanche events. Previous experiments have shown that the applied shear strain rate is sufficiently slow to separate individual avalanches and avoid avalanche overlap Denisov16 . This enables us to extract a wide range of predicted scaling exponents and scaling functions that identify the underlying slip statistics and dynamics Denisov16 ; Denisov17 . The second granular system, a rotary shear cell ( see Methods section), similarly shows strong intermittency: stress slowly increases, followed by sudden stress drops that correspond to acoustic emission (AE) events (Fig.1 middle). We also analyze numerical simulations of a two-dimensional emulsion using Lattice Boltzmann Equation (LBE) modelling LBE1 (see Methods section). The model aims at simulating two repelling fluids. Coarsening is strongly suppressed by using a frustration mechanism, which stabilizes the interface. The system exhibits a yield stress rheology with a non-Newtonian relation between stress and shear strain rate above the yield stress LBE2 . For a small imposed shear strain rate, such that the stress is below yield stress, a clear stick-slip behaviour is observed (Fig. 1 bottom).

Avalanche size distributions

The scaling properties of the avalanche sizes are similar for all our systems (figure 2). Results from two different packing fractions are shown: The distributions clearly collapse on top of each other, and all systems show a scaling with an exponent of $-1.33$ over 2-3 orders of magnitude in avalanche size, irrespective of their packing fraction.

Avalanche interevent time distributions

The scaling properties of interevent times, however, behave markedly differently depending on packing fractions. At high packing fraction, all data collapse onto a Gamma distribution given by

[TABLE]

with $\gamma=0.7$, $t_{m}=\langle t_{i}\rangle$ the average interevent time and $C$ a normalization constant. This is the same Gamma distribution as that observed for earthquakes Corral2004 . The quality of the fitting is shown in the inset of Fig. 3 where we plot the ratio $R(t_{i})\equiv P(t_{i})/G(t_{i})$ where $P(t_{i})$ are the different probability distributions shown in the figure and $G(t_{i})$ is given by Eq. (1). Furthermore, for short interevent times, most data display a clear scaling behavior $P(t_{i})\sim 1/t_{i}$ consistent with Omori’s law. The short time regime can also be collapsed onto the same Gamma distributon after proper rescaling Corral2004 , but this is not essential for our main findings.

While earthquake data do not allow us to explore a wide range of packing fractions, the laboratory experiments and numerical simulations can easily explore lower packing fractions. For data corresponding to low packing fractions, $P(t_{i})$ behaves markedly differently (Fig. 4). Two interesting features are seen in the probability distributions. For long $t_{i}/t_{m}$, the systems show an exponential behavior, corresponding to $\gamma=1$ in eq. (1), clearly distinct from the distributions displayed in Fig. 3. The short time behavior $P(t_{i})\sim 1/t_{i}$ is observed in all cases, except the rotary experiment r126, but no data collapse is obtained, i.e. the time scale separation between Omori’s law and the exponential behavior is expressed differently depending on the physical system.

Figs. (3 and 4) demonstrate the main point of our paper: $P(t_{i})$ changes upon varying the system properties, i.e.the packing fraction, although no relevant change is observed in the avalanche size distribution (Fig. 2). As the estimation of $t_{i}$ is done with the same methodology in all cases, the change of $P(t_{i})$ with packing fraction cannot be considered an artifact due to the chosen thresholds Santucci2016 . Rather, since the interevent times inform us on the relaxation time of the system, and their distribution on the relaxation process, a change of $P(t_{i})$ indicates a fundamental change in the nature of the underlying relaxation processes. In this context, it is noteworthy that at high packing fraction, all interevent time distributions are the same as those observed for earthquakes Corral2004 ; Corral2006 ; Saichev2006 ; Davidsen2 , suggesting a similar relaxation mechanism for earthquakes.

Discussion

From the above, we can draw some general and nontrivial conclusions.

The statistical properties of avalanche size, $S$, and duration, $t_{E}$ (Figs. 2 and 11) Denisov17 , are independent of the system stiffness (i.e. the packing fraction). Physically this means that, once the system starts to release elastic energy, it does this independently of the system’s packing fraction. This is one of the basic assumptions in many theoretical frameworks so far proposed to predict the (universal) scaling properties of avalanche size distribution.

The interevent time distributions, however, show a clear dependency on packing fractions (Figs. 3 and 4). They also display a clear signature of two different time scales shaping their probability distribution. At short time scales, we observe events clustered in time where the interevent time distribution is $P(t_{i})\sim 1/t_{i}$, regardless of packing fraction. This is consistent with Omori’s law observed for earthquakes, and is most likely due to smaller size events triggered by some master event (after-shocks). Note, however, that our analysis is independent of any definition of main- and after-shock. At longer time scales, the interevent time distribution is given by Eq. (1), where the exponent $\gamma$ depends clearly on packing fraction. At low packing fractions $\gamma$ is close to 1, and at high packing fractions it is about 0.7, close to the value found for earthquakes. It is also important to remark that different experiments don’t superimpose upon renormalization at low packing fractions, while the data collapse on the same functional form at high packing fractions. While we are not aware of any theoretical framework able to describe these features and/or explain how $P(t_{i})$ changes with material stiffness, it has been argued Corral2006 that for Eq. (1) to emerge with $\gamma\neq 1$, correlations have to be present in the system. The question thus is whether our systems display any correlations.

Fig. (2) demonstrates that our systems display strong correlations in space Denisov16 ; Barrat2018 . To investigate whether or not any correlations exist in time (”memory effect” in the avalanche dynamics), we follow Corral Corral2005 and consider the probability density of the interevent times conditioned on the avalanche size. We focus of the numerical simulation at high packing fraction ($R=0.03$) and define $P(t_{i}/t_{m},S_{th})$ as the interevent time probability density for a subset of events with size $S>S_{th}$. Upon increasing $S_{th}$, if $P(t_{i}/t_{m},S_{th})$ is not exponentially distributed and it does not change its shape, then time correlations (or memory effects) exist between the interevent times for different threshold values $S_{th}$. For a process with no correlations, one can show that the only scale invariant distribution is the exponential distribution Molchan . On the other hand, correlations introduce new invariant functions in the process. The robust shape of the distributions for different $S_{th}$ in Fig. (5) then clearly demonstrates that the interevent time is indeed correlated in a non-trivial way with the size of the previous event.

Motivated by this observation, we further looked for spatio-temporal correlations by analyzing sub-regions of the system. We consider two disjoint regions, say $A$ and $B$, and their set union, $A\cup B$ including both regions, over the full simulation time interval. If the regions are independent, the only possible interevent time distribution is exponential Molchan . If there are space-time correlations between the different regions, then we expect that the probability distributions $P(t_{i}/t_{m}(A))$, $P(t_{i}/t_{m}(B))$ and $P(t_{i}/t_{m}(A\cup B))$ are all the same and none of them is exponential Corral2005 ; Corral2006 . If at least in one of the regions $A$, $B$ and $A\cup B$ the interevent time is exponentially distributed, then we can argue that the Gamma distribution observed for the whole system is just by accident. This constitutes a severe test to uncover (although indirectly) space-time correlations in the system.

In the original picture of avalanche dynamics in amorphous systems, the usual assumption is that while the avalanches themselves are highly correlated events, they occur at random uncorrelated times i.e. with an exponential distribution for interevent times $P(t_{i}/t_{m})$. For such uncorrelated random events, we expect therefore that the interevent time distributions $t_{i}/t_{m}(A)$ and $t_{i}/t_{m}(B)$ are both also exponentially distributed. This is what happens for our LBE simulation at relatively low packing fraction ($R=0.09$) as shown in the upper panel of figure (6). All curves follow the exponential distribution. In the lower panel of the same figure we show a snapshot of the time series illustrating the avalanche events in the two regions (referred to as $box\,1$ and $box\,2$). In the upper panel we also show the interevent time distribution for the whole system, which, as we know, is also exponential (Fig. 4). Note that the curves do not collapse, as they apparently sense the short time scale differently.

The same analysis for the case of a large packing fraction $R=0.03$) gives a different picture (Fig. 7). In all cases, we indeed observe the same probability distribution for $P(t_{i}/t_{m})$, which is not exponential. As we have argued above, this can only be true if there are non-trivial correlations and/or memory effects between the two different regions. Note that in this case the curves collapse also for the short Omori time scales.

While it is hard to identify the exact nature of such correlations, a clear picture merges from our observations. Overall our results demonstrate that the interevent time distribution $P(t_{i})$ is able to disentangle the statistical properties of systems at different packing fractions, whereas this is not possible by looking only at the avalanche size distribution and duration. As the interevent times are directly related to the relaxation time of the system, their distribution should contain information about the relaxation mechanism. We thus conclude that the different interevent time distributions we observe, and their interpretation in terms of absence and presence of memory effects, indicate fundamentally different relaxation mechanisms of the systems at these different volume fractions or confining pressures.

Methods

Laboratory Experiments

We use two different granular systems, where we apply well-controlled normal forces and slow shear strain rates (Fig. 8).

The first system, a shear cell, consists of $3\cdot 10^{5}$ polymethyl methacrylate spheres with a diameter of $d=1.5$ mm and a polydispersity of $\sim 5\%$, and is deformed at a constant shear strain rate of $\dot{\gamma}=9.1\cdot 10^{-4}$ $s^{-1}$ and constant normal stress, using a shear cell with a confining top plate (Fig. 8a) Denisov17 . By placing weights onto the confining plate, we vary the normal stress between 4, 6.8 and 9.6 kPa, experiments P1, P2 and P3, respectively, resulting in a packing fraction $\phi$ between $55\%$ and $60\%$. We monitor shear-induced force fluctuations using force sensors included in the shearing walls.

The second granular system, a rotary shear cell, consists of glass bead layers that are sheared using a servo-controlled rotary shear apparatus (Fig. 8b). The average particle size is 0.5 mm and the standard deviation 0.1 mm. For each experiment, an approximately 4.5 mm thick layer of glass beads is deposited in an annular-shaped shear cell that consists of four independent steel rings: two of them act as pistons that provide vertical confinement, whereas the other two provide lateral support. The cell is then placed in the rotary shear apparatus, built inside an Instron 8862 testing machine. The servo-controlled Instron actuator is used to prescribe a constant normal stress condition (8 MPa for experiment r122 and 2 MPa for experiment r126), and a Parker MH-205 motor to maintain a constant angular velocity of $0.02~{}^{\circ}/s$ by rotating the bottom piston. An axially mounted load cell ($\pm 100$ kN range, 0.008 kN resolution) measures the normal stress, and a pair of laterally mounted load-cells (each 20 kN maximum load, 0.008 kN resolution) the traction. The experiments are started at a random close packing although there is an increase over time due to the breaking of particles. For the purpose of this paper, it is important to note that a confining pressure of 2 MPa is a proxy for a low packing fraction and 8 MPa for a high packing fraction experiment. AE activity is monitored via two arrays of 8 piezoelectric transducers each, mounted inside the two steel piston rings at $45^{\circ}$ intervals.

Numerical Simulations

The simulations use a model based on Lattice Boltzmann equations (LBE) for complex fluidsLBE1 , which are discussed in detail in LBE2 ; LBE3 ; LBE4 . The model simulates two repelling fluids, say $A$ and $B$, with the same density. Coarsening is strongly suppressed by using a frustration mechanism, which stabilizes the interface. The initial configuration is chosen such that $N$ droplets of the fluid $A$ are randomly created in space with small polydispersity. The interface is filled by fluid $B$. The system exhibits a yield stress rheology with a non-Newtonian relation between stress and shear strain rate above the yield stress LBE2 . The ratio $R$ between the interface area $A_{int}$ and the bubble area $A_{b}$ can be estimated as

[TABLE]

where $\delta$ is the interface thickness, $N$ is the overall number of bubbles and $L$ is the size of the system in LBE units. Note that, qualitatively, the quantity $1-R$ can be considered as the packing fraction of the system. Since the interfaces are not sharp in our system, the correct packing fraction should be $1-CR$, with $C$ a constant of order $1$. Since $\delta$ must be finite for the interface to be stable, the only way to decrease $R$ (i.e. to increase the packing fraction) is to decrease the ratio $\sqrt{N}/L$. By decreasing $R$, we also increase the value of the yield stress, i.e. we increase the rigidity of our system. Examples of initial configurations with $R=0.09$ and $R=0.03$ are shown in figure (9a) and (9b). The initial conditions $R=0.03$ together with $L=4096$ give a similar number of bubbles to that obtained for $R=0.09$ and $L=1024$. We will use results from both cases. The rheological properties of such systems are discussed in detail in LBE2 . We perform simulations with a small externally imposed shear strain rate, whose value is chosen such that the stress is below the yield stress transition. An example of the resulting stress $\sigma$ as a function of time is shown in figure (10) upper panel: a clear stick-slip behavior is observed. To perform a detailed statistical analysis of the dynamics, we used the method recently developed in GJI . We consider $n^{2}$ small squares of size $L/n$ and, for each of the squares, we compute the quantity $A_{i}\equiv\langle(\rho_{A}(x,y,t+\tau)-\rho_{A}(x,y,t))^{2}\rangle$, where $\langle...\rangle$ is the space average over the square $i$ and $i=1,2,\dots,n^{2}$. We chose $n=32$ and $\tau=1000$ LBE time steps. We checked that different choices do not change the results discussed in the rest of this section. The quantity $A_{i}$ is a measure of the relative number of points in square $i$ which move in time interval $\tau$, i.e. $A_{i}$ is the square of the displacement occurring in the square $i$. Plastic events are localized in space and correspond to the largest value of $A_{i}$ observed in the system at time $t$. Therefore the relevant quantity to consider is GJI :

[TABLE]

A large value of $A_{sup}(t)$ corresponds to a large stress drop in the system. In figure (10) we show in the upper panel $\sigma(t)$ and in the middle panel $A_{sup}(t)$ for the same time windows. To make a more quantitative analysis, we computed the quantity

[TABLE]

where $E_{r}$ represents the time averaged value of the energy release $-\sigma d\sigma/dt$ conditioned on a particular value of $A_{sup}$. In the lower left panel of figure (10) we show $E_{r}$ as a function of $A_{sup}$ for the two simulations at $R=0.09$ (red circles) and $R=0.03$ (blue triangles). A clear scaling law with a slope of $1/2$ is observed. Using this result, we can state that the following relation holds scaling wise

[TABLE]

Next we investigated the spatial correlation in the system and computed the variables

[TABLE]

where $B(r)$ is a box of side $r$ and $\langle...\rangle$ denotes the spacial average. Using $\psi_{r}$, we can construct the multi-fractal quantities $C_{q}(r)=E[{\langle\psi_{r}^{q}\rangle}]$, where $E[...]$ is the time average. For a multi-fractal system we should observe $C_{q}(r)\sim r^{(q-1)(D_{q}-d)}$, where the exponents $D_{q}$ are the generalized fractal dimensions and $d$ is the space dimension ($d=2$ in our case). Using the above definition of $C_{q}$, the space correlation of $A_{i}$ is associated with $D_{2}$ known as the correlation dimension, namely $E[\langle A_{i}A_{j}\rangle]\sim r^{D_{2}-2}$ where $A_{i}$ ans $A_{j}$ are separated by the distance $r$. The red circles in the lower right panel of figure (10) corresponds to the quantity $C_{2}(r)$ computed from $A_{i}$. A clear scaling is observed with exponent $\sim-0.7$ corresponding to $D_{2}\sim 1.3$. Therefore we can conclude that the system displays strong correlations in space. The same feature is observed using the overlap-overlap correlation function as discussed in LBE3 .

We then analyzed the probability distribution of $A_{sup}$. In Fig. (2), we show $P(A_{sup})$ as a function of $A_{sup}$ for $R=0.09$ and $R=0.03$. As already discussed in GJI , a clear scaling is observed for $A_{sup}>A_{th}$ over 2 decades, where $A_{th}$ is some threshold value (the same for both cases). The scaling exponent $P(A_{sup})\sim A_{sup}^{-\gamma}$ is given by $\gamma\sim 1.33$ . The quality of the scaling is best appreciated in the inset where we compensate for the power law, i.e we plot $P(A_{sup})/A_{sup}^{-1.33}$. The clear definition of $A_{th}$ appearing in the figure can be used to compute the scaling properties of the avalanche dynamics. We defined the time duration $t_{E}$ of the avalanche as the time interval when $A_{sup}>A_{th}$. Within each avalanche of duration $t_{E}$, we computed the size $S$ of the avalanche as the number of the $n^{2}$ region where $A_{i}>A_{th}$. The numerical simulations were performed for $40\times 10^{6}$ LBE steps in the case of $R=0.09$ and for $280\times 10^{6}$ LBE times steps for the case of $R=0.03$. A clear dynamical scaling $t_{E}\sim S^{z}$ is observed with $z\sim 1/2$, close to mean field predictions Fisher . From this dynamical scaling we therefore obtain

[TABLE]

It follows that the scaling exponent $\gamma$ previously defined is also the scaling exponent of the probability distribution of $S$, i,e, $P(S)\sim S^{-1.33}$ similar to what has been found in many numerical simulations, see the numerical values reported in Wyart .

Acknowledgements

This research was partly funded by the Shell-NWO/FOM programme ‘Computational sciences for energy research’ under project number 14CSER022, partially funded by ERC Starting Grant 335915 (SEISMIC) and partially funded by NWO VIDI grant 854.12.011. Numerical simulations for this work were carried out on the Dutch national e-infrastructure with the support of SURF Cooperative.

Author contributions statement

R.B., A.N., P.S., F. T. and J.T. conceived the study, P.K., E.K. and D.D. conducted the experiments, R.B. analysed the results. All authors reviewed the manuscript.