Theory of rigidity and numerical analysis of density of states of two-dimensional amorphous solids with dispersed frictional grains in the linear response regime

Amorphous materials consisting of dispersed grains such as powders, colloids, bubbles, and emulsions are ubiquitous in nature Jaeger96 ; Durian95 ; Wyart05 ; Baule18 . These materials behave like liquids at low densities and exhibit solid-like mechanical responses above their jamming point Liu98 . In systems consisting of frictionless grains, the rigidity changes continuously, but the coordination number of grains changes discontinuously at the jamming point as a function of density Durian95 ; Wyart05 ; Hern03 . The critical behavior near the jamming point is of interest to physicists as a non-equilibrium phase transition Hern01 ; Trappe01 ; Zhang05 ; Majmudar07 ; Ciamarra11 . Dispersed grains above the jamming point are fragile and exhibit softening and yielding transition under certain loads Nagamanasa14 ; Knowlton14 ; Kawasaki16 ; Leishangthem17 ; Bohy17 ; Ozawa18 ; Clark18 ; Boschan19 ; Singh20 ; Otsuki21 .

For amorphous solids consisting of frictionless grains, it is useful to analyze the dynamical matrix or the Hessian matrix, which is defined as the second derivative of the potential of a collection of grains with respect to the displacements from their stable configuration Wyart05 ; WyarLettt05 ; Ellenbroek06 ; Lerner16 ; Gartner16 ; Bonfanti19 ; Baule18 . For instance, the rigidity can be determined by eigenvalues and eigenvectors DeGiuli14 ; Mizuno16 ; Maloney04 ; Maloney06 ; Lemaitre06 ; Zaccone11 . It has been reported that the minimum nonzero eigenvalue of the Hessian matrix decreases with increasing strain and eventually becomes negative, where an irreversible stress drop takes place Maloney04 ; Maloney06 ; Manning11 ; Dasgupta12 ; Ebrahem20 .

It is known that amorphous solids have characteristic properties at low temperatures (e.g., thermal conductivity and specific heat) that are quite different from those of crystalline solids since a long-time ago Zeller71 . These days, we have recognized that amorphous solids consisting of dispersed grains exhibit unique elastic-plastic behavior as a mechanical response to an applied strain Nicolas18 . Because these properties are related to the density of states (DOS), there have been lots of studies on the DOS Hern03 ; Schirmacher07 ; Lerner16 ; Mizuno17 ; Wang19 . The DOS for systems composed of anisotropic grains, such as ellipses, dimers, deformable grains, and grains with rough surfaces, have been studied with the aid of the Hessian matrix Zeravcic09 ; Mailman09 ; Schreck10 ; Yunker11 ; Schreck12 ; Shiraishi19 ; Papanikolaou13 ; Britoa18 ; Treado21 ; Ikeda20 ; Ikeda21 . Because of the rotation of such anisotropic grains, there exists a rotational band in the DOS that is distinguishable from the translational band Mailman09 ; Yunker11 ; Schreck12 ; Britoa18 ; Ikeda20 ; Ikeda21 .

Even for systems of spherical grains that cannot be free from inter-particle friction, similar results are expected as a result of grain rotations. However, few studies have reported the existence of rotational bands in the DOS. Because the frictional force between the grains depends on the contact history, it cannot be expressed as a conservative force. Therefore, stability analysis for frictional grains based on the Hessian cannot be used. Nevertheless, the Hessian analysis using an effective potential for frictional grains has been performed Somfai07 ; Henkes10 . Recently, Liu et al. suggested that the Hessian analysis with another effective potential can be used even if slip processes exist Liu21 . The previous studies Somfai07 ; Henkes10 reported that friction between grains causes a continuous change in the functional form of the DOS from that of frictionless systems. However, there are only a few reports on whether an isolated band in the DOS originating from friction between grains is visible at lower frequencies.

Recently, Chattoraj et al. discussed the stability of the grain configuration under strain using the Jacobian matrix of frictional grains Chattoraj19 . They performed eigenvalue analysis under athermal quasi-static shear processes and determined the existence of oscillatory instability originating from inter-particle friction at a certain strain Chattoraj19 ; Chattoraj19E ; Charan20 . However, they did not discuss the rigidity or the DOS.

The theoretical determination of the rigidity of amorphous solids consisting of frictional grains is important for controlling amorphous solids. However, we do not know how to determine the rigidity from the Jacobian for the frictional grains.

The purpose of this study is to clarify the role of mutual friction between grains in terms of the rigidity and DOS. We focus on the response to an infinitesimal strain from a stable configuration of grains without any strain to obtain tangible results. In this study, we assume that there is no slip between grains because of an infinitesimal strain, and we then deal with friction as static friction.

The remainder of this paper is organized as follows. In the next section, we introduce the numerical method. In Sec. III, we introduce the Jacobian. Section IV consists of Sec IV.1, which deals with the theoretical prediction of rigidity in the linear response regime, and Sec. IV.2, which deals with the DOS. In the final section, we summarize the results of our study and discuss future work. The appendix consists of seven sections. In Appendix A, we summarize the method for preparing a stable grain configuration before applying shear. In Appendix B, we explain the implementation of the numerical integration method in the proposed system. In Appendix C, we summarize some properties of the Jacobian. In Appendix D, we present the explicit expressions of the Jacobian. In Appendix E, we investigate the effects of rattlers. In Appendix F, we write down the explicit results of Jacobian. In Appendix G, we derive the theoretical prediction of rigidity using the Jacobian. In Appendix H, we introduce the DOS using the Hessian analysis. In Appendix I, we investigate the system size dependence of the DOS. In Appendix J, we study the density dependence of the DOS.

Our system contains $N$ frictional spherical particles embedded in a monolayer configuration. We treat this system as a two-dimensional system (see Fig. 1).

To prevent the system from crystallizing Luding01 , we prepare an equal number of particles with diameters $d$ and $d/1.4$. We assume that the mass of particle $i$ is proportional to $d_{i}^{2}$, where $d_{i}$ is the diameter of $i$-th particle. We introduce $m$ as the mass of a particle with diameter $d$. In this study, $x_{i},y_{i}$, and $\theta_{i}$ denote $x,y$ coordinates, and the rotational angle of the $i$-th particle, respectively. We introduce the generalized coordinates of the $i$-th particle $\bm{q}_{i}:=(\bm{r}_{i}^{\textrm{T}},\ell_{i})^{\textrm{T}}$ with $\bm{q}_{i}:=(x_{i},y_{i})^{\textrm{T}}$ and $\ell_{i}:=d_{i}\theta_{i}/2$, where the superscript T denotes the transposition.

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

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

where $\bm{f}_{ij}$ and $T_{ij}$ are the force and $z$-component of the torque acting on the $i$-th particle from the $j$-th particle, respectively. Here, $T_{ij}$ is given by

where we have introduced the normal unit vector between $i$ and $j$ particles as $\bm{n}_{ij}:=\bm{r}_{ij}/|\bm{r}_{ij}|:=(\bm{r}_{i}-\bm{r}_{j})/|\bm{r}_{i}-\bm{r}_{j}|$. Here, $n_{ij}^{\zeta}$ and $f_{ij}^{\zeta}$ refer to $\zeta$-components of $\bm{n}_{ij}$ and $\bm{f}_{ij}$, respectively. Note that $\zeta$ expresses $\zeta=x$ or $y$ throughout this study. The force $\bm{f}_{ij}$ can be divided into normal $\bm{f}_{N,ij}$ and tangential $\bm{f}_{T,ij}$ parts as

where $d_{ij}:=d_{i}+d_{j}$ and $\Theta(x)$ is Heaviside’s step function, taking $\Theta(x)=1$ for $x>0$ and $\Theta(x)=0$ otherwise. We model the repulsive force between the contacted particles $i$ and $j$ as the Hertzian force in addition to the dissipative force proportional to the relative velocity with a damping constant $\eta_{D}$ Liu21 as follows:

where $k_{N}$ and $k_{T}$ are the stiffness parameters of normal and tangential contacts, respectively. For the normal compression length and its velocity are denoted as $\xi_{N,ij}:=d_{ij}/2-|\bm{r}_{ij}|$ and $\bm{v}_{N,ij}:=(\dot{\bm{r}}_{ij}\cdot\bm{n}_{ij})\bm{n}_{ij},$ respectively. For the tangential deformation, with the aid of $\bm{u}_{ij}:=(n_{ij}^{y},-n_{ij}^{x})^{\textrm{T}}$, the tangential velocity $\bm{v}_{T,ij}$ is defined as $\bm{v}_{T,ij}:=\dot{\bm{r}}_{ij}-\bm{v}_{N,ij}+\bm{u}_{ij}(d_{i}\omega_{i}+d_{j}\omega_{j})/2$, where we have introduced

with $\xi_{T,ij}:=|\bm{\xi}_{T,ij}|$ and $\bm{t}_{ij}:=-\bm{\xi}_{T,ij}/\xi_{T,ij}$. Here, $\dot{\bm{A}}:=d\bm{A}/dt$ and $\int_{C_{ij}}dt$ is the integration over the duration time of contact between $i$ and $j$ particles. Although the dissipative force between grains interacting with the Hertzian force is proportional to the product of the relative velocity and $\xi_{N,ij}^{1/2}$ Kuwabara87 ; Brilliantov96 ; Morgado97 , we adopt simple dissipative forces as in Eqs. (7) and (8) because we are not interested in the relaxation dynamics. We note that Eqs. (7) and (8) assume the Hertzian contact force for the static repulsion of contacting spheres, but all calculations in this study are those for two-dimensional systems. Here, we do not consider the effects of slips in the tangential equation of motion. This treatment can be justified if we restrict our interest in the linear response regime to a stable configuration of particles without any strain. This situation corresponds to frictional grains with infinitely large dynamical friction constant, in which the friction is only characterized by static friction. Therefore, our analysis does not apply to systems with finite strain Otsuki17 , where the effect of slip is important.

To generate a stable configuration of frictional particles, we prepare a stable configuration of frictionless particles in a square box of linear size $L$ using a fast inertial relaxation engine (FIRE) Bitzek06 . Subsequently, we turn on the tangential force using Eqs. (1) and (2) to achieve a stable configuration in the force-balanced (FB) state for frictional particle¹¹1 For simplicity, we prepare the configuration before applying shear for frictionless particles at first, and then considered the friction between particles. If we prepare a configuration before applying shear by compressing frictional particles, we confirm that the configuration had an oscillatory instability that resulted from the appearance of a pair of imaginary eigenvalues of the Jacobian divided by mass matrix: $\lambda^{\prime}=\lambda^{\prime}_{r}\pm i\lambda^{\prime}_{i}$ Chattoraj19 ; Chattoraj19E ; Charan20 , where $\lambda^{\prime},\lambda^{\prime}_{r}$, and $\lambda^{\prime}_{i}$ are the complex, real, and imaginary eigenvalues of $M^{-1}\mathcal{J}$, respectively. Here, $\mathcal{J}$ is the Jacobian defined in Eq. (18), and $M$ is the mass matrix whose explicit form is given by $M=\left[\begin{matrix}M_{1}&\\ &\ddots\\ &&M_{N}\end{matrix}\right]$, where $M_{i}:=\left[\begin{matrix}m_{i}&&\\ &m_{i}&\\ &&4I_{i}/d_{i}^{2}\end{matrix}\right]$. Because the linearized equation of motion is expressed as $d^{2}q_{i}^{\alpha}/dt^{2}=-\sum_{k,\kappa}\sum_{j,\beta}(M_{ik}^{\alpha\kappa})^{-1}\mathcal{J}_{kj}^{\kappa\beta}q_{j}^{\beta}$, there are four fundamental solutions $\bm{q}\propto e^{i\omega_{n}^{\prime}t}$, where $i\omega_{n}^{\prime}$ consists of $i\omega^{\prime}_{\pm 1}=\omega^{\prime}_{i}\pm i\omega^{\prime}_{r}$ and $i\omega^{\prime}_{\pm 2}=-\omega^{\prime}_{i}\pm i\omega^{\prime}_{r}$. Here, $\omega^{\prime}_{r}$ and $\omega^{\prime}_{i}$ satisfy the relation $\omega^{\prime}_{r}\pm i\omega^{\prime}_{i}=\sqrt{\lambda^{\prime}_{r}\pm i\lambda^{\prime}_{i}}$. Thus, to avoid the oscillatory instability of the configuration before applying shear, we adopt the protocol of creating the configuration with frictionless particles, and then let the system relax by adding static friction between particles. (see Appendix A for details). Here, the FB state satisfies the FB conditions $|F_{i}^{\zeta}|=0$ and $|T_{i}|=0$ for arbitrary particles. Note that we set $\theta_{i}=0$ when the tangential force is turned on.

We impose the Lees-Edwards boundary condition Lees72 ; Evans08 , where the direction parallel to the shear strain is $x$-direction. After applying a step strain $\Delta\gamma$ to all particles, $x$-coordinate of the position of the $i$-th particle is shifted by an affine displacement $\Delta x_{i}(\Delta\gamma):=\Delta\gamma y_{i}^{\textrm{FB}}(0)$, where the superscript FB denotes the FB state. The system is then relaxed to the FB state by the contact forces between the particles expressed in Eqs. (7) and (8). Here, $\zeta$-components of translational $\Delta\mathring{r}_{i}^{\zeta}(\Delta\gamma)$ and rotational $\Delta\mathring{\ell}_{i}(\Delta\gamma)$ nonaffine displacements of the $i$-th particle after the relaxation process are, respectively, defined as

Using Eqs. (10) and (11), we introduce the rate of nonaffine displacements as:

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 as $\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_{xy}(\gamma)$ at $\bm{q}^{\textrm{FB}}(\gamma)$ for one sample is given by:

The rigidity $g$ in the linear response regime for one sample is defined as

where the differentiation on the right-hand side (RHS) of Eq. (15) is defined as follows:

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

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

For the numerical FB condition, we use the condition $|\tilde{F}_{i}^{\alpha}|<F_{\textrm{Th}}$ for arbitrary $i$, where $F_{\textrm{Th}}$ is the threshold force for the simulation and $\tilde{\bm{F}}_{i}$ is the generalized force, defined as $\tilde{\bm{F}}_{i}:=(\tilde{F}_{i}^{x},\tilde{F}_{i}^{y},\tilde{F}_{i}^{\ell})^{\textrm{T}}:=(F_{i}^{x},F_{i}^{y},2T_{i}/d_{i})^{\textrm{T}}$.

In our simulation, we adopt $\eta_{D}=(mk_{N})^{1/2}d^{1/4}$ and $F_{\textrm{Th}}=1.0\times 10^{-14}k_{N}d^{3/2}$. The control parameters are the ratio of the tangential to normal stiffness $k_{T}/k_{N}$ and projected area fraction to two-dimensional space $\phi$. The operating ranges of $k_{T}/k_{N}$ and $\phi$ are $0.0$ to $10.0$ and $0.83$ to $0.90$, respectively. In this study, we mainly present the results for $N=4096$ and $\Delta\gamma=1.0\times 10^{-6}$ with the ensemble averages of $10$ samples for each $k_{T}/k_{N}$ and $\phi$. Some results are obtained with $N=1024$, $\Delta\gamma=1.0\times 10^{-6}$, and five samples for each $k_{T}/k_{N}$ and $\phi$. We verify that the results are independent of the choice of $\Delta\gamma$ for $1.0\times 10^{-6}\leq\Delta\gamma\leq 1.0\times 10^{-4}$. We ignore the effect of dissipation between particles because the velocity of each particle is sufficiently small for infinitesimal agitation from the FB state. The time step of the simulation, $\Delta t$, was set to $\Delta t=1.0\times 10^{-2}t_{0}$, and numerical integration was performed using the velocity Verlet method (see Appendix B), where $t_{0}:=(m/k_{N})^{1/2}d^{1/4}$. To obtain eigenvalues and eigenvectors of the Jacobian matrix, which will be introduced in detail in the next section, we have used the LAPACK, which provides a template library for linear algebra.

Figure 2 (a) shows an example of the affine displacements of particles, where the displacements exist only in the shear direction, and Fig. 2 (b) shows the nonaffine displacements.

In this section, we introduce the Jacobian, the DOS, and theoretical expressions of the linear rigidity. Here, we omit the effects of $\dot{\bm{q}}$ because the dissipative term proportional to $\dot{\bm{q}}$ vanishes under quasistatic shear.

In this section, we present the results of eigenvalue analysis and rigidity based on the formulation explained in the previous section. In Sec. IV.1, rigidity is determined using the eigenmodes of the Jacobian. Section IV.2 clarifies the effects of translational and rotational motions on the DOS.