Eigenvalue analysis of stress-strain curve of two-dimensional amorphous solids of dispersed frictional grains with finite shear strain

Amorphous materials consisting of dispersed repulsive grains, such as powders, colloids, bubbles, and emulsions, behave as fragile solids above jamming density Jaeger96 ; Durian95 ; Wyart05 ; Baule18 ; Liu98 . When we consider a response of such materials to an applied strain $\gamma$, the rigidity $G$ is independent of the strain in the linear response regime, whereas it exhibits softening in the nonlinear regime Nagamanasa14 ; Knowlton14 ; Kawasaki16 ; Leishangthem17 ; Bohy17 ; Ozawa18 ; Clark18 ; Boschan19 ; Singh20 ; Otsuki21 . Above the yielding points, there are some plastic events, such as stress avalanches in the collection of grains.

The Hessian matrix determined by the configuration of grains is commonly used for amorphous solids consisting of frictionless grains  Wyart05 ; WyarLettt05 ; Ellenbroek06 ; Lerner16 ; Gartner16 ; Bonfanti19 ; Baule18 ; Morse20 . To determine the rigidity, eigenvalue analysis of the Hessian matrix DeGiuli14 ; Mizuno16 ; Maloney04 ; Maloney06 ; Lemaitre06 ; Zaccone11 ; Palyulin18 , is commonly used, but we have only obtained semi-quantitative agreements between the theory and simulation so far Zaccone11 ; Palyulin18 . Recently, some researchers examined the applicability of the instantaneous normal mode analysis to systems that are not fully equilibrated and found qualitative agreement between the analysis and simulation Oyama21 ; Kriuchevskyi22 , although they could not get quantitatively accurate results. Some studies have suggested that the decrement of the non-zero smallest eigenvalue of the Hessian matrix with the strain is a precursor of an avalanche or stress drop near a critical strain $\gamma_{c}$  Maloney04 ; Maloney06 ; Manning11 ; Dasgupta12 ; Ebrahem20 . Correspondingly, some studies indicated that the rigidity $G$ and the stress $\sigma$ should behave as $G-G_{\rm reg}\propto-1/\sqrt{\gamma_{c}-\gamma}$ and $\sigma-\sigma_{\rm reg}\propto\sqrt{\gamma_{c}-\gamma}$ near $\gamma_{c}$, where $G_{\rm reg}$ and $\sigma_{\rm reg}$ are the regular parts of the rigidity and stress conversing to constants at $\gamma_{c}$, respectively Maloney04 ; Maloney06 ; Karmakar10PRL ; Richard20 .

In general, the frictional force between the grains cannot be ignored in physical situations. Because the frictional force generally depends on the contact history, Hessian analysis is not applicable to such systems. Thus, Chattoraj et al. adopted the Jacobian matrix instead of the Hessian matrix to discuss the stability of the configuration of frictional grains under strain Chattoraj19 . They performed eigenvalue analysis under athermal quasi-static shear processes and determined the existence of oscillatory instability originating from interparticle friction at a certain strain Chattoraj19 ; Chattoraj19E ; Charan20 . Moreover, some studies have performed an analysis of the Hessian matrix with the aid of an effective potential for frictional grains Somfai07 ; Henkes10 . Recently, Liu et al. suggested that Hessian matrix analysis with another effective potential can be used, even if slip processes exist Liu21 . Previous studies Somfai07 ; Henkes10 have reported that the friction between grains causes a continuous change in the functional form of the density of states (DOS), which differs from that of frictionless systems. So far, there have been few theoretical studies to determine the rigidity of the frictional grains.

In our previous study Ishima2022 , we developed an analysis of the Jacobian matrix to determine the rigidity of two-dimensional amorphous solids consisting of frictional grains interacting with the Hertzian force in the linear response to an infinitesimal strain. In the study, we ignore the dynamic friction caused by the slip processes of contact points. We found that there are two modes in the DOS for a sufficiently small tangential-to-normal stiffness ratio. Rotational modes exist in the region of low-frequency or small eigenvalues, whereas translational modes exist in the region of high-frequency or large eigenvalues. The rigidity determined by the translational modes is in good agreement with that obtained by the molecular dynamics simulations, whereas the contribution of the rotational modes is almost zero. Nevertheless, there are several shortcomings in the previous analysis. (i) The analysis can be used only in the linear response regime, where the base state is not influenced by the applied strain. (ii) As a result, we cannot discuss the behavior of plastic deformations or avalanches of grains. (iii) Even if we restrict our interest to the linear response regime, we cannot include the effect of tangential contact for the preparation of the initial configuration. (iv) We also ignored the effect of Coulomb’s slip between the contacted grains Ishima2022 .

The purpose of this study is to overcome the shortcomings of our previous study except for point (iv) Ishima2022 . Thus, we analyze a collection of two-dimensional grains interacting with repulsive harmonic potentials within the contact radius, without considering Coulomb’s slip between the grains. Owing to the special properties of the harmonic potential, the eigenvalue analysis of the Jacobian matrix becomes equivalent to that of the Hessian matrix. Subsequently, using the eigenvalue analysis of the Hessian matrix, we demonstrate that the theoretical rigidity under a large strain agrees with that obtained by the simulation.

The remainder of this paper is organized as follows. In Sec. II, we introduce the model to be analyzed in this study. In Sec. III, we summarize the theoretical framework for determining the rigidity of an amorphous solid consisting of frictional grains without considering Coulomb’s slip process. In Sec. IV, we present the results of the stress-strain relation obtained using the theory formulated in Sec. III. We also compare the theoretical results with the simulation results to demonstrate the relevancy of our theoretical analysis. In Sec. V, we summarize the obtained results and address future tasks to be solved. In Appendix A, we numerically demonstrate the absence of history-dependent tangential deformations in our system. In Appendix B, we explain the detailed behavior of the eigenvalue near the stress-drop points. In Appendix C, we present some detailed properties of the Hessian matrix in a harmonic system. In Appendix D, we present some properties of the Jacobian matrix in a harmonic system and demonstrate its equivalency to the Hessian matrix. Appendix E presents the detailed properties of rigidity.

Our system contains $N$ frictional circular disks embedded in a two-dimensional space. To prevent the system from crystallization, it contains an equal number of grains with diameters $d$ and $d/1.4$ Luding01 . We assume that the mass of grain $i$ is proportional to $d_{i}^{2}$, where $d_{i}$ is the diameter of $i$-th grain. For later convenience, we introduce $m$ as the mass of grain having a diameter $d$. In this study, $x_{i},y_{i}$, and $\theta_{i}$ denote $x$ and $y$ components of the position of $i$-th grain, and the rotational angle of the $i$-th grain, respectively. We introduce the generalized coordinates of the $i$-th grain as follow:

$$\bm{q}_{i}:=(\bm{r}_{i}^{\textrm{T}},\ell_{i})^{\textrm{T}},$$

(1)

where $\bm{r}_{i}:=(x_{i},y_{i})^{\textrm{T}}$, $\ell_{i}:=d_{i}\theta_{i}/2$, and the superscript T denotes the transposition.

Let the force, and $z$-component of the torque acting on the $i$-th grain be $\bm{F}_{i}:=(F_{i}^{x},F_{i}^{y})^{\textrm{T}}$ and $T_{i}$, respectively. Then, the equations of motion of $i$-th grain are expressed as

	$\displaystyle m_{i}\frac{d^{2}\bm{r}_{i}}{dt^{2}}$	$\displaystyle=\bm{F}_{i},$		(2)
	$\displaystyle I_{i}\frac{d^{2}\theta_{i}}{dt^{2}}$	$\displaystyle=T_{i},$		(3)

with mass $m_{i}$ and momentum of inertia $I_{i}:=m_{i}d_{i}^{2}/8$ of $i$-th grain. In a system without volume forces, such as gravity, we can write

	$\displaystyle\bm{F}_{i}$	$\displaystyle=\sum_{j\neq i}\bm{f}_{ij}-m_{i}\eta_{D}\dot{\bm{r}}_{i},$		(4)
	$\displaystyle T_{i}$	$\displaystyle=\sum_{j\neq i}T_{ij}-I_{i}\eta_{D}\dot{\theta}_{i},$		(5)

where we have adopted the notation $\dot{A}:=dA/dt$ for arbitrary variable $A$ such as $A=\bm{r}_{i}$ and ${\theta}_{i}$. Here, $\bm{f}_{ij}$ and $T_{ij}$ are the force and $z$-component of the torque acting on the $i$-th grain from the $j$-th grain, respectively. As a simplified description of the drag terms from the background fluid, $\eta_{D}$ is a damping constant uniformly acting on grains, and $T_{ij}$ is given by

$\displaystyle T_{ij}=-\frac{d_{i}}{2}({n}_{ij}^{x}{f}_{ij}^{y}-{n}_{ij}^{y}{f}_{ij}^{x}),$

(6)

where we have introduced the normal unit vector between $i$-th and $j$-th grains as $\bm{n}_{ij}:=\bm{r}_{ij}/|\bm{r}_{ij}|:=(\bm{r}_{i}-\bm{r}_{j})/|\bm{r}_{i}-\bm{r}_{j}|$. The force $\bm{f}_{ij}$ can be divided into the normal part $\bm{f}_{N,ij}$ and tangential part $\bm{f}_{T,ij}$ as

$\displaystyle\bm{f}_{ij}=(\bm{f}_{N,ij}+\bm{f}_{T,ij})H(d_{ij}/2-|\bm{r}_{ij}|),$

(7)

where $d_{ij}:=d_{i}+d_{j}$ and $H(x)$ is Heaviside’s step function, taking $H(x)=1$ for $x>0,$ and $H(x)=0$ otherwise. We assume that the contact force is expressed as

	$\displaystyle\bm{f}_{N,ij}:$	$\displaystyle=k_{N}\xi_{N,ij}\bm{n}_{ij},$		(8)
	$\displaystyle\bm{f}_{T,ij}:$	$\displaystyle=k_{T}\xi_{T,ij}\bm{t}_{ij},$		(9)

where $k_{N}$ and $k_{T}$ are the stiffness parameters of normal and tangential contacts, respectively. The contact force can be derived from the harmonic potential. In Eqs. (8) and  (9) we have introduced $\xi_{N,ij}:=d_{ij}/2-|\bm{r}_{ij}|$ and

	$\displaystyle\bm{\xi}_{T,ij}(t):$	$\displaystyle=\int_{C_{ij}(t^{\prime})}dt^{\prime}\bm{v}_{T,ij}(t^{\prime})$
		$\displaystyle\quad-\left[\left(\int_{C_{ij}(t^{\prime})}dt^{\prime}\bm{v}_{T,ij}(t^{\prime})\right)\cdot\bm{n}_{ij}(t)\right]\bm{n}_{ij}(t),$		(10)

where we have used

$\displaystyle\bm{v}_{T,ij}:=\dot{\bm{r}}_{ij}-\dot{\bm{\xi}}_{N,ij}+\bm{u}_{ij}(d_{i}\dot{\theta}_{i}+d_{j}\dot{\theta}_{j})/2$

(11)

with $\bm{u}_{ij}:=(n_{ij}^{y},-n_{ij}^{x})^{\textrm{T}}$, $\bm{t}_{ij}:=-\bm{\xi}_{T,ij}/|\xi_{T,ij}|$, $\xi_{T,ij}:=|\bm{\xi}_{T,ij}|$, and the integration over the duration time of contact between $i$-th and $j$-th grains $\int_{C_{ij}(t^{\prime})}dt^{\prime}$ with the trajectory $C_{ij}(t^{\prime})$ of the contact point between $i$-th and $j$-th grains at $t^{\prime}$. As shown in Appendix A, we have confirmed that the second term on the right-hand side (RHS) of Eq. (10) is zero for harmonic systems, although we do not have any mathematical proof for this statement thus far. For simplicity, we consider neither the effects of Coulomb’s slip in the tangential equation of motion nor the dissipative contact force, where the tangential forces $\bm{f}_{T,ij}^{C}$ including the effect of Coulomb’s slip process, satisfy

$\displaystyle\bm{f}_{T,ij}^{C}:=\left\{\begin{matrix}\bm{f}_{T,ij}&(\bm{f}_{T,ij}<\mu|\bm{f}_{N,ij}|)\\ \mu|\bm{f}_{N,ij}|\bm{t}_{ij}&(\textrm{otherwise})\end{matrix}\right.$

(12)

where $\mu$ is the friction constant.

We impose the Lees-Edwards boundary conditions Lees72 ; Evans08 , where the direction parallel to the shear strain is the $x$-direction. After generating a stable grain configuration via isotropic compression starting from a random configuration by using Eqs. (2)-(10) without strain (see detail in Ref. Ishima2022 ), we apply a step strain $\Delta\gamma$ to all grains, where $x$-coordinate of the position of the $i$-th grain is shifted by an affine displacement $\Delta x_{i}(\Delta\gamma):=\Delta\gamma y_{i}^{\textrm{FB}}(0)$. Here, the superscript FB denotes the force-balance (FB) state at which $\bm{F}_{i}=\bm{0}$ and $T_{i}=0$ for arbitrary $i$. As shown in Sec. III.1, the FB state is equivalent to the potential minimum for harmonic grains. Subsequently, the system is relaxed to an FB state. We further apply the step strain $\Delta\gamma$ associated with the subsequent relaxation process again to obtain the state at $2\Delta\gamma$. By repeating this process, we can reach a deformed state with the strain $\gamma$.

The plastic deformations for a large $\gamma$ depend on the choice of $\Delta\gamma$ Morse20 . Moreover, the theoretical formulation assumes $\Delta\gamma\to 0$. Therefore, we adopt the backtracking method Lerner09 ; Karmakar10PRE . If a plastic event is encountered under a fixed $\Delta\gamma_{\rm in}$, the system is restored to its original state without a plastic event. Subsequently, we apply a new strain, $0.1\Delta\gamma_{\rm in}$, to the system. Even if we encounter a plastic event with $0.1\Delta\gamma_{\rm in}$, we further examine the smaller step strain of $0.01\Delta\gamma_{\rm in}$. We repeat this procedure until we reach $\Delta\gamma<\Delta\gamma_{\rm Th}$ (see Appendix B).

We introduce the rate of nonaffine displacements for $r_{i}^{{\rm FB},\zeta}(\gamma)$ with $\zeta=x$ or $y$ and $\ell_{i}^{\rm FB}(\gamma)$ as:

	$\displaystyle\frac{d\mathring{r}_{i}^{\zeta}(\gamma)}{d\gamma}:$	$\displaystyle=\lim_{\Delta\gamma\to 0}\frac{r_{i}^{\textrm{FB},\zeta}(\gamma+\Delta\gamma)-r_{i}^{\textrm{FB},\zeta}(\gamma)}{\Delta\gamma}-\delta_{\zeta x}y_{i}^{\textrm{FB}}(\gamma),$		(13)
	$\displaystyle\frac{d\mathring{\ell}_{i}(\gamma)}{d\gamma}:$	$\displaystyle=\lim_{\Delta\gamma\to 0}\frac{\ell_{i}^{\textrm{FB}}(\gamma+\Delta\gamma)-\ell_{i}^{\textrm{FB}}(\gamma)}{\Delta\gamma}.$		(14)

Our system is characterized by the generalized coordinate $\bm{q}(\gamma):=(\bm{q}_{1}^{\textrm{T}}(\gamma),\bm{q}_{2}^{\textrm{T}}(\gamma),\cdots,\bm{q}_{N}^{\textrm{T}}(\gamma))^{\textrm{T}}$. The configuration in the FB state at strain $\gamma$ is denoted by $\bm{q}^{\textrm{FB}}(\gamma):=((\bm{q}_{1}^{\textrm{FB}}(\gamma))^{\textrm{T}},(\bm{q}_{2}^{\textrm{FB}}(\gamma))^{\textrm{T}},\cdots,(\bm{q}_{N}^{\textrm{FB}}(\gamma))^{\textrm{T}})^{\textrm{T}}$. The shear stress $\sigma(\gamma)$ at $\bm{q}^{\textrm{FB}}(\gamma)$ for one sample is given by:

	$\displaystyle\sigma(\bm{q}^{\textrm{FB}}(\gamma))=-\frac{1}{2L^{2}}\sum_{i}\sum_{j>i}$	$\displaystyle\left[f_{ij}^{x}(\bm{q}^{\textrm{FB}}(\gamma))r_{ij}^{y}(\bm{q}^{\textrm{FB}}(\gamma))\right.$
		$\displaystyle\ +\left.f_{ij}^{y}(\bm{q}^{\textrm{FB}}(\gamma))r_{ij}^{x}(\bm{q}^{\textrm{FB}}(\gamma))\right].$		(15)

The rigidity $g$ for one sample is defined as

$\displaystyle g$

$\displaystyle:=\left.\frac{d\sigma(\bm{q}(\gamma))}{d\gamma}\right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)},$

(16)

where the differentiation on the RHS of Eq. (16) is defined as follows:

	$\displaystyle\left.\frac{d\sigma(\bm{q}(\gamma))}{d\gamma}\right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}$
	$\displaystyle\quad\quad\quad:=\lim_{\Delta\gamma\to 0}\frac{\sigma(\bm{q}^{\textrm{FB}}(\gamma+\Delta\gamma))-\sigma(\bm{q}^{\textrm{FB}}(\gamma))}{\Delta\gamma}.$		(17)

In the numerical calculation, we use a non-zero but sufficiently small $\Delta\gamma$ for the evaluation of $g$. Then, the averaged rigidity $G$ is defined as

$\displaystyle G$

$\displaystyle:=\Braket{g},$

(18)

where $\langle\cdot\rangle$ is the ensemble average.

In this section, we introduce Hessian matrix in Sec. III.1 and theoretical expressions of rigidity in Sec. III.2.

We verify the validity of the shear modulus obtained by the eigenvalue analysis by comparing it with that obtained by the simulation. First, we have confirmed the quantitative accuracy of our analysis to obtain the rigidity in the linear response regime of our system (see Appendix E.2), as in Ref. Ishima2022 .

For the numerical FB condition, we use the condition $|\tilde{F}_{i}^{\alpha}|<F_{\textrm{Th}}$ for arbitrary $i$, where we adopt $F_{\textrm{Th}}=1.0\times 10^{-14}k_{N}d$ for the numerical calculation. In our simulation, we also adopt $\eta_{D}=\sqrt{k_{N}/m}$. In this study, we present the results for $k_{T}/k_{N}=1$, the area fraction $\phi=0.90$, and $\Delta\gamma_{\rm in}=1.0\times 10^{-4}$ with the ensemble averages of $30$ samples, except for Appendix E.2. Note that the we analyze only the systems with $\phi=0.90$ which is sufficiently larger than the jamming density for a two-dimensional system. We ignore the effect of dissipation in the eigenequation because the velocity of each grain is sufficiently small to incur infinitesimal agitation from the FB state. The time step used for the simulation, $\Delta t$, is set to $\Delta t=1.0\times 10^{-2}t_{0}$ with $t_{0}:=\sqrt{m/k_{N}}$, and numerical integration is performed using the velocity Verlet method.

Now, let us consider a nonlinear regime in which there are many plastic events caused by stress avalanches. Here, we regard an event as plastic if the condition (i) $\sigma(\gamma)-\sigma(\gamma-\Delta\gamma)<0$ or (ii) $G(\gamma-\Delta\gamma)-G(\gamma)>1.0\times 10^{-2}k_{N}$ is satisfied.

Figure 1 shows a typical example of the stress-strain curve obtained by one sample of the collection of grains based on both the simulation and eigenvalue analysis developed in the previous section under the condition $\Delta\gamma_{\rm Th}=\Delta\gamma_{\rm in}$. It should be noted that the difference between the theoretical and simulation results is almost invisible, even in the presence of avalanches. However, the eigenvalue analysis cannot be used immediately after plastic events, that is, for $\gamma\approx\gamma_{c}$ (see the inset of Fig. 1), because the stress is not determined by Eq. (51) immediately after a plastic event.

Figure 2 shows a comparison of the stress-strain curve for $\Delta\gamma_{\rm Th}=10^{-4}\Delta\gamma_{\rm in}$ (blue circles) with those for $\Delta\gamma_{\rm Th}=\Delta\gamma_{\rm in}$ (red triangles). From the figure, we cannot find any $\Delta\gamma_{\rm Th}$ dependence for $\gamma<0.08$, but some differences for larger $\gamma$ can be observed as a result of stress avalanches.

We might expect that some precursors of a stress-drop event can be detected from the behavior of the smallest non-zero eigenvalue. To verify this expectation, we plot the smallest non-zero eigenvalues in Fig. 3(a) near a critical strain with $\Delta\gamma_{\rm Th}=10^{-8}$. It can be observed that the eigenvalues changes discontinuously at the critical strain ($\gamma_{c}=0.34102862$ for $\gamma_{\rm Th}=10^{-8}$), where the critical strain converges if $\Delta\gamma_{\rm Th}<10^{-6}$ (see Appendix B for details). Notably, there is no precursor for the smallest eigenvalue below the critical strain, in contrast to Refs. Maloney04 ; Maloney06 ; Manning11 ; Dasgupta12 , where non-harmonic potentials are used. Correspondingly, we cannot find any singularity of the rigidity as $g-g_{\textrm{reg}}\sim-(\gamma_{c}-\gamma)^{-1/2}$ as in Fig. 3(b) for $\gamma\lesssim\gamma_{c}$ predicted in Refs. Maloney04 ; Maloney06 ; Karmakar10PRL ; Richard20 . The absence of the precursors and singularities in our model can be understood in the form of the Hessian matrix presented in Appendix D.3. In non-harmonic systems, some elements of the Hessian matrix become zero as $\xi_{N,ij}\to 0$ when the contact between the $i$-th and $j$-th grains disappears. This leads to the precursors and singularities Maloney04 ; Maloney06 . However, in the harmonic systems, the corresponding element approaches a non-zero constant in the limit $\xi_{N,ij}\to 0$ (see Appendix D.3), which results in the absence of the precursors.

Figure 4 shows a set of plots of the eigenvectors corresponding to the smallest eigenvalue at (a) $\gamma_{c-}$ and at (b)$\gamma_{c+}$, where $\gamma_{c+}$ is the strain immediately after the plastic event, and $\gamma_{c-}:=\gamma_{c+}-\Delta\gamma_{\rm Th}$ is the strain just before the event. As shown in Fig. 4, changes in eigenvectors owing to the stress drop event can be observed. Here, we find the existence of domains of grains of clockwise rotation and counter-clockwise rotation, and the collective motion of grains between two domains. We may observe the excitation of the quadrupole-like mode, although its structure is not sufficiently clear.

Figure 5 is the comparison of $\Ket{d\mathring{q}/d\gamma}$ obtained by the eigenvalue analysis (a) with that by the simulation (b) at $\gamma=\gamma_{c+}$, where $\Ket{d\mathring{q}/d\gamma}$ in the simulation is evaluated by $(\mathring{\bm{q}}(\gamma_{c+}+\Delta\gamma_{\textrm{Th}})-\mathring{\bm{q}}(\gamma_{c+}))/\Delta\gamma_{\textrm{Th}}$ with $\Delta\gamma_{\textrm{Th}}=1.0\times 10^{-8}$. It is obvious that the difference between the two figures is invisible, though the quadrupole-like structure cannot be clearly seen as in Fig. 4. Nevertheless, we can find the collective motion of grains in both figures.

Because we cannot use the eigenvalue analysis at the critical strain $\gamma_{c}$ for a plastic event, let us analyze the nonaffine displacement $\Delta\mathring{\bm{q}}$ between $\gamma_{c-}$ and $\gamma_{c+}$ caused by an avalanche using the simulation:

$\displaystyle\Delta\mathring{\bm{q}}|_{c}:=\left[\begin{matrix}\left.\Delta\mathring{\bm{q}}_{1}\right|_{c}\\ \left.\Delta\mathring{\bm{q}}_{2}\right|_{c}\\ \vdots\\ \left.\Delta\mathring{\bm{q}}_{N}\right|_{c}\end{matrix}\right],$

(52)

where

$\displaystyle\left.\Delta\mathring{\bm{q}}_{i}\right|_{c}:=\left[\begin{matrix}r_{i}^{\textrm{FB},x}(\gamma_{c+})-r_{i}^{\textrm{FB},x}(\gamma_{c-})-\Delta\gamma r_{i}^{\textrm{FB},y}(\gamma_{c-})\\ r_{i}^{\textrm{FB},y}(\gamma_{c+})-r_{i}^{\textrm{FB},y}(\gamma_{c-})\\ \ell_{i}^{\textrm{FB}}(\gamma_{c+})-\ell_{i}^{\textrm{FB}}(\gamma_{c-})\end{matrix}\right].$

(53)

Figure 6 shows a plot of the nonaffine displacement $\Delta\mathring{\bm{q}}|_{c}$ around a yielding point based on the simulation for $N=1024$ at $\Delta\gamma_{\rm Th}=1.0\times 10^{-8}$. This figure indicates that (i) grains move with rotations, which is one of the effects of mutual frictions between grains, and (ii) the existence of a quadrupole consisting of four domains of collective rotating grains exists. In particular, the rotations of the grains are sharply changed on the boundary between the domains. Unfortunately, $\Delta\mathring{\bm{q}}|_{c}$ cannot be described by the eigenvalue analysis, because $\Delta\mathring{\bm{q}}|_{c}$ expresses the configuration change, which is unstable for a small change in $\gamma$.

Figures 7 (a) and (b) show plots of the rigidity and smallest eigenvalue from $\gamma=0$ to $0.002$, which includes two plastic events based on the one-sample calculation of the collection of grains with $N=128$. One can find an almost perfect agreement of the rigidity between the eigenvalue analysis and simulation, except for the yielding points (see Fig. 7 (a)). We find discontinuous changes in the smallest eigenvalue at the yielding point, where the rigidity changes discontinuously (see Fig.  7 (b). As expected, the magnitude of the discontinuous change in rigidity at the yielding point in Fig. 7 (a) for $\gamma\approx 0.001$ is smaller than that for a point of a stress drop for $\gamma\approx 0.341$, as shown in Fig. 3.

Figure 8 shows the stress-strain curve corresponding to Fig. 7. It is difficult to find the plastic events in the main figure of Fig. 8, but we can find a small stress drop at this point if we use a close-up figure in the inset. We verify the creation and annihilation of contacting pairs at the stress drop points. The stress expression in Eq. (51) cannot be used at the yielding point; thus a disagreement exists between the eigenvalue analysis and simulation at the point in the inset of Fig. 8.

Figure 9 plots the rigidity of $G$ over 30 samples for $N=128$, where we have omitted the data if stress drop events take place. We verify that rigidity based on the eigenvalue analysis reproduces the results of the simulation. Note that non-monotonic changes in $G$ originate from changes in the contact points and configuration of grains.

In this paper, we have demonstrated that eigenvalue analysis of the Hessian matrix provides precise descriptions of the rigidity and stress of dispersed frictional grains in which the contact force is described by the harmonic potential, in spite of stress-drop events, such as stress avalanches. However, our model does not contain any slip processes between contacting grains. This success is a natural extension of the previous studies on frictionless grains Maloney04 ; Maloney06 to frictional grains and of our previous study on the linear response regime Ishima2022 to the nonlinear response regime. Two remarkable features of the contacting model are described by the harmonic potential. First, the tangential contact force in this model is no longer a history-dependent one. This leads to the significant simplification of the theoretical analysis. Second, unlike the naive expectation, the eigenvalues in our model do not indicate any precursors for the stress-drop events. In essence, stress-drop events take place suddenly by releasing contact points.

Some future tasks that need to be addressed are as follows. First, we need to consider the effect of slips, which causes a significant difference from our model because history-dependent contacts play important roles in the presence of slip events. Second, we must extend our analysis to nonlinear interacting models, such as the Hertzian contact model in a three-dimensional space. We plan on working on these points in the future.

The authors thank N. Oyama for the fruitful discussions and useful comments. This work was partially supported by a Grant-in-Aid from the MEXT for Scientific Research (Grant No. JP22K03459, JP21H01006, and JP19K03670) and a Grant-in-Aid from the Japan Society for Promotion of Science JSPS Research Fellow (Grant No. JP20J20292).