Noise Amplification in Frictional Systems: Oscillatory Instabilities

TL;DR

This paper investigates the oscillatory instabilities caused by complex eigenvalues in frictional granular systems, revealing how small perturbations can exponentially grow and lead to catastrophic events, with implications for earthquake physics.

Contribution

It explores the scaling laws, development scenarios, and effects of damping on oscillatory instabilities in frictional systems, advancing understanding of their dynamics and potential real-world relevance.

Findings

Giant amplification of perturbations can occur before catastrophic failure.

Instability onset follows specific scaling laws.

Damping influences the development and suppression of oscillatory instabilities.

Abstract

It was discovered recently that frictional granular materials can exhibit an important mechanism for instabilities, i.e the appearance of pairs of complex eigenvalues in their stability matrix. The consequence is an oscillatory exponential growth of small perturbations which are tamed by dynamical nonlinearities. The amplification can be giant, many orders of magnitude, and it ends with a catastrophic system-spanning plastic event. Here we follow up on this discovery, explore the scaling laws characterizing the onset of the instability, the scenarios of the development of the instability with and without damping, and the nature of the eventual system spanning events. The possible relevance to earthquake physics and to the transition from static to dynamic friction is discussed.

Equations

$$m_i\frac{d^2{\mathbf{r}}_i}{d{t}^2}={\mathbf{F}}_i({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots ,{\mathbf{r}}_N)\equiv -\frac{\partial U({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots {\mathbf{r}}_N)}{\partial {\mathbf{r}}_i}.$$

$${\mathbf{F}}_i\equiv -\frac{\partial U({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots {\mathbf{r}}_N)}{\partial {\mathbf{r}}_i}=0\text{in equilibrium}.$$

$$H_{ij}^{\alpha \beta }\equiv \frac{\partial F_i^{\alpha }}{\partial r_j^{\beta }}\equiv -\frac{{\partial }^{2}U({\mathbf{r}}_1,{\mathbf{r}}_2,\cdots {\mathbf{r}}_N)}{\partial r_i^{\alpha }\partial r_j^{\beta }}.$$

$${\mathbf{F}}_{ij}^{(n)}=k_n{\delta }_{ij}^{3/2}{\hat{r}}_{ij},{\hat{r}}_{ij}\equiv {\mathbf{r}}_{ij}/r_{ij}.$$

$$d{\mathbf{t}}_{ij}=d{\mathbf{r}}_{ij}-(d{\mathbf{r}}_{ij}\cdot {\mathbf{r}}_{ij}){\hat{r}}_{ij}+{\hat{r}}_{ij}\times ({\sigma }_{i}d{\mathbf{\theta }}_i+{\sigma }_{j}d{\mathbf{\theta }}_j).$$

$${\mathbf{F}}_{ij}^{(t)}=-k_t{\delta }_{ij}^{1/2}{t}_{ij}{\hat{t}}_{ij},$$

$${\mathbf{F}}_{ij}^{(t)}\le \mu {\mathbf{F}}_{ij}^{(n)},$$

$${\mathbf{F}}_{ij}^{(t)}=-k_t{\delta }_{ij}^{1/2}\left[1+\frac{t_{ij}}{t_{ij}^{\ast }}-{\left(\frac{t_{ij}}{t_{ij}^{\ast }}\right)}^2\right]t_{ij}{\hat{t}}_{ij},$$

$$t_{ij}^{\ast }\equiv \mu \frac{k_n}{k_t}{\delta }_{ij}.$$

$$m_i\frac{d^2{\mathbf{r}}_i}{d{t}^2}$$

$$I_i\frac{d^2{\mathbf{\theta }}_i}{d{t}^2}$$

$$J_{ij}^{\alpha \xi }\equiv \frac{\partial {\tilde{F}}_i^{\alpha }}{\partial q_j^{\xi }},{\tilde{\mathbf{F}}}_i\equiv \sum _j{\tilde{\mathbf{F}}}_{ij},$$

$$i{\omega }_{1,2}={\omega }_i\pm i{\omega }_r,i{\omega }_{3,4}=-{\omega }_i\pm i{\omega }_r,$$

$$\mathbf{q}(t)=\mathbf{q}(0)e^{{\omega }_{i}t}\sin ({\omega }_{r}t).$$

$$\begin{array}[]{l}{\ddot{x}_{1}+f_{1}(x_{1},x_{2},\gamma)=0}\\ {\ddot{x}_{2}+f_{2}(x_{1},x_{2},\gamma)=0}\end{array}$$

$$f_1({\tilde{x}}_1(\gamma ),{\tilde{x}}_2(\gamma ),\gamma )=0,$$

$$f_2({\tilde{x}}_1(\gamma ),{\tilde{x}}_2(\gamma ),\gamma )=0,\gamma \in ({\gamma }_1,{\gamma }_2).$$

$${\ddot{y}}_1+F_1(y_1,y_2,\gamma )=0$$

$${\ddot{y}}_2+F_2(y_1,y_2,\gamma )=0$$

$$F_k=f_k(y_1+{\tilde{x}}_1(\gamma ),y_2+{\tilde{x}}_2(\gamma ),\gamma ),k=1,2.$$

$$\begin{array}[]{l}{\ddot{y}_{1}+c_{11}(\gamma)y_{1}+c_{12}(\gamma)y_{2}+O(\left|y_{1}\right|^{2},\left|y_{2}\right|^{2},\left|y_{1}y_{2}\right|)=0}\\ {\ddot{y}_{2}+c_{21}(\gamma)y_{1}+c_{22}(\gamma)y_{2}+O(\left|y_{1}\right|^{2},\left|y_{2}\right|^{2},\left|y_{1}y_{2}\right|)=0}\\ {c_{kl}(\gamma)=\left.\frac{\partial F_{k}(y_{1},y_{2},\gamma)}{\partial y_{l}}\right|_{(y_{1},y_{2})=(0,0)}}\end{array}$$

$$(c_{11}(\gamma )-\lambda )(c_{22}(\gamma )-\lambda )-c_{12}(\gamma )c_{21}(\mu )=0,\lambda ={\omega }^2.$$

$${\lambda }_{1,2}=\frac{c_{11}(\gamma )+c_{22}(\gamma )\pm \sqrt{D(\gamma )}}{2},$$

$$D(\gamma )=(c_{11}(\gamma )-c_{22}(\gamma ){)}^2+4c_{12}(\gamma )c_{21}(\gamma ).$$

$${\lambda }_{1,2}=\frac{c_{11}({\gamma }_c)+c_{22}({\gamma }_c)\pm \sqrt{\alpha (\gamma -{\gamma }_c)}}{2}+O(\left|\gamma -{\gamma }_c\right|),$$

$$\displaystyle\alpha=D^{\prime}(\gamma)\Big{|}_{\gamma=\gamma_{c}}\ .$$

$${\lambda }_1-{\lambda }_2=\sqrt{\alpha (\gamma -{\gamma }_c)}+O(\left|\gamma -{\gamma }_c\right|),\gamma >{\gamma }_c$$

$$\Im ({\lambda }_1-{\lambda }_2)=\sqrt{\alpha |\gamma -{\gamma }_c|}+O(\left|\gamma -{\gamma }_c\right|),\gamma <{\gamma }_c.$$

$$|u^m⟩=|u_r^m⟩+i|u_i^m⟩.$$

$$\mathbf{J}|u^m⟩={\lambda }^m|u^m⟩$$

$$\mathbf{J}|u_r^m⟩$$

$$\mathbf{J}|u_i^m⟩$$

$$\mathbf{J}{\left[\begin{array}{cc}{\mathbf{u}}_r^m& {\mathbf{u}}_i^m\end{array}\right]}_{3N\times 2}={\left[\begin{array}{cc}{\mathbf{u}}_r^m& {\mathbf{u}}_i^m\end{array}\right]}_{3N\times 2}\times {\left[\begin{array}{cc}{\lambda }_r& {\lambda }_i\\ -{\lambda }_i& {\lambda }_r\end{array}\right]}_{2\times 2}.$$

$$\begin{gathered} {\left[\begin{array}{cc}{\ddot{\mathbf{u}}}^1& {\ddot{\mathbf{u}}}^2\end{array}\right]}_{3N\times 2}={\left[\begin{array}{cc}{\mathbf{u}}_r^m& {\mathbf{u}}_i^m\end{array}\right]}_{3N\times 2}\times {\left[\begin{array}{cc}{\lambda }_r& {\lambda }_i\\ -{\lambda }_i& {\lambda }_r\end{array}\right]}_{2\times 2} \\ . \end{gathered}$$

$$\ddot{\bf u}^{m}=\ddot{{\bf u}}^{1}+{\ddot{\bf u}}^{2}\$$

Full text

I Introduction

Understanding how granular materials fail is a problem of great geophysical importance due to its relevance to common phenomena like earthquakes and landslides. From a theoretical viewpoint, a micro-mechanical description of the failure process is difficult due to the absence of a Hamiltonian description of granular systems. To stress the importance of this, recall that starting from the pioneering works of Malandro and Lacks Malandro and Lacks (1998, 1999) plastic failures in athermal amorphous systems (with forces derivable from a Hamiltonian) has been identified as resulting from saddle node bifurcations occurring as a minimum of the energy landscape becomes a saddle point as the system deforms. In this approach, one considers a system of $N$ particles whose center of mass coordinates are ${\bm{r}}_{1},{\bm{r}}_{2},\cdots,{\bm{r}}_{N}$ at temperature $T=0$ which is endowed with a Hamiltonian $U({\bm{r}}_{1},{\bm{r}}_{2},\cdots{\bm{r}}_{N})$. The dynamics of the system is written as Newton’s equations of motion:

[TABLE]

When the system is in mechanical equilibrium the force ${\bm{F}}_{i}$ on each particle vanishes,

[TABLE]

The stability against perturbations of the equilibrium state is determined by the second derivative of the Hamiltonian which is the Hessian matrix:

[TABLE]

The Hessian matrix is evidently real and symmetric, and it has real eigenvalues which are all positive as long as the material is mechanically stable. Subjected to mechanical strain or stress the system may display a saddle node bifurcation with an eigenvalue going to zero. Generically this bifurcation results in the localization of the associated eigenfunction, in accordance with a local release of stress and energy Maloney and Lemaître (2004). The modes of the Hessian matrix which are associated with such instabilities are known as “plastic” or “soft” modes and their probability density function differs from the usual Debye density of stats in purely elastic materials Wyart, M. (2005); Karmakar et al. (2010). Work was devoted to understand the system size dependence of the eigenvalues of the Hessian Karmakar et al. (2010), their role in determining the mechanical characteristics like the elastic moduli Hentschel et al. (2011), the failure of nonlinear elasticity in such materials Hentschel et al. (2011); Procaccia et al. (2016); Dailidonis et al. (2017), and their relevance to shear banding and mechanical failure Dasgupta et al. (2012, 2013a, 2013b).

The dynamics of granular systems is governed by interparticle interaction forces, which have both normal and tangential components, obtained by coarse graining the microscopic degrees of freedom of the particles. A variety of different models exist Mindlin (1949); Walton and Braun (1986); Jacques Moreau (1994); Vu-Quoc and Zhang (1999); Vu-Quoc et al. (2000); Silbert et al. (2001), reflecting the different assumptions made in the derivation of the effective forces. In all cases, these forces cannot be expressed as derivatives of a Hamiltonian function with respect to the coarse grained degrees of freedom, which are the translational and the angular displacements of the particles. It is worth stressing that the lack of a Hamiltonian description is not exclusively related to the presence of viscous or hysteretic forces, but it is rather inherent in the coupling between the coarse grained tangential and normal forces. This implies that granular materials lack a Hamiltonian description even when dissipative forces may play a negligible role. As a consequence, the micro-mechanics of the failure process of granular systems cannot be understood using Hessian methods. Hence the question arises: how do granular materials fail?

Previous works investigating the loss of mechanical rigidity of granular assemblies have suggested that this could originate Nedderman (1992) from local events, related to the discontinuity of the frictional interaction between macroscopic objects. Indeed, the magnitude $|{\bm{F}}_{ij}^{(t)}|$ of the tangential force acting between two macroscopic objects $i$ and $j$ is bounded by $\mu|{\bm{F}}_{ij}^{(n)}|$, where ${\bm{F}}_{ij}^{(t)}$ is the normal interaction force and $\mu$ is Coulomb’s friction coefficient. Accordingly, contacts may reach their Coulomb threshold and start slipping as stress or stain is applied to a granular assembly. The failure of a contact could trigger subsequent rearrangements of other grains in a avalanche process Mcnamara et al. (2005); Welker and Mcnamara (2009), or lead to the emergence of a percolating cluster of unstable grains Cundall and Strack (1983); Staron et al. (2005, 2006), hence induce the macroscopic failure of the system. A complementary approach to rationalize the failure of granular systems has been devised building on the theory for the instability of elastic-plastic solids Hill (1958), and describes failure as a sharp increase of the kinetic energy of the system. When a granular system is slowly deformed, work is done on the system. Failure might occur when the system becomes unstable, meaning that the work done on the system, and possibly other work done by the system, is converted into kinetic energy. This so-called second-order work failure criterion Nicot et al. (2014) is formalized at the level of continuum equations, and it therefore not suitable to investigated the underlying microscopic features leading to the instability. Nevertheless, this criterion is indirectly related to the instability criterion developed for Hamiltonian systems, and hints towards the existence of collective failure modes. Support for the existence of a collective failure mechanism not initiated by the Coulomb failure of a contact also comes from detailed numerical simulations of model frictional granular particles, subject to an increasing in shear stress. These simulations have shown that a granular system may become unstable without any new contact reaching its Coulomb threshold Pica Ciamarra et al. (2010): only after the system becomes unstable and starts slipping new contacts are seen to reach their Coulomb limit. An indirect support towards the existence of collective failure mechanisms also stems from the existence of a linear regime in the response of granular systems to applied perturbations  Melosh (1996); Felzer and Brodsky (2006); Griffa et al. (2013); Giacco et al. (2015, 2018), which is not easily rationalized assuming failure to originate from an underlying discontinuous process such as the Coulomb failure of a contact.

In a recent Letter Chattoraj et al. (2019) we identified an instability mechanism in systems lacking a Hamiltonian description, such as frictional amorphous solids. This instability mechanisms involves an oscillating exponential growth of deviations from the state of mechanical equilibrium, and it is not related to the existence of a Coulomb threshold. It cannot exist in amorphous solids in which the inter-particle forces are derived from a Hamiltonian. These oscillatory instabilities furnish a micro-mechanical mechanism for a giant amplification of small perturbations that can lead to system spanning plastic events and mechanical failure. This physics was demonstrated in the context of amorphous assemblies of frictional disks, but the mechanism is generic for systems with friction. In this paper we follow up and discuss the phenomenon in greater detail.

The structure of the paper is as follows: in Section II we explicitly show that the coarse-grained description of the interparticle interaction of granular systems is not Hamiltonian, and introduce a modified interaction model that, while not-Hamiltonian, is at least differentiable. In Section III we review the novel instability mechanism for frictional systems Chattoraj et al. (2019), and then describe its critical features. Using the numerical model illustrated in Section IV, we test our theoretical predictions concerning the birth of the instability in Sec. V. We discuss the growth of the instability until nonlinear dynamics set in, culminating in a system spanning catastrophic event, in Sec. VI. Section VII provides information on the effect of damping. The upshot of the discussion is that a threshold damping frequency can be defined (based on the instability amplification frequency) below which the instability dynamics is unperturbed by damping. Sect. VIII raises up the important issue how generic is the instability discussed in this paper. Does its existence depend on the details of the coarse grained model used, or do we expect to appear generically in any models that employs a reduced set of coordinates like the positions of the centers of mass of the granules and their angular coordinates. We provide arguments for the generality of the phenomenon but propose that at this point in time experiments should be invoked as the final test. In Sec. IX we summarize the paper and discuss the possible connection of the present findings to remote triggering in earth quakes and to the transition from static to dynamic friction. The road ahead and future research are described.

II non Hamiltonian description of granular systems

Coarse grained descriptions of frictional amorphous solids, in which the effective degrees of freedom are the translational and the angular positions of the particles, do not admit a Hamiltonian description. To explicitly show that this is the case we recap here a popular and time honoured model used to describe the interaction between frictional granular particles, but the conclusion remains valid for other models.

The interaction between two particles, that we assume spherical for simplicity, has a normal and a tangential component. The normal force is determined by the overlap $\delta_{ij}\equiv\sigma_{i}+\sigma_{j}-r_{ij}$ between the particles, where ${\bm{r}}_{ij}\equiv{\bm{r}}_{i}-{\bm{r}}_{j}$ and $\sigma_{i}$ is the radius of particle $i$, and it is given by the Hertzian model,

[TABLE]

The tangential force is a function of the tangential displacement ${\bm{t}}_{ij}$ between the particles, a vector which is always orthogonal to $\hat{r}_{ij}$. Upon first contact between the particles, $t_{ij}=0$. Providing every particle with the angular coordinate ${\bm{\theta}}_{i}$ the change in tangential displacement is given by

[TABLE]

Accordingly, ${\bm{t}}_{ij}$ is obtained by integrating over time the relative velocity of the particles at the point of contact. In the Mindlin model, the tangential force depends on ${\bm{t}}_{ij}$ and on ${\bm{\delta}}_{ij}$ Mindlin (1949)

[TABLE]

and satisfy the Coulomb condition

[TABLE]

where $\mu$ is the friction coefficient. Usually this law is interpreted such that the tangential force reaches the limit abruptly, not analytically, thus not allowing to compute derivatives of the tangential force.

We can see now why the interaction is not Hamiltonian and the Hessian does not exist in this case. The first reason is somewhat trivial, stemming from the non-analyticity of the Coulomb law. This can be easily taken care of by smoothing out the Coulomb law such that the tangential force will have smooth derivatives; we choose:

[TABLE]

Now the derivative of the force with respect to $t_{ij}$ vanishes smoothly at $t_{ij}=t^{*}_{ij}$ and the Coulomb law Eq. (7) is fulfilled.

The second reason for the the loss of the Hessian matrix is not trivial at all Henkes et al. (2010). The time honored Hertz-Mindlin effective forces presented here, even for disks or balls, are not derivable from a potential, due to the coupling between the normal and the tangential displacement in the tangential force. Notice that this coupling is physical, as the tangential force depends on the normal force since it determines the contact area. This is the origin of the term $\delta_{ij}^{1/2}$ in Eqs. (6) and (8). It is easy to see that, because of this term, the derivative of the normal force with respect to $t$ does not equal the derivative of the tangential force with respect to $\delta$, which is what occurs when the forces are derived from a Hamiltonian $U(\delta,t)$.

III Oscillatory instabilities: theory

III.1 Stability matrix

While granular materials are not Hamiltonian, their dynamics is still Newtonian, with an extended set of coordinates ${\bm{q}}_{i}=\{{\bm{r}}_{i},{\bm{\theta}}_{i}\}$:

[TABLE]

where $m_{i}$ are masses for the center of mass coordinates, $I_{i}$ are moments of inertia for the angles, ${{\bm{F}}}_{i}$ are forces and ${{\bm{T}}}_{i}$ are torques, respectively. Using the smoothed out force Eq. (8), this allows to define the stability matrix which is an operator obtained from the derivatives of the force ${\bm{F}}_{i}$ and the torque ${\bm{T}}_{i}$ on each particle with respect to the coordinates. In other words

[TABLE]

where ${\bm{q}}_{j}$ stands for either a spatial position or a tangential coordinate, and $\tilde{{\bm{F}}}_{i}$ stands for either a force or a torque. Since in the usual case ${\bm{F}}_{i}=-\partial U/\partial{\bm{r}}_{i}$ we see that the operator ${\bm{J}}$ is an analog of the Hessian even when a Hamiltonian description is lacking. But with a huge difference: ${\bm{J}}$ is not a symmetric operator. Being real it can possess pairs of complex eigenvalues. When these appear, the system will exhibit oscillatory instabilities, since one of each complex pair will cause an oscillatory exponential divergence of any perturbation, and the other an oscillatory exponential decay. The actual calculation of the operator $J$ is somewhat cumbersome if conceptually straightforward. A detailed calculation for the present case of frictional disks interacting via Eqs. (4) and (8) is presented in Appendix A.

III.2 The oscillatory instability

When a pair of complex eigenvalues $\lambda_{1,2}=\lambda_{r}\pm i\lambda_{i}$ gets born, a novel instability mechanism develops. It should be stressed that the birth of a pair of complex eigenvalues is not a Hopf bifurcation. A pair of complex conjugate eigenvalues correspond to FOUR solutions $e^{i\omega t}$ to the linearized equation of motion with

[TABLE]

with $\omega_{r}\pm i\omega_{i}=\sqrt{\lambda_{r}\pm i\lambda_{i}}$. The first pair in Eq. (12) will induce an oscillatory motion with an exponential growth of any deviation ${\bm{q}}(0)$ from a state of mechanical equilibrium,

[TABLE]

The second pair represents an exponentially decaying oscillatory solution. The actual spatial dynamics that sets in due to this instability will be discussed below in VI.1.

III.3 Universal Scaling Laws

Noticing that the instability under discussion is somewhat unusual, we address the question how the coalescence of two real eigenvalues and the bifurcation of two imaginary parts take place as a function of the imposed strain. As noticed above, we need at least four degrees of freedom to have this instability, which means either four first order equations or two second order equations. Since we are solving Newton’s equations of motion we will discuss the instability in the latter form. Since the bifurcation is of co-dimension 1, it is enough to take into account a single parameter which we denote as $\gamma$ for obvious reasons. So, we consider the pair of equations:

[TABLE]

We assume that both functions $f_{1},{\rm\;}f_{2}$ are analytic with respect to all their arguments. Let us also assume that the system (14) has a family of equilibrium points $(\tilde{x}_{1}(\gamma),\tilde{x}_{2}(\gamma))$ for some interval of $\gamma$, so that

[TABLE]

Next define the new variables $y_{1}=x_{1}-\tilde{x}_{1}(\gamma),y_{2}=x_{2}-\tilde{x}_{2}(\gamma)$. Then, the system (14) is rewritten as:

[TABLE]

By construction Eqs. (16) possesses equilibria $(y_{1},y_{2})=(0,0)$ for $\gamma\in(\gamma_{1},\gamma_{2})$. Therefore, using analyticity, we rewrite equations (16) as follows:

[TABLE]

Removing the nonlinearities, and taking $y_{1}=y_{1}^{(0)}\exp(i\omega t),y_{2}=y_{2}^{(0)}\exp(i\omega t)$, one obtains the following equation for the eigenvalues:

[TABLE]

The solution is obvious:

[TABLE]

Finally we identify the critical point where the real eigenvalues coincide as $\gamma=\gamma_{c},\gamma_{c}\in(\gamma_{1},\gamma_{2})$. Then, it must be $D(\gamma_{c})=0$. In the vicinity of $\gamma_{c}$ we assume generic dependence of all functions on the parameters, and get the following:

[TABLE]

Without loss of generality, we can assume $\alpha>0$. Finally, for the eigenvalues near the bifurcation point we have:

[TABLE]

We thus conclude that the bifurcation is characterized by a square-root singularity for both the coalescence of the real eigenvalues and for the bifurcation of the imaginary parts. We verify this prediction in Sec.V.2.

IV Numerical model

To demonstrate the predictions presented above (which are based on a linearized stability matrix) we employ the following numerical simulations. In addition, we will use numerical simulations to investigate the evolution of the instability away from the linear regime, where no theoretical predictions are available, as well as the effect of damping.

We focus on a specific example of a binary assembly of $N$ frictional disks half of which with radius $\sigma_{1}=0.5$ and the other half with $\sigma_{2}=0.7$, unit mass $m=1$ and moment of inertia $I_{i}=0.5m_{i}\sigma_{i}^{2}$. The normal interaction between the grains is given by Eq. (4), while the tangential one is given by Eq. (8), with $k_{t}=2k_{n}/7$. We use $m$, $2\sigma_{1}$ and $\sqrt{m(2\sigma_{1})^{-1/2}k_{n}^{-1}}$ as our units of mass, length and time, respectively. We consider different values for the friction coefficient $\mu$ and the system size $N$.

The equation of motion are solved using two types of algorithms: “Newtonian” and “Overdamped”. The first is simply a solution of the Newton equations of motion with the given forces Eqs. (4) and (8). The second algorithm is solving the same equations of motion but with a damping force that is proportional to the velocities of the disks with a coefficient of proportionality $\eta_{v}=m\eta_{0}$. If not otherwise mentioned we use $\eta_{0}=2.2\times 10^{-2}$ expressed in reduced units. This value of $\eta_{0}$ ensures that the dynamics is overdamped as the damping timescale $\eta_{v}^{-1}$ is of the order of the time that sounds needs to travel one particle diameter. We use LAMMPS Plimpton (1995) to perform the numerical integration for these two algorithms, with integration timestep $10^{-5}\sqrt{(2\sigma_{1})^{1/2}k_{n}m^{-1}}$.

V Critical feature of the instability

V.1 Birth of complex conjugate pairs

To demonstrate the novel instability mechanism and describe its critical features, we start recapping the results of Ref. Chattoraj et al. (2019) with regard to the emergence of a pair of complex eigenvalues. An initial configuration is prepared by arranging binary particles randomly in a two dimensional box and then perform two consecutive runs of overdamped dynamics to bring the configuration at mechanical equilibrium. The initial configuration is prepared focussing on a frictionless system (i.e. $\mu=0$), and hence has no complex eigenvalues. Afterwards, we switch on friction, and perform athermal quasi static (AQS) simulations: starting from the initial stable configuration we shear the simulation box along the horizontal direction ($x$) by the amount $\delta\gamma$ and then we run the overdamped dynamics until the system reaches mechanical equilibrium. Operatively, we consider the system to be in mechanical equilibrium when the net force on each particle is less than $5\times 10^{-14}$. Here $\delta\gamma$ varies in the range $10^{-4}$ to $10^{-6}$ depending on the precision needed for the identification of the instability. After every AQS step we diagonalize the matrix ${\bm{J}}$ to find its eigenvalues. At some value of $\gamma$ we find for the first time the birth of conjugate pair of complex eigenvalues as seen in Fig. 1.

If we continue to increase the strain using the same protocol, we see the continuous emergence of other complex pairs, as well as the death of existing ones. We remark that we observe complex eigenvalues, which correspond to an unstable system, as we are using in this investigation an overdamped dynamics, which kills the growth of the instability. The effect of damping is discussed in Sec. VII

V.2 Numerical tests of the universal square-root singularity

The calculation in Sec. III.3 shows that the approach to the instability has critical features. Here we examine the numerical ramifications of this criticality. We note that when the strain is increased two real eigenvalues can collide to produce a complex conjugate pair, but the opposite can also happen as long as we increase the strain subject to overdamping; a pair of complex conjugate pair can give rise to two real eigenvalues. The square root singularity applies to both transitions.

A good indicator for the birth and presence of two complex conjugate eigenfunctions over some range of strain values is provided by the scalar product of the two colliding modes. The scalar product of two eigenvectors is defined as $\langle{u}^{m}|{u}^{n}\rangle=|\sum_{j=1}^{3N}(u_{r;j}^{m}+iu_{i;j}^{m})\cdot(u_{r;j}^{n}+iu_{i;j}^{n})|$. $\langle{u}^{m}|{u}^{n}\rangle$ varies in the interval $[0,1]$. At the critical value $\gamma_{c}$ the scalar product reaches unity. This is quite obvious since at the critical point also the eigenfunctions become complex conjugates. We show this first in Fig. 2 for a system of 10 disks and in Fig. 3 for 500 disks. Note that the eigenvectors of a real non-symmetric matrix are linearly independent, but in general not orthogonal.

A direct test of the universality predicted in the last subsection is provided by the difference between the two real eigenvalues. The data are consistent with $|\lambda_{r}^{1}-\lambda_{r}^{2}|\propto(\gamma_{c}-\gamma)^{x}$, and $x\simeq 0.5$. This is shown in Fig. 4.

VI The development of the instability

VI.1 Linear response regime

To investigate the development of an instability, we consider the Newtonian (no damping) evolution of a configuration having a pair of complex conjugate eigenvalues. As long as the deviations from mechanical equilibrium are sufficiently small, the eigenvalues and the eigenfunctions of ${\bm{J}}$ furnish an accurate prediction for the developing instability. To see this denote the complex mode $m$ of ${\bf J}$ whose complex eigenvalue $\lambda^{m}=\lambda_{r}+i\lambda_{i}$ as

[TABLE]

Since the dimension of ${\bf J}$ is $3N\times 3N$ then $|u^{m}_{r}\rangle$ and similarly $|u^{m}_{i}\rangle$ have $3N$ components. As

[TABLE]

we can rewrite the above equation in the form:

[TABLE]

In matrix notation this can be expressed as follows:

[TABLE]

The right hand side of the above equation produces a dynamic matrix of dimension $3N\times 2$:

[TABLE]

In particular, under the operation of this matrix, the resultant of $\ddot{{\bf u}}^{1}$ and $\ddot{{\bf u}}^{2}$, i.e.

[TABLE]

is subject to a rotation-scaling operation, and is therefore expected to predict a spiral trajectory with exponentially increasing speed. To demonstrate this, we consider a system in equilibrium and set the velocities of the particles (both translational and rotational) along the direction fixed by $\ddot{\bf u}^{m}$, and then follow the evolution of the system. Fig. 5 clearly shows that particles display the expected spiral motion.

The corresponding evolution of the mean square displacement of the center of mass $M(t)$ and the mean square change in the angular coordinate $A(t)$ can also be determined. Denoting

[TABLE]

we define

[TABLE]

According to Eq. (13) the system should behave as $M(t)\propto A(t)\propto e^{2\omega_{i}t}\sin^{2}(\omega_{r}t)$. As shown in Fig. 6 both quantities display sinusoidal motion since the beginning of the Newtonian dynamics with the expected frequency $\omega_{r}$ and the expected exponential growth $e^{2\omega_{i}t}$.

VI.2 The Nonlinear regime

When the instability develops sufficiently, the system exits from the linear response regime. It is of interest to consider how this happens, and what is the long time fate of the instability. We address this point following the evolution of a number of unstable configurations in which the linear instability is associated with complex eigenvalues with very different frequencies. It is noteworthy that while these initial configurations usually have more than a single pair of complex eigenfunctions the dynamics in the linear regime is dominated by the complex mode with the largest $\omega_{i}$, that defines the growth rate of the instability. We illustrate the evolution of these different systems in Fig. 7, which reports the time evolution of the mean square displacement (a), of the pressure (b), and of the shear stress (c). From these figure, we identify in the dynamics the following regimes:

at short times, the dominant effect of the oscillatory instabilities results in an exponential growth in $M(t)\sim\sin^{2}(\omega_{r}t)\exp(2\omega_{i}t)$. Similarly, as illustrated in the inset of Fig. 7b, the pressure $P$ oscillates around its initial value with the oscillation amplitude increasing exponentially as $P(t)-P(0)\sim\sin(\omega_{r}t)\exp(\omega_{i}t)$. The dashed lines in the inset illustrate the envelope $\pm\exp(\omega_{i}t)$, where $\omega_{i}\sim 2.3\times 10^{-5}$. The shear stress is also oscillating as the pressure around a constant value in the linear response regime, already shown explicitly in Chattoraj et al. (2019). 2. 2.

the exponential growth is interrupted when $M(t)\sim 10^{-7}$, which occurs at a characteristic time proportional to $\omega_{i}^{-1}$. At this time the two unstable complex conjugate modes annihilate, as a consequence of the breaking of few contacts (see below) and $M(t)$ enters a plateau regime. In this plateau regime, $P$ and $\sigma$ exhibit enhanced fluctuations, but no noticeable change in their mean values. 3. 3.

the plateau regime is interrupted by a fast growth of $M(t)$ driven by a second instability, which is triggered by the emergence of a number of complex eigenvalues. Afterwards, the system enters a diffusive regime. This second instability leads to a sharp increase in the pressure, and to a dramatic drop in shear stress $\sigma_{xy}$.

Fig. 8 illustrates how the development of the instability is related to the change in the spectrum. The inset follows the evolution of the imaginary part of the complex eigenvalues dominating at short time, and show that this vanishes when the exponential grow phase ends, for the considered $\omega_{i}\sim 2.3\times 10^{-5}$. The main panel follows the evolution of the imaginary part of all eigenvalues, and demonstrates that the secondary instability is triggered by the emergence of several unstable modes. The two consecutive vertical black dashed lines in Fig. 8 represent the end of the exponential growth and the beginning of the second instability.

The comparison of Fig. 8 and Figs. 7b,c clarifies that the effect of a dying complex mode might be eventful from a mechanical viewpoint. On the one hand it activates many new complex modes. On the other hand the shear stress $\sigma_{xy}$ drops almost to zero. Regarding the structure of the system, we observe in Fig. 9(a) that the end of the exponential growth is associated to the destruction of few contacts. In addition a few contacts achieve the Coulomb threshold (inset). A significant restructuring of the contact network only occurs during the second instability, which is associated with a drop in the number of contacts between grains. This drop leads to a significant change in the probability distribution of the number of contacts per particle, $P(z_{i})$ due to the emergence of particles with few contacts, as in Fig. 9(b).

We finally clarify that the sharp increase in the pressure observed during the second instability is associated to the fluidization of the system. To this end, we introduce a per-particle pressure defined as $P_{i}=\frac{1}{4L^{2}}\sum_{j=1}^{N}{\bf r}_{ij}\cdot{\bf F}_{ij}$ such that $P=\sum_{i=1}^{N}P_{i}$. Fig. 10a illustrates the distribution of $P_{ij}={\bf r}_{ij}\cdot{\bf F}_{ij}$ at time $t=0$, and at a later time after the secondary instability. Clearly, in the final state the pressure distribution has a longer tail, consistent with the observation of a larger average pressure value. Fig. 10b illustrates the average value of $P_{i}$ of the particles having $z$ contacts. In the final configuration, the pressure is higher for all values of $z$, as a consequence of the existence of large particle deformation associated to the flow of the system.

VII The effects of damping

While we have investigated the development of the oscillatory instability solving Newton’s equations with no dissipation, in most real system applications there is some form of damping. It can be due to the water and grit in a fault, to the interparticle interaction, or to additional friction with a confining substrate. It is therefore important to ask if and how damping affects the observed instabilities. In this section we show that damping tames but does not kill the instabilities. In particular, the instabilities associated with a small exponential growth rate (and hence emerging from the collision of low frequency modes) are more affected by the damping.

For our investigation we consider a starting configuration obtained from our athermal quasistatic simulations, which has a complex eigenpair with eigenfrequency $\omega_{r}\pm i\omega_{i}=0.3900643\pm i0.0000228$. We then evolve the system according to a Newtonian dynamics, where in Eq. (9) we add a damping force $-\tilde{m}_{i}\eta_{0}\dot{{\bm{r}}}_{i}$ on the right hand side.

Fig. 11 shows that the evolution of this damped dynamics depends on how the damping rate $\eta_{0}$ compares with the instability growth rate $\omega_{i}$, focussing on the time dependence of the mean square displacement (a) and of the pressure (P). For $\eta_{0}\gtrsim\omega_{i}$, the exponential growth is suppressed. For $\eta_{0}<\omega_{i}$, the dynamics has all the same features we have found in Fig. 7. As an example, in the inset of the lower panel we show that the pressure $P$ oscillates around its initial value while the oscillation amplitude increases exponentially as $P(t)-P(0)\sim\sin(\omega_{r}t)\exp(\omega_{i}t)$, as observed in the absence of damping in Fig. 7b. Overall, this analysis clarifies that damping suppresses the instabilities with a small growth rate $\omega_{i}$, not those corresponding to a large growth rate $\omega_{i}$.

VIII How generic is the oscillatory instability?

Having presented the oscillatory instability and its ability to self-amplify small perturbations, it is important to ask how generic is this instability. In particular, is the oscillatory instability a consequence of the structure of the Hertz-Mindlin model as interpreted for example by Cundall and Strack Cundall and Strack (1979), or is it expected to appear in any dynamical model of frictional disks. The answer to this question depends on what do we mean by “a model” of frictional disks. During the last few decades the granular community has accepted that assemblies of disks can be described by coarse-grained coordinates, for example the coordinates ${\bm{r}}_{i}$ of the centers of mass and $\theta_{i}$ for the angle of the disk (with reference to an initial angle $\theta_{i}(t=0)$). From the knowledge of these coordinates one writes down a physical model of the normal and tangential forces at each contact between neighboring disks. The other possibility is to track the exact micro-dynamics at each contact, paying attention to its leading and trailing edges, and taking into account the roughness of the contact and the plastic and elastic events that take place. This approach is not realistic for large assemblies of frictional disks, and even less so for arbitrarily shaped granules. Thus, if one follows the approach of coarse-grained coordinates, the options are somewhat limited. Imagine for example that we decide to use the harmonic-type normal interaction and delete $\sqrt{\delta_{ij}}$ term from the tangential force in Eq. (6). Then on the face of it, the forces could be derived from a Hamiltonian of the form

[TABLE]

Will this suffice to eliminate the oscillatory instability? The answer is no, as long as we do not remove the Coulomb law and the requirement of smoothness of forces to allow derivatives. Any smoothing of the type done in Eq. (8) brings back the dependence of the tangential force on the normal force and makes the system non-Hamiltonian.

To demonstrate this we have rerun our dynamics after removing the term $\sqrt{\delta_{ij}}$ from the tangential force in Eq. (6), but keeping the smoothing of this force at the Coulomb threshold. The resulting oscillatory instability is displayed in Fig. 12. As explained, the system still does not have a Hamiltonian structure, and the oscillatory instability remains generic. To drive home the message even further, we now remove the requirement of a smooth tangential force at the Coulomb threshold, allowing the tangential force to reach the limit in an abrupt way. Now the oscillatory instability disappears, since the Hamiltonian (31) remains now valid all the way to the Coulomb threshold. The price is of course that the Jacobian matrix ${\bm{J}}$ does not exist since the tangential force is not derivable at the threshold.

The above findings lead us to an interesting juncture for studying frictional properties in granular materials. The Hamiltonian models which oversimplify the role of deformation between two particles on static friction do not possess the oscillatory instabilities. In reality, the deformation between contacts is complex as found by various experiments and atomistic simulations, and a direct intrusion of overlap distance in the tangential force is the culmination of this fact. The overall dynamics becomes non-Hamiltonian that further displays oscillatory instabilities with a cost of no guarantee in energy conservation. At this point one depends on experiments to verify or refute the existence of this micromechanics in granular materials.

IX Summary and the Road Ahead

It appears from our discussion that oscillatory instabilities should be expected as a generic feature present in physical system whose dynamics are not derivable from a Hamiltonian. This novel instability was recently reported for frictional assemblies of disks Chattoraj et al. (2019). In this article, we have explored further the characteristics and the consequences of such instabilities on similar frictional granular systems.

We have shown that the birth and death of an oscillatory instability upon athermal quasi-static shear are a critical phenomenon. A pair of complex eigenpair gets born when two real eigenmodes align to one another, the relative difference between the two real eigenvalues vanishes with shear strain following a power law function with exponent 1/2. Similarly, the death of the complex eigenpair takes place when the alignment between the two eigenvectors is lost. The relative difference between the two real eigenvalues gradually increases following the same power law as a function of the strain.

With Newtonian dynamics the presence of a complex eigenpair generates instabilities where the oscillatory amplitudes increase exponentially. During this linear response regime particles move in spiral trajectories whose directions are determined by the two complex eigenvectors of the corresponding eigenpair. We have presented a formalism describing the direction of the spiral trajectories in the real plane; it is determined by the two complex eigenvectors. The linear response regime comes to an end when the oscillatory amplitude becomes so large that it breaks or makes a few pair contacts. This has a dramatic consequences; this triggers a second instability due to which the system goes through a structural change. The number of contacts of a particle becomes broadly distributed resulting in an increase of pressure. The dynamics becomes diffusive, helping the system to dissipate the shear stress.

The above phenomena are generic for a whole range of complex eigenvalues and even in the presence of damping force as long as the damping frequency is smaller than the imaginary eigenfrequency.

Probably the most important conclusion is that we have uncovered here a mechanism for a giant amplification of small perturbations. This finding may have further consequences in the context of a number of physics problems that involve frictional interactions. One fundamental question is how static friction turns into dynamical friction Urbakh et al. (2004); Braun et al. (2009); Trømborg et al. (2011). So far it is not quite clear what is the precise instability that allows two compressed interfaces to start moving with respect to one another. It is worthwhile in the coming future to examine whether the type of instability discussed above may be responsible for setting in this interesting transition. Another context of interest is that of remote triggering. It is known that one earthquake can induce another earthquake far away Felzer and Brodsky (2006); Brodsky et al. (2000); Brodsky and van der Elst (2014). Since strong seismic waves are strongly damped, only weak perturbation can travel a long distance. So the phenomenon of remote triggering requires a mechanism for self amplification. While we are fully aware that geological faults are very much more complex than assemblies of frictional disks, it appears highly worthwhile to explore the relevance of the kind of self amplification of small perturbation discussed above also to the problem of remote triggering. It is certainly on our agenda for the near future.

Acknowledgements.

This work had been supported in part by the ISF-Singapore exchange program and by the US-Israel Binational Science Foundation. We thank Jacques Zylberg and Yoav Pollack for useful discussions and exchanges at the early stages of this project. JC and MPC acknowledge NSCC Singapore for granting the computational facility under project 12000621.

Appendix A Calculation of the operator ${\bm{J}}$

The Jacobian operator ${{\bm{J}}}$, which is the dynamical response of the system, represents the derivative of the forces and of the torques acting on the particles with respect to all the degrees of freedom. The interaction forces used in this study is recalled in Sec. A.1, then the expression of the tangential displacement and its derivative (Sec. A.2) and the derivatives of forces and torques with respect to generalized coordinates (Sec. A.3) are stated. Finally, we show the expressions for all the components of ${{\bm{J}}}$ (Sec. A.4) and how these components are arranged as a matrix (Sec. A.5).

A.1 Interaction force

In our simulation, a pair of granular particles interacts when they overlap. The overlap distance $\delta_{ij}$ is measured as

[TABLE]

where $r_{ij}$ is the center-to-center distance of a pair-$i$ and $j$, and $\sigma_{i}$ is the radius of particle-$i$. The pair vector ${{\bm{r}}}_{ij}$ is defined as

[TABLE]

The pair-interaction force ${{\bm{F}}}_{ij}$ has two contributions. ${{\bm{F}}}^{(n)}_{ij}$ is the force acting along the normal direction of the pair $\hat{r}_{ij}$, and ${{\bm{F}}}^{(t)}_{ij}$ is the force acting along the tangential direction of the pair ${\hat{t}}_{ij}$. The normal force is Hertzian:

[TABLE]

where $k_{n}$ is the force constant with dimension: Force per length3/2. The tangential force ${{\bm{F}}}^{(t)}_{ij}$ is a function of both the overlap distance $\delta_{ij}$ and the tangential displacement ${{\bm{t}}}_{ij}$. We have modified the standard expression for ${{\bm{F}}}^{(t)}_{ij}$ and included a few higher order terms of $t_{ij}$ (i.e., $|{{\bm{t}}}_{ij}|$) such that the derivative of the force function $F^{(t)}_{ij}$ with respect to tangential distance $t_{ij}$ becomes continuous and it goes to zero smoothly. We use the following form:

[TABLE]

where $k_{t}$ is the tangential force constant. Its dimension is force per length3/2. $t^{*}_{ij}$ is the threshold tangential distance:

[TABLE]

where $\mu$ is the friction coefficient, a scalar quantity, which essentially determines the maximum strength of the tangential force with respect to the normal force at a fixed overlap $\delta_{ij}$. The derivative of $F^{(t)}_{ij}$ with respect to $t_{ij}$ vanishes at $t^{*}_{ij}$, as it turns out

[TABLE]

We stress here that the above forces imply a non Hamiltonian dynamics. That is, there is not a function $U(\delta,t)$ such that $F^{(n)}=-\frac{\partial U}{\partial\delta}$ and $F^{(t)}=-\frac{\partial U}{\partial t}$.

A.2 Tangential displacement

The tangential force is a function of both ${\bm{t}}_{ij}$ and ${{\bm{r}}_{ij}}$. The derivative of this force thus includes the derivative of the two latter quantities. Here we evaluate these derivatives using the chain rule.

The derivative of tangential displacement ${{\bm{t}}}_{ij}$ with respect to time $t$ is

[TABLE]

where ${{\bm{v}}}_{ij}={{\bm{v}}}_{i}-{{\bm{v}}}_{j}$ is the relative velocity of pair-$i$ and $j$. ${{\bm{v}}}^{n}_{ij}$ is the projection of ${{\bm{v}}}_{ij}$ along the normal direction $\hat{r}_{ij}$. ${{\bm{v}}}_{ij}-{{\bm{v}}}^{n}_{ij}$ is the tangential component of the relative velocity. ${{\bm{\omega}}}_{i}$ and ${{\bm{\omega}}}_{j}$ are the angular velocity of $i$ and $j$, respectively. In differential form, the above equation reads:

[TABLE]

where $\mathrm{d}{{\bm{\theta}}}_{i}$ is the angular displacement of $i$ which follows the relation: $\mathrm{d}{{\bm{\omega}}}_{i}=\frac{\mathrm{d}{{\bm{\theta}}}_{i}}{\mathrm{d}t}$.

Here on, we assume the two-dimensional ($\bf 2D$) system. Therefore, $\omega_{i}$, and so $\theta_{i}$, only have one component along $\hat{z}$, perpendicular to the xy plane, and $\hat{r}_{ij}\times\mathrm{d}{{\bm{\theta}}}_{i}=\mathrm{d}\theta_{i}(y_{ij}\hat{x}-x_{ij}\hat{y})/r_{ij}$. This allows to write Eq. (39) as

[TABLE]

where $\alpha$ and $\beta$ can take value 0 and 1 which correspond to x and y components, respectively. Now if particle-$i$ changes its position the angular displacement remains unaffected, i.e. $\frac{\mathrm{d}\theta_{i}}{\mathrm{d}r^{\alpha}_{i}}=0$. Thus, the change in tangential displacement along $\beta$ due to the change in position of particle-$i$ along $\alpha$ only contributes in translations, and it can be written as (using (40))

[TABLE]

where $\Delta_{\alpha\beta}$ is the Kronecker delta which is one when $\alpha=\beta$, or else zero. Similarly, a change in rotational coordinates does not modify the particles relative distance, i.e. $\frac{\mathrm{d}r^{\beta}_{ij}}{\mathrm{d}\theta_{i}}=0$. Thus, the change in tangential displacement along $\beta$ due to the change in $\theta_{i}$ is (from (40))

[TABLE]

In the above equation $\alpha$ and $\beta$ are always different. Now the magnitude of tangential distance $t_{ij}$ can be obtained from the relation $t^{2}_{ij}=\sum_{\alpha}{t^{\alpha}_{ij}}^{2}$. Its differential follows $\mathrm{d}t_{ij}=\sum_{\alpha}\frac{t^{\alpha}_{ij}}{t_{ij}}\mathrm{d}t^{\alpha}_{ij}$. The derivatives of tangential distance $t_{ij}$ with respect to $r^{\alpha}_{i}$ and $\theta_{i}$ can be expressed as

[TABLE]

With the help of equations (41) and (42) we can solve the above two differential equations. As the tangential threshold is a linear function of overlap distance $\delta_{ij}$ (see (36)), it also gets modified due to a change in $r^{\alpha}_{i}$ as

[TABLE]

and it is unaffected by the change in rotation, i.e. $\frac{\mathrm{d}t^{*}_{ij}}{\mathrm{d}\theta_{i}}=0$.

A.3 Evaluation of ${{\bm{J}}}$

The derivative of tangential force (equation (35)) with respect to $r^{\alpha}_{i}$:

[TABLE]

Here we use the notation $\tilde{t}$ to represent the ratio $t_{ij}/t^{*}_{ij}$, and the notation ${\tilde{t}}^{\beta}$ for ${t_{ij}}^{\beta}/t^{*}_{ij}$. The expressions for all the three partial differentiation in (LABEL:Eq:dFsbdra) are already shown in (42), (43), and (45).

Similarly, the derivative of tangential force with respect to $\theta_{i}$ (using the same notation as above) can be found as

[TABLE]

From the above two equations it is then understood that if ${{\bm{r}}}_{ij}$ and ${{\bm{t}}}_{ij}$ are known the differential equations can be solved easily. When ${\tilde{t}}^{\beta}$ is negligible for all $\beta$, then ${\tilde{t}}\approx 0$. This translates to $\frac{\partial{F^{(t)}_{ij}}^{\beta}}{\partial\theta_{i}}=-(-1)^{\beta}k_{t}\sigma_{i}\delta^{1/2}_{ij}\frac{r_{ij}^{\alpha}}{r_{ij}}$ with $\alpha\neq\beta$, implying that even in the case of zero tangential displacement and therefore, zero tangential force, the above derivative can be finite.

The derivative of normal force (equation (34)) with respect to $r^{\alpha}_{i}$:

[TABLE]

where $\Delta_{\alpha\beta}$ is the Kronecker delta. The derivative of total force which reads:

[TABLE]

can be solved using (48), (LABEL:Eq:dFsbdra), and (47).

The torque of particle-$j$ due to tangential force ${{\bm{F}}^{(t)}}_{ij}$ is ${{\bm{T}}}_{j}=-\sigma_{j}\left({\hat{r}}_{ij}\times{{\bm{F}}^{(t)}}_{ij}\right)\equiv\sigma_{j}\tilde{{\bm{T}}}_{ij}$. In 2D, $\tilde{{\bm{T}}}_{ij}$ has only z-component:

[TABLE]

The derivative of ${\tilde{T}}_{ij}^{z}$ then becomes:

[TABLE]

where $\delta_{\alpha x}$ (similarly, $\delta_{\alpha y}$) is the Kronecker delta, such that $\delta_{xx}=1$ and $\delta_{yx}=0$, and

[TABLE]

The above two differential equations can be solved using (LABEL:Eq:dFsbdra), and (47). If the tangential displacement $t_{ij}^{\beta}$ is negligible compared to the threshold $t_{ij}^{*}$, i.e., ${\tilde{t}}^{\beta}\approx 0$ for all $\beta$. This results in ${\tilde{t}}\approx 0$. Therefore, $\frac{\partial{\tilde{T}}_{ij}^{z}}{\partial\theta_{i}}=k_{t}\sigma_{i}\delta_{ij}^{1/2}$.

A.4 Jacobian

The dimension of Jacobian operator $J$ is force over length. To be consistent with the dimension we redefine the torque $T$ and rotational coordinate $\theta$ as

[TABLE]

In addition, the dynamic matrix has a contribution from the moment of inertia $I_{i}=I_{0}m_{i}\sigma_{i}^{2}$ as $\Delta{{\bm{\omega}}}_{i}={{\bm{T}}}_{i}/I_{i}\Delta t$. In our calculation, we assume that mass $m_{i}$ and $I_{0}$ both are one. The remaining contribution of $I_{i}$, i.e. $\sigma_{i}^{2}$, is taken care of by rescaling the torque and the angular displacement as ${\tilde{T}}_{i}$ and ${\tilde{\theta}_{i}}$ (54). For $I_{0}\neq 1$, the contribution of $I_{0}$ can be correctly anticipated if we rewrite (38) as below:

[TABLE]

$J$ essentially contains four different derivatives:

• •

First type: Derivative of force with respect to the position of particles:

[TABLE]

where $N$ is the total number of particles. $J^{\alpha\beta}_{ij}$ is symmetric if we change pairs, i.e.: $J^{\alpha\beta}_{ij}=J^{\alpha\beta}_{ji}$, however the symmetry is not guaranteed with the interchange of $\alpha$ and $\beta$.

• •

Second type: Derivative of force with respect to rotational coordinate:

[TABLE]

The negative sign makes sure that in stable systems all the eigenvalues are positive. $J^{\beta}_{ij}$ is asymmetric: $J^{\beta}_{ij}=-J^{\beta}_{ji}$.

• •

Third type: Derivative of torque with respect to position:

[TABLE]

$J^{\alpha}_{ij}$ is also asymmetric: $J^{\alpha}_{ij}=-J^{\alpha}_{ji}$.

• •

Fourth type: Derivative of torque with respect to rotational coordinate:

[TABLE]

The negative sign makes sure that in stable systems all the eigenvalues are positive. $J_{ij}$ is symmetric: $J_{ij}=J_{ji}$.

A.5 Arrangement of Jacobian matrix

In two dimension $D=2$, for $N$ particles the total number of elements in $J$ is $(D+1)N\times(D+1)N$. In the matrix, first $DN\times DN$ elements contain the first type of force derivative, i.e. $J^{\alpha\beta}_{ij}$. Here the row-index $ro$ and column-index $co$ of $J$ runs in the range $0\leq ro<DN$ and $0\leq co<DN$. Rows from $DN\leq ro<(D+1)N$ and columns $0\leq co<DN$ of $J$ contain $J^{\beta}_{ij}$, i.e., the second type of derivative. Rows from $0\leq ro<DN$ and columns $DN\leq co<(D+1)N$ of $J$ contain the third type $J^{\alpha}_{ij}$. Finally, rows from $DN\leq ro<(D+1)N$ and columns $DN\leq co<(D+1)N$ of $J$ hold $J_{ij}$, i.e., the fourth type of derivative. For a fixed type of derivative, at a fixed row, the column-index first runs over $j$ starting from 0 to $N-1$. Then $\beta$ is incremented, if it exists for that particular derivative type. Similarly, at a fixed column, row-index first runs over $i\in[0,N)$ and then $\alpha$ is incremented.