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

HAL
open science

# Stability and fluctuations in granular materials during the critical state regime

Xiaoxiao Wang, Yang Liu, A. Wautier, François Nicot

► To cite this version:

Xiaoxiao Wang, Yang Liu, A. Wautier, François Nicot. Stability and fluctuations in granular materials during the critical state regime. Computers and Geotechnics, 2025, 183, pp.107201. 10.1016/j.compgeo.2025.107201. hal-04996493

HAL Id: hal-04996493

https://hal.inrae.fr/hal-04996493v1

Submitted on 19 Mar 2025

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Stability and fluctuations in granular materials during the critical state regime

Xiaoxiao WANG¹, Yang LIU¹*, Antoine WAUTIER², François NICOT³

(1) Department of Civil Engineering, University of Science and Technology Beijing, Beijing, 100083, China

(2) INRAE, Aix-Marseille Université, RECOVER, Aix-en-Provence, France

(3) ISTerre, Université Savoie Mont Blanc, Le Bourget du Lac, France

Corresponding author: Yang LIU (yangliu@ustb.edu.cn)

## ABSTRACT:

Referring to small strain continuum mechanics, a material is unstable if a spontaneous release in kinetic energy can occur following a small perturbation, consistently to the second-order work Hill's instability criterion. Based on a 2-D discrete element method (DEM), this manuscript investigates the stability of granular assemblies in relation with the elastic energy processes along a biaxial loading. The evolution of relevant stability indexes shows that the strong degree of instability very close to the stress peak is linked with the limited capacity of storing elastic energy, resulting in the bifurcation from unstable state to stable state with the formation of one or multiple shear bands. In the post-bifurcation stage, the significant release in net elastic energy within a dense specimen brings the system in a more stable state where a fluctuation regime develops. On the other hand, the loose specimen undergoes stability fluctuation from the very beginning up to the critical state of biaxial loading. From the perspective of configurational mechanisms, the stability fluctuation stems from loops transformations entailing repeated energy storing and release. Such changes in microstructural topology continuously adjust the mesostructures under external loading to approach a stable bifurcated branch in the form of stability fluctuations.

**KEYWORDS:** Granular materials, Critical state, Second-order work theory, Stability and fluctuation, Elastic energy, Grain loops

## 1 Introduction

Granular materials exhibit a wide spectrum of emergent properties under a given loading even though the interaction between contacting granules can be described by very simple contact
laws. This drives the complexity of granular assemblies, and has attracted a marked interest in the past years (Jin et al., 2024; Kruyt, 1999; Kruyt and Rothenburg, 2006; Liu and Wang, 2024; Nicot et al., 2023; Prévost and Höeg, 1975; Yan et al., 2024; Yu and Liu, 2022). Due to the discrete nature of granular materials, continuous rearrangement occurs within the grain assembly under shearing, which is the source of complexity. For a given microstructural configuration, the system is reputed stable if any small incremental loading is accommodated by small incremental responses at every point of the system. In some configurations, a small perturbation may be amplified and result in finite release of kinetic energy. It is essential for many engineering purposes to characterize the stability of granular systems as they are constantly subjected to diverse perturbations.

From a theoretical point of view, the so-called second-order work criterion as the pioneering work of Hill (1958) applies for stability analysis of rate-independent materials (Lopera Perez et al., 2016; Nicot et al., 2017; Yan et al., 2024; Zhao et al., 2018). A sufficient condition of stability for a given mechanical state of a material system refers to a strictly positive second-order work for any incremental loading directions. On the other hand, the vanishing of the second-order work is a necessary condition for the existence of instability. Using for example a discrete element method (DEM) (Cundall and Strack, 1979), a directional analysis, where a series of probes are imposed along all the directions of the loading space, can be carried to analyze the stability of a given mechanical state. Prior results (Nicot et al., 2011; Wang et al., 2021) have illustrated the possible existence of instability cones, for a given mechanical state, gathering all loading directions (i.e., the class of directional perturbations) for which the second-order work is negative. More recently, the second-order work theory was formalized further by linking the negative values of second-order work to a burst of kinetic energy, corresponding to a transition from a quasi-static
regime toward a dynamic regime (Clerc et al., 2021; Daouadji et al., 2011; Nguyen et al., 2016; Nicot et al., 2007; Nicot et al., 2012; Wan et al., 2013; Wautier et al., 2018). Furthermore, the set of all potentially unstable states for a given material constitutes the bifurcation domain inside the stress space (Darve et al., 2007; Neilsen and Schreyer, 1993; Veylon et al., 2018; Wan, 2011), in which an effective failure of the material can occur according to the loading direction applied, and to the loading control adopted.

It is well-known that granular materials evolve towards a critical state following a strain-softening or hardening regime, in which the shear strain continuously progresses without any change in both stress and volume (Li and Wang, 1998; Roscoe, 1958; Schofield and Wroth, 1968). Using photo elastic grains (Drescher and De Jong, 2006; Tordesillas et al., 2009; Zhang et al., 2008) and numerical simulations (Bardenhagen et al., 2000; Liu et al., 2018; Nguyen et al., 2009; Liu et al., 2023), the meso-structures such as force chains and meso-loops have proven to be relevant to characterize the force transmission and geometrical rearrangement within the granular assemblies under shearing. Using a multiscale approach, prior results (Liu and Li, 2018; Liu et al., 2020b; Nguyen et al., 2009; Zhu et al., 2016a) have shown that the statistical descriptors of meso-structures converge to a stationary state corresponding to the critical state regime. Deng et al. (2022) and Wang et al. (2024) unraveled an underlying dynamic equilibrium behind this stationary regime, which involves contact-generating and contact-vanishing processes entailing energy storage and energy release at contact scale. The dynamic processes at the critical state regime give rise to the configurational reactions that can be regarded as local fluctuations (Liu et al., 2018; Liu et al., 2020a; Nicot et al., 2024). As a result, the impact of local fluctuations on the mechanical stability at critical state remains an open issue.

From a microscopic energy point of view, it has been shown that dense materials hold a
74 maximum capacity for storing elastic energy at peak stress state (Nicot et al., 2023; Wang et al., 2024), giving rise to a bifurcation commonly accompanied by the development of shear bands (Chen et al., 2022; Guo and Zhao, 2016; Oda and Kazama, 1998; Rattez et al., 2022). Then, for a dense assembly to accommodate the external loading, the material releases energy which is dissipated by contact sliding. As a result, the bifurcated response mode of dense systems is characterized by a stationary elastic-to-plastic energy transfer at macroscopic critical state. This stationary regime is regarded as an optimal dissipative regime where specific structures develop by concentrating all the dissipation processes inside the ultimate shear band (Wang et al., 2024).
81 82 Although the stability of granular assemblies has been extensively studied in the last decades, the following open issues remain to be investigated: (i) Does the system evolve to a more ‘stable’ or ‘unstable’ state driven by the elastic energy processes? (ii) What is the contribution of the configurational fluctuations to the overall stability of granular assemblies?

86 This manuscript tackles these open issues in the framework of second-order work theory.
87 The DEM simulation of the macroscopic response of granular assemblies with different initial densities subjected to a 2D biaxial loading is detailed in Section 2. The formation of the shear band 88 along a biaxial loading, in relation with the evolution of elastic energy processes, is discussed in 89 Section 3. After revisiting the classical strain controlled directional analysis procedure in the 90 context of the second-order work theory, the stability of granular assemblies is investigated in 91 Section 4 by introducing two stability indexes. Finally, Section 5 discusses the underlying 92 configurational reorganization mechanisms. 93

## 2 DEM simulation of a biaxial loading path

94 In this study, the software PFC2D (2D Particle Flow Code; Itasca Consulting Group Inc., 95 2008) is used to simulate the drained biaxial loadings under constant lateral pressure based on the
discrete element method (DEM) proposed by Cundall and Strack (1979). Down to the microscale, the simple linear elasto-frictional contact law as illustrated in Fig. 1 is adopted to model the interaction between two particles. The normal and tangential contact forces $F_n$ and $F_t$ in the contact model considered are computed as:

$$ \left\{ \begin{array}{l} F_n = k_n u_n \\ dF_t = k_t du_t, \quad F_t \le F_n \tan \phi \end{array} \right. \qquad (1) $$

where $k_n$ and $k_t$ are the normal and tangential stiffness, respectively. $u_n$ and $du_t$ are the total normal relative displacement and incremental tangential displacement between two contacting particles, respectively. $\phi$ is the friction angle involved in the Mohr–Coulomb criterion to limit the tangential contact force $F_t$.

<figure><img src="image_3.png" /><figcaption>Fig. 1. Elasto-frictional contact law used in DEM simulations.</figcaption></figure>

Two specimens are generated by containing 18,000 idealized disc particles with the default thickness of 1.0 m. The density of particles is 2600 kg/m³. The particles with an average diameter $d_{50} = 0.009$ m and $d_{max} / d_{min} = 2$ are installed within a 1.0 m (width)×1.5 m (height) rectangle domain, as illustrated in Fig. 2. For the parameters of the contact model considered, the friction coefficient of 0.0 and to 0.5 are used for the two specimens during the preparation process, respectively, in order to obtain different initial densities. $k_n$ and $k_t$ are set to $5 \times 10^8$ N/m for
particles and walls. Note that, a very small value of the damping coefficient is considered in order to limit numerical dissipation in the system, which is set to 0.01 in this study.

<figure><img src="image_4.png" /><figcaption>Fig. 2. DEM specimen under a biaxial loading and 2D loop tessellation in the granular assemblies.</figcaption></figure>

<table><caption>Table 1. Numerical parameters used for the simulation.</caption><thead><tr><td>Parameters</td><td>Numerical values</td></tr></thead><tbody><tr><td>Initial size of specimen (m)</td><td>1.5 (height) × 1.0 (width)</td></tr><tr><td>Average radii of particles D<sub>50</sub> (m)</td><td>0.009</td></tr><tr><td>Density of particles (kg/m<sup>3</sup>)</td><td>2600</td></tr><tr><td>Contact normal stiffness k<sub>n</sub> (N/m)</td><td>5 × 10<sup>8</sup></td></tr><tr><td>Contact tangential stiffness k<sub>t</sub> (N/m)</td><td>5 × 10<sup>8</sup></td></tr><tr><td>Interparticle friction coefficient</td><td>0.5</td></tr><tr><td>Damping coefficient</td><td>0.01</td></tr><tr><td>Confining pressure (kPa)</td><td>400</td></tr><tr><td>Loading rate (s<sup>-1</sup>)</td><td>0.01</td></tr></tbody></table>

A confining pressure of 400 kPa is considered to consolidate the specimens in an isotropic state. After the confining stage, a relative dense state with an initial porosity of 0.15 is achieved for the numerical specimen with zero friction, while another specimen with a friction coefficient of 0.5 reaches a relative loose state, with an initial porosity of 0.19. Then, a biaxial loading is applied to specimens after resetting the coefficient of friction to 0.5, i.e., imposing a compression in the vertical direction ($\sigma_1$ and $\epsilon_1$) with a constant loading rate of 0.01 s$^{-1}$, and keeping the stress
in horizontal direction ($\sigma_2$ and $\varepsilon_2$) constant to 400 kPa, as shown in Fig. 2. Note that, the Inertial Number $I$ of the granular system ranges from $5.5 \times 10^{-6}$ to $7.3 \times 10^{-6}$ as the pressure varies from 400 kPa to 700 kPa, suggesting that the system remains in a quasi-static state (Da Cruz et al., 2005; Deng et al., 2021). All the numerical parameters used for the biaxial loading simulation are summarized in Table 1.

The macroscopic responses involving the evolutions of deviatoric stress $q$ ($q = \sigma_1 - \sigma_2$) and volumetric strain $\varepsilon_v$ ($\varepsilon_v = \varepsilon_1 + \varepsilon_2$) along the biaxial loading for dense and loose specimens are illustrated in Fig. 3 and Fig. 4, respectively, where $\sigma_1$ and $\sigma_2$ represents principal stress components, and $\varepsilon_1$ and $\varepsilon_2$ are principal strain components. As observed in Fig. 3(a) and Fig. 4(a), the deviatoric stress $q$ exhibits typical strain-softening behavior for dense specimen and strain-hardening behavior for loose specimen. In terms of the volumetric strain $\varepsilon_v$, the dense specimen first contracts and then dilates as shown in Fig 3(b), while the loose specimen contracts continually as observed in Fig 4(b). It is worth noting that, the critical state is obtained at large axial stains where both the evolutions of deviatoric stress $q$ and volumetric strain $\varepsilon_v$ fluctuate around residual value, no matter the initial densities of the specimen considered.

According to the macroscopic mechanical response of dense specimen, four characteristic states are emphasized, named state A-D referring to the maximal contractancy state, peak stress state, end of softening stage, and critical state, respectively. Additionally, with reference to these characteristic states, the whole biaxial loading path is divided into four stages, corresponding to the initial loading stage I, the strain softening stage II and the critical state stage III-IV for dense specimen (in Fig. 3), and the strain hardening stage I-II and the critical state stage III-IV for loose specimen. In order to analyze the stability evolution of granular assemblies, a sufficient number of
mechanical states are saved along the biaxial loading paths as marked by black crosses in Fig. 3 and Fig. 4. For both dense and loose specimens, these mechanical states are selected at intervals of 0.5% axial strain along the biaxial loading path, which ensures that there are enough mechanical states considered in each loading stage. Note that, several additional mechanical states in the proximity of the peak stress state are saved to analyze the specific changes in mechanical stability.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Evolution of deviatoric stress $q$ (a) and volumetric strain $\varepsilon_v$ (b) with respect to axial strain $\varepsilon_1$ along a biaxial loading path for dense specimen.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 4. Evolution of deviatoric stress $q$ (a) and volumetric strain $\varepsilon_v$ (b) with respect to axial strain $\varepsilon_1$ along a biaxial loading path for loose specimen.</figcaption></figure>
161

# 3 Shear band initiation and the elastic energy evolution

162

## 3.1 Shear band formation

163

2D granular assemblies can be tessellated into loops as shown in Fig. 2. Note that, particles with only one contact or no contact are excluded from the tessellation, as they are not involved in supporting external loading. In this study, meso-loops are divided into different categories according to the number of grains, which plays an important role in describing the configurational evolution within the granular assemblies. The loop structure with $n$ branches is denoted as “Ln”, where $n$ is the order (or valence) of the loop. According to Zhu et al. (Zhu et al., 2016a; 2016b), loops with higher order ($n \ge 6$) hold similar structural evolution tendency along biaxial loading path, and will be categorized as the L6+ in this study.

164

165

166

167

168

169

170

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

171

172

<caption>Fig. 5. Mesoscopic strain definition of a loop L.</caption>

173

174

The definition of mesoscopic strain proposed by Liu et al. (2020a; 2022) is adopted, and the incremental strain tensor of loop L is expressed as follows:

175

$$d\epsilon_{ij}^{L} = \frac{1}{V^L} \sum_{k=1}^{c} \frac{(du_{i}^{0k} + du_{i}^{1k})}{2} n_{j}^{k} l^{k} \quad (2)$$

176

177

where $du_{i}^{0k}$ and $du_{i}^{1k}$ are incremental displacement in the direction $i$ of the two particles connected by the $k$th edge; $l^k$ is the length of the $k$th edge; $n^k$ is an outward normal vector to the $k$th edge;
178

$V^L$ is the area of the loop $L$ referring to the polygon enclosed by the blue boundary in Fig. 5.

179

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

180

<caption>Fig. 6. Incremental deviatoric strain field for dense specimen at different loading states. (a) State A: $\epsilon_1=0.41\%$; (b) State B: $\epsilon_1=0.95\%$; (c) State C: $\epsilon_1=4.0\%$; (d) State D: $\epsilon_1=8.0\%$.</caption>

181

182

183

The incremental deviatoric strain field for each meso-loop, referring to a small axial strain increment of 0.01% at the four key mechanical states considered (A, B, C and D), are shown in Fig. 6. The development of local deformation patterns in the dense specimen is clearly observed. Along the biaxial loading path, the deformation of the dense specimen transitions from initial homogeneity (Fig. 6a) to increasing heterogeneity (Fig. 6b–d). A single shear band begins to emerge at peak stress state and becomes more evident during the critical state regime (Wan et al., 2017).

184

185

186

187

188

189

190

## 3.2 Elastic energy processes

191

As a biaxial loading is incrementally applied, the incremental external work $dW^{ext}$ ($dW^{ext} = V\sigma_i d\epsilon_i$, where $V$ is the current volume of specimen) provided to the granular system at each numerical time step remains positive. The evolution of the external work $W^{ext}$ with respect to the axial strain for dense and loose specimens is shown in Fig. 7(a) and Fig. 8(a), respectively. During the biaxial loading, a fraction of $W^{ext}$ is stored at contacts in the form of elastic energy.

192

193

194

195
Within the granular assembly, the total elastic energy $E^{el}$ at a certain loading state is computed by summing up the elastic energy stored at each contact $c$, which includes a normal part and a tangential counterpart as follows:

$$E^{el} = \frac{1}{2} \sum_{c=1}^{N_c} \left( (F_n^c)^2 / k_n + (F_t^c)^2 / k_t \right) \quad (3)$$
where $N_c$ is the total number of contacts.

The evolution of total elastic energy $E^{el}$ with respect to the axial strain for dense specimen is shown in Fig. 7(b). It is observed that the dense specimen holds a limited capacity of storing elastic energy after experiencing a sharp increase in $E^{el}$ at the beginning of the loading. And interestingly, the limit in elastic energy storage is reached very close to the peak stress (state B). As broadly discussed over the past years, the peak stress state is unstable, according to the second-order work theory (Nicot et al., 2012; Ord et al., 2007; Tordesillas, 2007; Wan et al., 2013). Simultaneously, the shear band appears, corresponding to a break in symmetry within the dense specimen. This break in symmetry is characterized by the spatial bifurcation (or phase transition) from a homogeneous strain field to a heterogeneous one. Then, when the loading program is pursued beyond state B, the system begins to release a fraction of the stored elastic energy, corresponding to the macroscopic strain softening and the development of the shear band. At state C (end of macroscopic strain-softening phase), the elastic energy in the dense specimen approaches a stationary regime with fluctuations around a non-zero value. The specimen reaches the macroscopic critical state as the same time as the shear band is steadily structured. In addition, the level of elastic energy $E^{el}$ stored in the loose specimen, as plotted in Fig. 8(b), is much lower than that for the dense specimen. The capacity of storing elastic energy increases continuously, and also
approaches to a stationary state regime in a very smooth way, accompanied with the macroscopic

strain-hardening regime and no occurrence of strain localization.

<figure><img src="image_9.png" /><figcaption>Fig. 7. Evolution of the external work (a) and of the elastic energy stored within dense specimen (b) along the biaxial loading path.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 8. Evolution of the external work (a) and of the elastic energy stored within loose specimen (b) along the biaxial loading path.</figcaption></figure>

# 4 Stability analysis for granular assemblies along a biaxial loading path

In this study, the standard strain-controlled directional analysis procedure is used to investigate the stability of granular assemblies in the framework of the second-order work criterion.
The numerical procedure is composed of two main steps including a pre-stabilization step and an incremental strain probe, which are introduced in detail hereafter. Based on the transient mechanical response during the incremental strain probe procedure, two instability indexes are defined, namely the probability of instability and the degree of stability. Then, using these indexes, the stability variations with the development of elastic energy process for dense and loose specimens are analyzed.

## 4.1 Strain controlled directional analysis

### 4.1.1 Pre-stabilization step

Along the biaxial loading path, the mechanical response of specimen is inherently dynamic as it depends on particle rearrangement obtained from the integration of Newton's second law of motion for each particle. Bursts of fast reorganization always occur, no matter how small is set the axial loading rate to ensure quasi-static loading conditions. As a result, a pre-stabilization step is first achieved by keeping the current lateral pressure and axial stress constant and equal to the stress level of the saved mechanical state considered. The equilibrium state is assumed to be reached once the unbalanced force ratio verifies the condition $F_{unb} < 10^{-5}$, where $F_{unb}$ is the ratio between the mean summation force on the $N_p$ particles and the mean contact force magnitude on $N_p$ particles,

$$F_{unb} = \frac{\frac{1}{N_p} \sum_{p=1}^{N_p} \left\| \sum_{F/p} F \right\|}{\frac{1}{N_p} \sum_{p=1}^{N_p} \left\| F \right\|} \quad (4)$$

where, $N_p$ is the number of particles within the granular assembly considered; $(F/p)$ represents the force $F$ acting on particle $p$.
During the pre-stabilization step, the particle rearrangement is still possible to result in volumetric deformation even if the stress state is kept constant. The stress-strain states at the end of the pre-stabilization step are presented in Fig. 3 and Fig. 4 by blue hollow diamonds. Compared with the virgin stress-strain states marked by black crosses in Fig. 3 and Fig. 4, negligible delayed deformations are observed for most of the mechanical states considered. However, as shown in the zoomed insets of the volumetric strain curve in Fig. 3(b) and Fig.4(b), one can observe that non-negligible delayed deformations larger than 0.1% occur at $\epsilon_i \in \{4.0, 5.0, 8.0, 13.5\}\%$ for the dense specimen and at $\epsilon_i \in \{4.5, 11.5\}\%$ for the loose specimen. For these mechanical states, the quasi-static assumption becomes questionable. Grain rearrangements continue to occur during the pre-stabilization step which result in delayed permanent deformations at the specimen scale. A parallel can be drawn here with the concept of “delayed plasticity” introduced by Di Prisco (1996), which clearly shows that the time dependency of strain response is associated with a delayed microstructural evolution, as also mentioned in Wautier et al. (2018). The delayed deformation will affect the mechanical stability since the state variables obtained after the pre-stabilization step are different from the initial ones.

### 4.1.2 Directional analysis

After the pre-stabilization step, a strain increment $d\epsilon$ with the orientation angle $\theta_e$ is imposed to the specimen. The directional analysis is achieved by varying the orientation angle $\theta_e$ from 0° to 360° with a step of 10°, while keeping the magnitude of $d\epsilon$ constant for each incremental strain loading as shown in Fig. 9(a). Note that, all the strain loading increments are applied to the specimen at the end of pre-stabilization step with the same strain loading rate.
<figure><img src="image_11.png" /><figcaption>Fig. 9. Strain probe (a) and stress response (b) in the axisymmetric planes of strain and stress increments respectively.</figcaption></figure>

271

272

273

274

275

Along each direction, the strain increments in axial direction $d\epsilon_1$ and lateral direction $d\epsilon_2$

276

are expressed as:

$$ \left\{ \begin{array}{l} d\epsilon_1 = \|d\epsilon\| \sin \theta_{\epsilon} \\ d\epsilon_2 = \|d\epsilon\| \cos \theta_{\epsilon} \\ \|d\epsilon\| = \sqrt{d\epsilon_1^2 + d\epsilon_2^2} \end{array} \right. \qquad (5) $$

277

278

279

280

At the end of each imposed strain increment $d\epsilon$, a stress increment $d\sigma$ is obtained with an orientation angle $\theta_{\sigma}$ within the incremental stress space, as illustrated in Fig. 9(b). Then, the normalized second-order work $W_2^{norm}$ can be defined at the global specimen scale as:

$$ W_2^{norm} = \frac{d\epsilon : d\sigma}{\|d\epsilon\|\|d\sigma\|} \qquad (6) $$

281

282

where, the range of the normalized second-order work is $-1.0 \le W_2^{norm} \le 1.0$.

283

284

285

286

Based on the second-order work criterion, and excluding flutter instabilities, the mechanical state is reputed stable if the second-order work is strictly positive along all loading directions (Nicot et al., 2011). When the second-order work vanishes or becomes negative along some incremental loading directions, the mechanical state is reputed unstable.
The magnitude of incremental strain $||d\varepsilon||$ is an important factor that strongly affects the incremental response of the specimen. Indeed, slight boundary perturbation characterized by a too small value of $||d\varepsilon||$ cannot spread to the whole specimen. Indeed, the specimen has a finite characteristic length that requires a large enough magnitude of $||d\varepsilon||$ to make the microstructure significantly reorganize in response to boundary incremental loading. To account for this scale effect, the influence of incremental strain magnitude is discussed specifically in Appendix A, by examining the evolution of the second-order work. According to this discussion, setting the size of incremental strain $||d\varepsilon||$ as 0.05% is thought to be appropriate: the incremental strain is appropriately large to propagate the incremental perturbation throughout the whole specimen.

## 4.2 Instability cone and degree of stability

Figure 10 presents the variation of second-order work $W_2$ in terms of the strain probe direction $\theta_e$ for the dense specimen at mechanical states A and B. Note that, state A is a stable state, as the second-order work takes strictly positive values along all directions, while state B is an unstable state where the second-order work vanishes along some loading directions (i.e., $W_2 \le 0$). For the unstable state B, an instability cone emerges which gathers all potentially unstable loading directions with $W_2 \le 0$, as marked by the pink zone in Fig. 10. The size $\alpha$ of the instability cone denotes the amplitude of the directional interval in which the specimen is likely to destabilize when such an incremental loading is applied. Therefore, the ratio of $\alpha$ over the total amplitude of the incremental loading direction interval ($2\pi$) is defined as the probability of instability $P^{\text{uns}}$ to observe a potential instability triggered off when random perturbations are likely to occur.
<figure><img src="image_12.png" /><figcaption>Fig. 10. Variation of second-order work in terms of the strain probe direction for dense specimen at mechanical states A and B. The ratio between the amplitude of the instability cone (pink zone) and the whole directional range (360°) corresponds to the probability of instability (P<sup>uns</sup>), while the minimum value (star point) corresponds to the degree of stability (D<sup>s</sup>).</figcaption></figure>

In addition, even though the values of P<sup>uns</sup> for the specimens at two different mechanical states are the same, the amplitude of the second-order work inside the instability cone may be different. The physical meaning of the second-order work vanishing is clarified in terms of a second-order energy balance (Nicot et al., 2007; Nicot et al., 2011; Wan et al., 2017). Considering a given mechanical system originally at equilibrium at time t (the kinetic energy E<sub>c</sub>(t) = 0), the variation of kinetic energy (ΔE<sub>c</sub> = E<sub>c</sub>(t + Δt) - E<sub>c</sub>(t)) over a small time range [t, t + Δt] equals to the difference between the external second-order work W<sub>2</sub><sup>ext</sup> and the internal second-order work W<sub>2</sub><sup>int</sup>:

$$ \Delta E_c = \frac{1}{2} (W_2^{\text{ext}} - W_2^{\text{int}}) = \frac{1}{2} \left( \int_{\partial\Omega} df \cdot du \, dS - \int_{\Omega} d\sigma : d\varepsilon dV \right) \quad (7) $$

where the external second-order work W<sub>2</sub><sup>ext</sup> is computed from the conjugate incremental force df and incremental displacement du acting on the boundary ∂Ω of the system Ω; the internal
second-order work corresponds to the standard second-order work $W_2$ as introduced by Hill (1958), and recalled in Eq. (6) where no incremental geometrical changes are accounted for.

According to Eq. (7), an abruptly increase in kinetic energy ($\Delta E_c > 0$) from an equilibrium configuration is related to the competition between the external loading applied to the boundary of the system and the constitutive behavior of the material. In the case of a nil external second-order work (e.g., the loading applied to the system is assigned to remain constant), the vanishing of the second-order work automatically results in an outburst in the kinetic energy of the system, corresponding to a transition from a quasi-static regime to an inertial regime (Clerc et al., 2021; Wautier et al., 2018).

Based on Eq. (7), it transpires that a larger amplitude of the second-order work along a direction inside the instability cone entails a more significant outburst in kinetic energy. Then, another complementary index $D^s$ (degree of stability) can be defined by referring to the minimum value of second-order work:

$$D^s = (W_2)_{\min} \qquad (8)$$

where $(W_2)_{\min}$ is the minimum of the second-order work, as marked by the black star point in Fig. 10. For the unstable state, the smaller the value of $D^s$, the stronger the degree of instability of the granular specimen. Conversely, for the stable state, the larger the value of $D^s$, the stronger the degree of stability of the granular specimen.

## 4.3 Fluctuating stability regime along biaxial loading

The variations of instability probability $P^{uns}$ and degree of stability $D^s$ for the dense specimen at different mechanical states along the biaxial loading path are plotted in Fig. 11. Considering the role of the elastic energy in the occurrence of a spatial bifurcation, as discussed in
Section 3.2, the elastic energy $E^{el}$ stored within the granular specimens along the biaxial loading is recalled in Fig. 11.

As illustrated in Fig. 11(a), in the initial loading stage I (before the peak stress state B), the probability of instability $P^{uns}$ of dense specimen at these considered mechanical states is mostly zero (no instability cones). However, with the sharp increase in elastic energy storage $E^{el}$, the degree of stability $D^s$ gradually decreases as shown in Fig. 11(b). For the considered mechanical states, the second-order work vanishing characterized by the non-zero $P^{uns}$ and the negative value of $D^s$ first occurs before the stress peak. Note that, the first vanishing of the second-order work corresponds to the boundary of bifurcation domain. Prior studies (Darve et al., 2007; Neilsen and Schreyer, 1993; Prunier et al., 2009; Wan, 2011) have repeatedly shown that this boundary is clearly within the plastic limit surface. As a limited capacity of elastic energy storage is achieved after but very close to the peak state B, the probability of instability $P^{uns}$ increases to a local maximum and the degree of stability $D^s$ decreases to a local minimum. This demonstrates a resulting high instability level at the maximum of elastic energy storage within the dense specimen. After state B, with the release of elastic energy, the unstable dense assembly breaks down, resulting in macroscopic strain-softening and the development of a shear band. Correspondingly, the probability of instability $P^{uns}$ of dense specimen decreases to zero at $\epsilon_1 = 2\%$, combined with the increase in $D^s$, which suggests that the system bifurcates from an unstable mechanical state to a stable one. It is worth noting that the elastic energy evolves in the form of fluctuation with a gradual increase in the magnitude of the fluctuations, especially at the critical state. Accompanying the observation is the fact that the degree of stability $D^s$ displays a fluctuating trend around zero value during stages II~IV. Noticeably, this fluctuating trend of $D^s$ gradually increases after state C, in the critical state regime. As a result, the dense specimen undergoes a fluctuation regime of stability
368 immediately after state B, which corresponds exactly to the starting point of the fluctuation

369 evolution of elastic energy storage.

<figure><img src="image_13.png" /><figcaption>Fig. 11. Evolution of probability of instability $P^{\text{uns}}$ (a) and degree of stability $D^s$ (b) for dense specimen at different mechanical states along the biaxial loading path.</figcaption></figure>

The variations of instability probability $P^{\text{uns}}$ and degree of stability $D^s$ for the loose specimen at different mechanical states along the biaxial loading path are illustrated in Fig. 12.
Unlike the dense specimen, the instability probability repeatedly varies between $P^{\text{uns}} = 0$ and $P^{\text{uns}} > 0$, while a continuous increase in elastic energy storage can be observed. Likewise, the
degree of stability $D^s$ fluctuates around zero value from the very beginning to the end of biaxial loading. This suggests that the loose specimen undergoes a fluctuation regime of stability throughout the loading process. Similar with what was observed for the dense specimen, the fluctuation of stability tends to be more significant after state C as the loose specimen approaches to the critical state regime. Correspondingly, the fluctuating evolution of elastic energy for loose specimen starts from the very beginning of loading, and the magnitude of fluctuations also gradually increases after state C.

<figure><img src="image_14.png" /><figcaption>Fig. 12. Evolution of probability of instability $P^{\text{uns}}$ (a) and degree of stability $D^s$ (b) for loose specimen at different mechanical states along the biaxial loading path.</figcaption></figure>
At the light of these results, two types of stability evolution can be summarized for the
granular assemblies along a biaxial loading path:

(i) A stability stage, which arises mainly in the initial loading stage of dense specimen. In this stage, the mechanical states considered are stable with gradually decreasing stability degree, governed by the smoothly increasing elastic energy storage within the dense granular assembly.

(ii) A stability fluctuation regime is related to the fluctuating evolution of elastic energy for the dense granular assembly after peak stress state and from the very beginning of the loading for the loose granular assembly. The magnitude of the stability fluctuation tends to be more significant with the increase of the magnitude of elastic energy fluctuation, whatever the initial density of the specimen considered.

It is worth noting that the fluctuating characteristic of the elastic energy processes shows a close correlation with the fluctuating evolution of macroscopic deviatoric stress *q* curve (in Fig. 3(a) and 4(a)). Fig. 13 zooms up an example of elastic energy fluctuations (red curve) and stress fluctuations (black curve) for the dense specimen at several mechanical states belonging to the critical state. It can be observed that the local fluctuations between stress and elastic energy are almost synchronized, e.g., they synchronously reach the corresponding local peaks or valleys. Similar results can be observed at different mechanical states considered, as shown in Fig. B1 and B2 in Appendix B. The local fluctuations of elastic energy indicate that the elastic energy within the system is constantly stored and released, although the elastic energy is released as a whole during the softening stage and reaches a stationary state at the critical state regime. The storage and release of elastic energy cause the increase and decrease of stress, i.e., the stress fluctuation. This fluctuation regime starts from the peak stress state for the dense specimen and from the initial loading for the loose specimen and becomes well marked in the critical state regime.
<figure><img src="image_15.png" /><figcaption>Fig. 13. Zoomed picture of elastic energy fluctuations (the red curve) and stress fluctuations (the black curve) for dense specimen at mechanical states belonging to the critical state regime.</figcaption></figure>

Considering the synchronization between the local increase/decrease in the stress and the storage/release of elastic energy, the correlation between the location of a mechanical state on the stress curve and the elastic energy regime (storage or release) at this considered state is worth discussing. The locations of these mechanical states considered in this study on the zoomed *q*-curve, as well as the probability of instability *P*<sup>uns</sup>, are reported in Appendix B. Contrary to elastic energy fluctuation, it should be highlighted that the stability of the mechanical state considered is not sensitive to its position along the *q*-curve. The mechanical states corresponding to the local stress peaks where the elastic energy storage increases to a local maximum, are not always detected as unstable, as observed for ε₁ ∈ {3.5%, 4.0%} in Fig. B1 and ε₁ ∈ {1.5%, 3.0%, 8.5%, 13.5%} in Fig. B2. Likewise, the mechanical states corresponding to the local stress valleys, where the elastic energy storage decreases to a local minimum, are not always detected as stable, as observed for ε₁ ∈ {1.1%, 14.5%} in Fig. B1. This finding holds even when no delayed deformations (responsible for important microstructure modify) are observed.
5 Meso-loops reorganization mechanisms for stability fluctuation

Prior studies (Clerc et al., 2021; Liu et al., 2020a; Tordesillas et al., 2011; Wanjura et al., 2020) have shown that the critical state regime of granular systems under shear loading undergoes repeated dynamic configurational transitions characterized by the destruction-construction of meso-loops. This was shown to be a key ingredient responsible for the fluctuating evolutions of macroscopic responses (Nicot et al., 2024). In this sense, this stability fluctuation phenomenon probably roots into the dynamic rearrangement of microstructure within the granular assemblies.

Dynamic configurational transitions are accompanied by processes of contact generation and disappearance, entailing dynamic elastic energy processes. In this sense, when focusing on a loading increment from time t to time t + Δt, the contribution of a contact c storing or releasing energy can be quantified by the positive or negative increment of elastic energy $\Delta E_c^{el}$, $\Delta E_c^{el} = (E_c^{el})^{t+\Delta t} - (E_c^{el})^t$. Then, the amount of stored elastic energy ($\Delta E_c^{el}$)$^{store}$ and released elastic energy ($\Delta E_c^{el}$)$^{release}$ during a loading increment can be obtained by counting all contacts storing energy and all contacts releasing energy. Correspondingly, the stored power $P^{store}$ and released power $P^{release}$ can be computed as the stored energy and released energy per unit time, e.g.,

$$P^{store} = \frac{(dE^{el})^{store}}{dt} = \frac{(dE^{el})^{store}}{d\epsilon_1} \dot{\epsilon}_1 \qquad (9)$$

Figure 14 plots the evolutions of stored and released elastic powers along the biaxial loading path, where a small incremental axial strain $d\epsilon_1$ of 0.01% is considered. Note that, the difference in stored and released elastic powers equals to the net elastic power, which represents the incremental rate of net elastic energy in Fig. 7(b). It can be observed in Fig. 14(a) that the evolutions of $P^{store}$ and $P^{release}$ are smooth before state B with a much larger amount of energy
storage than energy release. The dominant role of energy storage representing the increase of net elastic storage gives rise to the decrease of the degree of stability for the dense specimen until the system is unable to store more elastic energy at state B, where a transient balance between $P^{store}$ and $P^{release}$ is obtained. As a result, the system is unstable with the emergence of an instability cone together with a negative degree of stability. After state B, an overall decreasing of net elastic energy occurs, while the energy-storing fluctuates with a weaker amplitude. This suggests that the contact-opening process is activated significantly together with the contact-closing process, giving rise to the microstructural reorganization to accommodate the shearing loading. As a result, a new state with fluctuating stability is achieved for the dense specimen, where the microstructural reorganization allows the shear band to develop. For the loose specimen, as illustrated in Fig. 14(b), the fluctuations in the energy-storing $P^{store}$ and energy-releasing $P^{release}$ increase simultaneously with high initial values. This indicates that the contact-opening and closing processes are active from the very beginning of biaxial loading, which is likely the reason for the stability fluctuation from the beginning of the shear loading.

<figure><img src="image_16.png" /><figcaption>Fig. 14. Evolutions of stored and released power along the biaxial loading path (a) for dense specimen and (b) for loose specimen.</figcaption></figure>
Moreover, contact closing and opening over an incremental loading cause a change in configuration with no frictional dissipation. From a mesoscopic perspective, the destructed meso-loop transform by contact closing (contact c) into two new meso-loops with lower orders. The opening of contact c induces the reverse of this process. The meso-configurational change deserves further discussion in this section to unravel the physics underlying the stability fluctuation for the granular assemblies. The ratios of the number of disappearing loops $Ln$ or generated loops $Ln$ within an incremental axial strain of 0.1%, to the number of loops $Ln$ at the current state are denoted as $R_{Ln}^{-}$ or $R_{Ln}^{+}$.

Figure 15 plots the evolution of ratios $R_{Ln}^{-}$ or $R_{Ln}^{+}$ for different loop categories within the dense specimen. For each category of loops, we can observe a significant difference between the ratios of disappearing and generating loops before state B, where the larger ratios of $R_{L3}^{-}$ and $R_{L6}^{+}$ are consistent with the results in Zhu et al. (2016b), suggesting the mutual conversion between L3 and L6+ in the dense specimen. A loop with a larger size (e.g., L6+) is more deformable, and therefore less stable. On the contrary, a small loop (e.g., L3) is prone to hold its original structure without any geometrical transformation and thus it keeps a structural stability. In this sense, the disappearing of L3 and generation of L6+ before state B indicate that the global stability of the topological structures is degrading, consistent with the decrease in the degree of stability $D^s$. After state B, the ratios of $R_{Ln}^{-}$ and $R_{Ln}^{+}$ for each loop category converge in a fluctuation regime around a constant average value. Correspondingly, the dense specimen undergoes stability fluctuation at the post-bifurcation stage (after state B). For the loose specimen, as shown in Fig. 16, the ratios $R_{Ln}^{-}$ and $R_{Ln}^{+}$ for each loop category converge and fluctuate around a constant average value from
493 the very beginning of the loading to the critical state regime, highlighting the stability fluctuation

494 starting from the initial loading.

<figure><img src="image_17.png" /><figcaption>Fig. 15. Evolutions of the number ratios of disappearing loops and generated loops along the biaxial loading path (a) for L3; (b) for L4; (c) for L5; and (d) for L6+ within the dense specimen.</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Fig. 16. Evolutions of the number ratios of disappearing loops and generated loops along the biaxial loading path (a) for L3; (b) for L4; (c) for L5; and (d) for L6+ within the loose specimen.</figcaption></figure>

It should be noted that when the microstructure within the specimen is composed of a robust topology that fully accommodates perturbation in all directions, the specimen is reputed stable. On the contrary, if the microstructure within the specimen is disrupted without being able to accommodate perturbations in all directions, the likelihood of potential instability raises. In this sense, this phenomenon of fluctuating stability at critical state regime can be regarded as stemming mainly from the configurational transformation of meso-loops that continues to occur within the granular assembly. This entails the dynamic elastic energy processes corresponding to repeated
energy storage and release, together with the continuous increase in frictional dissipation (Wang et al., 2024). Thus, at the critical state regime, a currently stable assembly can be destabilized by local fluctuations likely to occur during any subsequent loading (as represented by the non-zero ratios of disappearing and generating loops). As the shear loading proceeds, the meso-loops within the granular assembly continue to adjust and reorganize towards a new robust topology which can accommodate perturbation in all directions. A stable state is achieved again with $P^{\text{uns}} = 0$. In this sense, the fluctuating stability regime could be regarded as a stable bifurcated branch, where the configurational transitions continuously adjust the topological structures to resist against the external loading.

## 6 Conclusions

Assembling the second-order work theory with configurational arguments, this paper discusses the stability of granular assemblies with different initial densities along a biaxial loading using DEM simulations. The formation of the shear band is linked with the evolution of elastic energy processes occurring during the loading. Based on a strain controlled directional analysis procedure, the stability of sheared granular assemblies at the critical state regime is discussed by elucidating some of the meso-loop transformation mechanisms.

The dense specimen holds a limited capacity of storing elastic energy very close to the stress peak resulting in the spatial bifurcation with the formation of a shear band. After that, the system releases a fraction of the whole elastic energy to reach a new state where the shear band continues to develop and stabilize in the critical state regime. On the other hand, the capacity of storing elastic energy for the loose specimen increases continuously, and also approaches a stationary state regime in a very smooth way, following a macroscopic strain-hardening regime with no occurrence of strain localization. The dependence of spatial bifurcation on the limited
capacity of storing elastic energy suggests the intricate relation between the elastic energy storage and the stability of granular assemblies.

Based on the second-order work theory, a particular attention has been paid to the magnitude of the strain increment required to trigger off potential instabilities. A threshold value exits for strain increments below which the sign of the normalized second-order work is very sensitive to the magnitude of strain increment. This threshold can be regarded as an appropriate value of strain increment to detect the stability of specimens at different mechanical states. By defining two stability indexes, i.e., the probability of instability and the degree of stability, it is shown that the dense specimen is stable at the beginning of biaxial loading, while the degree of stability gradually decreases due to the smoothly increase in elastic energy storage. The vanishing of the second-order work is firstly observed before the stress peak. As the maximum capacity of storing elastic energy is achieved, a resulting high instability level is confirmed by observing a local maximum value of probability of instability together with a local minimum value of degree of instability. In the post-bifurcation stage, the significant release in net elastic energy in the dense specimen brings the system to a stable stage on average, but with fluctuating stability. The loose specimen undergoes stability fluctuations from the very beginning of the biaxial loading. Noticeably, the fluctuation of stability tends to be more significant as the granular specimen approaches the critical state regime, governed by the increase in the magnitude of elastic energy fluctuation. Moreover, considering the synchronization between the local increase/decrease in the stress and the storage/release of elastic energy, the stability fluctuation closely relates to the stress fluctuations which start from the stress peak state for dense specimen and from the very beginning of biaxial loading for loose specimen. Contrary to elastic energy fluctuation, the stability of the mechanical state considered is not sensitive to its location along the *q*-curve.
Considering the dynamic microstructural rearrangement within the granular assembly under shearing, the abundant generation of more deformable L6+ loops, together with the disappearing of L3 loops which hold a structural stability, contribute together to the decreasing stability degree before the peak stress state for the dense specimen. When the stability fluctuation regime takes place, the rates of creation and destruction of loops of the same order converge around a non-zero constant average value, entailing repeated energy storing and releasing. That is, the configurational transition is continuously adjusting the topological structures under external loading to approach a stable bifurcated branch in the form of stability fluctuations.

Finally, insofar as 2D specimens are considered, this manuscript sheds a new light on the mechanical stability of granular assemblies in relation with the elastic energy processes taking place in the configurational transformations which operate within the assembly. It turns out that the stability of a granular assembly is strongly governed by the transformations that occur between the loops of different orders. It is worth emphasizing that the critical state regime can be interpreted within the statistical mechanics framework. This is indeed a stationary regime, where fluctuations occur, as a response of a continuous input of external energy. The loop transformations that permanently operate, associated with important elastic energy mechanisms, reveal that the critical state regime is not at statistical equilibrium. The microstructure constantly rearranges and remains adapted against the external loading conditions. Extending this investigation in 3D conditions, where grain loops should be replaced with fully 3D clusters, is now the challenge at hand.

## 579 Acknowledgements

This work was supported by the financial support provided by National Natural Science Foundation of China (No.52278327) and China Postdoctoral Science Foundation (No.2024M760208). We also gratefully acknowledge the International Research Network
GeoMech for having offered the opportunity to kick start this project through the organization of

many workshops with topics connected to the concepts developed in this paper (http://gdr-mege.univ-lr.fr/).

## Appendix A

In this study, as shown in Fig. A1, the specific influence of the magnitude of incremental strain $||d\varepsilon||$, ($||d\varepsilon|| \in \{0.01\%, 0.025\%, 0.05\%, 0.075\%, 0.1\%$), on the normalized second-order work envelops $W_2^{norm}$ is discussed for the pre-stabilized dense specimen at mechanical state A (maximal contractancy state) and state B (peak stress state). Note that, the strain probe procedure remains the same for each strain increment as described in section 4.1. In particular, the same strain loading rate is used. A circular representation is given in Fig. A1 as the second-order work is positive outside the grey dashed circle and negative inside. For mechanical state A, no vanishing of the second-order work is observed in Fig. A1(a), whatever the magnitude of incremental strain $||d\varepsilon||$ considered. However, for mechanical state B as illustrated in Fig. A1(b), no vanishing of the second-order work $W_2^{norm}$ is observed when applying a small incremental strain of $||d\varepsilon|| = 0.01\%$, while the second-order work comes to vanish when $||d\varepsilon|| \ge 0.025\%$. This suggests that for the specimen size considered, $||d\varepsilon|| \ge 0.025\%$ is the minimum perturbation able to trigger significant particle reorganization.
<figure><img src="image_19.png" /><figcaption>Fig. A1. The normalized second-order work envelops for mechanical state A (a) and B (b) under different strain increment. $W_2^{norm} = 0$ is indicated by the grey dashed circle.</figcaption></figure>

600

601

602

603

604

605

606

607

608

609

610

611

612

613

614

615

In addition, the negative value of $W_2^{norm}$ and the instability cone (gathering all strain probe directions such that $W_2^{norm} \le 0$) are sensitive to the variation of $||d\varepsilon||$ as illustrated in Fig. A1(b).

The loading direction of $\theta_e = 152°$ for the dense specimen, which corresponds to the incremental loading direction with negative values of $W_2^{norm}$ for state B in Fig. A1, is considered to illustrate the evolution of $W_2^{norm}$ with the numerical loading time step (numerical time iterations), as shown in Fig. A2. It should be noted that, the corresponding end of incremental strain loading for different magnitudes considered ($||d\varepsilon|| \in \{0.01\%, 0.025\%, 0.05\%, 0.075\%, 0.1\%\}$) are highlighted in Fig. A2.

For these mechanical states belonging to the initial loading stage in Fig. A2(a), a slight increase of $W_2^{norm}$ is observed at the initial stage of the incremental loading. As the perturbation continues to be applied, $W_2^{norm}$ decreases sharply and reaches a stationary plateau when $||d\varepsilon|| \ge 0.025\%$ (around 4000 numerical time steps). Following the initially slight increase, $W_2^{norm}$ decreases sharply to a negative minimum for $\varepsilon_1 \in \{1.1\%, 1.5\%, 2.5\%, 4.5\%, 8.5\%\}$, and takes negative values when
$\|d\varepsilon\| \ge 0.05\%$ (around 7000 numerical time steps) or continues to evolve to positive values. While $W_2^{\text{norm}}$ decreases slowly for $\varepsilon_1 \in \{2.0\%, 3.0\%, 3.5\%, 6.5\%, 7.5\%, 10.5\%, 12\%, 15\%\}$, and then it either approaches a stationary plateau with positive value when $\|d\varepsilon\| \ge 0.05\%$ or tends to negative value when $\|d\varepsilon\| > 0.075\%$ (around 10000 numerical time steps).

<figure><img src="image_20.png" /><figcaption>Fig. A2. Evolution of normalized second-order work with numerical loading timestep (numerical time iterations) for dense specimen at different mechanical state under loading direction $\theta_g = 152^\circ$.</figcaption></figure>

These different evolutions of second-order work probably indicate that the microstructure of the granular assembly reorganizes continuously in response to the external perturbation.
However, a too small perturbation magnitude is not enough to produce reorganization of the whole specimen. For example, $W_2^{norm}$ is positive in all mechanical states when $||d\varepsilon|| = 0.01\%$.

Contact sliding as one of the important factors associated with microstructural rearrangement of granular assemblies is considered to highlight the specific influence of strain increment on the assembly response. By considering five magnitudes of incremental strain loading $||d\varepsilon||$ applied to the dense specimen at state B, the ratio ($N_s^f / N_s^i$) of the sliding number at each end of strain probe procedure $N_s^f$ to the one at the initial pre-stabilized state $N_s^i$ is presented in Fig. A3 for $\theta_e \in \{47^\circ, 152^\circ\}$. Note that $\theta_e = 47^\circ$ corresponds to the incremental loading direction with positive second-order work as shown in Fig. A1(b). The case of $N_s^f / N_s^i > 1.0$ suggests that a more active structural rearrangement occurs within the granular assembly under incremental perturbation, while the other situations correspond to a weaker structural rearrangement.

Different results for the two directions considered can be observed in Fig. A3. For $\theta_e = 47^\circ$, the ratio of $N_s^f / N_s^i$ decreases with the increase of $||d\varepsilon||$ and keeps constant when $||d\varepsilon|| \ge 0.025\%$. It suggests that the continuous perturbation (incremental loading) does not cause strong contact rearrangement within specimens, which corresponds to positive values of $W_2^{norm}$ in Fig. A1(b). On the other hand, the sliding contact number $N_s^f$ increases significantly when $||d\varepsilon|| < 0.025\%$ for $\theta_e = 152^\circ$, which indicates that a small incremental loading can cause a strong microstructural rearrangement. Then, for $||d\varepsilon|| \ge 0.025\%$, the ratio of $N_s^f / N_s^i$ decreases to a constant level whereas it is always greater than 1.0, accompanied by the occurrence of negative values of $W_2^{norm}$ in Fig. A1(b). The constant ratio greater than 1.0 probably suggests that the microstructure continues reorganizing to accommodate incremental strain perturbation.
As a result, the influence of $||d\varepsilon||$ on stability analysis seems to admit a threshold, which is consistent with the numerical results in Wautier et al. (2018). The magnitude of $||d\varepsilon||$ should be selected larger than this threshold to ensure that the incremental perturbation propagates to the whole specimen. According to these observations in Fig. A2 and A3, the threshold magnitude corresponding to different mechanical states along the biaxial loading path is different. However, a magnitude of $||d\varepsilon||=0.05\%$ is adopted in this study, which is appropriately large to propagate the incremental perturbation throughout the whole specimen.

<figure><img src="image_21.png" /><figcaption>Fig. A3. Evolution of the sliding ratio with the size of incremental strain for the dense specimen at state B, for direction $\theta_c = 152^\circ$ and direction $\theta_c = 47^\circ$. The sliding ratio is the ratio of the sliding contact number at the end of strain probe procedure to that at the initial pre-stabilized state.</figcaption></figure>

# Appendix B

The zoomed pictures of elastic energy fluctuations (red curve) and stress fluctuations (black curve) for dense and loose specimens are illustrated in Figure B1 and B2, where the mechanical states considered in the manuscript are marked by blue crosses.
<figure><img src="image_22.png" /><figcaption>Fig. B1. Zoomed pictures of elastic energy fluctuations (red curve) and stress fluctuations (black curve) for dense specimen.</figcaption></figure>

664

665

666
<figure><img src="image_23.png" /><figcaption>Fig. B2. Zoomed pictures of elastic energy fluctuations (red curve) and stress fluctuations (black curve) for loose specimen.</figcaption></figure>

References

Bardenhagen, S., Brackbill, J., Sulsky, D., 2000. Numerical study of stress distribution in sheared granular material in two dimensions. Phys. Rev. E 62, 3882.
https://doi.org/10.1103/PhysRevE.62.3882

Chen, Q., Gao, Y., Yuan, Q., Wang, Y.-H., 2022. The correlation of macro deformation and micro kinematics for dense granular material subjected to shearing. Comput. Geotech. 141, 104523.
https://doi.org/10.1016/j.compgeo.2021.104523

Clerc, A., Wautier, A., Bonelli, S., Nicot, F., 2021. Meso-scale signatures of inertial transitions in granular materials. Granul. Matter 23, 28. https://doi.org/10.1007/s10035-021-01087-5
679 Cundall, P.A., Strack, O.D., 1979. A discrete numerical model for granular assemblies.
680 Géotechnique 29, 47-65. https://doi.org/10.1680/geot.1979.29.1.47

681 Da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N., Chevoir, F., 2005. Rheophysics of dense granular materials: Discrete simulation of plane shear flows. Phys. Rev. E 72, 021309.
682 https://doi.org/10.1103/PhysRevE.72.021309

683 Daouadji, A., Darve, F., Al Gali, H., Hicher, P.-Y., Laouafa, F., Lignon, S., Nicot, F., Nova, R., Pinheiro, M., Prunier, F., 2011. Diffuse failure in geomaterials: experiments, theory and modelling. Int. J. Numer. Anal. Methods Geomech. 35, 1731-1773.
687 https://doi.org/10.1002/nag.975

688 Darve, F., Sibille, L., Daouadji, A., Nicot, F., 2007. Bifurcations in granular media: macro- and micro-mechanics approaches. C. R. Mecanique 335, 496-515.
689 https://doi.org/10.1016/j.crme.2007.08.005
690

691 Deng, N., Wautier, A., Thiery, Y., Yin, Z., Hicher, P.Y., Nicot, F., 2021. On the attraction power of critical state in granular materials. J. Mech. Phys. Solids 149, 104300.
692 https://doi.org/10.1016/j.jmps.2021.104300

694 Deng, N., Wautier, A., Tordesillas, A., Thiery, Y., Yin, Z., Hicher, P.Y., Nicot, F., 2022. Lifespan dynamics of cluster conformations in stationary regimes in granular materials. Phys. Rev. E 105, 014902. https://doi.org/10.1103/PhysRevE.105.014902

697 di Prisco, C., Imposimato, S., 1996. Time dependent mechanical behaviour of loose sands. Mechanics of Cohesive-frictional Materials: An International Journal on Experiments, Modelling and Computation of Materials and Structures 1, 45-73.

700 Drescher, A., De Jong, G.D.J., 2006. Photoelastic verification of a mechanical model for the flow of a granular material. Soil Mechanics and Transport in Porous Media: Selected Works of G. de Josselin de Jong, 28-43.

703 Guo, N., Zhao, J., 2016. 3D multiscale modeling of strain localization in granular media. Comput. Geotech. 80, 360-372. https://doi.org/10.1016/j.compgeo.2016.01.020

705 Hill, R., 1958. A general theory of uniqueness and stability in elastic-plastic solids. J. Mech. Phys. Solids 6, 236-249. https://doi.org/10.1016/0022-5096(58)90029-2

707 Inc, I.C.G., 2008. Particle Flow Code in 2 Dimensions, PFC2D Version 4.0. Minneapolis: Itasca Consulting Group Inc.

709 Jin, Z., Liu, J., Ma, G., Hu, C., Yang, Q., Shi, X., Wang, X., 2024. How Does the Largest Cluster
710 in the Strong Network Rule Granular Soil Mechanics? A DEM Study. Int. J. Numer. Anal. Geomech. 0,1-21. https://doi.org/10.1002/nag.3903

711 Methods Geomech. 0,1-21. https://doi.org/10.1002/nag.3903

712 Kruyt, N.P., 1999. Structured deformation in granular materials. Mech. Mater. 31(6), 407-429.

713 https://doi.org/10.1016/S0167-6636(99)00010-1

714 Kruyt, N.P., Rothenburg, L., 2006. Shear strength, dilatancy, energy and dissipation in quasi-static

715 deformation of granular materials. J. Stat. Mech. P07021. https://doi.org/10.1088/1742-

716 5468/2006/07/P07021

717 Li, X.S., Wang, Y., 1998. Linear Representation of Steady-State Line for Sand. J. Geotech.

718 Geoenviron. 124, 1215-1217. https://doi.org/10.1061/(ASCE)1090-0241(1998)124:12(1215)

719 Liu, J., Nicot, F., Zhou, W., 2018. Sustainability of internal structures during shear band forming

720 in 2D granular materials. Powder Technol. 338, 458-470.

721 https://doi.org/10.1016/j.powtec.2018.07.001

722 Liu, J., Wautier, A., Bonelli, S., Nicot, F., Darve, F., 2020a. Macroscopic softening in granular

723 materials from a mesoscale perspective. Int. J. Solids Struct. 193-194, 222-238.

724 https://doi.org/10.1016/j.ijsolstr.2020.02.022

725 Liu, J., Wautier, A., Zhou, W., Nicot, F., Darve, F., 2022. Incremental shear strain chain: a

726 mesoscale concept for slip lines in 2D granular materials. Granul. Matter 24(4), 119.

727 https://doi.org/10.1007/s10035-022-01258-y

728 Liu, J., Wu, X., Jiang, J., Ding, Z., Lü, C., Shi, X., 2023. A network-based investigation on the

729 strong contact system of granular materials under isotropic and deviatoric stress states.

730 Comput. Geotech. 153, 105077. https://doi.org/10.1016/j.compgeo.2022.105077

731 Liu, Y., Li, S., 2018. Numerical simulation and analysis of meso-mechanical structure

732 characteristic at critical state for granular media. Rock and Soil Mechanics 39, 2237-2248.

733 (in Chinese)

734 Liu, Y., Wang, X., 2024. Micro-mechanism of stress-dilatancy anisotropy in granular materials:

735 Affine and nonaffine deformation. Int. J. Numer. Anal. Methods Geomech. 48, 1460-1482.

736 https://doi.org/10.1002/nag.3693

737 Liu, Y., Wang, X., Yu, P., 2020b. Critical state of granular materials at mesoscale: Transition from

738 local to global. Int. J. Numer. Anal. Methods Geomech. 44, 1676-1694.

739 https://doi.org/10.1002/nag.3084

740 Lopera Perez, J., Kwok, F., Huang, X., Hanley, K., 2016. Exploring the micro-mechanics of triaxial
741 instability in granular materials. Géotechnique 66, 725-740.

742 https://doi.org/10.1680/jgeot.15.P.206

743 Neilsen, M.K., Schreyer, H.L., 1993. Bifurcations in elastic-plastic materials. Int. J. Solids Struct.

744 30, 521-544. https://doi.org/10.1016/0020-7683(93)90185-A

745 Nguyen, H.N.G., Prunier, F., Djeran-Maigre, I., Nicot, F., 2016. Kinetic energy and collapse of

746 granular materials. Granul. Matter 18, 5. https://doi.org/10.1007/s10035-016-0609-1

747 Nguyen, N.-S., Magoariec, H., Cambou, B., Danescu, A., 2009. Analysis of structure and strain at

748 the meso-scale in 2D granular materials. Int. J. Solids Struct. 46, 3257-3271.

749 https://doi.org/10.1016/j.ijsolstr.2009.04.019

750 Nicot, F., Darve, F., Khoa, H.D.V., 2007. Bifurcation and second-order work in geomaterials. Int.

751 J. Numer. Anal. Methods Geomech., 31(8), 1007-1032. https://doi.org/10.1002/nag.573

752 Nicot, F., Hadda, N., Bourrier, F., Sibille, L., Darve, F., 2011. Failure mechanisms in granular

753 media: a discrete element analysis. Granul. Matter 13, 255-260.

754 https://doi.org/10.1007/s10035-010-0242-3

755 Nicot, F., Lerbet, J., Darve, F., 2017. Second-order work criterion: from material point to boundary

756 value problems. Acta Mech 228, 2483-2498. https://doi.org/10.1007/s00707-017-1844-1

757 Nicot, F., Lin, M., Wautier, A., Wan, R., Darve, F., 2024. Configurational mechanics in granular

758 media. Granul. Matter 26, 74. https://doi.org/10.1007/s10035-024-01443-1

759 Nicot, F., Sibille, L., Darve, F., 2012. Failure in rate-independent granular materials as a bifurcation

760 toward a dynamic regime. Int. J. Plasticity 29, 136-154.

761 https://doi.org/10.1016/j.ijplas.2011.08.002

762 Nicot, F., Wang, X., Wautier, A., Wan, R., Darve, F., 2023. Shear banding as a dissipative structure

763 from a thermodynamic viewpoint. J. Mech. Phys. Solids 179, 105394.

764 https://doi.org/10.1016/j.jmps.2023.105394

765 Oda, M., Kazama, H., 1998. Microstructure of shear bands and its relation to the mechanisms of

766 dilatancy and failure of dense granular soils. Géotechnique 48, 465-481.

767 https://doi.org/10.1680/geot.1998.48.4.465

768 Ord, A., Hobbs, B., Regenauer-Lieb, K., 2007. Shear band emergence in granular materials - A

769 numerical study. Int. J. Numer. Anal. Methods Geomech. 31, 373-393.

770 https://doi.org/10.1002/nag.590

771 Prévost, J.-H., Höeg, K., 1975. Soil mechanics and plasticity analysis of strain softening.
772 Géotechnique 25, 279-297. https://doi.org/10.1680/geot.1975.25.2.279

773 Prunier, F., Nicot, F., Darve, F., Laouafa, F., Lignon, S., 2009. Three-Dimensional Multiscale

774 Bifurcation Analysis of Granular Media. J. Eng. Mech. 135, 493-509.

775 https://doi.org/10.1061/(ASCE)EM.1943-7889.0000003

776 Rattez, H., Shi, Y., Sac-Morane, A., Klaeyle, T., Mielniczuk, B., Veveakis, M., 2022. Effect of

777 grain size distribution on the shear band thickness evolution in sand. Géotechnique 72, 350-

778 363. https://doi.org/10.1680/jgeot.20.P.120

779 Roscoe, K.H., 1958. On yielding of soils. Géotechnique 8(1), 22-53.

780 https://doi.org/10.1680/geot.1958.8.1.22

781 Schofield, A., Wroth, P., 1968. Critical State Soil Mechanics. McGraw-Hill, London.

782 Tordesillas, A., 2007. Force chain buckling, unjamming transitions and shear banding in dense

783 granular assemblies. Philos. Mag. 87, 4987-5016.

784 https://doi.org/10.1080/14786430701594848

785 Tordesillas, A., Lin, Q., Zhang, J., Behringer, R., Shi, J., 2011. Structural stability and jamming of

786 self-organized cluster conformations in dense granular materials. J. Mech. Phys. Solids 59,

787 265-296. https://doi.org/10.1016/j.jmps.2010.10.007

788 Tordesillas, A., Zhang, J., Behringer, R., 2009. Buckling force chains in dense granular assemblies:

789 physical and numerical experiments. Geomech. Geoengin. 4, 3-16.

790 https://doi.org/10.1080/17486020902767347

791 Veylon, G., Nicot, F., Zhu, H., Darve, F., 2018. Microstructure incidence on the bifurcation domain

792 topology in granular materials. J. Eng. Mech. 144, 04018031.

793 https://doi.org/10.1061/(ASCE)EM.1943-7889.0001446

794 Wan, R., Nicot, F., Darve, F., 2017. Failure in geomaterials: a contemporary treatise. Elsevier.

795 Wan, R., Pinheiro, M., Daouadji, A., Jrad, M., Darve, F., 2013. Diffuse instabilities with transition

796 to localization in loose granular materials. Int. J. Numer. Anal. Methods Geomech. 37(10),

797 1292-1311. https://doi.org/1292-1311.10.1002/nag.2085

798 Wan, R.G., 2011. Instabilities and bifurcations in geomechanics. Int. J. Numer. Anal. Methods

799 Geomech. 35(2), 111-111. https://doi.org/111-111.10.1002/nag.1007

800 Wang, T., Wautier, A., Liu, S., Nicot, F., 2021. How fines content affects granular plasticity of

801 under-filled binary mixtures. Acta Geotech. 17, 2449-2463. https://doi.org/10.1007/s11440-

802 021-01430-3
803 Wang, X., Liu, Y., Nicot, F., 2024. Energy processes and phase transition in granular assemblies.
804 Int. J. Solids Struct. 289, 112634. https://doi.org/10.1016/j.ijsolstr.2023.112634

805 Wanjura, C.C., Gago, P., Matsushima, T., Blumenfeld, R., 2020. Structural evolution of granular
806 systems: theory. Granul. Matter 22, 91. https://doi.org/10.1007/s10035-020-01056-4

807 Wautier, A., Bonelli, S., Nicot, F., 2018. Micro-inertia origin of instabilities in granular materials.
808 Int. J. Numer. Anal. Methods Geomech. 42, https://doi.org/1037-1056_10.1002/nag.2777

809 Yan, Z., Liu, Y., Zhao, D., 2024. Instability of binary mixtures subjected to constant shear drained
810 stress path: Insight from macro and micro perspective. Int. J. Numer. Anal. Methods Geomech.
811 48(16), 1-17. http://dx.doi.org/10.1002/nag.3823

812 Yu, P., Liu, Y., 2022. Micromorphic modeling of wave propagation in granular materials:
813 Considering particle nonaffine micro-deformation and kinematics decomposition. Powder
814 Technol. 411, 117919. http://dx.doi.org/10.1016/j.powtec.2022.117919

815 Zhang, J., Majmudar, T., Behringer, R., 2008. Force chains in a two-dimensional granular pure
816 shear experiment. Chaos 18, 1079. http://dx.doi.org/10.1063/1.2997139

817 Zhao, C., Yin, Z.-Y., Hicher, P.-Y., 2018. A multiscale approach for investigating the effect of
818 microstructural instability on global failure in granular materials. Int. J. Numer. Anal.
819 Methods Geomech. 42(17), 2065-2094. http://dx.doi.org/10.1002/nag.2842

820 Zhu, H., Nguyen, H.N.G., Nicot, F., Darve, F., 2016a. On a common critical state in localized and
diffuse failure modes. J. Mech. Phys. Solids 95, 112-131.
821 http://dx.doi.org/10.1016/j.jmps.2016.05.026

823 Zhu, H., Nicot, F., Darve, F., 2016b. Meso-structure evolution in a 2D granular material during
824 biaxial loading. Granul. Matter 18, 3. http://dx.doi.org/10.1007/s10035-016-0608-2

825