# Lifespan dynamics of cluster conformations in stationary regimes in granular materials

Na Deng, Antoine Wautier, Antoinette Tordesillas, Yannick Thiery, Zhen-Yu

Yin, Pierre-Yves Hicher, François Nicot

► To cite this version:

Na Deng, Antoine Wautier, Antoinette Tordesillas, Yannick Thiery, Zhen-Yu Yin, et al., Lifespan dynamics of cluster conformations in stationary regimes in granular materials. Physical Review E, 2022, 105 (1), pp.1-16. 10.1103/PhysRevE.105.014902. hal-03572105

HAL Id: hal-03572105

https://hal.inrae.fr/hal-03572105

Submitted on 14 Feb 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Lifespan dynamics of cluster conformations in stationary regimes in granular materials

Na DENG

Grenoble Alps University, INRAE, UR ETNA, 2 rue de la Papeterie-BP 76, F-38402 St-Martin-d'Hères, France.

Antoine WAUTIER*

Aix-Marseille University, INRAE, UMR RECOVER,
3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence Cedex 5, France

Antoinette TORDESILLAS

School of Mathematics and Statistics, The University of Melbourne, Melbourne, Victoria, 3010, Australia

Yannick THIERY

BRGM (French Geological Survey), Risk and Prevention Division,
3 Av. Claude Guillemin, 45100 Orléans, France

Zhen-Yu YIN

Department of Civil and Environmental Engineering,
The Hong Kong Polytechnic University, Hong Kong, China

Pierre-Yves HICHER

Research Institute in Civil Engineering and Mechanics (GeM), UMR CNRS 6183,
Ecole Centrale de Nantes, 1 Rue de la Noë, 44300 Nantes, France

François NICOT

Grenoble Alps University, INRAE, UR ETNA, 2 rue de la Papeterie-BP 76, 38402 St-Martin-d'Hères, France.
Université Savoie Mont Blanc, Laboratoire EDYTEM, UMR 5204,
5 bd. de la Mer Caspienne, 73376 Le Bourget-du-Lac, France

(Dated: January 5, 2022)

We examine stationary regimes in granular materials from a dynamical systems theory (DST) perspective. The aim is to enrich the classical view of the critical state regime in granular materials, and more broadly, improve fundamental understanding of the underlying mesoscale mechanisms responsible for macroscopic stationary states in complex systems. This study is based on a series of discrete element method (DEM) simulations, in which 2D assemblies of non-uniformly sized circular particles are subjected to biaxial compression under constant lateral confining pressure. The lifespan and life expectancy of specific cluster conformations, comprising particles in force chains and grain loops, are tracked and quantified. Results suggest that these conformational clusters reorganize at similar rates in the critical state regime, depending on strain magnitudes and confining pressure levels. We quantified this rate of reorganization and found that the material memory rapidly fades, with an entirely new generation of force chains and grain loops replacing the old within a few percent strain.

**Keywords:** Lifespan, Life expectancy, Birth-death, Detailed balance, Grain loop, Force chain, Granular materials, Critical state, Dynamical systems

# I. INTRODUCTION

Stationary regimes exist in a great variety of complex systems, which comprise large numbers of interacting units [1–3]. At the macroscopic scale, such systems may be characterized by statistical descriptors that remain constant over time, while underlying local scale interactions between components continually change. In the case of granular materials, where the internal degrees of freedom

are large, a multitude of conformations (arrangements of particle structures) and conformational transitions [4–12] exist for which the macroscopic stress and porosity remain constant under continuous shearing. Such macroscopic states, named critical states (CS) in the mechanics community, form the focus of Critical State Theory (CST) [13–15] which stands as a cornerstone of constitutive theory for granular materials. Constitutive laws built on CST are however phenomenological and cannot account for the underlying conformational transitions which are responsible for the emergence of critical states. To this respect, only multiscale models may

<footnote>* Corresponding author: antoine.wautier@inrae.fr</footnote>
capture the microscale detailed balance responsible for¹¹⁸ critical state. But, many aspects of these transitions re-¹¹⁷ main poorly understood, and no multiscale model exists¹¹⁸ in which i) the minimal set of conformational transitions¹¹⁹ responsible for critical state and ii) the rate of the reor-¹²⁰ ganisation are embedded. Indeed, novel structures, tran-¹²¹ sitions and regimes are still being uncovered in granu-¹²² lar materials under shear, mostly thanks to the use of¹²³ Discrete Element Modelling (DEM), as well as advanced¹²⁴ experimental equipment [16].¹²⁵

One line of advance, which is pursued in this study, is to¹²⁶ adopt a dynamical systems theory (DST) approach to im-¹²⁷ prove fundamental understanding of the critical state by¹²⁸ focusing on conformational transitions from the perspec-¹²⁹ tive of two important classes of mesostructures, namely¹³⁰ force chains and grain loops [4–6] (see the definitions in¹³¹ Appendix B). Here we similarly employ a DST approach¹³² to uncover insights on these structures’ lifespan dynam-¹³³ ics. Our approach can help enrich multiscale constitutive¹³⁴ models, such as the promising H-model [17] which can in-¹³⁵ corporate the physics of grain loops and force chains to¹³⁶ the extent that the *emergent* critical states can be re-¹³⁷ produced. More broadly, a proper understanding of the¹³⁸ birth-death dynamics of these structures benefits a gamut¹³⁹ of applications where granular materials are subjected to¹⁴⁰ shear, including: the reduction of energy consumption¹⁴¹ and greenhouse gas emissions in processes like comminu-¹⁴² tion [18, 19], discovery and design of granular materials¹⁴³ with enhanced properties, and geohazard forecasting [20],¹⁴⁴ to name a few examples.¹⁴⁵

Prior studies have shown that the emergent mesostruc-¹⁴⁶ tures of force chains and grain loops are important build-¹⁴⁷ ing blocks for self-organization in 2D and 3D systems¹⁴⁸ [5, 6, 21–24]. Prior work has highlighted that some statis-¹⁴⁹ tical descriptors of the microstructure (e.g. fabric tensor, ¹⁵⁰ proportions of chained particles, proportions of the differ- ¹⁵¹ ent categories of loops) reach constant values at critical¹⁵⁰ state [22, 25–30]. Indeed birth-death dynamics and re-¹⁵² lated concepts like “lifespan”, is a promising way to think¹⁵¹ of the critical state regime. Granular shear simulations¹⁵² and experiments, including those focusing on stick-slip¹⁵³ dynamics, have shown that this regime is governed by¹⁵⁴ the competing mechanisms of growth/birth versus col-¹⁵⁵ lapse/death of force chains and these events drive evo-¹⁵⁶ lution towards, and during, stationary states of systems¹⁵⁷ under shear [4, 5, 31–33]. Experimental evidence of de-¹⁵⁸ tailed balance in sheared dense planar granular systems¹⁵⁹ is recently reported in Sun et al., 2021 [34]. On the¹⁶⁰ other hand, a DST analysis of DEM simulations for a¹⁶¹ 3D system in the critical state regime uncovered bistable¹⁶² dynamics [5], which is governed by broken detailed bal-¹⁶³ ance and underlying conformational transitions having¹⁶⁴ connections to the rattler dynamics reported in Wautier¹⁶⁵ et al [35]. Despite the significant attention paid to the¹⁶⁶ topic, the open question(s) that arise from all these stud-¹⁶⁷ ies remains: Why do granular systems evolve towards a¹⁶⁸ stationary limiting state under continued shear? What¹⁶⁹ mechanisms attract the system towards the critical state?¹⁷⁰

To quantify whether jamming dominates unjamming in evolving systems, the biology-inspired concepts of *lifespan* and *life expectancy* are particularly useful. Starting from a reference microstructure conformation, the lifespan of a mesostructure corresponds to the “time” it has continuously existed since its formation. Its life expectancy corresponds to the remaining time it will survive before disappearing. By comparing the distribution of lifespans and life expectancies, one can assess whether the birth rate of a particular class of mesostructure dominates its death rate. These concepts have been applied to granular materials for instance by [30, 36, 37]. Indeed, these concepts enable the quantification of the renewal rate of the microstructure in both 2D and 3D conditions under an external forcing time parameter. Specifically, it measures how fast the system is losing the memory of its past state under an external forcing. Consequently, they are particularly well adapted to treat granular materials from the perspective of dynamical systems with memory. In all these contexts, the time may correspond to physical time or any relevant forcing parameter, depending on the system of interest (for biaxial tests considered in the present study, the axial strain is the time parameter). The paper is organized as follows. In section II, the mechanical responses of granular samples under biaxial tests with discrete element modelling (DEM) are analyzed. This is followed in Section III by a discussion on critical state from a stationarity perspective. To this respect, the concepts of lifespan and life expectancy of mesostructures are applied to force chains and grain loops. Birth and death rates of chained particles and grain loops are analyzed in Section IV. The relations between microstructure dynamics and loading conditions are presented in Section V.

## II. NUMERICAL SETUP

The discrete element method (DEM) proposed by [38] is widely used to simulate granular materials. The open source DEM software YADE [39] is used in this study. The sample considered is a quasi-2D granular sample subjected to biaxial loadings under constant lateral pressure. The assembly consists of a single layer of 20,000 spherical particles contained to planar displacements within a domain of 1 meter wide per 1.5 meters high, as shown in Fig. 1. The particle radii are uniformly distributed with an average diameter $d_{50} = 0.008$ m and a size ratio $d_{max}/d_{min} = 2$. The contact law is described by a standard cohesionless elasto-frictional law. The contact parameters between two grains contain normal and tangential linear springs of respective stiffness $k_n$ and $k_t$, as well as a sliding limit characterized by a friction angle $\phi = 35^{\circ}$. $E = k_n/D_s$, given as 300 MPa, where $D_s = 2R_1R_2/(R_1 + R_2)$ and $R_1, R_2$ are the radii of particles.

For sample preparation, two schemes can be used to obtain an isotropic compression of the sample in Yade
<figure><img src="image_2.png" /><figcaption>FIG. 1. Quasi-2D DEM sample for biaxial tests where loading parameters are $\varepsilon_2$ and $\sigma_1$. The direction 2 is always the vertical loading control direction at the representative element volume (REV) scale.</figcaption></figure>

DEM: the boundary moving scheme and the internal compacting scheme. In this study, particles are homothetically enlarged, at first, while keeping the boundary walls fixed to generate a preliminary sample with a certain level of confining pressure. Then, to get closer to an experimental set up, the sample undergoes an isotropic consolidation by imposing an equivalent incremental strain on the boundaries in the vertical and lateral directions up to a specific consolidation state. Three confining levels are considered, characterized by a dimensionless quantity $P_0$ representing the ratio $p_0/E = \frac{\sigma_1^0+\sigma_2^0}{2E}$ of 206 the confining pressure over the stiffness, as shown in Table I. They are labeled as A, B and C from the lowest to 209 the highest. For preparing dense and loose samples, the 210 contact friction angle is set to 2° and 35°, respectively, 211 during the confining stage. Dense and loose samples are 212 prepared and named as indicated in Table II. The initial 213 void ratio $e_0$ and the ratio $<u_n>/d_{50}$ of the average 214 contact overlap and particle size of the samples prepared are presented. More details about sample preparation can be found in [8, 25].

The stress and strain states are described using dimensionless quantities, the normalized deviatoric stress $Q$ be-217 ing defined as $Q = (\sigma_2 - \sigma_1)/E$, and the volumetric strain 218 being defined as $\varepsilon_v = \varepsilon_1 + \varepsilon_2$, where $\sigma_1$ and $\sigma_2$ are the 219 principal stresses and $\varepsilon_1$ and $\varepsilon_2$ are the principal strains, 220 as shown in Fig. 1. Axes 1 and 2 refer to lateral and vertical directions, respectively. Note that stresses $\sigma_1$ and 221 $\sigma_2$ are computed in quasi-2D DEM simulations by introducing an out of plane dimension equal to the largest 222 grain diameter ($d_{max}$). If full 2D frameworks were pre-223 ferred, the stress $\sigma$ would be expressed in N/m and the 224 dimensionless stresses $P_0$ and $Q$ should be normalized by 225 $Ed_{max}$ instead of $E$.

<figure><img src="image_3.png" /><figcaption>FIG. 2. Evolution of normalized deviatoric stress $Q$ and volumetric strain $\varepsilon_v$ with respect to the axial strain $\varepsilon_2$ along biaxial tests. Note that soil mechanics convention is adopted with positive compression and positive contraction.</figcaption></figure>

### III. THE SO-CALLED CRITICAL STATE

In the literature, it is common to report fluctuations of the deviatoric stress and the volumetric strain around steady-state values rather than a smoothly constant evolution under biaxial loading. The evolution of micro and meso indices like fabric tensor and ratio of loops population also exhibits fluctuations at critical state. In this section, macroscopic stress and strain are put in relation with the evolution of force chains and grain loops populations, respectively.

#### A. Stress-strain analysis

Macroscopic responses of the tested samples are illustrated in Fig. 2, in which normalized deviatoric stress $Q$ and volumetric strain $\varepsilon_v$ evolve with respect to the axial strain $\varepsilon_2$. The figure shows typical stress and strain responses of dense and loose samples under the biaxial loadings.

In the dense cases, the normalized deviatoric stress $Q$ increases at first then decreases to a stationary value, while the sample contracts first and then dilates before reaching a steady state. With an increase in the confining pressure there is an increase in the shear strength and a decrease in the dilation. On the other hand, the normalized deviatoric stress and volumetric strain in the loose samples
<table><caption>TABLE I. Different initial confining levels</caption><thead><tr><th>Normalized initial confining pressure<br>P<sub>0</sub></th><th>A<br>1.3 × 10<sup>-4</sup></th><th>B<br>3.3 × 10<sup>-4</sup></th><th>C<br>1.3 × 10<sup>-3</sup></th></tr></thead></table>

<table><caption>TABLE II. Initial void ratios $e_0$ and average overlaps in prepared samples.</caption><thead><tr><th>Sample</th><th>A-loose</th><th>B-loose</th><th>C-loose</th><th>A-dense</th><th>B-dense</th><th>C-dense</th></tr></thead><tbody><tr><td>e<sub>0</sub></td><td>0.267</td><td>0.261</td><td>0.234</td><td>0.197</td><td>0.191</td><td>0.164</td></tr><tr><td>&lt; u<sub>n</sub> &gt; /d<sub>50</sub></td><td>9.2 × 10<sup>-4</sup></td><td>2.1 × 10<sup>-3</sup></td><td>7.6 × 10<sup>-3</sup></td><td>6.6 × 10<sup>-4</sup></td><td>1.6 × 10<sup>-3</sup></td><td>5.6 × 10<sup>-3</sup></td></tr></tbody></table>

<figure><img src="image_4.png" /><figcaption>FIG. 3. Incremental deviatoric strain (IDS) maps at the critical states, estimated for axial strain increments of 0.1%. Two samples are considered: (a) B-dense sample ($\epsilon_2 = 6.1\%$) and (b) B-loose sample ($\epsilon_2 = 8.02\%$).</figcaption></figure>

evolve monotonously in the course of loading. It can be observed that, at critical state, both normalized deviatoric stress $Q$ and volumetric strain $\epsilon_v$ fluctuate around a constant value. Based on these mechanical responses, both the stationarity and dynamical properties of CS will be investigated in a micro-mechanical framework. The incremental deviatoric strain distribution can be analyzed to characterize kinematic patterns, as introduced in literature [30, 40]. The maps are computed based on the interpolation of particle incremental displacements over meshed sample [30, 41]. When dense samples are considered, shear bands may develop, traversing the whole sample as shown in Fig. 3. It is the spatial domain within shear bands for dense cases and the whole sample for loose cases that are considered in this study [26, 30]. The spatial extension of the shear band was determined based on the method proposed in [41].

## B. Force chains and normalized deviatoric stress

Force chains have attracted considerable interest owing to their significant role in force transmissions in granular systems [23, 42]. To consider force chains as individual mesostructures (and not as a network of strong contacts), the following three conditions were proposed by Peters et al. (2005) [43]:

- The number of contacting particles in a force chain

is not less than 3.

- The major principal stress in each particle belonging to a force chain is larger than the averaged assembly particle major principal stress.

- The major principal stress direction is aligned with the contact direction with the deviation angle less than 45°.

In this paper, the particles belonging to force chains are called *chained particles*. The proportions of chained particles is given by $r_{cha} = N_{cha}/N_{tot}$, where $N_{cha}$ and $N_{tot}$ are the number of chained particles and the total number of particles within the investigated domain, respectively. The ratio within the shear band $r_{cha}^*$ and the whole sample $r_{cha}$ are compared in Fig. 4.

It can be seen that the evolution with strain of both $r_{cha}^*$ and $r_{cha}$ tends to reach a constant value (with fluctuations) when the normalized deviatoric stress and the volumetric strain become constant at CS. $r_{cha}^*$ is larger than $r_{cha}$ after the generation of the shear band, with greater fluctuations. This may be results in a stronger rearrangement of chained particles within the shear band. It is worth noting that the ratio of chained particles fluctuating around a mean constant value highlights that important force chain reorganizations take place at the critical state.

## C. Grain loops and volumetric strain

In a 2D granular material, grain loops, enclosed by contact branches, formed by tessellating the material area. It is usually categorized according to the edge number (from L3 with three particles to L6p for loops containing six or more particles), a factor which has a significant influence on the behavior of a loop, especially its deformability. The percentage of each category reaches a steady state that corresponds to the critical state regime at the sample scale [30]. It is worth noting that the reason for categorizing $L_{6p}$ is that the percentage of $L_{7p}$ (percentage of categories with a side number equal or greater than 7) evolves in a similar way as that of $L_6$ [30]. In this section, we link the grain loop populations (L3, L4, L5, L6p) to the volumetric evolution along the biaxial test by introducing the volumetric strain of the different
<figure><img src="image_5.png" /><figcaption>FIG. 4. Proportions of chained particles. Two samples are considered: (a) B-dense sample and (b) B-loose sample. Note that $r_{cha}^*$ represents the value when the spatial domain of shear band are considered, and $r_{cha}$ when the whole sample is considered. The normalized deviatoric stresses are recalled.</figcaption></figure>

grain loops. The incremental strain tensor of each type³²⁴ of loop $i$ can be expressed as follows:

$$d\bar{\epsilon}_{Li} = \frac{1}{|\Omega|} \sum_{k} |\Omega_{Li}^{k}| d\bar{\epsilon}_{Li}^{k} \quad i \in (3, 4, 5, 6p) \quad (1) \quad 325$$

where $\Omega$ is the volume of the global domain considered³²⁵ and $\Omega_{Li}^k$ is the volume of a given loop $k$ belonging to the³³⁰ category $L_i$. The incremental strain $d\bar{\epsilon}_{Li}^k$ of each loop 'k'³³¹ can be computed assuming a uniform deformation within³³² that loop according to the strain definition proposed by Liu et al (2018) [41].

The contribution of local incremental strains to the global³³³ incremental strain can be expressed by equation (2):

$$< d\bar{\epsilon} > = \frac{1}{|\Omega|} \int_{\Omega} d\bar{\epsilon} dS = d\bar{\epsilon}_{L3} + d\bar{\epsilon}_{L4} + d\bar{\epsilon}_{L5} + d\bar{\epsilon}_{L6p} \quad (2) \quad 335$$

Figure 5 presents the volumetric strain of each category³³⁸ of loops and their integration during the biaxial loading³³⁹ as well as the evolution of the global volumetric strain.³⁴⁰ Figure 5 shows a good agreement between the macroscop-³⁴¹ ically computed volumetric strain and the mesoscopically³⁴² averaged one. Despite reaching a zero volume change on³⁴³ average, it can be seen that each loop category reaches³⁴⁴ constant dilative or contractive rates at critical state.³⁴⁵ The dilatancy from $L_3$, $L_4$ and $L_5$, and the contractancy³⁴⁶ from L6p should be related with the inherent property of³⁴⁷ each loop: the grains belonging to $L_i$ ($i \le 5$) have limited³⁴⁸ space to move inward to induce contractancy, whereas³⁴⁹ L6p loops have larger void space to move inward. As³⁵⁰ result, critical state can only be achieved from a collec-³⁵¹ tive process by loops transforming from one category to³⁵² another.

## IV. HIDDEN DYNAMICS AT CRITICAL STATE

As highlighted in the previous sections, the jamming and unjamming behavior at CS includes dynamic mechanisms at the mesoscale [4, 31, 33, 34]. To this respect, the dynamical rate are investigated in this section through computing the lifespan and life expectancy of chained particles and grain loops. The incremental evolutions of the sets of chained particles and grain loops are checked for small axial strain increments of $\delta\epsilon_2 = 0.1\%$.

### A. Generating and vanishing process of chained particles

All chained particles are tracked from birth, i.e. the first time they become part of a given force chain, to death, i.e. the first time they no longer belong to any force chain. This enables to construct probability density functions (PDF) for the lifespan and the life expectancy of chained particles at any time of the simulation. The lifespan and life expectancy PDFs are displayed in Fig. 6 for two arbitrary microstructure conformations (M and N) in the critical state regime. Because the sample is known to experience shear band localization (see Fig. 3), the analysis is restricted to particles belonging to the shear band in Fig. 6.

At the reference conformation M (resp. N) in Fig. 6 (a), the PDFs for the lifespan and life expectancy demonstrate the generating and vanishing processes of chained particles within the shear band domain. Few chained particles existing at M or N comes from the initial conformation, which proves that the initial arrangement of particles has been erased when the sample reaches CS.
<figure><img src="image_6.png" /><figcaption>FIG. 5. Volumetric strain of different types of loops and the global assembly during the biaxial tests in B-dense sample (a) and B-loose sample (b).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6. Lifespan and life expectancy PDFs (probability density functions) of chained particles in two samples: (a) B-dense sample where only the spatial domain within the shear band is considered and (b) B-loose sample with the whole domain considered. For each sample, two reference conformations belonging to the stationary regime are adopted, labeled as M and N. The normalized deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled.</figcaption></figure>

The symmetric shape of the lifespan and life expectancy$^{364}$ PDFs shows that the rearrangement of chained particles$^{365}$ reaches an equilibrium at CS, which is expected for a$^{366}$ stationary regime. More interestingly the shape of the$^{367}$ PDFs for the reference conformations M and N is almost$^{368}$ the same. Force chains are always rearranging at crit-$^{369}$ ical state and the rate of renewal of the set of chained$^{370}$ particles is constant. These features can be observed in$^{371}$ B-loose with lower renewal rates from Fig. 6 (b).$^{372}$ The PDFs have been averaged to get rid of statistical$^{373}$

noise coming from the finite size of the chained particle population at different strain levels belonging to the critical regime ($\epsilon_2$ = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the B-dense sample and $\epsilon_2$ = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the B-loose sample). Since force chains have a limited life duration at critical state, the six microstructure conformations can be considered as statistically independent (not correlated), which allows us to build a single PDF from the six data sets. The corresponding lifespan and life ex-
<figure><img src="image_8.png" /></figure>

FIG. 7. The averaged PDFs of chained particles with respect to the incremental strain evolution ($\Delta\varepsilon_2$), over six microstructure conformations at different strain levels belonging to the critical regime ($\varepsilon_2$ = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the B-dense sample and $\varepsilon_2$ = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the B-loose sample).

pectancy PDFs are presented in Figure 7 with respect to incremental evolution of the axial strain ($\Delta\varepsilon_2$). Thus, $\Delta\varepsilon_2$ at current state is 0. $\Delta\varepsilon_2 > 0$ and $\Delta\varepsilon_2 < 0$ correspond to the lifespan and life expectancy of chained particles, respectively.

The renewal rate in dense and loose samples appears to be different at first glance, which would contradict the existence of a unique critical state, whatever the initial sample density. To find the origin of the apparent faster renewal rate within the shear band in the dense sample, it is meaningful to compare the incremental strain magnitude $d\varepsilon_m = \sqrt{d\varepsilon_1^2 + d\varepsilon_2^2}$ in the shear band domain for the dense sample, to the incremental strain magnitude in the whole sample domain for the loose sample. Similar to the incremental strain of each type of loops computed in Subsection 3.3, the incremental strain within the shear band can be computed by averaging the incremental strains of loops within the shear band over the area of the shear band. When averaged over the whole sample domain $d\varepsilon_m$ is the same in dense and loose samples at critical state, as shown in Fig. 8. However, because of the strain localization in the dense sample, a larger $d\varepsilon_m^*$ develops within the shear band. The ratio between $<d\varepsilon_m^*>$ and $<d\varepsilon_m>$ at CS inside the shear band is around 2.8. In other words, an axial strain increment of 1% from the boundaries of the sample imposes a strain increment of 2.8% in the shear band domain because of strain localization. This means that the actual forcing parameter (the "time") inside the shear band is 2.8 time faster than for the whole sample. As a result, the axial strain applied at sample scale is not the right forcing parameter in the dense case as the strain is 2.8 times more intense in the

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

FIG. 8. The strain magnitudes within the shear band domain (B-dense*), whole dense sample (B-dense) and the whole loose sample (B-loose). The strain magnitude ratio between B-dense* and B-dense is around 2.8.

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

FIG. 9. PDFs for chained particles from B-loose sample and from B-dense sample stretched by the strain magnitude ratio 2.8.

shear band domain. In order to account for this observation, we need to stretch the curves of B-dense* by the ratio 2.8 to make the lifespan and life expectancy PDF comparable between the dense and loose sample cases. Then, it can be observed in Fig. 9 that the PDF in dense and loose samples are coinciding. Similar results can be observed for different confining pressure levels, as shown in Fig. 20 and 21 in Appendix A. It can be concluded that the strain magnitude (observed in homogeneous domain of the samples) is the major factor affecting the rate of rearrangement of chained particles. The fact that the PDFs coincide also demonstrates that the strain rate does not affect the reorganization process (at least for
the order of magnitude considered). We can also conclude that the uniqueness of the critical state (namely whatever the initial density, statistical descriptors converge towards the same value in biaxial tests of similar confining pressure) that holds for statistical descriptors, also holds for the underlying dynamical processes.

## B. Generating and vanishing process of grain loops

To further understand the critical state dynamics, the generating and vanishing of grain loops are analyzed for loops composed of 3, 4, 5 and 6 or more particles. Lifespan and life expectancy PDFs for grain loops are computed. Similar observations as reported in the previous subsections are visible in Fig.10 - 12: (1) PDFs show symmetrical shapes at CS which suggests birth and death rates of grain loops are equal at CS; (2) PDFs within the shear band from B-dense coincide with those obtained in the loose case after being stretched by the appropriate ratio of 2.8. These results extends the conclusion drawn by Zhu et al. (2016) [30], that the whole loose sample and the inner part of the shear band of the dense sample share the same statistical distribution of mesostructures once the critical state regime is reached.

Based on Figures 6 - 12, a novel way for understanding CS from the perspective of a mesoscale dynamics can be inferred: the critical state results from the generating and vanishing of meso-clusters until a dynamic equilibrium is reached under continuous shearing. Even though dynamic rearrangements permanently occur, at any time, the material admits the same fabric characterized by the fabric tensor, the distribution of loop orders etc.. The particle rearrangements at critical state can be quantitatively characterized by the lifespan and life expectancy of force chains and grain loops. These observations still hold when other confining pressure levels (A and C) are considered as shown in Appendix A.

## V. MICROSTRUCTURE REORGANIZATION
DYNAMICS UNDER DIFFERENT CONFINING
PRESSURE LEVELS

In this section, the influence of the confining pressure on the dynamics in the three loose samples is analyzed. As shown in Figs.13 and 14, the renewal rate decreases with the increase in the confining pressure. The increase in confining pressure helps the mesostructures live longer (with respect to the forcing parameter of the system, i.e. the incremental axial strain). This is consistent with the fact that larger confining pressure leads to larger elastic energy storage and larger overlap at contacts, which delays contact opening and contact sliding.

In order to quantify the rate of microstructure reorganization with respect to the axial strain, it is interesting to propose a fit for the PDFs presented so far. As a first guess, an exponential fit with one characteristic strain

<table><caption>TABLE III. The fitting parameters for different confining pressure levels.</caption><thead><tr><th>Sample</th><th>Parameters</th><th>Chained particles</th><th>L3</th><th>L4</th><th>L5</th><th>L6p</th></tr></thead><tbody><tr><td rowspan="4">A-loose</td><td>ε<sub>c1</sub></td><td>0.25%</td><td>0.30%</td><td>0.33%</td><td>0.30%</td><td>0.25%</td></tr><tr><td>ε<sub>c2</sub></td><td>2.29%</td><td>1.61%</td><td>1.92%</td><td>1.81%</td><td>1.59%</td></tr><tr><td>P<sub>c</sub></td><td>0.165</td><td>0.229</td><td>0.228</td><td>0.288</td><td>0.290</td></tr><tr><td>ε<sub>c1</sub>/ε<sub>c2</sub></td><td>9.02</td><td>5.44</td><td>5.78</td><td>5.97</td><td>6.41</td></tr><tr><td rowspan="4">B-loose</td><td>ε<sub>c1</sub></td><td>0.35%</td><td>0.43%</td><td>0.41%</td><td>0.34%</td><td>0.34%</td></tr><tr><td>ε<sub>c2</sub></td><td>3.30%</td><td>2.57%</td><td>2.79%</td><td>2.40%</td><td>2.25%</td></tr><tr><td>P<sub>c</sub></td><td>0.146</td><td>0.205</td><td>0.162</td><td>0.209</td><td>0.271</td></tr><tr><td>ε<sub>c1</sub>/ε<sub>c2</sub></td><td>9.35</td><td>5.92</td><td>6.89</td><td>7.03</td><td>6.71</td></tr><tr><td rowspan="4">C-loose</td><td>ε<sub>c1</sub></td><td>0.61%</td><td>0.77%</td><td>0.64%</td><td>0.59%</td><td>0.50%</td></tr><tr><td>ε<sub>c2</sub></td><td>5.33%</td><td>5.56%</td><td>5.60%</td><td>4.35%</td><td>3.73%</td></tr><tr><td>P<sub>c</sub></td><td>0.129</td><td>0.124</td><td>0.113</td><td>0.169</td><td>0.209</td></tr><tr><td>ε<sub>c1</sub>/ε<sub>c2</sub></td><td>8.76</td><td>7.25</td><td>8.74</td><td>7.41</td><td>7.52</td></tr></tbody></table>

may be proposed as $PDF(Δε_2) = \frac{a}{ε_c}e^{-|Δε_2|/ε_c}$. $ε_c$ accounts for the typical size of the axial strain increment $Δε_2$ needed to renew the whole mesostructures. The smaller $ε_c$, the faster the microstructure renewal under the external forcing $Δε_2$. However, such a fit cannot capture well the shape of the whole PDFs. Depending on $ε_c$, the fit is able to account either for the beginning or for the tail of the PDFs. This observation motivated the use of a double exponential fit in the following form:

$$PDF(Δε_2) = P_c \frac{1}{ε_{c1}} e^{-\frac{|Δε_2|}{ε_{c1}}} + (1 - P_c) \frac{1}{ε_{c2}} e^{-\frac{|Δε_2|}{ε_{c2}}} \quad (3)$$

where $ε_{c1}$ and $ε_{c2}$ ($ε_{c1} < ε_{c2}$) correspond to two characteristic strains accounting for the dynamic evolution of the system. $ε_{c1}$ and $ε_{c2}$ account for the typical amplitude of the axial strain increments $Δε_2$ needed to renew the mesostructures on the short term and long term, respectively. The characteristic probability $P_c \in [0, 1]$ can show the relative contribution of the short and long term mechanisms. When $P_c < 0.5$, it is the long term mechanism that dominates; the short term mechanism will be more important when $P_c > 0.5$; the two mechanisms are evenly matched when $P_c = 0.5$. Figs. 15 and 16 show that double exponential fit is able to account for the shape of all the lifespan and life expectancy PDFs. The values of $ε_{c1}$ and $ε_{c2}$ for different confining pressure levels are presented in Table III and shown in Fig. 17 with respect to the normalized critical mean stress $P = \frac{σ_1+σ_2}{2E}$.

As shown in Table III and Fig. 17, large loops (L5 and L6p) live shorter (faster reorganization) than small loops (L3 and L4) and force chains. It is consistent with the assumption that loops connected to force chains open prior to force chain bending [44]. When it comes to the effect of confining pressures, the renewal rate of mesostructures which can be characterized by $\frac{1}{ε_{c1}}$ and $\frac{1}{ε_{c2}}$ decreases with the increase in the confining pressure. This effect differs according to the type of mesostructures. L3 and L4 seem to be more sensitive to the confining pressures, especially in the long term mechanism, which could be due to the fact that these small, less compressible loops can store larger (stabilizing) elastic energy before breaking. On
<figure><img src="image_11.png" /><figcaption>FIG. 10. Lifespan and life expectancy PDFs (probability density functions) of grain loops in two samples with respect to the axial strain $\epsilon_2$: (a) B-dense sample where only the spatial domain within the shear band is considered and (b) B-loose sample with the whole domain considered. For each sample, two reference conformations belonging to the stationary regime are adopted. The normalized deviatoric stress curve (dot-dashed line) and the volumetric strain (dashed line) curves are recalled in the subfigure of Loop 4.</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>FIG. 11. The averaged PDFs of grain loops with respect to the incremental strain evolution ($\Delta\epsilon_2$), over six microstructure conformations at different strain levels belonging to the critical regime ($\epsilon_2$ = 0.0610, 0.0708, 0.0806, 0.0904, 0.1002 and 0.1100 for the B-dense sample and $\epsilon_2$ = 0.0605, 0.0708, 0.0802, 0.0906, 0.1005 and 0.1106 for the B-loose sample).</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIG. 12. PDFs for grain loops from B-loose sample and from B-dense sample stretched by the strain magnitude ratio 2.8.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIG. 13. Lifespan and life expectancy PDFs for chained particles under the different confining pressure levels with samples A-loose, B-loose and C-loose considered.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 14. Lifespan and life expectancy PDFs for grain loops with respect to the incremental strain evolution ($\Delta\epsilon_2$) under the different confining pressure levels with samples A-loose, B-loose and C-loose considered.</figcaption></figure>

the other hand, large loops (L5 and L6p) are more compressible, then less sensitive to confining pressure, as the increase in $p$ does not have as much stabilization effect as for smaller loops. As for chained particles, even though force chains can share particles with large loops, their stability (and thus lifespan) is strongly dependent on the existence of small order loops around them. Thus, the effect of confining pressure on force chains lifespan stem from the different responses of loops under the confining pressure according to their category (small vs. large
<figure><img src="image_16.png" /><figcaption>FIG. 15. Double-exponential fits given for the chained particle PDFs with respect to the absolute value of incremental evolution $|\Delta\epsilon_2|$ under the different confining pressure levels.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>FIG. 16. Double-exponential fits given for the loop PDFs with respect to the absolute value of incremental evolution $|\Delta\epsilon_2|$ under the different confining pressure levels.</figcaption></figure>

ble exponential PDF is the signature of the existence of two mechanisms responsible for microstructure reorganization acting over two different “time” scales. A first mechanism provokes microstructure reorganizations over small axial strain scales of typical magnitude $\epsilon_{c1}$, while a second mechanism induces microstructure reorganizations over much larger axial strain scales of typical magnitude $\epsilon_{c2}$. The existence of two reorganization mechanisms is indeed consistent with the visualization of the incremental deviatoric strain map in Fig. 3 for the loose sample. The map is far from being homogeneous over the whole sample domain, and zones with large increments of deviatoric strain are visible (in the form of red zones in Fig. 3.). Contrary to shear bands, these zones are not persisting when the axial strain increases as shown. Indeed, these zones may be seen as shear transformation zones (STZ) [45] and characterized by the recent concept of shear chain [46]. As a result, we conjecture that $\epsilon_{c1}$ relates to the axial strain needed for a shear transition zone to develop, while $\epsilon_{c2} > \epsilon_{c1}$ corresponds to the axial strain needed to develop shear transition zones everywhere in the sample domain. This conjecture can be initially verified by the color maps of lifespan and life expectancy of loops in the sample B-loose, as shown in Figs. 18 and 19, respectively.

The short and long threshold are set to $4\delta\epsilon_2$ and $30\delta\epsilon_2$, respectively. In Fig.18 (a), red loops are newly-generated at the current conformation. The loops that have lived longer than $4\delta\epsilon_2$ are filled in grey. Fig.18 (a) shows that the microstructure is experiencing a series of localized reorganizations (visible as clusters of loops of the same colors) affecting a limited proportion of the sample domain. In the long term, these series of reorganizations will affect the whole sample domain as shown in Fig.18 (b), where loops that have been existing not less than $30\delta\epsilon_2$ are very limited. In Fig.19 (a), grey loops are those that will live longer than four steps. It can also be observed in Fig.19 (b) that beyond the long threshold, only a small part of the loops will still be alive.

## VI. CONCLUSION

We revisit the critical state which has been emphasized as a stationary state in geomechanics and characterize it using a dynamical systems analysis of grain structures and their rearrangements. We have shown that the critical state results from a balance between the birth and death of meso-clusters. We have established that force chains and grain loops have a regular, short life duration (0.4-3% compared to the whole 15% loading) at critical state, corresponding to a rapid fading of the material memory along its loading history. Critical state "forgot" not just the initial microstructure conformation but its recent past in a sample under continuous shearing. We have shown that the uniqueness of the critical state that holds for statistical descriptors also holds for the underlying dynamical processes. By fitting the data, we have

ones), which may explain why the slope of the $\epsilon_{c1}$ and $\epsilon_{c2}$ for chained particles is smaller than the one for small loops but larger than the one for large loops. In addition, the long term mechanism increasingly has dominate effect, since $P_c$, always being less than 0.5, decreases with the growth of critical pressures. The form of the proposed fit may appear arbitrary at first glance; but a double
<figure><img src="image_18.png" /><figcaption>FIG. 17. Two characteristic strains ($\epsilon_{c1}$ and $\epsilon_{c2}$), the characteristic probability ($P_c$) and the ratio $\epsilon_{c1}/\epsilon_{c2}$ with respect to the normalized mean stress ($P$) at the critical state. The normalized critical mean stress under the different confining pressure levels (A, B and C) are $1.97 \times 10^{-4}$, $4.77 \times 10^{-4}$ and $18.63 \times 10^{-4}$, respectively.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>FIG. 18. Color maps of the lifespan of grains loops for the reference axial strain $\epsilon_2 = 8.02\%$. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%.</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>FIG. 19. Color maps of the life expectancy of grains loops for the reference axial strain $\varepsilon_2$ = 8.02%. Two thresholds are considered corresponding to the two characteristic strains: (a) 0.4% and (b) 3%.</figcaption></figure>

shown that force chains and loops of different sizes have 618 different lifespan and life expectancy which increases with 619 the normalized mean stress level. Larger loops have a 620 shorter life due to their deformability and force chain 621 life duration is similar to the life duration of small or-622 der loops which is consistent with the known relationship 623 between force chains and their small order supporting 624 loops. Moreover, the microstructure reorganization has been proved to rely on two mechanisms that acts over 625 very different axial strain increments. A local mecha- 625 nism related to sheared mesostructures provoke a rapid 626 (i.e. for small axial strain increments) reorganization of 627 the microstructure, while the characteristic strain that 628 relates to the nucleation of new sheared mesostructures, 629 account for the reorganizations over larger axial strain 630 increments. The recent concept of shear chain [46] may 631 also shed further light on the elastic and plastic mecha- 632 nisms responsible for the observed dynamics of the criti- 633 cal state. The velocity of strain waves in the sample mer- 634 its future investigation, since the dynamic load transfer 635 characteristics can be also interpreted within the context 636 of wave propagation theory [47]. In particular, such work 637 may help elucidate whether the unjammed microstruc- 638 ture retains a partial memory of its previous jammed, 639 state. The findings presented in this paper also provide 640 further clues on how to explicitly incorporate CS fea- 641 tures in micro-mechanical constitutive models, as well as 642 to limit the number of different loading paths needed to 643 construct the data base in data-driven models because 644 of limited hysteresis effect. In relation to a particular 645 micromechanical model, it would be possible to consider 646 the rearrangement of the hexagonal loop mesostructures 647

relies on within the H-model[17, 48] to achieve a better description of the critical state, based on the renewal rate of mesostructures as quantified in this study. Although the present study deals with the critical state dynam- ics in granular materials, the concepts and methodology proposed here are generic and may be applied to uncover the hidden dynamics of many other complex systems.

## ACKNOWLEDGEMENT

The authors express their sincere thanks to the IRN (International Research Network) GeoMech for promoting positive and convivial interactions and discussions among researchers. The support from CSC (China Scholarship Council) under the Grant CSC Number 201801810030 and from the BRGM (French Geological Survey), France under the contract RP19DRP023 is gratefully acknowledged. The constructive comments by the reviewers are appreciated.

## Appendix A: Complementary results on critical state dynamics

In addition to the results presented in the main part of the paper, figures 20 and 21 provides complementary results on a) the detailed balance achieved at critical state for varying confining pressure levels and b) the uniqueness of the PDF for lifespan and life expectancy in dense and loose samples provided the strain rate is rescaled in the dense case to account for strain
<figure><img src="image_21.png" /><figcaption>FIG. 20. PDFs for chained particles under two confining pressure levels: (a) A, PDFs from the A-dense are stretched by strain magnitude ratio 1.9; (b) C, the C-dense stretched by 3.2.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIG. 21. PDFs for grain loops under two confining pressure levels: (a) A, PDFs from the A-dense are stretched by strain magnitude ratio 1.9; (b) C, the C-dense stretched by 3.2.</figcaption></figure>

## Appendix B: Glossary of terms

### Lifespan and life expectancy

### Grain loop
rials. The term grain loop refers here to the concept⁶⁷⁹ of minimal cycles in graph theory. In 2D, grain loops⁶⁸⁰ enclose polygons that forms a partition of the sample⁶⁸¹ domain able to characterize the dilatant/contractive⁶⁸² behavior of granular materials. Even if cycles exist in⁶⁸³ 3D, the bijective link with a partition of the material⁶⁸⁴ volume is lost which restrict the definition of grain loops⁶⁸⁵ to 2D conditions. Current work aims at generalizing the⁶⁸⁶ definition of such kind of mesostructures to both 2D and⁶⁸⁷ 3D conditions [49]. Its difference to the minimal cycles⁶⁸⁸ in literature [22] is that grain loops have the connotation⁶⁸⁹ of volume.

## Force chain

Force chain consists of a set of particles within a⁶⁹⁴ compressed granular material that are held together and⁶⁹⁵ jammed into place by a network of mutual compressive⁶⁹⁶ forces [43]. To consider force chains as individual⁶⁹⁷ mesostructures (and not as a network of strong contacts⁶⁹⁸ [42]), the definition used in the paper includes three⁶⁹⁹

thresholds as detailed in Subsection III B.

## Critical state (CS)

Critical state in soil mechanics is a macroscopic state where stress and porosity of a material remain constant under continuous shearing. In standard soil mechanics, a set of critical state points form a line in the $(e, p, q)$ space $(e$ being the void ratio, $p$ the mean stress and $q$ the deviatoric stress), known as the critical state line (CSL). From a statistical physics viewpoint, critical state is a steady state relying on equivalent microstructure conformations. Macroscopic variables (such as stress or void ratio) computed on these different microstructure conformations are equal.

## Conformation

Conformation is used to describe any spatial arrangement of grains in a given structure, while conformational transition is a change in the conformation [4].

[1] M. R. Gardner and W. R. Ashby. Connectance of large⁷³⁶ dynamic (cybernetic) systems: critical values for stabil-⁷³⁷ ity. *Nature*, 228(5273):784–784, 1970.

[2] R. M. May. Will a large complex system be stable? *Na-739 ture*, 238(5364):413–414, 1972.

[3] L. Papadopoulos, M. A. Porter, K. E. Daniels, and D. S.⁷⁴¹ Bassett. Network analysis of particles and grains. *Journal of Complex Networks*, 6(4):485–565, 2018.

[4] A. Tordesillas, D. M. Walker, G. Froyland, J. Zhang, and⁷⁴⁴ R. P. Behringer. Transition dynamics and magic-number-⁷⁴⁵ like behavior of frictional granular clusters. *Physical Review E*, 86(1):011306, 2012.

[5] S. Pucilowski, A. Tordesillas. Rattler wedging and force⁷⁴⁸ chain buckling: metastable attractor dynamics of lo-⁷⁴⁹ cal grain rearrangements underlie globally bistable shear⁷⁵⁰ banding regime. *Granular Matter*, 22(1):1, 2019.

[6] D. M. Walker, A. Tordesillas, N. Brodu, J. A. Dijksman,⁷⁵² R. P. Behringer, and G. Froyland. Self-assembly in a⁷⁵³ near-frictionless granular material: conformational struc-⁷⁵⁴ tures and transitions in uniaxial cyclic compression of hy-⁷⁵⁵ drogel spheres. *Royal Society of Chemistry*, 11(11):2157,⁷⁵⁶ 2015.

[7] R. P. Behringer and B. Chakraborty. The physics of⁷⁵⁸ jamming for granular materials: a review. *Reports on Progress in Physics*, 82(1):012601, 2018.

[8] J. Liu, A. Wautier, S. Bonelli, F. Nicot, and F. Darve.⁷⁶¹ Macroscopic softening in granular materials from a⁷⁶² mesoscale perspective. *International Journal of Solids and Structures*, 193:222–238, 2020.

[9] T. Matsushima and R. Blumenfeld. Universal structural⁷⁶⁵ characteristics of planar granular packs. *Physical review letters*, 112(9):098003, 2014.

[10] R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan,⁷⁶⁸ D. Chklovskii, and U. Alon. Network motifs: simple building blocks of complex networks. *Science*,⁷⁷⁰ 298(5594):824–827, 2002.

[11] A. Tordesillas. Force chain buckling, unjamming transitions and shear banding in dense granular assemblies. *Philosophical Magazine*, 87(32):4987–5016, 2007.

[12] D. M. Walker, A. Tordesillas, M. Small, R. P. Behringer, and C. K. Tse. A complex systems analysis of stick-slip dynamics of a laboratory fault. *Chaos: An Interdisciplinary Journal of Nonlinear Science*, 24(1):013132, 2014.

[13] K. H. Roscoe, A. Schofield, and a. P. Wroth. On the yielding of soils. *Geotechnique*, 8(1):22–53, 1958.

[14] K. Been, M. Jefferies, and J. Hachey. The critical state of sands. *Geotechnique*, 41(3):365–381, 1991.

[15] A. Casagrande. Characteristics of cohesionless soils affecting the stability of slopes and earth fills. *J. Boston Society of Civil Engineers*, 23(1):13–32, 1936.

[16] J. W. Rocks, A. J. Liu, and E. Katifori. Hidden topological structure of flow network functionality. *Physical Review Letters*, 126(2):028102, 2021.

[17] F. Nicot and F. Darve. The h-microdirectional model: accounting for a mesoscopic scale. *Mechanics of materials*, 43(12):918–929, 2011.

[18] D. M. Walker, A. Tordesillas, I. Einav, and M. Small. Complex networks in confined comminution. *Physical Review E*, 84(2):021301, 2011.

[19] O. Ben-Nun, I. Einav, and A. Tordesillas. Force attractor in confined comminution of granular materials. *Physical review letters*, 104(10):108001, 2010.

[20] A. Tordesillas, S. Kahagalage, L. Campbell, P. Bellett, E. Intrieri, and R. Batterham. Spatiotemporal slope stability analytics for failure estimation (sssafe): linking radar data to the fundamental dynamics of granular failure. *Scientific Reports*, 11(1):1–18, 2021.

[21] S. Antony. Link between single-particle properties and macroscopic properties in particulate assemblies: role of structures within structures. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and En-
gineering Sciences, 365(1861):2879-2891, 2007. 821

[22] A. Tordesillas, D. M. Walker, and Q. Lin. Force cycles and force chains. *Physical Review E*, 81(1):011302, 2010. 823

[23] A. Drescher and G. D. J. De Jong. Photoelastic verification of a mechanical model for the flow of a granular material. *Journal of the Mechanics and Physics of Solids*, 82(5):337-340, 1972. 827

[24] T. Majmudar, M. Sperl, S. Luding, and R. P. Behringer. Jamming transition in granular systems. *Physical review letters*, 98(5):058001, 2007. 830

[25] N. Deng, A. Wautier, Y. Thiery, Z.-Y. Yin, P.-Y. Hicher, 831 and F. Nicot. On the attraction power of critical state in granular materials. *Journal of the Mechanics and Physics of Solids*, 83(3):033-041, 2007. 834

[26] P. Fu and Y. F. Dafalias. Fabric evolution within shear bands of granular materials and its relation to critical state theory. *International Journal for numerical and analytical methods in geomechanics*, 35(18):1918-1948, 838, 2011. 839

[27] N. P. Kruyt and L. Rothenburg. On micromechanical characteristics of the critical state of two-dimensional granular materials. *Acta Mechanica*, 225(8):2301-2318, 842, 2014. 843

[28] L. Rothenburg and N. P. Kruyt. Critical state and evolution of coordination number in simulated granular materials. *International Journal of Solids and Structures*, 84(6):1117-1134, 2001. 846

[29] R. Wang, Y. F. Dafalias, P. Fu, and J.-M. Zhang. Fabric evolution and dilatancy within anisotropic critical state theory guided and validated by dem. *International Journal of Solids and Structures*, 188:210-222, 2020. 851

[30] H. Zhu, H. N. Nguyen, F. Nicot, and F. Darve. On a common critical state in localized and diffuse failure modes. *Journal of the Mechanics and Physics of Solids*, 95:112-131, 2016. 855

[31] M. R. Kuhn. The critical state of granular media: convergence, stationarity and disorder. *Géotechnique*, 66(11):902-909, 2016. 858

[32] M. R. Kuhn. Maximum disorder model for dense steady state flow of granular materials. *Mechanics of Material*, 93:63-80, 2016. 861

[33] A. Rechenmacher, S. Abedi, and O. Chupin. Evolution of force chains in shear bands in sands. *Géotechnique*, 60(5):343-351, 2010. 864

[34] X. Sun, Y. Wang, Y. Wang, R. Blumenfeld, and J. Zhang. Experimental evidence of detailed balance in granular systems. *arXiv preprint arXiv:2105.01355*, 2021. 866

[35] A. Wautier, S. Bonelli and F. Nicot Rattlers' contribution to granular plasticity and mechanical stability. *Interna* 869

tional Journal of Plasticity, 2019.

[36] N. Deng, A. Wautier, Y. Thiery, Z.-Y. Yin, P.-Y. Hicher, and F. Nicot. Dynamic view of critical state regime in granular materials: a mesoscale perspective. In *EPJ Web of Conferences*, volume 249, page 11011. EDP Sciences, 2021.

[37] A. Wautier, S. Bonelli, and F. Nicot. Flow impact on granular force chains and induced instability. *Physical Review E*, 98(4):042909, 2018.

[38] P. A. Cundall and O. D. Strack. A discrete numerical model for granular assemblies. *geotechnique*, 29(1):47-65, 1979.

[39] V. Šmilauer, E. Catalano, B. Chareyre, S. Dorofeenko, J. Duriez, A. Gladky, J. Kozicki, C. Modenese, L. Scholtès, L. Sibille, et al. Yade documentation 2nd ed. the yade project. *Transp. Porous Media*, 2015.

[40] C. O'Sullivan, J. D. Bray, and S. Li. A new approach for calculating strain for particulate media. *International Journal for numerical and analytical methods in geomechanics*, 27(10):859-877, 2003.

[41] J. Liu, F. Nicot, and W. Zhou. Sustainability of internal structures during shear band forming in 2d granular materials. *Powder technology*, 338:458-470, 2018.

[42] F. Radjai, M. Jean, J.-J. Moreau, and S. Roux. Force distributions in dense two-dimensional granular systems. *Physical review letters*, 77(2):274, 1996.

[43] J. Peters, M. Muthuswamy, J. Wibowo, and A. Tordesillas. Characterization of force chains in granular material. *Physical review E*, 72(4):041307, 2005.

[44] A. Wautier, S. Bonelli, and F. Nicot. Micro-inertia origin of instabilities in granular materials. *International Journal for Numerical and Analytical Methods in Geomechanics*, 42(9):1037-1056, 2018.

[45] P. Schall, D. A. Weitz, and F. Spaepen. Structural rearrangements that govern flow in colloidal glasses. *Science*, 318(5858):1895-1899, 2007.

[46] F. Darve, F. Nicot, A. Wautier, and J. Liu. Slip lines versus shear bands: two competing localization modes. *Mechanics Research Communications*, page 103603, 2020.

[47] M. H. Sadd, G. Adhikari, and F. Cardoso. Dem simulation of wave propagation in granular materials. *Powder Technology*, 109(1-3):222-233, 2000.

[48] A. Wautier, G. Veylon, M. Miot, M. Pouragha, F. Nicot, R. Wan, and F. Darve. Multiscale modelling of granular materials in boundary value problems accounting for mesoscale mechanisms. *Computers and Geotechnics*, 134:104143, 2021.

[49] J. Chueire, A. Wautier, F. Nicot, and A. Daouadji. On the extension of the grain loop concept from 2d to 3d granular assemblies. (Under review)