#  Quasistatic kinetic avalanches and self-organized criticality in
deviatorically loaded granular media

Jordi Baró<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Calgary. 2500 University Drive NW Calgary, Alberta T2N 1N4, Canada  Centre for
Mathematical Research, Campus de Bellaterra, Edifici C 08193 Bellaterra,
Barcelona, Spain  Mehdi Pouragha  Civil Engineering Department, University of
Calgary. 2500 University Drive NW Calgary, Alberta T2N 1N4, Canada  Department
of Civil and Environmental Engineering, Carleton University. 1125 Colonel By
Drive, Ottawa, ON K1S 5B6 Canada  Richard Wan  Civil Engineering Department,
University of Calgary. 2500 University Drive NW Calgary, Alberta T2N 1N4,
Canada  Jörn Davidsen  Department of Physics and Astronomy, University of
Calgary. 2500 University Drive NW Calgary, Alberta T2N 1N4, Canada  Hotchkiss
Brain Institute, University of Calgary. 3330 Hospital Drive NW, Calgary,
Alberta T2N 4N1, Canada

###### Abstract

The behavior of granular media under quasi-static loading has recently been
shown to attain a stable evolution state corresponding to a manifold in the
space of micromechanical variables. This state is characterized by sudden
transitions between metastable jammed states, involving the partial
micromechanical rearrangement of the granular medium. Using numerical
simulations of two-dimensional granular media under quasistatic biaxial
compression, we show that the dynamics in the stable evolution state is
characterized by scale-free avalanches well before the macromechanical
stationary flow regime traditionally linked to a self-organized critical
state. This, together with the non-uniqueness and the non-monotony of
macroscopic deformation curves, suggests that the statistical avalanche
properties and the susceptibilities of the system cannot be reduced to a
function of the macromechanical state. The associated scaling exponents are
non-universal and depend on the interactions between particles. For stiffer
particles (or samples at low confining pressure) we find distributions of
avalanche properties compatible with the predictions of mean-field theory. The
scaling exponents decrease below the mean-field values for softer interactions
between particles. These lower exponents are consistent with observations for
amorphous solids at their critical point. We specifically discuss the
relationship between microscopic and macroscopic variables, including the
relation between the external stress drop and the internal potential energy
released during kinetic avalanches.

##  I Introduction

Mechanical systems with evolving internal structure often exhibit
intermittency or avalanche dynamics under quasistatic mechanical driving (see
Sethna  _et al._ (2001); Salje and Dahmen (2014); Nicolas  _et al._ (2018) for
reviews). Specifically, intermittent athermal dynamics dominated by slip
dislocations are observed in the plastic deformation of single crystals Miguel
_et al._ (2001); Zaiser (2006); Friedman  _et al._ (2012) as well as amorphous
materials such as bulk metallic Wang  _et al._ (2009); Antonaglia  _et al._
(2014) and colloidal Schall  _et al._ (2007) glasses, high entropy alloys Hu
_et al._ (2018); Rizzardi  _et al._ (2018), and granular assemblies Dahmen
_et al._ (2011); Denisov  _et al._ (2016); Berthier  _et al._ (2019a), with
their prevalent dynamics being successfully captured through discrete element
(DEM) simulations in idealized athermal quasistatic (AQS) conditions
considering Lennard-Jones (LJ) particles — simulating molecular dynamics in
glassy materials Maloney and Lemaître (2004); Salerno  _et al._ (2012);
Salerno and Robbins (2013); Parisi  _et al._ (2017); Zhang  _et al._ (2017);
Shang  _et al._ (2020) — as well as soft spheres Lerner and Procaccia (2009);
Karmakar  _et al._ (2010) — simulating assemblies of larger spherical
aggregates — and elastic interactions — simulating granular matter Gimbert
_et al._ (2013); Pouragha and Wan (2016); Karimi and Barrat (2018). In
granular materials, where potentials are exclusively repulsive, intermittency
occurs upon the activation of fragile jammed states Cates  _et al._ (1998);
Liu and Nagel (1998), where the system is consolidated in marginally stable
solid-like configurations that become unstable under infinitesimal stress
variations. Jamming states can only exist within certain conditions Liu and
Nagel (1998). Under external mechanical driving, and in absence of thermal
fluctuations, jamming is observed to be limited below a macroscopic shear
stress threshold, $\sigma_{c}$, often referred to as a critical point, which
becomes stationary after an initial transient regime Bi  _et al._ (2011). When
the system is pushed above $\sigma_{c}$, jammed states are unlikely and the
excess of free energy (corresponding to the excess of shear stress) is
dissipated through particle rearrangements. Under strain driving, such
dissipation reestablishes the mechanical stability at stress values below
$\sigma_{c}$. Such a sequence of quasistatic metastable states leads to an
elastic/plastic stress-deformation response and the macromechanical state of
the system eventually hovers around $\sigma_{c}$, which can be considered a
case of self-organized criticality (SOC) Bak and Chen (1991); Chen  _et al._
(1991); Zapperi  _et al._ (1995); Watkins  _et al._ (2016) characterized by
critical avalanche dynamics exhibiting scale-free behavior.

The hypothesis of SOC in the deformation of amorphous materials inspired the
development of the theory of elasto-plasticity, bridging the gap between
macromechanical rheology models and SOC (see Nicolas  _et al._ (2018) for an
extensive review). Elastoplastic models (EPM), as defined in soft-matter
physics Nicolas  _et al._ (2018), represent amorphous materials as lattices of
bistable – either elastic or plastic – representative volume elements or
sites. Under strain or stress driving, each site experiences an affine elastic
deformation. When the local stress at a site overcomes a local threshold,
representing a local mechanical instability, the element undergoes a local
non-affine plastic-deformation event involving a predefined stress propagator
interacting with other elements and triggering an avalanche as a cascading
event. Conceptual mean field models interpret the residual stresses to local
instabilities as scalars undergoing a stochastic process Hébraud and Lequeux
(1998). Within this later category, the simplified and mathematically solvable
slip mean field theory (SMFT) Ben-Zion and Rice (1993); Dahmen  _et al._
(2009, 2011); Dahmen (2017), considering all-to-all homogeneous site
interactions, provides a remarkably apt conceptual model for a wide spectrum
of processes including the macromechanics and the statistical features of
externally measured avalanche dynamics in granular media Dahmen  _et al._
(1998); Mehta  _et al._ (2006); Dahmen  _et al._ (2011); Geller  _et al._
(2015); Denisov  _et al._ (2016); Zhang  _et al._ (2017). Such results
reinforced the idea of a predominant mean-field universality class in the
avalanche statistics Corwin  _et al._ (2005); Dahmen  _et al._ (2011); Geller
_et al._ (2015); Denisov  _et al._ (2016); Zhang  _et al._ (2017); Parisi  _et
al._ (2017); Nicolas  _et al._ (2018). In SMFT, SOC appears in the
thermodynamic limit as an asymptotic approach to a regime at stress value
$\sigma_{c}$ characterized by a divergence of scales in avalanche statistics.
Following scaling relations Dahmen (2017), the macromechanical evolutionary
path of strain-driven ($\epsilon$) amorphous solids leading to the steady-
state flow regime is uniquely represented in the SMFT by the stress-strain
relation:

$\sigma(\epsilon)\approx\sigma_{c}\left({1-\exp\left({-\epsilon}\right)}\right)$
(1)

Such a generalized macromechanical perspective is by construction insufficient
to describe the dynamics and statistics of slip avalanches as a function of
state Nicolas  _et al._ (2018). Both experiments and DEM simulations display a
non-monotonous evolution towards $\sigma_{c}$, implying that the relation
$\sigma(\epsilon)$ does not represent a function of state Karmakar  _et al._
(2010). Furthermore, initial configurations determine the $\sigma(\epsilon)$
response during such a transient regime and do not converge to $\sigma_{c}$
through the same evolutionary path (e.g. Fig. 1). A widespread scientific
endevour aims to understand the limitations and discrepancies of atomistic
DEM, coarse-grained full tensorial Lin  _et al._ (2014); Budrikis  _et al._
(2017); Karimi and Barrat (2018); Korchinski  _et al._ (2021), scalar Chen
_et al._ (1991); Baret  _et al._ (2002); Picard  _et al._ (2005), and mean-
field Hébraud and Lequeux (1998); Picard  _et al._ (2005); Agoritsas  _et al._
(2015); Lin and Wyart (2016); Ferrero and Jagla (2019) EPM under the lens of
the micromechanics of avalanche dynamics Rodney  _et al._ (2011); Budrikis
_et al._ (2017); Nicolas  _et al._ (2018). Other studies focusing on the
internal topology of the granular material suggest that the state of compact
granular media is better parameterized by a micromechanical description of the
particle configuration and the network of contact forces, which can be only
partially represented in terms of its macromechanics. Specifically, the role
of internal structures to the mechanics of granular media has been discussed
in early conceptual models Cates  _et al._ (1998); Liu and Nagel (1998),
experiments Miller  _et al._ (1996); Mueth  _et al._ (1998); Liu  _et al._
(1995); Cates  _et al._ (1998); Howell  _et al._ (1999); Andò  _et al._
(2012); Gendelman  _et al._ (2016); Berthier  _et al._ (2019a, b), DEM
Tordesillas (2007); Guo and Zhao (2013); Kuhn (1999); Pouragha and Wan (2017a)
and new mechanical-statistical theories of granular matter Edwards and
Oakeshott (1989); Blumenfeld and Edwards (2009); Blumenfeld  _et al._ (2012);
Bi  _et al._ (2015); Baule  _et al._ (2018); Bililign  _et al._ (2019); Ball
(2019); Bi  _et al._ (2011). Indeed, the definition of fragile matter is
intrinsically linked to the existence of force-chain structures Cates  _et
al._ (1998) and experiments reveal how the nucleation of failure events is
directly linked to _local_ topological measurements such as the betweenness
centrality Berthier  _et al._ (2019b). Parallel DEM studies identified a
stable evolution state (SES) in a space of independent structural variables.
The SES is reached after an adaptation transient, reached considerably before
the stationary flow regime observed in the macromechanical variables Pouragha
and Wan (2016). As an extended micromechanical representation of the jamming
transition, the SES describes a subset of micromechanical states upon which
the internal structure of granular media starts to evolve leading to stick-
slip like events or avalanches. In particular, the same spontaneous shear
localization leading to the non-monotonous macromechanical state evolution can
degenerate into different evolutionary paths towards the flow regime even
under similar intial conditions Masson and Martinez (2001); Desrues and
Viggiani (2004); Desrues and Hammad (1989); Pouragha and Wan (2017b), reaching
the same topological SES at different values of stress. The parametric study
in Pouragha and Wan (2016) shows that once SES is reached, regardless of the
initial condition and material parameters, the granular assembly starts to
exhibit irreversible deformations that are predominantly governed by micro-
avalanches.

Here, we analyze the statistics of kinetic avalanches identified during AQS
biaxial compression in DEM simulations of granular assemblies. We interpret
the results in the context of different hypothesis regarding the coarse-
grained description of amorphous materials. We show how the distribution of
avalanche sizes is independent of the macroscopic state of the system once SES
is reached, even when the macromechanical paths from different initial
conditions (e.g. loose vs. dense assemblies) differ significantly. More
importantly, the SES itself exhibits the properties of a self-organized
critical state. Thus, while this critical state includes the final steady-
state regime at the customary critical stress value $\sigma_{c}$, it is
reached significantly earlier. Because of that, rather than using the unique
term critical point to designate both, we will distinguish between the final
steady state regime at $\sigma_{c}$ and criticality as the statistical
mechanics concept refering to a divergence of susceptabilities. We also show
that while the statistical avalanche behavior is stationary within SES, each
individual avalanche modifies noticeably the state of the system, introducing
a characteristic delay, a.k.a. a pseudogap, in the interevent times between
consecutives avalanches. This leads to stronger correlations between the
states at which avalanches are triggered than between the states at which
avalanches stop and suggests that the SES defines an evolving stability limit
leading to the steady-state regime at the macroscopic scale rather than a
stability basin towards which the system is relaxed after each avalanche.
Moreover, we find non-universality in the power-law exponents defining the
statistics of the energy transfers involved in avalanches. In this case, we
expose a dependence of the exponents on the softness of the interactions
between particles. The external stress drop is found to be approximately
proportional to the drop of internal elastic energy caused by configurational
changes, but, due to the presence of soft lateral walls, this is not the case
for the energy dissipated by frictional forces and damping. The relation
between internal energy drops and the peaks of kinetic energy are consistent
with a specific set of assumptions within the SMFT formulation. This
consistency allows us to relate the non-universality of our results with a
potential transition between the expectations of SMFT and elasto-plasticity
based on the estimated critical exponents.

The manuscript is organized as follows: In section II we introduce our
numerical DEM simulations of granular assemblies, review the definition of SES
and explain how the definition of kinetic avalanches is used to implement the
quastistic biaxial compression driving. In section III we present the
statistical analysis of the avalanche properties, focusing on the average
relations between avalanche-size measurements, the existing correlations
between events and their full distributions with the corresponding fitted
exponents. In section IV we discuss the conditions and properties of avalanche
criticality and compare them with the mean-field solution of elastoplastic
models. Final conclusions are presented in section V.

##  II Background and Methods

###  II.1 Computational Model

The software package PFC2D Inc (2008) has been used here for discrete element
(DEM) simulation of 2d granular assemblies under AQS compression. A linear
contact model with constant stiffness is considered along normal and
tangential directions with their respective stiffnesses, $k_{n}$ and $k_{t}$
being the same and equal to $10^{9}N/m$. The tangential force is restricted by
a Coulomb friction law $f_{t}\leq\mu f_{n}$ with $\mu=0.5$ for all
simulations. A numerical damping of $5\%$ is also considered which acts as a
drag-like force on particles by artificially ramping changes in velocities
(see Inc (2008) for exact definition). The radii of circles are uniformly
distributed over $r_{\min}$ and $r_{\max}$ with $r_{\max}/r_{\min}=1.5$,
around an average radius of $1~{}mm$. To avoid overcrowding the figures and
for a better readability, the results in this paper are presented in their
dimensional forms with SI units. The material and loading parameters adopted
for numerical simulations are summarized in Table 1.

name |  num. of particles $N$ |  confining pressure $\sigma_{x}(N/m)$ |  driving rate $\dot{\epsilon}_{y}(\times 10^{-9}s^{-1})$ |  initial porosity $\phi_{0}$
---|---|---|---|---
D20kSc5 | 19520 | $10^{5}$ | 2.3 | 0.156
L20kSc5 | 19353 | $10^{5}$ | 2.3 | 0.190
D5kSc5 | 6374 | $10^{5}$ | 2.4 | 0.159
L5kSc5 | 5504 | $10^{5}$ | 2.4 | 0.192
D2kSc5 | 1593 | $10^{5}$ | 1.3 | 0.165
D5kSc6 | 6374 | $10^{6}$ | 2.4 | 0.154
D5kFc5 | 6374 | $10^{5}$ | 7.0 | 0.159
D5kFc6 | 6374 | $10^{6}$ | 7.0 | 0.154
D5kFc7 | 6374 | $10^{7}$ | 7.0 | 0.120
Table 1:  Simulation details for the nine numerical samples. The confining
pressure is kept constant, the driving rate is halted during the propagation
of avalanches. The initial porosity is measured at the beginning of the
simulation.

The rectangular sample has an initial aspect ratio of 2-to-1 which is confined
in between frictionless walls. During the initial so-called compaction stage,
the particles are slowly inflated, as a common numerical technique, to obtain
the target isotropic pressure on the walls. During this stage the
interparticle friction is artificially modified to achieve the desired initial
density. The friction is set to $0.5$ at the end of compaction. Next, the
deviatoric load is applied by moving the top and bottom walls towards each
other with a straining rate of $\dot{\epsilon_{y}}$ while the stress on the
lateral walls, $\sigma_{x}$, is kept constant through a servo-control
mechanism. The wall velocity is sufficiently small to keep the inertia number
well below the quasi-static limit Da Cruz  _et al._ (2005). Moreover, to
ensure a quasi-static loading, the driving is halted during surges in kinetic
energy, as is explained in detail later. The typical stress-strain responses
and internal displacement maps are exemplified by a pair of initially loose
and dense assemblies in Fig. 1. A relatively monotonic and homogeneous average
trend is seen for loose assembly (L5kSc5) while the dense assembly (D5kSc5)
experiences a peak of stress before dropping to the steady-state regime at
larger axial strains. This previous transient regime is linked to the
formation of shear localization bands (see right-bottom map in Fig. 1)

Figure 1:  (color online) Stress-strain response from simulation results
(sim.) of samples D5kSc5 and L5kSc5. Avalanche catalogs (cat.) only consider
the intervals represented in wider and darker lines. The black dotted line
represents a fitted SMFT prediction (1) to L5kSc5. The statistical analysis
represented in Fig. 5 correspond to the subcatalogs constrained within strain
intervals $A$ and $B$ of the evolution. Right figures represent the changes in
deviatoric strain ($\Delta\epsilon$) of both samples during the whole
simulation, showing strong localization in D5kFc5.

###  II.2 Internal dynamics

At any given instance, the energy of the granular system comprises the
following parts: The kinetic energy of particles, $E_{K}$, the internal
elastic energy stored at contact points, $U$, the dissipation due to
frictional and drag forces, $D=D_{f}+D_{d}$, and the external work by vertical
and lateral stresses, $W_{y}$ and $W_{x}$, respectively. The balance of
energies is given by $W_{y}+W_{x}=U+D+E_{K}$. For a system of $N$ particles,
the total instantaneous kinetic energy can be expressed as
$E_{K}(t)=\frac{1}{2}\sum_{i=1}^{N}(m_{i}v_{i}^{2}+I_{i}\omega_{i}^{2})$,
where $m_{i}$ and $I_{i}$ are the mass and moment of inertia of $i^{th}$
particle, respectively, and $v_{i}$ and $\omega_{i}$ are translational and
angular velocities. Considering that the velocity of the lateral walls is, on
average, proportional to the vertical driving rate, a characteristic value for
total kinetic energy of a homogeneous sample can be defined as:
$K_{D}:=\frac{1}{6}N\langle m\rangle(h\dot{\epsilon_{y}})^{2}$, where $h$ is
the sample’s height. We use $K_{D}$ as a reference parameter to define the
energy thresholds for avalanche detection.

The microscopic configuration of contact points creates a heterogeneous
landscape of elastic potentials. Even under monotonic driving the system
exhibits a non-linear mechanical response involving sharp local accelerations,
complex internal velocity fields and global bursts of internal kinetic energy
that cannot be simplified in terms of the external measures such as
$\dot{\epsilon_{y}}^{2}$ Radjai and Roux (2002). Considering a quastistically
slow driving, the system is episodically trapped in a series of metastable
jammed configurations well represented by the potential-energy Bonfanti  _et
al._ (2019). The external AQS driving pushes the system through local
mechanical instabilities that trigger local dynamics in the form of frictional
sliding between particles at contact points, or the collective buckling of
chained contact forces Cates  _et al._ (1998); Tordesillas (2007); Nicot  _et
al._ (2017). These initial instabilities can develop into collective internal
avalanches whose propagation is almost independent of the driving speed and is
controlled, instead, by particle-scale dissipative forces due to the drag-like
forces and inter-particle friction.

###  II.3 Stable evolution state

The mechanical evolution of the assembly is determined by a sequence of
episodes of quiescent deformation and sudden avalanches involving
configurational rearrangements. Recent studies Pouragha and Wan (2016, 2017a,
2017b, 2018a) suggested that mechanically stable configurations are
constrained by a stable evolution state surface (SES) defined as a manifold
observed in the space of three micromechanical variables: coordination number
(average number of contacts per particle), rigidity ratio (average contact
deformation relative to the particle size, which scales with
$(\sigma_{x}+\sigma_{y})/k_{n}$), and fabric anisotropy (a measure for
directional distribution of contacts), as represented in Fig. 2. The concept
of SES has been proven to be effective in developing new micromechanical
constitutive models for granular materials Pouragha and Wan (2018a).

Considering a dense sample, the initial state of the assembly is located on
the stable side (right side in Fig. 2) of the SES. As the sample is loaded
deviatorically, the microstructure quickly evolves through loss and gain of
contacts approaching the SES limit Pouragha and Wan (2017a). Upon reaching the
SES, the assembly appears to have just enough constraints to retain metastable
static equilibrium Pouragha and Wan (2016); Bonfanti  _et al._ (2019). As
such, the additional deviatoric straining introduced by the monotonic driving
acts as an initial perturbation pushing the assembly beyond the stability
limit and triggering the avalanche as a sudden increase in the internal
kinetic energy. Kinematic momentum is known to modulate the jamming transition
in shear processes Lemaître and Caroli (2009); LU  _et al._ (2007); Nicot  _et
al._ (2012); Gimbert  _et al._ (2013) and alters the stability conditions
during the transition. The kinematic momentum caused by avalanches is rapidly
dissipated, hence acting as a transient overshooting mechanism Maimon and
Schwarz (2004) pushing the system away from the fragile SES towards more
stable (non-fragile) configurations under quasistatic conditions, as discussed
in more detail later in Section III.4.

The resulting sequence of avalanches involves rearrangements adjusting the
number and orientation of contacts bringing gradually the micromechanical
state of the granular assembly towards the steady state at $\sigma_{c}$, where
the deviatoric straining is sustained with no appreciable change in the
micromechanical state of the assembly. The evolution within the SES towards
the final steady state is not necessarily unique and depends on the original
configuration and loading path. This is most clearly observed for the
initially loose cases whose initial state falls on SES, see Fig. 2.

Figure 2: (color online) Schematic visualization of SES surface and the
material state path for initially loose and dense samples. The micromechanical
state of the loose sample (red) is situated on the SES surface all along the
deviatoric loading history. The response of the dense sample (blue) starts
outside of SES but rapidly approaches the surface and remains on it
afterwards. The meta-stable, stick-slip like response in stress is observed
only when the micromechanical state of granular assembly reaches the SES.

###  II.4 Quasistatic driving and kinetic avalanches

Figure 3: (color online) Close-up example (D5kFc6) of the macromechanical
evolution of the system during the span of two consecutive avalanches as a
function of time. a) Evolution of the total kinetic energy ($E_{K}(t)$) used
to define the avalanches and the driving protocol, which halts the driving
above $K_{h}$ (upper dashed line) and restarts it when $E_{K}(t)>K_{r}$ (lower
dashed line). The middle dashed line indicates the expected average kinetic
energy during driving. The intervals of the avalanche relaxation time without
driving is indicated by the gray area. b) Energy transfer during driving and
avalanche propagation. The supplied work ($W_{y}$) is represented as a green
line. The green area represents the part of the work invested in lateral
expansion of the system against the confining pressure ($W_{x}$). The blue
line represents the rest of the work transformed into heat (dissipation $D$
and transient kinetic energy $E_{K}$) and the final variation of the internal
energy ($U$). Energy is dissipated by a medium drag-like force ($D_{d}$) and
friction between particles ($D_{f}$) represented as red and grey areas
respectively. c) Avalanches are always accompanied by a drop in the internal
energy that matches a drop in the external vertical stress ($\sigma_{y}$). We
use the derivative of this stress drop to measure the duration of the
avalanche. d) Shows the driving stopping during the relaxation time of the
avalanche. The relation between $\sigma_{y}$ (c) and $\epsilon_{y}$ (d) is
represented in Fig. 1 for two other simulations.

To guarantee a quasitatic condition and in particular, a clean separation
between independent avalanches, in our simulations, the driving is halted as
soon as an avalanche is detected. Such a detection occurs when the total
kinetic energy ($E_{K}(t)$) exceeds a selected threshold value $K_{h}=3K_{D}$.
A similar strategy was used by Salerno  _et al._ (2012); Salerno and Robbins
(2013) for the DEM simulations of amorphous materials. As an example, Fig. 3
shows the assembly magnitudes in a typical simulation (D5kFc6) during the time
interval containing 2 avalanches. The jumps in $E_{K}(t)$, up to several
decades (Fig. 3.a), signals the onset of avalanches triggered by the sudden
release of elastic potential energy (Fig. 3.b). Such a release of elastic
energy is signified also by the drop in the boundary stress, $\sigma_{y}$
(Fig. 3.c), while the lateral strain is adjusted to keep $\sigma_{x}$
constant. The surplus of kinetic energy is dissipated over time by frictional
sliding and drag-like forces (Fig. 3.a,b). Since the simulations are athermal,
the heat transfer due to dissipation compensates the loss in internal energy.
The numerical results exhibit dissipation intervals (grey areas in Fig. 3.a)
considerably longer than the duration $T$ of the first rebound ($\dot{U}\geq
0$) of the potential energy release (Fig. 3.c). Once $E_{K}(t)$ is damped
below $K_{r}=0.15K_{D}$, the external driving is resumed. The factors defining
$K_{h}$ and $K_{r}$ from the base $K_{D}$ have been selected to optimize the
resolution and the computational speed. The values of $K_{h}$ are just above
noise fluctuations in all of our simulations, and $K_{r}$ is low enough to
represent a precise measure of the dissipated energy guaranteeing the
quasistatic limit.

We obtain an avalanche catalog by labeling individual avalanches with a sorted
index $i$ and represent them as a marked point process. Instead of the instant
in simulation time $t_{0}$ when the $i$-th avalanche initiates at
$E_{K}(t)\geq K_{h}$, we select $\epsilon_{y}(t_{0})$ as our temporal scale,
which is stationary during the propagation of the avalanche because of the
quasistatic driving. We also record the “size” of each avalanche as measured
by different means. In particular, we focus in this study on the duration $T$
of the potential energy release $\Delta U:=U(t_{0})-U(t_{0}+T)$, the stress
drop $\Delta\sigma_{y}:=\sigma_{y}(t_{0})-\sigma_{y}(t_{0}+T)$, the peak in
kinetic energy $K=E_{K}^{\max}-K_{D}$, the total dissipated energy until the
quasi-static driving resumes, $D=D_{f}+D_{d}$, which accounts for the friction
at the contact points ($D_{f}$) and the drag-like forces ($D_{d}$), the work
against the lateral walls
$W_{x}=\int_{t_{0}}^{t_{0}+T}V(t)\sigma_{x}\dot{\epsilon}_{x}dt$, where
$V(t):=V_{0}\epsilon_{x}(t)\epsilon_{y}(t)=V_{0}(1+\epsilon_{v})$ is the 2d
volume of the assembly whose variation during a given avalanche can be safely
neglected. The term kinetic avalanches is used here to specify the avalanche
detection method and driving protocol, based on the internal kinetic energy.
However, following the observations presented below, we will argue that most
avalanche detection protocols based on other physical measures such as
external stress, or internal elastic energy drop would render similar results.

In the following sections, we study avalanches occurring within the SES,
including avalanches before reaching the steady state (see intervals A in Fig.
1). The behavioural characteristics of all assemblies at SES are argued to be
the same Pouragha and Wan (2016). We will verify that avalanche statistics are
stationary within the SES regime from the analysis of subcatalogs, such as the
intervals A and B represented in Fig. 1.

##  III Results

###  III.1 Energy partition and avalanche scaling relations

As illustrated for the interval represented in Fig. 3.b, the energy supplied
to the system as external work along the y-axis ($W_{y}=\int
V\sigma_{y}d\epsilon_{y}=\int V_{0}(1+\epsilon_{v})\sigma_{y}d\epsilon_{y}$)
is mostly dissipated during the avalanche processes ($D$) or transferred to
the lateral walls ($W_{x}$). The difference corresponds to variations in the
elastic potential energy $\Delta U$. Instead of following a random sequence of
drops and loads, $U$ self-organizes towards a stationary value — determined by
the flow regime in terms of elastoplasticity, or the SES in terms of granular
microstructure — and fluctuates around that value through a sequence of
loading and avalanche cycles.

Figure 4: (color online) Relation between normalized stress released
($\Delta\sigma_{y}/\sigma_{x}$) and (a) average elastic potential released
($\langle\Delta U\rangle(\Delta\sigma_{y})$) or (b) dissipated energy
($\langle D\rangle(\Delta\sigma_{y})$) for stress released. Error-bars
correspond to one standard deviation. Dashed lines represent an asymptotic
relation between the energy magnitudes and $N\sigma_{x}^{0.85}$ which fits the
results as discussed in section III.1.

This metastable cascade appears also in the external macromechanics of the
assembly, i.e. in the stress drop ($\Delta\sigma_{y}$). For kinetic
avalanches, Fig. 4 tests the linearity between the external $\Delta\sigma_{y}$
and the internal elastic potential drop averaged at $\Delta\sigma_{y}$
windows, $\langle\Delta U\rangle(\Delta\sigma_{y})$ (left), as well as the
average dissipated energy, $\langle D\rangle(\Delta\sigma_{y})$ (right).
Recalling that the lateral stress is kept constant, $\sigma_{x}$ is considered
as the characteristic stress measure and, as such, the normalized stress
release ($\Delta\sigma_{y}/\sigma_{x}$) is used for describing the avalanche
size. On average, and for sufficiently large avalanches, $\Delta U$ is
approximately proportional to $\Delta\sigma_{y}$ in all simulations. This is
different from the behavior of $\langle D\rangle$: For most simulations,
$\langle D\rangle$ is not proportional to $\Delta\sigma_{y}$, though such a
proportionality has been observed in other works Salerno  _et al._ (2012). The
difference between $\langle D\rangle(\Delta\sigma_{y})$ and $\langle\Delta
U\rangle(\Delta\sigma_{y})$ reflects the contribution of the lateral walls to
the energy balance ($W_{x}$), which appears to be more prominent for large
avalanches but disappears at high confining pressures. In addition, we observe
a scaling factor, matching an apparent relation $\Delta U\propto
N\sigma_{x}^{1.85(5)}\Delta\sigma_{y}/\sigma_{x}$. This relation is
represented as horizontal dashed lines in Fig. 4. This scaling relation is
reasonably close to the expectations from a dimensional analysis: Following a
linear elastic relation, the change in elastic energy can be written in terms
of the wall stresses and the normalized stress released:

$\Delta
U=V\left({\frac{\sigma_{y}}{\sigma_{x}}\frac{1}{E_{y}}-\frac{\nu_{xy}}{E_{x}}}\right)\sigma_{x}^{2}\frac{\Delta\sigma_{y}}{\sigma_{x}},$
(2)

where $V:=l_{x}l_{y}\propto N\pi\langle r^{2}\rangle/(1-\phi)$ is the system’s
volume with $\phi$ being the porosity, $E_{x}$ and $E_{y}$ the macroscopic
Young’s moduli, and $\nu_{xy}$ the associated Poisson ratio. The macroscopic
Young’s moduli also depends on coordination number ($Z$) and the total volume,
$E\propto Z\/V^{-1}\propto Z(1-\phi)$ Chang  _et al._ (1995); Pouragha and Wan
(2018b). Taking into account that the stress ratio
$(\frac{\sigma_{y}}{\sigma_{x}})$ fluctuates in no more than 1/3 of its
average value within the catalog (Fig. 1), Eq. 2 can be simplified to give the
following scaling:

$\Delta U\propto
Z^{-1}(1-\phi)^{-2}\sigma_{x}^{2}\frac{\Delta\sigma_{y}}{\sigma_{x}}.$ (3)

The remaining mismatch between the exponent ${1.85(5)}$ used in the scaling of
Fig. 4 and the exponent 2 in Eq. (3) is found in an effective relation between
$Z(1-\phi)^{-2}$ and $\sigma_{x}$ which renders $\Delta
U\propto\sigma_{x}^{1.91(4)}{\Delta\sigma_{y}}/{\sigma_{x}}$ (see Fig. S2 in
supplementary material).

###  III.2 Distribution of avalanche sizes

Figure 5:  (color online) Distributions of (a) kinetic energy peaks ($K$) (b)
internal energy drops ($\Delta U$), for L5kSc5 and D5kSc5 separated into
stages A and B (see Fig. 1) compared with the total sequences in the SES
regime. Exponent values ($P(K)\sim K^{-4/3}$ and $P(\Delta U)\sim\Delta
U^{-3/2}$) expected from SMFT are shown for comparison (black lines).

In addition to the scaling relations between different measures of avalanche
sizes, the distribution of these sizes close to criticality also satisfy
scaling relations of interest for the comparison with other avalanche
phenomena and modelling approaches (see section IV). Fig. 5 shows the binned
avalanche distributions for $K$ and $\Delta U$ on a log-log scale for a loose
and a dense assembly of $N\sim 5000$ particles under similar external and
internal conditions during different time intervals. In particular, the
external stress in the assemblies is considerably different during interval
$A$, and evolves towards almost indistinguishable stress values ($B$) as both
systems approach the steady-state regime (see Fig. 1). Note that in a
macromechanical formulation, the state of the material is different during
strain hardening (interval $A$ for the loose case), strain softening (interval
$A$ for the dense case), and close to the steady state (interval $B$).
However, despite such differences in stress state and the softening/hardening
rate, no pronounced differences between the distributions are present in Fig.
5. Both distributions of $\Delta U$ and $K$ appear to be compatible with the
scale-free hypothesis, i.e., the distributions follow a power law decay,
$P(\Delta U)\sim\Delta U^{-\kappa}$ and $P(K)\sim K^{-\varepsilon}$, to a good
approximation over a wide range of scales, although finite size effects and
limited statistics might play a role. The collapse of data with different
initial densities and stress values further confirms that the state of the
material associated with meta-stable avalanches is not sufficiently
represented by macromechanical state variables, such as stress ratio and
porosity, as commonly used in coarse-grained models. Considering that
avalanche statistics can be represented in terms of state variables one needs
a higher-dimensional space capable of representing all states with coinciding
statistical properties, such as intervals $A$ and $B$ of dense and loose
assemblies, in a unique manifold. In this regard, here we propose the SES as a
viable candidate for such representation, since all evaluated intervals fall
inside the SES manifold.

Figure 6:  (color online) Distributions of (a,b) internal energy drops
($\Delta U$); (c,d) kinetic energy peaks ($K$) for all simulations. All
distributions are shown on a regular scale (a,c) and rescaled to fit the
apparent scaling relations in the x-axes (b,d). The y-axes are scaled to
collapse the distributions to a single curve in the region of large values.
Vertical dashed lines denote the intervals used to estimate the power-law
exponents (dotted lines correspond to the lower limit for $\sigma_{x}\geq
10^{6}$).

To further investigate the collapse of the avalanche data, we focus on the
standard scaling ansatz for self-similar behavior:

$P(x)dx=x^{-\tau_{x}}\Phi_{x}(x/x^{*})dx,$ (4)

where $x$ is a measure of the avalanche size, $\tau_{x}$ is the associated
critical exponent, $\Phi_{x}$ is a universal scaling function and $x^{*}$ is a
characteristic avalanche size. This characteristic $x^{*}$ should be
controlled by the finite system size as well as by other parameters. We
attempt the following general ansatz motivated by the usual relations found in
finite size scaling of scale-free phenomena Privman (1990)

$x^{*}(N,\sigma_{x},\dot{\epsilon_{y}})=\tilde{x}^{*}N^{\gamma_{N}^{x}}\sigma_{x}^{\gamma_{\sigma}^{x}}\dot{\epsilon}_{y}^{\gamma_{\epsilon}^{x}}.$
(5)

Figure 7: (color online) Average relations between the potential energy drop
($\widetilde{\Delta U}$), kinetic peak ($\widetilde{K}$) and duration of the
energy drop ($T$). Units have been rescaled according to the exponents
represented in Fig. 6.b,d, see main text for details. Black dashed lines:
fitted power-law in the intervals represented in Fig. 6.b,d. (a) Mean kinetic
peak ($\langle\widetilde{K}|\widetilde{\Delta U}\rangle$) for all simulations.
Averages scale to a single relation with exponents close to an average
$\gamma=1.71(9)$. (b) Mean duration ($\langle T|\widetilde{\Delta U}\rangle$)
for simulations D5kFc5, D5kSc6, D5kFc6, D5kFc7, averages scale as $\langle
T|\widetilde{\Delta U}\rangle\sim\widetilde{\Delta U}^{0.177(9)}$.

If the scale-free hypothesis is valid, $\Phi_{x}$ should be universal and all
the distributions of $x$ obtained with different parameters should collapse
when rescaled by the characteristic event size, i.e. the distribution of
$\tilde{x}:=x/x^{*}$ should be independent of system size and parameters. Fig.
6 shows the distribution of $\Delta U$ (a,b) and $K$ (c,d) without scaling
(a,c) and scaled (b,d) according to Eq. 5 using the best exponents found to
collapse the exponential cutoffs into a single curve on visual inspection.
Note that the y-axes are normalized independently, since the normalization of
the distributions are non-trivially affected by the lower cutoffs. We find a
reasonable collapse for ${\gamma_{\sigma}^{\Delta U}}=1.85$ for potential
energy drops and ${\gamma_{N}^{\Delta U}}=0.50$ corresponding to a fractal
dimension $d_{f}^{\Delta U}:=d\times\gamma_{N}^{\Delta U}=1$, where $d$
corresponds to the dimension of the specific system, which is two for our
model. The dependence on $\sigma_{x}$ is reasonable from the analysis in
section III.1, if we assume that the intrinsic measure
$\Delta\sigma_{y}/\sigma_{x}$ is also characteristic of the avalanche
propagation at a given simulation size Lerner and Procaccia (2009). For the
peaks of kinetic energy ($K$) we find a good agreement with
${\gamma_{\sigma}^{K}}=1.5$, and ${\gamma_{N}^{K}}=0.0$. Therefore, $K$
appears to be an intensive quantity. We argue that the instantaneous $E_{k}$
is probably dominated by a few fast-moving particles at a size-independent
characteristic speed. In agreement with the quastatistic hypothesis, we do not
observe any apparent relation with the speed of the external driving
($\dot{\epsilon}_{y}$), such that ${\gamma_{\epsilon}^{x}}=0$, except for a
shortening of the scale-free regime. The latter is caused by the lower bound
of the power-law behavior controlled by the applied threshold $K_{h}$ in the
definition of the avalanches, which results in a sharp cutoff in $K$ which in
turn translates into a smoother transition in $\Delta U$. The distribution of
avalanche durations are presented in the Fig. S1 of the Supplementary
Material. Such distributions appear to be compatible with a power-law
distribution $P(T)\sim T^{-\alpha}$ with exponent $\alpha\approx 2$. Yet, the
power-law range is too short to provide a reliable estimate of the exponent.

| stiffness |  $\\#$ events | $\kappa$ | $\varepsilon$ | $\gamma$
---|---|---|---|---|---
D2kSc5 | (stiff limit) | 1684 | 1.62(10) | 1.32(10) | 1.95(5)
L5kSc5 | ” ” | 979 | 1.60(7) | 1.34(5) | 1.85(10)
D5kSc5 | ” ” | 788 | 1.71(8) | 1.33(6) | 1.83(4)
L20kSc5 | ” ” | 130 | 1.77(17) | 1.45(6) | 1.85(15)
D20kSc5 | ” ” | 236 | 1.49(11) | 1.36(4) | 1.69(6)
D5kFc5 | ” ” | 1215 | 1.46(6) | 1.36(4) | 1.71(5)
D5kSc6 |  | 396 | 1.41(8) | 1.14(11) | 1.65(8)
D5kFc6 |  | 851 | 1.32(5) | 1.14(6) | 1.71(4)
D5kFc7 | (soft limit) | 633 | 1.08(3) | 1.02(8) | 1.48(7)
SMFT(1) |  |  | 1.5 |  $1+\frac{\kappa-1}{2-\varsigma\nu z}=$ 1.33 |  $2-\varsigma\nu z=$1.5
SMFT(2) |  |  | 1.5 |  $(\mu+1)/2=$1.5 |  $2\varsigma\rho=$1
2D EPM |  |  |  1.25–1.28 |  $\sim$1.2 Budrikis  _et al._ (2017) |  $\sim$1.45Budrikis  _et al._ (2017)
Table 2:  Summary of the number of kinetic avalanches and critical exponents
estimated from the simulations. Errors in brackets represent the std.
deviation estimate from the likelihood function given a Gaussian
approximation. For comparison, the predictions of SMFT theory considering
$K=E$ (SMFT(1)) and $K=E_{m}$ (SMFT(2)) Dahmen (2017) (see Section IV) are
also shown together with published findings of 2D full tensorial FEM
simulations Budrikis  _et al._ (2017) and EPM simulations Nicolas  _et al._
(2018). Note that similar exponents are also found in some DEM simulations of
shear processes Liu  _et al._ (2016); Zapperi  _et al._ (1997) while other DEM
simulations return values as low as $\kappa=0.98(1)$ Shang  _et al._ (2020).

The relation between ${\Delta U}$ and $K$ of individual avalanches agrees with
the scaling relation found in Fig. 6. In particular, Fig. 7 shows the average
duration $\langle T|\widetilde{\Delta U}\rangle$ (a) and the rescaled kinetic
peaks $\langle\widetilde{K}|\widetilde{\Delta U}\rangle$ (b) expressed in bins
of $\widetilde{\Delta U}$. Over a range of a factor 25 in $\Delta U$ (x-axis
in Fig. 7) and 150 in $K$ (y-axis in Fig. 7.b), the conditional averages for
all simulations collapse and agree with a power-law relation between avalanche
energies: $\langle\widetilde{K}|\widetilde{\Delta
U}\rangle\propto\widetilde{\Delta U}^{\gamma}$ with $\gamma=1.71(9)$. As a
particular case, a clear lower exponent is found at high confining pressures
(D5kFc7). In a similar way, we can assess the average relation between
$\widetilde{\Delta U}$ and the duration of the avalanche $T$, rendering a
power-law relation $\langle T|\widetilde{\Delta
U}\rangle\propto\widetilde{\Delta U}^{0.20(3)}$. Notice that the values of $T$
have not been scaled since no clear parameter nor size dependencies have been
observed in its distribution.

Visually, the exponent values seem to be in reasonable agreement with the
expected values of the slip mean field theory (SMFT), $\varepsilon=4/3$,
$\kappa=3/2$, represented as dashed lines in Fig. 6 and the $\gamma=3/2$, see
Section IV for an extended discussion. However, a closer inspection reveals a
clear dependency on confining pressure. In order to estimate their actual
values we evaluate the avalanche sizes in their rescaled form considering the
exponents used in Fig. 6 and the ranges of scale invariance shown as black
segments in Fig.7, to mitigate the effect of $K_{h}$ on $\Phi_{x}$.

Figure 8: (color online) Estimated exponents $\kappa$ (blue triangles) and
$\varepsilon$ (red circles) as function of the confining pressure within the
intervals represented in Figs. 6 and 7. Errorbars represent one standard
deviation of the normal approximation to the likelihood function. Smaller
points represent the estimation of each of the exponents with the help of Eq.
(6) by using the other exponent and the estimated $\hat{\gamma}$ (see Table 2)
to check the consistency of the results. Dashed lines represent the expected
values of $\kappa$ (blue) and $\varepsilon$ (red) from SMFT.

The exponents $\kappa$ and $\varepsilon$ characterizing the distributions of
$\Delta U$ and $K$, respectively, have been estimated using the maximum
likelihood method following the recipe in (Baró and Vives, 2012) within the
intervals $-13.05<\log_{10}({\widetilde{\Delta U}})<-12.05$ and
$-11.3<\log_{10}({\tilde{K}})<-9.8$ represented in Fig. 7. For simulations at
high confining pressures ($\sigma_{x}\geq 10^{6}$) we increase the fitting
ranges to $-13.55<\log_{10}({\widetilde{\Delta U}})<-12.05$ and
$-11.8<\log_{10}({\widetilde{K}})<-9.8$ to take advantage of the extended
power-law regime to lower values (see individual distributions and selected
estimation ranges in Fig. S4, S5 of the supplementary material). All estimated
exponents are summarized in Table 2 and represented in Fig. 8 as a function of
the confining pressure $\sigma_{x}$. Table 2 includes the least square fit of
$\gamma$ for each experiment and Fig. 8 represents as well the expected
exponents considering the approximated scaling relation
$\widetilde{K}(\widetilde{\Delta
U})\approx\langle\widetilde{K}|\widetilde{\Delta
U}\rangle\propto\widetilde{\Delta U}^{\gamma}$, which renders the scaling
relation between exponents:

$\gamma=\frac{1-\kappa}{1-\varepsilon}.$ (6)

This relation can be deduced from Eq. (4). The maximum likelihood estimation
method reveals that both exponents $\hat{\kappa}$ and $\hat{\varepsilon}$
decrease with the confining pressure, from $\hat{\kappa}=1.58(11)$ and
$\hat{\varepsilon}=1.36(4)$ for $\sigma_{x}=10^{5}$ to $\hat{\kappa}=1.14(6)$
and $\hat{\varepsilon}=1.02(8)$ for $\sigma_{x}=10^{7}$. The exponent
$\hat{\gamma}$ also decreases from $1.82(10)$, close to the value $1.71(9)$
found by averaging all the conditional average dependencies of Fig. 7, to the
value $1.48(7)$ for $\sigma_{x}=10^{7}$ (golden lines in Fig. 7). In general,
the maximum likelihood method provides a more robust estimation of the
exponent $\varepsilon$ than $\kappa$. This is probably a consequence of the
extended scale-free regime in $K$, given the lower exponent values and the
sharper cutoff at $K_{h}$. The scaling relation captured by Eq. (6) is
consistent with the estimated exponents when used to determine
$\kappa(\hat{\varepsilon},\hat{\gamma})$ and
$\varepsilon(\hat{\kappa},\hat{\gamma})$, especially at low confining
pressures. The consistency between the two estimation methods asserts the
robustness of the results and the existence of an actual dependence beyond
statistical artifacts that could alter the exponent estimation within the
fitting intervals. We notice as well that this agreement survives at high
confining pressures, close to the singular undetermined point at
$\varepsilon\approx\kappa\approx 1$. Such values are close to the values
reported in DEM simulations with Lennard-Jones (LJ) potentials Shang  _et al._
(2020). Note that lower scaling exponents ($\tau_{x}<1$) are not expected in
the framework of scale invariance since they would imply a divergence of
energy released for large system sizes. However, it is possible that
_estimated_ exponents ($\hat{\tau}_{x}$) in _finite_ systems are lower even
though the true scaling exponent is $\tau_{x}=1$. This difference corresponds,
in that case, to the form of the scaling function $\Phi_{x}(x/x^{*})$
Christensen  _et al._ (2008). Variations of exponents as function of input
parameters or conditions are fairly common in DEM Salerno and Robbins (2013);
Karimi and Barrat (2016) and FEM Budrikis  _et al._ (2017) simulations of
amorphous materials exhibiting SOC, but they are seldom observed in conceptual
models and they are usually linked to the superposition of different
subcritical processes Amitrano (2003); Lord-May  _et al._ (2020), the sweeping
through criticality Didier Sornette (1994), or rare regions, a.k.a Griffiths
phases Vazquez  _et al._ (2011).

###  III.3 Temporal correlations

Figure 9: (color online) (a,b) Distribution of interevent times as measured by
the strain difference $\Delta\epsilon_{y}$ for different $K$-thresholds
(corresponding to $100\%$, $\sim 50\%$, and $\sim 10\%$ of the events) and
scaled by average $\langle\Delta\epsilon_{y}\rangle(K)$ in (a) log-log and (b)
linear-log scales, exemplified by three different simulations. For comparison,
we show the Weibull renewal model (black lines) predicted by EPM (see Eq. (9))
with $\theta=1.0$, $\theta=0.5$ and $\theta=0.0$, the latter corresponding to
a homogeneous Poisson process. (c) Kolmogorov–Smirnov test over the $H$-values
obtained by the Bi-test for Poisson rejection (see text) where $n$ is the
number of interevent times. Color bands represent confidence intervals for
rejection as indicated in the labels.

Considering the selected AQS driving protocol, the time between avalanches is
best defined in terms of the strain component $\epsilon_{y}$, since
$\dot{\epsilon_{y}}$ is constant during external driving and $0$ during the
avalanche propagation intervals, i.e., avalanches are instantaneous in terms
of $\epsilon_{y}$. Thus, interevent times are represented as
$\Delta\epsilon_{y}$. Correlations between avalanches are often observed in
natural mechanical processes. In particular, we are interested in the presence
of temporal clustering, linked to triggering phenomena, which is often
observed in natural systems Utsu  _et al._ (1995); Weiss and Carmen Miguel
(2004); Crane and Sornette (2008); Sornette and Utkin (2009); Baró  _et al._
(2013); Hainzl  _et al._ (2014); Baró  _et al._ (2014); Davidsen  _et al._
(2017); Kumar  _et al._ (2020) but seldom reproduced in conceptual
constitutive models unless explicitly introduced Dieterich (1972); Hainzl  _et
al._ (1999); Mehta  _et al._ (2006); Zöller  _et al._ (2006); Ben-Zion and
Lyakhovsky (2006); Jagla (2014); Zhang and Shcherbakov (2016); Baró and
Davidsen (2018). Conversely, conceptual and numerical elastoplastic models
identify a power-law relation in the distribution of the local residual-
stresses or stability limits close to the steady state, a.k.a. pseudogap
(Hébraud and Lequeux, 1998; Picard  _et al._ , 2005; Karmakar  _et al._ ,
2010; Lin  _et al._ , 2014; Agoritsas  _et al._ , 2015) leading to
quasiperiodicity in the strain increments between consecutive avalanches Lin
and Wyart (2016). Recent studies revealed that the pseudogap is more prominent
for large system sizes Ruscher and Rottler (2020); Ferrero and Jagla (2019);
Korchinski  _et al._ (2021), contrasting with the naive expectation of
quasiperiodicity being caused solely by finite-size effects. We address this
question in section IV.2 through the analysis of $\Delta\epsilon_{y}$
presented here.

Fig. 9 shows the statistical analysis of the time between avalanches as
measured by the axial strain increment $\Delta\epsilon_{y}$ for selected
samples L5kSc5, D5kSc5, D2kSc5 used here as examples. Each of the lines
represent the strain increments between consecutive events with energy $K$
above different thresholds values. The distribution of interevent times in
linear-log scale (Fig. 9.b) follows an exponential distribution for long
interevent times, suggesting that avalanches are independent of one another
and follow a Poisson process with a stationary characteristic avalanche rate
Touati  _et al._ (2009); Baró  _et al._ (2013); Baró and Davidsen (2017).
Similar results have been recently reported by Kuhn et al. (2019) (Kuhn and
Daouadji, 2019) for non-spherical particles. In contrast, self-excitations
usually involve strong correlations between large events Utsu  _et al._
(1995); Baró  _et al._ (2013) which would be easily identified as an excess of
short interevent times. No such an excess is observed in Fig. 9(a). Therefore
we can rule out the presence of self-excitations in our simulations. On the
contrary, most thresholds values reveal a deficit of short interevent times
compared to a homogeneous Poisson process in all simulations. This implies
regularity, i.e. the existence of characteristic (quasiperiodic) interevent
times, which might be linked to finite size effects or the aforementioned
pseudogap.

This deviation from Poisson behavior in the avalanche process can also be
verified by the rejection method first proposed in (Bi  _et al._ , 1989), so-
called Bi test, consisting of the evaluation of the statistical variable
$h_{i}=\delta^{(1)}_{i}/(\delta^{(1)}_{i}+0.5\delta^{(2)}_{i})$, where
$\delta^{(1)}_{i}$ is the shortest of the intervals preceding or following
event $i$ (i.e.,
$\delta^{(1)}_{i}:=\min(\epsilon_{y,i+1}-\epsilon_{y,i},\epsilon_{y,i}-\epsilon_{y,i-1})$)
in terms of strain and $\delta^{(2)}_{i}$ is the next interevent time in the
same temporal direction of $\delta^{(1)}_{i}$ (i.e.
$\delta^{(2)}_{i}:=\epsilon_{y,i-1}-\epsilon_{y,i-2}$ if
$\delta^{(1)}_{i}=\epsilon_{y,i}-\epsilon_{y,i-1}$ and
$\delta^{(2)}_{i}:=\epsilon_{y,i+2}-\epsilon_{y,i+1}$ if
$\delta^{(1)}_{i}=\epsilon_{y,i+1}-\epsilon_{y,i}$ ). The Kolmogorov-Smirnov
test over the $h_{i}$ values shown in Fig. 9.c displays a rotated ’S’ profile
characteristic of regular processes, in opposition to the a rotated ’Z’
profile, which is associated with clustering or self-excitation mechanisms
Baró  _et al._ (2014). The ’S’ profile is clear in all simulations (see Fig.
S3 in the supplementary material). The Poisson hypothesis is rejected at 95%
confidence interval in all simulations and at 99.95% in most of them,
including D5kSc5 and L5kSc5. In the case of D2kSc5 the $F_{n}(H)$ distribution
falls short to reject the Poisson hypothesis at the 99.95% confidence
interval. This discrepancy is in line with the better agreement with a Poisson
process observed for the distribution of waiting times of D2kSc5 (Fig. 9). We
verified numerically that the deficit of short interevent times is more
pronounced in larger simulations. In particular, the fraction of rescaled
interevent times below the expected first quartile of a normalized exponential
distribution is $\sim 0.2$ for $\sim 2000$ particles, $\sim 0.13$ for $\sim
5000$ particles and $\sim 0.08$ for$\sim 20000$ particles with respect to the
expected $0.25$ if the process were Poisson (see Fig. S6 in the supplementary
material). Unfortunately, since the threshold selection to define avalanches
is size dependent, we cannot straightforwardly infer the relation
$\langle\delta(N)\rangle$ which is expected to go as $N^{d/d_{f}-1}$ Lerner
and Procaccia (2009).

Figure 10: (color online) (a) The macromechanical evolution of stress
($\sigma_{y}$) vs strain ($\epsilon_{y}$) (blue line). The dashed and dotted
lines connect the values of stress on the onset ($t^{-}$) and at the end
($t^{+}$) of each avalanche $i$, respectively, for visual purposes. (b)
Autocorrelation functions for the sequence of $\sigma_{y}(t^{-})$ and
$\sigma_{y}(t^{+})$, respectively. Lines represent best fits to a stretched
exponential decay over the interval $j\in[1,30]$ (c,d) Scatter plots of the
relations between the size of the (c) previous and (d) next avalanche
($\Delta\sigma_{y}$) and the interevent strain $\Delta\epsilon_{y}$ for three
different simulations. Conditional averages represented in colored lines and
linear fitting in black. Dashed line in (c) represents an empirical minimal
loading period $\Delta\epsilon\sim\Delta\sigma_{y}^{0.36(8)}$.

###  III.4 Stress correlations

Fig. 1 shows that the macromechanical state of the system is constrained
within an evolving narrow stripe following a well-defined path at the large
scale. However, the commonly used smooth average stress-strain relation is
clearly a broad approximation. At the scale displayed in Fig. 3 the intrinsic
metastable dynamics arising from the competition between the slow driving and
the fast relaxation in the form of avalanches exhibits its non-smooth
behaviour. As an example, Fig. 10.a shows a close-up to the macromechanics of
one of the simulations. We mark the stress ($\sigma_{y}$) at the onset of the
avalanche ($\sigma_{y}(t^{-})$) and at the end of the avalanche
($\sigma_{y}(t^{+})$). Fig. 10.a also shows an example of interevent loading
time $\Delta\epsilon_{y}^{k}=\epsilon_{y,k+1}-\epsilon_{y,k}$, preceded by an
avalanche $k$ with stress drop $\Delta\sigma_{y,k}$ and followed by an
avalanche $k+1$ with stress drop $\Delta\sigma_{y,k+1}$. The regularity in the
distribution of interevent times exposed in Fig. 9 could be a consequence of a
characteristic release or reloading of stress correlated to the interevent
times, similar to predictions of the elastic rebound theory Reid (1910) and
related earthquake models Zöller  _et al._ (2006). Fig. 10.b shows the
autocorrelation function (ACF) of the sequences of $\sigma_{y}(t^{+})$ and
$\sigma_{y}(t^{-})$, respectively, for a pair of loose and dense assemblies
with similar simulation parameters. Since the macromechanics of the system are
relatively smooth at the long range, as shown in Fig. 1, both autocorrelations
persist for long times and remain at significant positive values over the
shown range. In particular, correlations decay to a good approximation as a
stretched exponential $C\approx A\exp{(-j^{\xi})}$ for $j<30$ with
$\xi=0.598(16)$ for $\sigma(t^{-})$ and a faster $\xi=0.317(18)$ for
$\sigma(t^{+})$. Therefore, $\sigma(t^{-})$ values are significantly more
persistent than $\sigma(t^{+})$, suggesting that the state triggering the
avalanches is more regular than the state arresting it. This hypothesis is
confirmed by a comparison between the scatter plots of $\Delta\sigma_{y,k}$
and $\Delta\sigma_{y,k+1}$ as a function of $\Delta\epsilon_{y}$ (Fig. 10.c
and d). The duration of an interevent loading time does not affect the size of
the next avalanche (Fig. 10.d). Instead, the size of an avalanche correlates
with the next interevent loading time (Fig. 10.c). Specifically, after an
event with stress drop $\Delta\sigma_{y}$, the system needs to load a
characteristic minimum amount of strain
$\Delta\epsilon^{\min}\propto\Delta\sigma_{y}^{0.36}$ to trigger the next
avalanche. This so called time-predictability Shimazaki and Nakata (1980) is
consistent with the lack of short interevent times (Fig. 9.a). In terms of the
SES representation, and at least for the testing procedure adopted in this
study, the upper values of stress, i.e. $\sigma(t^{-})$, represent a better
candidate for an objective reference state, and can be associated with SES.
The fact that the interevent times are correlated with their preceding events
implies that typically the avalanches can be considered as excursions bringing
the system away from the critical SES limit towards more stable
configurations. The following loading has then a tendency to return the system
to the SES.

We have shown that both distributions of avalanche sizes and interevent times
are stationary once the SES is reached. On the other hand, the effective
macroscopic Young’s moduli $E_{y}:=\dot{\sigma}_{y}(\epsilon_{y})$ within the
SES is non-stationary and evolves towards the steady state in non-unique
evolutionary paths (Fig. 1). This implies that a time dependence exists in the
average stress reloading amount between consecutive avalanches
($\langle\Delta\sigma_{y}^{r}\rangle:=\langle\sigma_{y,k}(t^{-})-\sigma_{y,k-1}(t^{+})\rangle$)
or the avalanche rates ($r:=N_{\mathrm{aval}}/\Delta\epsilon_{y}$), even
within the SES. To verify this hypothesis we evaluate such average reloading
sizes in addition to the the effective elastic modulus
($\hat{E}_{y}:=\Delta\sigma_{y}/\Delta\epsilon_{y}$) and average avalanche
sizes
($\langle\Delta\sigma_{y}\rangle:=\langle\sigma_{y,k}(t^{+})-\sigma_{y,k}(t^{-})\rangle$)
in relatively broad intervals of strain, i.e. driving, of
$\Delta\epsilon=0.005$ containing around 70-100 avalanches each. We build a
catalog considering all strain intervals in L5kSc5 and D5kSc5 within SES
together, since the input parameters are similar and include a large
variability in $E_{y}$, mostly positive in L5kSc5 and mostly negative in
D5kSc5 (see Fig. 1). We report the following Pearson correlation coefficients
with respect to $\hat{E}_{y}$:

$\rho_{\hat{E}_{y},r}=0.37\;,\;\;\rho_{\hat{E}_{y},\langle\Delta\sigma_{y}^{r}\rangle}=0.30\;,\;\;\rho_{\hat{E}_{y},\langle\Delta\sigma_{y}\rangle}=0.058$
(7)

the latter being effectively null. Therefore, we argue that the effective
macroscopic Young’s modulus is determined by the rate of avalanches and the
reloading between avalanches in a way that low rates and high stress reloading
correlates with positive effective moduli, while the characteristic size of
avalanches is effectively constant, as already argued in section III.2. Some
of the issues pertaining to time averaging such meta-stable quasistatic
behaviours are discussed in Combe and Roux (2000).

##  IV Discussion

Non-linear mechanics and avalanches in amorphous solids are often studied
under the lens of AQS shear driving. Both atomistic and coarse-grained
simulations usually consider overdamped systems of soft particles, usually
modeled with attractive LJ potentials in DEM models, returning short-range or
far-field homogeneous stress transfer functions in a coarse-grained
representation, either isotropic when considering interactions along the
maximum shear plane Ben-Zion and Rice (1993); Dahmen  _et al._ (2011), or
polar Eshelby inclusions, when considering off-plane interactions Hébraud and
Lequeux (1998); Picard  _et al._ (2004); Albaret  _et al._ (2016). The DEM
simulations of compact granular assemblies presented here differ from some of
the previous considerations. Most notably, the system is driven by biaxial
compression, instead of shear, particles have inertia and the stress
distribution is channeled in complex force-chain networks resulting in highly
disordered long-range interactions.

###  IV.1 Biaxial compression

Biaxial and triaxial compression experiments Lockner (1993); Baró  _et al._
(2013); Lennartz-Sassinek  _et al._ (2014); Davidsen  _et al._ (2017); Baró
_et al._ (2018), and numerical simulations Amitrano (2003); Kun  _et al._
(2014) such as the ones presented in this work, serve as analogues for
geomechanical processes. Mechanics at the mesoscale are dominated by internal
shear processes introduced by residual deviatoric stresses at contact points.
The correspondence with macroscopic shear is especially clear when activity is
localized in a shear band Lennartz-Sassinek  _et al._ (2014); Kwiatek  _et
al._ (2014); Davidsen  _et al._ (2017) (see deformation of D5kSc5 in Fig. 1)
but not so much when deformation is homogeneously distributed within the
sample (L5kSc5 in Fig. 1). Avalanche statistics shown in previous sections
appear to be independent of internal details once the system reaches SES,
whether system exhibits softening behaviour associated with localized failure
(e.g. in the case of D5kSc5), or monotonic hardening leading to diffuse
failure (e.g. in the case of L5kSc5). This suggests that the similarities in
avalanche statistics is linked to particle scale mechanisms, such as
frictionally arrested stick-slip events and force-chain buckling, rather than
the coarse-grained elasto-plastic macromechanics Tordesillas  _et al._ (2009).
On the other hand, the introduction of soft, instead of fixed, lateral walls
affects the energy partition of the system dynamics in general (Fig. 3) and
the avalanche energy budget in particular. For purely strain-controlled
loading programmes, such as in constant-volume tests, an average
proportionality is reported between $D$ and $\Delta\sigma$ Salerno  _et al._
(2012); Salerno and Robbins (2013). In our study presented here, the
equivalence between $D$ and $\Delta U$, and thus, the proportionality between
the stress drop and the energies, does not stand, because of the free
production of lateral work $W_{x}$. Specifically, Fig. 4 shows that such a
proportionality is maintained in $\Delta U$ and, therefore, lost for $D$ in
general. Recalling the different scaling exponents in $\Delta
U^{*}\propto\sigma_{x}^{1.85}$ and $W_{x}\propto\sigma_{x}$, the equivalence
$\Delta U\approx D$ is asymptotically recovered at high confining pressures,
given $D=\Delta U+W_{x}\approx\Delta U$.

###  IV.2 Effect of inertia and the pseudogap

Recent experiments in granular aggregates reported interevent time
distributions consistent with triggering Kumar  _et al._ (2020). Inertia is
known to affect statistics of both avalanche sizes and temporal correlations
and could potentially constitute a mechanism of triggering explaining such
observation. The relatively slow dissipation of the kinetic energy with
respect to the stress drops we observe (see Fig. 3) might trigger secondary
avalanches in the form of new spikes in $\dot{U}(t)$ and $E_{K}(t)$. Similar
phenomena have been reported in the signal profiles of acoustic emission
events during the uniaxial compression of concrete samples Vu and Weiss
(2020), We do observe, in rare occasions, a double peak in the kinetic energy
within the same avalanche profile (not shown). Nevertheless, given the
selected AQS driving protocol and low values of $K_{r}$, we can be certain
that such double peaks are absorbed in a single avalanche event, recovering a
practical independence between avalanches under the definition suggested in
Ref. Janićević  _et al._ (2016). No evidence of self-excitation is found in
either the Poisson-rejection test (Figs. 9.c) or the distribution of
interevent times shown in Fig. 9.b, which exhibits an exponential decay. On
the contrary, we report a deficit of short interevent times linked to
regularity in the macromechanical state of the system at the onset of an
avalanche (Fig. 10). This indicates that the process has a certain time-
predictability Shimazaki and Nakata (1980), well identified as a
characteristic minimum recharging strain as represented in Fig. 10.c.

This regularity of interevent times is predicted by finite size effects as
considered in some earthquake hazard assessment frameworks Salditch  _et al._
(2020). In small systems, all avalanches dissipate macroscopic fractions of
elastic energy, causing longer reloading times after large avalanches. This
effect disapears in the thermodynamic limit, where all subcritical and
critical avalanches are massless. Alternatively, a similar quasiperiodic
behavior is predicted in elastoplasticity by the presence of a so-called
’pseudogap’ in the density of residual-stresses, i.e. the local stress
differences to failure through different mechanisms, usually represented on a
coarse-grained scale Hébraud and Lequeux (1998); Picard  _et al._ (2004);
Karmakar  _et al._ (2010); Lin  _et al._ (2014); Lin and Wyart (2016); Ferrero
and Jagla (2019). Rather than a gap i.e., absence of small residual stresses
($x$), the pseudogap is characterized by a power-law decay in the density of
low values attributed to the non-monotonous evolution of residual-stresses
towards the stability limit which corresponds to an absorbing boundary in a
stochastic process representations Lin and Wyart (2016):

$N(x)\sim{x}^{\theta}\;\;\mathrm{with}\;\;\theta>0.$ (8)

The pseudogap introduces a characteristic delay in the interevent times
between consecutive avalanches, defined here as the strain $\epsilon_{y}$
needed to activate the lowest $x$ value in Eq. (8), and rendering the Weibull
renewal behavior Karmakar  _et al._ (2010):

$P(\Delta\epsilon)=\frac{1+\theta}{\langle\Delta\epsilon\rangle}\left({\frac{\Delta\epsilon}{\langle\Delta\epsilon\rangle}}\right)^{\theta}\exp\left({-\left({\dfrac{\Delta\epsilon}{\langle\Delta\epsilon\rangle}}\right)^{\theta+1}}\right).$
(9)

Such a behavior is shown in Fig. 9(a,b) as the black lines. Recent studies in
DEM Ruscher and Rottler (2020), FEM Korchinski  _et al._ (2021) and conceptual
Ferrero and Jagla (2019) models identified an additional system-size dependent
plateau overlapping with the pseudogap, which vanishes in the thermodynamic
limit. According to these results, quasi-periodicity would be more prominent
in large systems while interevent times in small systems would be better
described by a Poisson process. This system-size dependence, opposite to the
expectations from finite size effects, qualitatively agrees with our results
presented in Fig. 9 and is numerically verified in the supplementary material
(Fig. S3, S6). However, the analysis presented in Fig. 9 is insufficient to
determine the system-size dependency of the plateau or the exponent $\theta$
defining the pseudogap following the analysis presented in Lin  _et al._
(2014). Instead one can use hyperscaling relations to address the latter as we
discuss in section IV.5. Yet, Fig. 9 shows that the distribution of short
interevent times are close to the expectations of the Weibull renewal model
with an exponent between 0.5 and 1, with a generally higher exponent,
suggesting a more prominent pseudogap in the large samples (Fig. 9.a). At the
region of long interevent times, the distributions are better reproduced by
the Poisson hypothesis instead (Fig. 9.b). This might indicate a change in the
distribution of $x$, marking the end of the power-law regime described by Eq.
(8).

Finally, inertial effects were found to affect the critical exponents
determinining bulk avalanche statistics in DEM simulations of LJ particles
Salerno  _et al._ (2012); Salerno and Robbins (2013); Karimi and Barrat
(2016); Karimi  _et al._ (2017). This feature could as well be related to the
observation that the states triggering avalanches, herein considered to be
constrained to a SES manifold, are more regular than the arrested states (see
section III.4). The kinetic energy introduced by ‘fast’ external driving is
known to affect mechanical equilibrium Anthony and Marone (2005); Lemaître and
Caroli (2009); LU  _et al._ (2007); Nicot  _et al._ (2012). We propose that
the SES represents a collection of metastable configurational limits to the
mechanical equilibrium that are marginally reachable under quasistatic
conditions Bonfanti  _et al._ (2019). When an avalanche is initiated at SES,
the excess kinetic energy prevents the relaxation of the system at a fragile
state, e.g. the same configurational state at the onset of the avalanche. The
configurational state arresting the avalanche, which is reached at $t_{0}+T$,
needs to be stable while the system still retains a relatively high amount of
kinetic energy (see Fig. 3). Since the excess of kinetic energy is dampened
below $K_{r}$ before resuming the driving, we believe that the particle
inertia during the avalanche has to affect the distribution of residual
stresses, and the properties of the pseudogap. This idea of stability
overshooting was conceptually captured in a deppining model with transient
effects predicting a system-size dependence on the time-predictability of
transitions Maimon and Schwarz (2004) qualitatively similar to the one
observed in the pseudogap, verified at the stiff particles limit in our
simulations. A similar effect is reproduced by the addition of dynamic
weakening in some EPM, and the SMFT in particular Dahmen  _et al._ (2009).
However, the effect on avalanche sizes is far from being trivial. Numerical
studies with LJ particles found non-monotonous relations between critical
exponents and inertia, and regimes with characteristic and time-predictable
avalanches Salerno and Robbins (2013), being the latter predicted in SMFT as
well Dahmen  _et al._ (2009).

###  IV.3 Stiff granular particles

Apart from particle inertia, most coarse-grained EPM also overlook the non-
linearity of the internal stress fields Liu  _et al._ (1995); Mueth  _et al._
(1998); Howell  _et al._ (1999). The hypothesis of the uniformity in the
stress propagation enabling the coarse-graining in elastoplastic models is
known to fail in the presence of stiff particles Nicolas  _et al._ (2018). As
mentioned, stress distribution in granular media is governed by the structural
topology of contact points, generating complex force-chain networks, highly
susceptible to small variations of stress. A cross-over between fragile and
amorphous materials has long been suggested to exist as a function of softness
or deformability Cates  _et al._ (1998), even in the absence of attractive
particle potential Liu and Nagel (1998); Lerner and Procaccia (2009). The
variation in the power-law exponents in Fig. 8 is in agreement with such a
transition (see section IV.5 for a more detailed discussion). Recent studies
have suggested the introduction of internal structural variables such as
compativity Edwards and Oakeshott (1989); Frenkel  _et al._ (2008); Blumenfeld
_et al._ (2012), angoricity Henkes  _et al._ (2007); Blumenfeld and Edwards
(2009); Bi  _et al._ (2015); Baule  _et al._ (2018), or the more recent
concept of keramicity Bililign  _et al._ (2019); Ball (2019); Baule  _et al._
(2018) complementing the stress-strain formulation traditionally considered in
EPM. Based on our results, we propose the SES concept as an alternative
micromechanical representation of the state of granular materials determined
by the structure of internal contact network rather than the macromechanics
Pouragha and Wan (2018a).

###  IV.4 Interpretation of avalanche sizes in terms of elastoplasticity

The EPM formulation separates plastic avalanche events from the elastic
deformation caused by the AQS driving. Although the separation of temporal
scales is not as straightforward in DEM simulations, one can still find a
correspondence between the variables defining avalanches in the different
modeling approaches. In this sense, we are interested in the conceptual
picture provided by the SMFT, which is arguably the simplest EPM framework to
characterize avalanche statistics, and has shown good agreement with the
statistics of avalanches in granular assemblies sharing in particular the same
universal critical exponents Dahmen  _et al._ (2011); Denisov  _et al._
(2016); Geller  _et al._ (2015). In the SMFT, the avalanches are measured in
terms of a global signal ($v(t)$ in short intervals) which accounts for the
collective velocity of displacement of the constituting elements (a.k.a.
cells) of the model, and is proportional to the rate of external stress drop
under strain driving. Henceforth, the size $S:=\int_{t_{0}}^{t_{0}+T}v(t)dt$
of an avalanche of duration $T$ starting at time $t_{0}$ is a measure
proportional to stress drop. Similarly, $E:=\int_{t_{0}}^{t_{0}+T}v^{2}(t)dt$
is proportional to the total kinetic energy dissipated, and
$E_{m}:=v_{\max}^{2}(t)dt$ is proportional to the maximum peak of the
instantaneous kinetic energy within the avalanche. In order to compare DEM
results with the expectations from elastoplasticity, we first need to
establish reliable links between the avalanche sizes defined within the
conceptual framework proposed in the SMFT and the measures defining avalanche
sizes in DEM.

Considering EPM formulations, it is natural to assume an approximate
proportionality between elastic energy and small deformations. Given its mean-
field nature, the SMFT assumes an exact proportionlity ($\dot{\sigma}\propto
v(t)$) during avalanches, and predicts an approximate proportionality with the
elastic energy release. This is consistent with the observed proportionality
between internal elastic energy $\Delta U$ and $\Delta\sigma_{y}$ illustrated
in Fig. 4. We adopt the strain drop as the avalanche size
$S:=\Delta\sigma_{y}$ following the convention in Refs. Salerno  _et al._
(2012); Salerno and Robbins (2013). Moreover, to remain consistent with SMFT,
the avalanche duration $T$ has been defined in this work to be the time
interval of the potential energy release, which is virtually indistinguishable
from the time interval of the stress drop.

Although $K:=E_{K}^{\max}-K_{D}$ represents the peak in our simulations and,
hence, it is reminiscent of $E_{m}$, the dissipation process is significantly
slower compared to the rise of kinetic energy. In contrast, SMFT considers a
strict overdamped limit where all inertial effects are instantly dissipated.
Therefore, we argue that the pick of kinetic energy $K$ reached during the
initial avalanche growth corresponds to the collective cumulative energy $E$
in terms of SMFT. In Fig. 3.a we observe a separation of temporal scales
between the kinetic energy rise time, roughly matching the temporal interval
($t_{0}$, $t_{0}+T$), and the characteristic dissipation time (see the grey
areas in Fig. 3.a). Almost no energy has been dissipated at $T$. We therefore
adopt $K\propto E$ in terms of the variables in SMFT. This assumption can be
verified by comparing the values of $K$ with an auxiliary avalanche variable
$\Delta(U^{2}):=\int_{t_{0}}^{t_{0}+T}dt(dU/dt)^{2}$ (see Fig. S7 in
supplementary material). If we assume a temporal proportionality between
$\dot{U}(t)$ and $v(t):=dS/dt$, for consistency $(dU/dt)^{2}(t)\propto
v(t)^{2}$. Therefore the proportionality $K\propto E$ implies
$K\propto\Delta(U^{2})$, as verified in Fig. S7. However, for the sake of
completeness, in the next section we compare our results considering both
$K\propto E$ (herein referred as model SMFT(1)) and $K\propto E_{m}$ (SMFT(2))
to validate our hypothesis in light of the exponent values.

###  IV.5 Non-universal critical exponents

In the absence of self-excitation the spatial extension of avalanches serves
as a proxy of local spatial correlations and, therefore, one can assess system
susceptibilities from avalanche statistics. Since no divergence of activity
rate is observed in the catalogs, criticality has to be linked to a divergence
of characteristic scales in avalanche distributions. Hence the power-law
distributions $P(\Delta U)$ and $P(K)$ reported in section III.2 are a
trademark of criticality. Close to criticality such distributions obey a
scaling form:

$\begin{array}[]{c}P(\Delta U)\>d\Delta U=\Delta U^{-\kappa}\>\Phi_{\Delta
U}(\Delta U/\Delta U^{*})\>d\Delta U,\\\
P(K)\>dK=K^{-\varepsilon}\>\Phi_{K}(K/K^{*})\>dK,\end{array}$ (10)

where $\kappa$ and $\varepsilon$ are the so-called critical exponents (Table
2), $\Phi_{\Delta U},\Phi_{K}$ are universal scaling functions (often in the
form of an exponential cut-off) and $\Delta U^{*},K^{*}$ are characteristic
avalanche sizes, related to the correlation length. When approaching
criticality, the correlation length diverges such that for any finite system
it eventually takes on the system size value. Therefore, $\Delta U^{*}\sim
L^{d_{f}}$, where $L$ is the lateral size of the system. Considering that the
system stays close to criticality, at the stiff particle limit
($\sigma_{x}=10^{5}$) our simulation results are consistent with $d_{f}\approx
1$ (see section III.2). Values $d_{f}\sim 1$ have been commonly found in
previous studies Maloney and Lemaître (2004) and DEM simulations Salerno and
Robbins (2013); Shang  _et al._ (2020). Other FEM simulations Karimi  _et al._
(2017) reported a dependence of exponents $d_{f}$ and $\theta$ on a damping
coefficient $\Gamma$ that could, potentially, include the contribution of
confining pressure to the avalanche propagation. Here for the sake of the
discussion we will consider $d_{f}$ to be invariant to $\sigma_{x}$.

In SMFT, the relations between the average sizes, energies, and durations take
the following power-law form (see Dahmen (2017)):

$\begin{array}[]{lcl}\langle T|S\rangle\sim S^{\varsigma\nu
z}&\;\mathrm{where}&\varsigma\nu z=1/2\\\ \langle E|S\rangle\sim
S^{2-\varsigma\nu z}&\;\mathrm{where}&2-\varsigma\nu z=3/2\\\ \langle
E_{m}|S\rangle\sim S^{2\varsigma\rho}&\;\mathrm{where}&2\varsigma\rho=1\;.\\\
\end{array}$ (11)

Considering the proportionality between differential deformation and avalanche
sizes, one can prove that the monotonic driving brings the system towards the
aforementioned stationary SOC flow regime around $\sigma_{c}$ with scale-free
statistics, leading to Eq. 1. Considering the exponents above, and the
proportionalty between $\Delta U$ and $S$, criticality imposes power-law
distributions with the following exponent relations (see Dahmen (2017)):

$\begin{array}[]{lcl}P(S)\sim S^{-\kappa}&\;\mathrm{where}&\kappa=3/2\\\
P(E)\sim E^{-1-\frac{\kappa-1}{2-\varsigma\nu
z}}&\;\mathrm{where}&1+\frac{\kappa-1}{2-\varsigma\nu z}=4/3\\\ P(E_{m})\sim
E_{m}^{-\frac{1+\mu}{2}}&\;\mathrm{where}&\frac{1+\mu}{2}=3/2\;.\\\
\end{array}$ (12)

In Table 2 we compare our results with the predictions of SMFT in both
interpretations of $K$ (models SMFT(1) and SMFT(2)) and the typical range of
exponents reported in numerical studies considering similar relation between
observables Budrikis  _et al._ (2017); Liu  _et al._ (2016); Nicolas  _et al._
(2018). We notice that studies investigating the effect of dimensionality
reported no significant differences between 2D and 3D systems Budrikis  _et
al._ (2017); Salerno and Robbins (2013). The most remarkable observation is
that the exponents in our simulations are non-universal. When the confining
pressure is low, or contacts are relatively stiff ($\sigma_{x}=10^{5}$), the
granular compact assemblies in all three exponents $\kappa$, $\epsilon$,
$\gamma$ is consistent with SMFT(1), i.e. considering $K=E$, although
exponents $\kappa$ and $\gamma$ are slightly higher than expected by SMFT.
Notice that this agreement holds for both loose and dense assemblies, even
when the macromechanical evolutionary path is different from SMFT in the
latter. The correspondence with the SMFT predictions is progressively lost for
higher confining pressures. In the extreme case of $\sigma_{x}=10^{7}$, the
exponent $\gamma=1.48(7)$ is still close to the mean field prediction, and in
agreement with recent results in 2D EPM simulations (e.g. $\gamma\approx 1.45$
in Budrikis  _et al._ (2017)). However, $\kappa$ and $\epsilon$ decrease
towards values close to 1.0. Note that $\gamma$ plays a minor role in the
relation between $\kappa$ and $\varepsilon$ when both exponents are close to
$1$, see Eq. (6). The exponents $\kappa$ and $\varepsilon$ at high confining
pressures, which depart from SMFT predictions, are closer to the predictions
from EPM within 2D lattices, where lower exponents are reported: $\kappa\sim
1.28$ and $\varepsilon\sim 1.20$ Budrikis  _et al._ (2017); Liu  _et al._
(2016), $\kappa\approx 1.1-1.2$ Jagla (2015) or even down to $\kappa=0.98(1)$
Shang  _et al._ (2020). The ranges found in $\kappa$ are also consistent with
previous reports comparing experiments and DEM simulations Barés  _et al._
(2017), and simulations with variable heterogeneity Ma  _et al._ (2020). It
has been argued that this change of behavior can be interpreted as a
transition between different universality classes Budrikis  _et al._ (2017).
In that case, our results might indicate a transition between conceptual
limits.

To sum up, the statistics of avalanches in non-cohesive assemblies of stiff
particles are relatively close to the expectations of the SMFT, specially for
the most robust exponent $\varepsilon$, whilst the avalanches for soft
particles are close to the predictions of elastoplasticity with short range
interactions or geometric propagators. This might signal a dependence to
softness of the propagation of interactions, as already predicted Cates  _et
al._ (1998) and observed in impact experiments Clark  _et al._ (2015a). A
plausible explanation is that the force chains formations are more prominent
at the stiff-particle limit, and as such, the effective interactions are long-
ranged Clark  _et al._ (2015b) and avalanches can be represented in terms of
isotropic homogeneous and independent branching processes. On the contrary, at
the soft-particle limit, the deformability of particles accommodates the
stress variation in finite structures well modelled by Eshelby-like
propagators in the coarse-grained mesoscale. Therefore, we argue that the non-
equilibrium statistical physics of fragile matter do not significantly differ
from those of amorphous materials with attractive potentials as long as
correlations, or influence of force-chain structures, are damped at short
length scales. Similar parameter-dependent exponents have been found in other
simulations and linked to changes in the effective dimension of avalanche
location Amitrano (2003); Niemann  _et al._ (2012). From standard
renormalization theory it is known that, in general, effective system
dimensions and the range of interactions condition the universality class,
i.e. the critical exponents. The behavior at intermediate softness could
therefore represent a cross-over phenomenon, giving rise to an “effective” set
of exponents that would depend on the system size.

Finally, following the hyperscaling relation between avalanche-size
distribution and the pseudogap Lin and Wyart (2016):

$\kappa=2-\frac{\theta}{\theta+1}\frac{d}{d_{f}}.$ (13)

Assuming $d_{f}$ as constant, the exponent $\theta$, determining the
distribution of residual stresses, Eq. (8), and interevent times, transitions
from $\theta\approx 0.5$ for stiff particles to $\theta\approx 1$ for soft
particles if we use the fitted values of $\kappa$ from Table 2. These two
extreme cases (represented in Fig. 9) are in agreement, respectively, with the
expected $\theta=0.5$ found in mean field models Lin and Wyart (2016) and
$\theta\approx 0.9-1$ found in DEM simulations considering LJ potentials Liu
_et al._ (2016); Jagla (2015); Shang  _et al._ (2020). Overall, we observe a
general correspondence between the exponents fitted at low confining pressure
with the predictions of SMFT, and we recover the usual ranges of values found
in EPM and other simulations of amorphous materials at higher confining
pressures.

###  IV.6 Avalanche statistics on the SES are stationary and scale free

The interpretation of the SES manifold as a stability limit solves the problem
of the non-monotonous evolution from a macroscopic point of view, once
structural configurational aspects are taken into account. Beyond that, and as
far as we can tell, the SES defines the ’state’ of the assembly where
avalanche statistics is non-degenerated. Upon reaching the SES at different
stress values, and once the variations in the avalanche activity rate are
accounted for, no difference in the statistics of avalanches is observed (Fig.
5). The distributions of avalanche sizes are virtually identical,
independently of the $\sigma_{y}$ values, the evolutionary path, and the
initial configuration. Focusing on the result at the stiff-particle limit, and
despite the complexity of the internal processes involved, the macromechanical
evolution of loose assemblies, such as L5kSc5 in Fig. 1, is in agreement with
the predictions of the SMFT (Eq. (1)). In particular, the SMFT predictions
extend way beyond the proximity-limit of $\sigma_{c}$, and apply to the entire
observation interval. In contrast, initially dense assemblies, such as D5kSc5
in Fig. 1, display the typical peak in stress and post-peak strain softening
which is not captured by Eq. (1) nor the scaling relations in SMFT. However,
once SES is reached, the distributions in Fig. 5 are already compatible with
the scale-free behavior with SMFT universality, besides the aforementioned
finite size effects, way before the steady state flow regime. Therefore, we
argue that the SES is not only an attractor manifold to the internal state of
the system, but might as well correspond to the self-organized critical (SOC)
state Bak and Chen (1991); Chen  _et al._ (1991); Zapperi  _et al._ (1995);
Watkins  _et al._ (2016) from the perspective of out-of-equilibrium
statistical physics. This is also demonstrated in Ref. Pouragha and Wan
(2018a) where a preliminary measure for a distance to criticality is suggested
based on the concept of SES. Further studies are nevertheless necessary to
determine whether the SES is indeed an exact representation of avalanche
criticality or yet another approximation. This question can be addressed in
principle through a finite size scaling analysis of avalanche distributions in
terms of the parameters representing the distance to SES. Yet, it is a
comptuationally cumbersome task if one wants to obtain statistically reliable
conclusions Budrikis  _et al._ (2017), which remains a challenge for the
future.

##  V Conclusions:

We presented a statistical analysis of kinetic avalanches recorded during DEM
simulations of granular assemblies with elastic potentials under AQS
conditions. It showed that criticality, in the form of power-law statistics,
is degenerated in the space of the commonly used macroscopic state variables
such as stress and strains. Instead, internal measurements, such as the fabric
anisotropy, the rigidity ratio, and average coordination number, in which the
manifold herein referred to as stable evolutionary state (SES) is defined,
appear to provide a more appropriate representation of criticality, and
closeness to criticality, in terms of avalanche statistics and
susceptibilities. In particular, we have presented evidence that points
towards the existence of a link between (i) the SES, identified as an
attractor of the evolution of granular systems in the space of internal
measurements Pouragha and Wan (2018a) and (ii) criticality characterized by a
divergence of susceptibilities and scale-free avalanches. Therefore, we argue
that the SES is itself a self-organized critical (SOC) state, representing a
basin of attraction to the state of the system with scale-free independent
avalanches before reaching the stationary flow regime.

We interpret the observed time predictability as evidence that the SES
corresponds to a stability limit under AQS conditions. In terms of the
residual-stress formulation, we conjecture that the SES, determined by the
potential energy landscape, dominates the dynamics during driving as
considered in EPM, although the high susceptibility of the system might
introduce non-affine stress rearrangements through the force-chain network.
Conversely, based on the higher correlation between states at the onset of
avalanches than between states arresting avalanches, we believe that the
avalanche propagation is strongly affected by the transient kinetic energy in
underdamped inertial systems. Therefore, the stability of the assembly during
the compression process is determined by the combined effect of equilibrium
constraints governing the quasi-static stability and the transient addition of
kinetic energy during avalanche propagation. Testing this hypothesis further
remains a challenge for the future, potentially involving simulations with
alternative control protocols, the inclusion of thermal effects and/or the
tuning of frictional and damping parameters. On the other hand, the reported
avalanche interevent times are compatible with the power-law pseudogap
hypothesis, rendering a Weibull renewal process for large system sizes, while
the Poisson hypothesis is still valid for small systems. This is in agreement
with the possible existence of a system-size dependent plateau in the
distribution of residual stresses recently discussed in the literature Ruscher
and Rottler (2020); Korchinski  _et al._ (2021); Ferrero and Jagla (2019).

Despite the particular aspects of our simulations, the reported avalanche
statistics do not differ significantly from the expectations of EPM and can,
in fact, bridge the gap between the depinning universality class associated
with the slip mean field theory (SMFT) and the extended non-universal
exponents in yielding Lin  _et al._ (2014) that could represent a cross-over
towards a new universality class of elastoplasticity and, therefore, be
invariant with respect to the dimensionality of the system Budrikis  _et al._
(2017). We found that the external stress drop ($\Delta\sigma$) is
proportional to the internal energy drop $\Delta U$ and a good proxy for the
internal displacement, or $S$ in terms of the slip mean field theory (SMFT)
Dahmen  _et al._ (2011); Dahmen (2017). The peak of kinetic energy ($K$)
serves as a good proxy for the avalanche energy $E$. The observed exponents
characterizing the critical state are non-universal and display a transition
from (i) the expected values for SMFT, for stiff particles or low confining
pressures in our simulations, to (ii) the lower exponents observed in 2d
simulations of amorphous materials Budrikis  _et al._ (2017); Liu  _et al._
(2016); Jagla (2015); Shang  _et al._ (2020), for soft particles or high
confining pressures in our simulations. This could be linked to the known
change in spatial distribution and anisotropy of correlations in the two
limits, which depends on the structure and robustness of the internal force-
chain network connecting the contact points Cates  _et al._ (1998); Clark  _et
al._ (2015a). In particular, for stiff particles, small displacements in a
single contact can propagate through the force-chain and span the whole
system, resulting in a homogeneous and isotropic long-ranged response at the
coarse-grained scale, similar to the interactions in the SMFT. In contrast,
soft-particles accommodate variations of stress at a short-range scale,
reminiscent of the typical Eshelby inclusion modeled in 2d EPM Albaret  _et
al._ (2016). On the other hand, intertial effects are found to depend on
stiffness and are known to affect avalanche statistics as well Salerno  _et
al._ (2012); Salerno and Robbins (2013); Karimi and Barrat (2016); Karimi  _et
al._ (2017). Finally, the relation between stiffness and the universality
class is further reinforced by the indirect estimation of the exponent
$\theta$ defining the distribution of residual stresses, which satisfies the
expectations of mean-field models Lin and Wyart (2016) for stiff particles and
2D EPM for soft particles Liu  _et al._ (2016); Jagla (2015); Shang  _et al._
(2020).

To summarize, based on our results we propose a static stability limit (SES)
as a candidate for a SOC state reached before the traditional flow regime, and
a link between the collection of critical exponents defining the critical
state and the stiffness of the particles. These two novel conjectures should
serve as a starting point to extend our understanding of the discussed non-
universality in the deformation of amorphous materials, and the role of the
internal structure and its dynamics during the deformation. Indeed, we
strongly believe that this might be a fundamental aspect to understand the fit
of fragile matter within elastoplastic theory.

###### Acknowledgements.

Daniel Korchinski for fruitful discussions. J.B. thanks financial support from
the Banco de Santander - María de Maetzu fellowship program, and the AXA
Research Fund through the project RheMechFail. This work was financially
supported by the Natural Sciences and Engineering Research Council of Canada
(NSERC).

## References

  * Sethna  _et al._ (2001) J. P. Sethna, K. A. Dahmen, and C. R. Myers, Nature 410, 242 (2001).
  * Salje and Dahmen (2014) E. K. Salje and K. A. Dahmen, Annual Review of Condensed Matter Physics 5, 233 (2014),
  * Nicolas  _et al._ (2018) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Rev. Mod. Phys. 90, 045006 (2018).
  * Miguel  _et al._ (2001) M. C. Miguel, A. Vespignani, S. Zapperi, J. Weiss, and J.-R. Grasso, Nature 410, 667 (2001).
  * Zaiser (2006) M. Zaiser, Advances in Physics 55, 185 (2006), .
  * Friedman  _et al._ (2012) N. Friedman, A. T. Jennings, G. Tsekenis, J.-Y. Kim, M. Tao, J. T. Uhl, J. R. Greer, and K. A. Dahmen, Phys. Rev. Lett. 109, 095507 (2012).
  * Wang  _et al._ (2009) G. Wang, K. Chan, L. Xia, P. Yu, J. Shen, and W. Wang, Acta Materialia 57, 6146 (2009).
  * Antonaglia  _et al._ (2014) J. Antonaglia, W. J. Wright, X. Gu, R. R. Byer, T. C. Hufnagel, M. LeBlanc, J. T. Uhl, and K. A. Dahmen, Phys. Rev. Lett. 112, 155501 (2014).
  * Schall  _et al._ (2007) P. Schall, D. A. Weitz, and F. Spaepen, Science 318, 1895 (2007), .
  * Hu  _et al._ (2018) Y. Hu, L. Shu, Q. Yang, W. Guo, P. K. Liaw, K. A. Dahmen, and J.-M. Zuo, Communications Physics 1, 61 (2018).
  * Rizzardi  _et al._ (2018) Q. Rizzardi, G. Sparks, and R. Maaß, JOM 70, 1088 (2018).
  * Dahmen  _et al._ (2011) K. a. Dahmen, Y. Ben-Zion, and J. T. Uhl, Nature Physics 7, 554 (2011).
  * Denisov  _et al._ (2016) D. V. Denisov, K. A. Lörincz, J. T. Uhl, K. A. Dahmen, and P. Schall, Nature Communications 7, 10641 (2016).
  * Berthier  _et al._ (2019a) E. Berthier, J. E. Kollmer, S. E. Henkes, K. Liu, J. M. Schwarz, and K. E. Daniels, Phys. Rev. Materials 3, 075602 (2019a).
  * Maloney and Lemaître (2004) C. Maloney and A. Lemaître, Phys. Rev. Lett. 93, 016001 (2004).
  * Salerno  _et al._ (2012) K. M. Salerno, C. E. Maloney, and M. O. Robbins, Phys. Rev. Lett. 109, 105703 (2012).
  * Salerno and Robbins (2013) K. M. Salerno and M. O. Robbins, Phys. Rev. E 88, 062206 (2013).
  * Parisi  _et al._ (2017) G. Parisi, I. Procaccia, C. Rainone, and M. Singh, Proceedings of the National Academy of Sciences 114, 5577 (2017),
  * Zhang  _et al._ (2017) D. Zhang, K. A. Dahmen, and M. Ostoja-Starzewski, Phys. Rev. E 95, 032902 (2017).
  * Shang  _et al._ (2020) B. Shang, P. Guan, and J.-L. Barrat, Proceedings of the National Academy of Sciences 117, 86 (2020),
  * Lerner and Procaccia (2009) E. Lerner and I. Procaccia, Phys. Rev. E 79, 066109 (2009).
  * Karmakar  _et al._ (2010) S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E 82, 055103(R) (2010).
  * Gimbert  _et al._ (2013) F. Gimbert, D. Amitrano, and J. Weiss, EPL (Europhysics Letters) 104, 46001 (2013).
  * Pouragha and Wan (2016) M. Pouragha and R. Wan, Granular Matter 18, 38 (2016).
  * Karimi and Barrat (2018) K. Karimi and J.-L. Barrat, Scientific Reports 8, 4021 (2018).
  * Cates  _et al._ (1998) M. E. Cates, J. P. Wittmer, J.-P. Bouchaud, and P. Claudin, Phys. Rev. Lett. 81, 1841 (1998).
  * Liu and Nagel (1998) A. J. Liu and S. R. Nagel, Nature 396, 21 (1998).
  * Bi  _et al._ (2011) D. Bi, J. Zhang, B. Chakraborty, and R. P. Behringer, Nature 480, 355 (2011).
  * Bak and Chen (1991) P. Bak and K. Chen, Scientific American 264, 46 (1991).
  * Chen  _et al._ (1991) K. Chen, P. Bak, and S. P. Obukhov, Phys. Rev. A 43, 625 (1991).
  * Zapperi  _et al._ (1995) S. Zapperi, K. B. Lauritsen, and H. E. Stanley, Phys. Rev. Lett. 75, 4071 (1995).
  * Watkins  _et al._ (2016) N. W. Watkins, G. Pruessner, S. C. Chapman, N. B. Crosby, and H. J. Jensen, Space Science Reviews 198, 3 (2016).
  * Hébraud and Lequeux (1998) P. Hébraud and F. Lequeux, Phys. Rev. Lett. 81, 2934 (1998).
  * Ben-Zion and Rice (1993) Y. Ben-Zion and J. R. Rice, Journal of Geophysical Research: Solid Earth 98, 14109 (1993),
  * Dahmen  _et al._ (2009) K. A. Dahmen, Y. Ben-Zion, and J. T. Uhl, Phys. Rev. Lett. 102, 175501 (2009).
  * Dahmen (2017) K. A. Dahmen, Mean field theory of slip statistics, in  _Avalanches in Functional Materials and Geophysics_, edited by E. K. Salje, A. Saxena, and A. Planes (Springer International Publishing, Cham, 2017) pp. 19–30.
  * Dahmen  _et al._ (1998) K. Dahmen, D. Ertaş, and Y. Ben-Zion, Phys. Rev. E 58, 1494 (1998).
  * Mehta  _et al._ (2006) A. P. Mehta, K. A. Dahmen, and Y. Ben-Zion, Phys. Rev. E 73, 056104 (2006).
  * Geller  _et al._ (2015) D. A. Geller, R. E. Ecke, K. A. Dahmen, and S. Backhaus, Phys. Rev. E 92, 060201(R) (2015).
  * Corwin  _et al._ (2005) E. I. Corwin, H. M. Jaeger, and S. R. Nagel, Nature 435, 1075 (2005).
  * Lin  _et al._ (2014) J. Lin, A. Saade, E. Lerner, A. Rosso, and M. Wyart, EPL (Europhysics Letters) 105, 26003 (2014).
  * Budrikis  _et al._ (2017) Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser, and S. Zapperi, Nature Communications 8, 15928 (2017).
  * Korchinski  _et al._ (2021) D. Korchinski, C. Ruscher, and J. Rottler, arXiv preprint arXiv:2101.10688 (2021).
  * Baret  _et al._ (2002) J.-C. Baret, D. Vandembroucq, and S. Roux, Phys. Rev. Lett. 89, 195506 (2002).
  * Picard  _et al._ (2005) G. Picard, A. Ajdari, and F. Lequeux, L. Bocquet, Phys. Rev. E 71, 010501(R) (2005).
  * Agoritsas  _et al._ (2015) E. Agoritsas, E. Bertin, K. Martens, and J.-L. Barrat, The European Physical Journal E 38, 71 (2015).
  * Lin and Wyart (2016) J. Lin and M. Wyart, Phys. Rev. X 6, 011005 (2016).
  * Ferrero and Jagla (2019) E. E. Ferrero and E. A. Jagla, Soft Matter 15, 9041 (2019).
  * Rodney  _et al._ (2011) D. Rodney, A. Tanguy, and D. Vandembroucq, Modelling and Simulation in Materials Science and Engineering 19, 083001 (2011).
  * Miller  _et al._ (1996) B. Miller, C. O’Hern, and R. P. Behringer, Phys. Rev. Lett. 77, 3110 (1996).
  * Mueth  _et al._ (1998) D. M. Mueth, H. M. Jaeger, and S. R. Nagel, Physical Review E 57, 3164 (1998).
  * Liu  _et al._ (1995) C. h. Liu, S. R. Nagel, D. A. Schecter, S. N. Coppersmith, S. Majumdar, O. Narayan, and T. A. Witten, Science 269, 513 (1995),
  * Howell  _et al._ (1999) D. Howell, R. P. Behringer, and C. Veje, Phys. Rev. Lett. 82, 5241 (1999).
  * Andò  _et al._ (2012) E. Andò, S. Hall, G. Viggiani, J. Desrues, and P. Bésuelle, Géotechnique Letters 2, 107 (2012).
  * Gendelman  _et al._ (2016) O. Gendelman, Y. G. Pollack, I. Procaccia, S. Sengupta, and J. Zylberg, Physical Review Letters 116, 078001 (2016),
  * Berthier  _et al._ (2019b) E. Berthier, M. A. Porter, and K. E. Daniels, Proceedings of the National Academy of Sciences 116, 16742 (2019b),
  * Tordesillas (2007) A. Tordesillas, Philosophical Magazine 87, 4987 (2007),
  * Guo and Zhao (2013) N. Guo and J. Zhao, Computers and Geotechnics 47, 1 (2013).
  * Kuhn (1999) M. R. Kuhn, Mechanics of materials 31, 407 (1999).
  * Pouragha and Wan (2017a) M. Pouragha and R. Wan, International Journal of Solids and Structures 110, 94 (2017a).
  * Edwards and Oakeshott (1989) S. Edwards and R. Oakeshott, Physica A: Statistical Mechanics and its Applications 157, 1080 (1989).
  * Blumenfeld and Edwards (2009) R. Blumenfeld and S. F. Edwards, The Journal of Physical Chemistry B 113, 3981 (2009),
  * Blumenfeld  _et al._ (2012) R. Blumenfeld, J. F. Jordan, and S. F. Edwards, Phys. Rev. Lett. 109, 238001 (2012).
  * Bi  _et al._ (2015) D. Bi, S. Henkes, K. E. Daniels, and B. Chakraborty, Annual Review of Condensed Matter Physics 6, 63 (2015),
  * Baule  _et al._ (2018) A. Baule, F. Morone, H. J. Herrmann, and H. A. Makse, Rev. Mod. Phys. 90, 015006 (2018).
  * Bililign  _et al._ (2019) E. S. Bililign, J. E. Kollmer, and K. E. Daniels, Phys. Rev. Lett. 122, 038001 (2019).
  * Ball (2019) P. Ball, Nature Materials 18, 196 (2019).
  * Masson and Martinez (2001) S. Masson and J. Martinez, Journal of engineering mechanics 127, 1007 (2001).
  * Desrues and Viggiani (2004) J. Desrues and G. Viggiani, International Journal for Numerical and Analytical Methods in Geomechanics 28, 279 (2004).
  * Desrues and Hammad (1989) J. Desrues and W. Hammad, in _Proc. of the Int. Workshop on Numerical Methods for Localization and Bifurcation of Granular Bodies_ (1989) pp. 57–67.
  * Pouragha and Wan (2017b) M. Pouragha and R. Wan, Journal of Engineering Mechanics 143, C4016002 (2017b).
  * Inc (2008) I. C. G. Inc, Minneapolis: Itasca Consulting Group Inc. (2008).
  * Da Cruz  _et al._ (2005) F. da Cruz, S. Emam, M. Prochnow, J.-N. Roux, and F. Chevoir, Physical Review E 72, 021309 (2005).
  * Radjai and Roux (2002) F. Radjai and S. Roux, Phys. Rev. Lett. 89, 064302 (2002).
  * Bonfanti  _et al._ (2019) S. Bonfanti, R. Guerra, C. Mondal, I. Procaccia, and S. Zapperi, Phys. Rev. E 100, 060602(R) (2019).
  * Nicot  _et al._ (2017) F. Nicot, H. Xiong, A. Wautier, J. Lerbet, and F. Darve, Granular Matter 19, 18 (2017).
  * Pouragha and Wan (2018a) M. Pouragha and R. Wan, Mechanics of Materials 126, 57 (2018a).
  * Lemaître and Caroli (2009) A. Lemaître and C. Caroli, Phys. Rev. Lett. 103, 065501 (2009).
  * LU  _et al._ (2007) K. LU, E. E. Brodsky, and H. P. Kavehpour, Journal of Fluid Mechanics 587, 347–372 (2007).
  * Nicot  _et al._ (2012) F. Nicot, N. Hadda, F. Bourrier, L. Sibille, R. Wan, and F. Darve, International Journal of solids and Structures 49, 1252 (2012).
  * Maimon and Schwarz (2004) R. Maimon and J. M. Schwarz, Phys. Rev. Lett. 92, 255502 (2004).
  * Chang  _et al._ (1995) C. S. Chang, S. J. Chao, and Y. Chang, International journal of solids and structures 32, 1989 (1995).
  * Pouragha and Wan (2018b) M. Pouragha and R. Wan, International journal of solids and structures 138, 97 (2018b).
  * Privman (1990) V. Privman, _Finite size scaling and numerical simulation of statistical systems_ (World Scientific Singapore, 1990).
  * Liu  _et al._ (2016) C. Liu, E. E. Ferrero, F. Puosi, J.-L. Barrat, and K. Martens, Phys. Rev. Lett. 116, 065501 (2016).
  * Zapperi  _et al._ (1997) S. Zapperi, A. Vespignani, and H. E. Stanley, Nature 388, 658 (1997).
  * Baró and Vives (2012) J. Baró and E. Vives, Phys. Rev. E 85, 066121 (2012).
  * Christensen  _et al._ (2008) K. Christensen, N. Farid, G. Pruessner, and M. Stapleton, The European Physical Journal B 62, 331 (2008).
  * Karimi and Barrat (2016) K. Karimi and J.-L. Barrat, Phys. Rev. E 93, 022904 (2016).
  * Amitrano (2003) D. Amitrano, Journal of Geophysical Research: Solid Earth 108 (2003), 10.1029/2001JB000680,
  * Lord-May  _et al._ (2020) C. Lord-May, J. Baró, D. W. Eaton, and J. Davidsen, Phys. Rev. Research 2, 043324 (2020).
  * Didier Sornette (1994) Didier Sornette, J. Phys. I France 4, 209 (1994).
  * Vazquez  _et al._ (2011) F. Vazquez, J. A. Bonachela, C. López, and M. A. Muñoz, Phys. Rev. Lett. 106, 235702 (2011).
  * Utsu  _et al._ (1995) T. Utsu, Y. Ogata, _et al._ , Journal of Physics of the Earth 43, 1 (1995).
  * Weiss and Carmen Miguel (2004) J. Weiss and M. Carmen Miguel, Materials Science and Engineering: A 387-389, 292 (2004), 13th International Conference on the Strength of Materials.
  * Crane and Sornette (2008) R. Crane and D. Sornette, Proceedings of the National Academy of Sciences 105, 15649 (2008),
  * Sornette and Utkin (2009) D. Sornette and S. Utkin, Phys. Rev. E 79, 061110 (2009).
  * Baró  _et al._ (2013) J. Baró, A. Corral, X. Illa, A. Planes, E. K. H. Salje, W. Schranz, D. E. Soto-Parra, and E. Vives, Phys. Rev. Lett. 110, 088702 (2013).
  * Hainzl  _et al._ (2014) S. Hainzl, J. Moradpour, and J. Davidsen, Geophysical Research Letters 41, 8818 (2014),
  * Baró  _et al._ (2014) J. Baró, J.-M. Martín-Olalla, F. J. Romero, M. C. Gallardo, E. K. H. Salje, E. Vives, and A. Planes, Journal of Physics: Condensed Matter 26, 125401 (2014).
  * Davidsen  _et al._ (2017) J. Davidsen, G. Kwiatek, E.-M. Charalampidou, T. Goebel, S. Stanchits, M. Rück, and G. Dresen, Phys. Rev. Lett. 119, 068501 (2017).
  * Kumar  _et al._ (2020) P. Kumar, E. Korkolis, R. Benzi, D. Denisov, A. Niemeijer, P. Schall, F. Toschi, and J. Trampert, Scientific Reports 10, 626 (2020).
  * Dieterich (1972) J. H. Dieterich, Journal of Geophysical Research (1896-1977) 77, 3771 (1972),
  * Hainzl  _et al._ (1999) S. Hainzl, G. Zöller, and J. Kurths, Journal of Geophysical Research: Solid Earth 104, 7243 (1999),
  * Zöller  _et al._ (2006) G. Zöller, S. Hainzl, Y. Ben-Zion, and M. Holschneider, Tectonophysics 423, 137 (2006), spatiotemporal Models of Seismicity and Earthquake Occurrence.
  * Ben-Zion and Lyakhovsky (2006) Y. Ben-Zion and V. Lyakhovsky, Geophysical Journal International 165, 197 (2006),
  * Jagla (2014) E. A. Jagla, Phys. Rev. E 90, 042129 (2014).
  * Zhang and Shcherbakov (2016) X. Zhang and R. Shcherbakov, Scientific Reports 6, 36668 EP (2016), article.
  * Baró and Davidsen (2018) J. Baró and J. Davidsen, Phys. Rev. E 97, 033002 (2018).
  * Ruscher and Rottler (2020) C. Ruscher and J. Rottler, Soft Matter 16, 8940 (2020).
  * Touati  _et al._ (2009) S. Touati, M. Naylor, and I. G. Main, Phys. Rev. Lett. 102, 168501 (2009).
  * Baró and Davidsen (2017) J. Baró and J. Davidsen, European Physical Journal: Special Topics 226, 3211 (2017).
  * Kuhn and Daouadji (2019) M. R. Kuhn and A. Daouadji, Granular Matter 21, 10 (2019).
  * Bi  _et al._ (1989) H. Bi, G. Börner, and Y. Chu, Astronomy and Astrophysics 218, 19 (1989).
  * Reid (1910) H. F. Reid, The California Earthquake of April 18, 1906, Report of the State Earthquake Investigation Commission (1910).
  * Shimazaki and Nakata (1980) K. Shimazaki and T. Nakata, Geophysical Research Letters 7, 279 (1980),
  * Combe and Roux (2000) G. Combe and J.-N. Roux, Physical Review Letters 85, 3628 (2000).
  * Picard  _et al._ (2004) G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, The European Physical Journal E 15, 371 (2004).
  * Albaret  _et al._ (2016) T. Albaret, A. Tanguy, F. Boioli, and D. Rodney, Phys. Rev. E 93, 053002 (2016).
  * Lockner (1993) D. Lockner, International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts 30, 883 (1993).
  * Lennartz-Sassinek  _et al._ (2014) S. Lennartz-Sassinek, I. G. Main, M. Zaiser, and C. C. Graham, Phys. Rev. E 90, 052401 (2014).
  * Baró  _et al._ (2018) J. Baró, K. A. Dahmen, J. Davidsen, A. Planes, P. O. Castillo, G. F. Nataf, E. K. H. Salje, and E. Vives, Phys. Rev. Lett. 120, 245501 (2018).
  * Kun  _et al._ (2014) F. Kun, I. Varga, S. Lennartz-Sassinek, and I. G. Main, Phys. Rev. Lett. 112, 065501 (2014).
  * Kwiatek  _et al._ (2014) G. Kwiatek, T. H. W. Goebel, and G. Dresen, Geophysical Research Letters 41, 5838 (2014),
  * Tordesillas  _et al._ (2009) A. Tordesillas, J. Zhang, and R. Behringer, Geomechanics and Geoengineering: An International Journal 4, 3 (2009).
  * Vu and Weiss (2020) C.-C. Vu and J. Weiss, Phys. Rev. Lett. 125, 105502 (2020).
  * Janićević  _et al._ (2016) S. Janićević, L. Laurson, K. J. Måløy, S. Santucci, and M. J. Alava, Phys. Rev. Lett. 117, 230601 (2016).
  * Salditch  _et al._ (2020) L. Salditch, S. Stein, J. Neely, B. D. Spencer, E. M. Brooks, A. Agnon, and M. Liu, Tectonophysics 774, 228289 (2020).
  * Karimi  _et al._ (2017) K. Karimi, E. E. Ferrero, and J.-L. Barrat, Phys. Rev. E 95, 013003 (2017).
  * Anthony and Marone (2005) J. L. Anthony and C. Marone, Journal of Geophysical Research: Solid Earth 110 (2005), 10.1029/2004JB003399,
  * Frenkel  _et al._ (2008) G. Frenkel, R. Blumenfeld, Z. Grof, and P. R. King, Phys. Rev. E 77, 041304 (2008).
  * Henkes  _et al._ (2007) S. Henkes, C. S. O’Hern, and B. Chakraborty, Phys. Rev. Lett. 99, 038002 (2007).
  * Jagla (2015) E. A. Jagla, Phys. Rev. E 92, 042135 (2015).
  * Barés  _et al._ (2017) J. Barés, D. Wang, D. Wang, T. Bertrand, C. S. O’Hern, and R. P. Behringer, Phys. Rev. E 96, 052902 (2017).
  * Ma  _et al._ (2020) G. Ma, Y. Zou, K. Gao, J. Zhao, and W. Zhou, Geophysical Research Letters 47, e2020GL090458 (2020), e2020GL090458 2020GL090458,
  * Clark  _et al._ (2015a) A. H. Clark, A. J. Petersen, L. Kondic, and R. P. Behringer, Phys. Rev. Lett. 114, 144502 (2015a).
  * Clark  _et al._ (2015b) A. H. Clark, A. J. Petersen, L. Kondic, and R. P. Behringer, Physical review letters 114, 144502 (2015b).
  * Niemann  _et al._ (2012) R. Niemann, J. Baró, O. Heczko, L. Schultz, S. Fähler, E. Vives, L. Mañosa, and A. Planes, Phys. Rev. B 86, 214101 (2012).