# Non-backtracking eigenvector delocalization for random regular graphs

Xiangyi Zhu  Department of Mathematics, University of California Irvine,
Irvine, CA 92697<EMAIL_ADDRESS>and  Yizhe Zhu  Department of Mathematics,
University of California Irvine, Irvine, CA 92697<EMAIL_ADDRESS>
###### Abstract.

The non-backtracking operator of a graph is a powerful tool in spectral graph
theory and random matrix theory. Most existing results for the non-
backtracking operator of a random graph concern only eigenvalues or top
eigenvectors. In this paper, we take the first step in analyzing its bulk
eigenvector behaviors. We demonstrate that for the non-backtracking operator
$B$ of a random $d$-regular graph, its eigenvectors corresponding to
nontrivial eigenvalues are completely delocalized with high probability.
Additionally, we show complete delocalization for a reduced $2n\times 2n$ non-
backtracking matrix $\tilde{B}$. By projecting all eigenvalues of $\tilde{B}$
onto the real line, we obtain an empirical measure that converges weakly in
probability to the Kesten-McKay law for fixed $d\geq 3$ and to a semicircle
law as $d\to\infty$ with $n\to\infty$.

##  1\. Introduction

The non-backtracking operator is an important object in the study of spectral
graph theory [31, 4, 5]. It has recently been used as a powerful tool for
studying random matrices [12, 9, 8] and for designing efficient algorithms in
community detection and matrix completion [22, 12, 28, 11, 29, 30]. Many
recent results on the spectrum of the non-backtracking operator have been
established in various graph models, such as Erdős-Rényi graphs [33],
stochastic block models [12, 14], inhomogeneous random graphs [8, 18], random
regular graphs [9], and bipartite biregular graphs [13]. Very recently, some
of these results have also been extended to hypergraphs [17, 29]. However, so
far, these results mainly focus on top eigenvalues and eigenvectors and global
spectral distribution, while not much is known about bulk eigenvector
behavior.

This paper represents a first step towards understanding the bulk eigenvectors
of the non-backtracking operator in random graphs. Eigenvector delocalization
is an important topic in the study of random matrices, demonstrating that many
matrix models exhibit behavior similar to their Gaussian analog, a phenomenon
known as universality. A unit eigenvector is considered delocalized if its
infinity norm is close to that of a uniformly distributed random vector on the
unit sphere up to a polylog factor. Eigenvector delocalization has been shown
for many random matrix models, including non-Hermitian ones. However, for
sparse and non-Hermitian models, the literature includes only a few results,
such as those on eigenvector delocalization in random regular digraphs [23]
and Erdős-Rényi digraphs [19].

The non-backtracking operator $B$ for a random $d$-regular graph is a sparse
and non-Hermitian operator, where each row and column of $B$ has exactly
$(d-1)$ many nonzero entries. It can also be viewed as the adjacency matrix of
a random $(d-1)$-regular digraph with $(nd)$ many vertices. However, not all
$d$-regular digraphs can be obtained in this way (see [26, Theorem 2.3]), and
the eigenvalue distribution of $B$ is very different from the conjectured
oriented Kesten-McKay law for a uniformly chosen random $d$-regular digraph
[10]. A simulation of the non-backtracking spectrum of a random $d$-regular
graph is shown in Figure 1.

Figure 1. Simulation for eigenvalues of $\tilde{B}$ and $B$ for a random
$5$-regular graph with 200 vertices.

In Theorem 3.3, we show that with high probability, all eigenvectors
corresponding to nontrivial non-backtracking eigenvalues are completely
delocalized. The main idea behind the proof is to use the delocalization [21]
and spectral gap [21, 9] results for the adjacency matrix of a random
$d$-regular graph to understand eigenvector delocalization for the non-
backtracking matrix with the help of the Ihara-Bass formula [5]. The algebraic
structure of the non-backtracking operator enables the translation of
eigenvalue and eigenvector information from the corresponding adjacency
matrix. However, this precise algebraic connection between these two operators
is only available for regular graphs. Exploring beyond $d$-regular graphs is
an interesting future direction.

In the course of proving non-backtracking eigenvector delocalization, we show
in Theorem 3.1 the convergence in probability to the Kesten-McKay law and
semicircle law for fixed $d$ and growing $d$ when projecting the eigenvalues
of the non-backtracking operator onto the real line. From the Ihara-Bass
formula, another $2n\times 2n$ non-Hermitian block matrix $\tilde{B}$ closely
related to $B$ emerges, which we call a reduced non-backtracking operator [4,
33, 14, 27]. We also demonstrate that the eigenvectors of $\tilde{B}$ are
completely delocalized with high probability in Theorem 3.2.

### Organization of the paper

The rest of the paper is organized as follows. In Section 2, we provide some
basic notations and definitions for regular graphs and the non-backtracking
operator. In Section 3, we state the main results. Sections 4 to 6 contain all
the proofs. In Section 7, we discuss the open problem of generalizing the
results to Erdős-Rényi graphs.

Figure 2.  Simulation of a $4$-regular random graph with $n=200$

##  2\. Preliminaries

Let $G=(V,E)$ be a graph. $G$ is $d$-regular of size $n$ if each vertex has
degree $d$ and $|V|=n$. For given $n$ and $d$, we say $G$ is a random
$d$-regular graph if it is uniformly chosen from all $d$-regular graphs with
$n$ vertices, see Figure 2 for a simulation of a random regular graph.

The $(i,j)$-th entry of the adjacency matrix $A$ of a graph $G$ is defined as

$A_{ij}=\begin{cases}1&\text{if }\\{i,j\\}\in E,\\\
0&\text{otherwise.}\end{cases}$

The degree matrix $D$ of a graph $G$ is a diagonal matrix where
$D_{ii}=\sum_{j\in V}A_{ij}$. Define the oriented edge set $\vec{E}$ for $G$
as

(1) $\displaystyle\vec{E}=\\{(i,j):\\{i,j\\}\in E\\}.$

Each edge yields two oriented edges; therefore, $|\vec{E}|=2|E|$. The non-
backtracking operator $B$ of $G$ is a non-Hermitian operator of size
$|\vec{E}|\times|\vec{E}|$. For any $(u,v),(x,y)\in\vec{E}$, $B$ is defined as

$B_{(u,v),(x,y)}=\begin{cases}1&\text{$v=x$, $u\neq y$},\\\
0&\text{otherwise.}\end{cases}$

In particular, for a $d$-regular graph with $n$ vertices, the corresponding
$B$ is of size $nd\times nd$. A useful identity we will use in this paper is
the following Ihara-Bass formula [5].

######  Lemma 2.1 (Ihara-Bass formula).

For any graph $G=(V,E)$, and any $z\in\mathbb{C}$, the following identity
holds:

(2) $\displaystyle\det(B-zI)=(z^{2}-1)^{|E|-n}\det(z^{2}I-zA+D-I).$

Define a block matrix

(3) $\displaystyle\tilde{B}=\begin{bmatrix}0&D-I\\\
-I&A\end{bmatrix}\in\mathbb{R}^{2n\times 2n}.$

Then from (2), we have

(4) $\displaystyle\det(B-zI)=(z^{2}-1)^{|E|-n}\det(\tilde{B}-zI).$

The identity (4) implies that $B$ and $\tilde{B}$ share the same spectrum, up
to the multiplicity of trivial eigenvalues $\pm 1$. We also call $\tilde{B}$
the reduced non-backtracking matrix of $G$. When $G$ is a $d$-regular graph
with $n$ vertices, (2) can be further simplified to

(5) $\displaystyle\det(B-zI)=(z^{2}-1)^{nd/2-n}\det(z^{2}I-zA+(d-1)I).$

##  3\. Main results

Let $\mu_{1},\dots,\mu_{2n}$ be the eigenvalues of $\tilde{B}$. Define $x_{i}$
to be the real part of $\mu_{i}$. The empirical measure of $\\{x_{i},1\leq
i\leq 2n\\}$ is denoted as

$\mu=\frac{1}{2n}\sum_{i=1}^{2n}\delta_{x_{i}}.$

The following theorem characterizes the eigenvalue distribution of $\tilde{B}$
after projecting all eigenvalues to the real line.

######  Theorem 3.1 (Projecting eigenvalues of $\tilde{B}$ to the real line).

The following holds for the reduced non-backtracking matrix $\tilde{B}$ of a
uniformly chosen random $d$-regular graph:

  1. (1)

When $d\geq 3$ is a fixed integer, $\mu$ converges weakly in probability to a
rescaled Kesten-McKay distribution $\mu_{\mathrm{KM}}$ supported on
$[-\sqrt{d-1},\sqrt{d-1}]$, where

$\mu_{\mathrm{KM}}(x)=\frac{2d\sqrt{(d-1)-x^{2}}}{\pi(d^{2}-4x^{2})}~{}\mathbf{1}\left\\{|x|\leq\sqrt{d-1}\right\\}.$

  2. (2)

When $d\to\infty$ as $n\to\infty$, the empirical measure of
$\left\\{\frac{2x_{i}}{\sqrt{d-1}},i\in[2n]\right\\}$ converges weakly in
probability to a semicircle law $\mu_{\mathrm{SC}}$, where

$\mu_{\mathrm{SC}}(x)=\frac{1}{2\pi}\sqrt{4-x^{2}}~{}\mathbf{1}\\{|x|\leq
2\\}.$

Figure 3 is a simulation for the projected eigenvalues on the real line of a
random regular graph. Similar results have been established for Erdős-Rényi
graphs $G(n,p)$ with $np=\omega(\log n)$ [33] and for stochastic block models
[14] in the same regime.

Figure 3. Simulation for the empirical measure of the eigenvalue real parts of
$\tilde{B}$ for a random $5$-regular graph with $2000$ vertices, excluding the
deterministic eigenvalue $\mu_{1}=4$.

We now move on to study the eigenvectors of $B$ and $\tilde{B}$. Note that
each unit eigenvector of $\tilde{B}$ is in $\mathbb{C}^{2n}$. We show they are
completely delocalized based on the eigenvector delocalization results in [21]
for the adjacency matrix of a random $d$-regular graph.

######  Theorem 3.2 (Eigenvector delocalization for $\tilde{B}$).

Let $d\geq 3$ be fixed, and $\tilde{B}$ be the reduced non-backtracking matrix
of a random $d$-regular graph. Let $u_{i},i\in[2n]$ be the
$\ell_{2}$-normalized eigenvector associated with $\mu_{i}$. Then there exist
absolute constants $C_{1},C_{2}>0$ such that with probability at least
$1-n^{-C_{1}}$, for all $i\in[2n]$,

(6) $\displaystyle\|u_{i}\|_{\infty}\leq\frac{\log^{C_{2}}(n)}{\sqrt{n}}.$

By a different argument with the help of the Ihara-Bass formula, we can also
show that unit eigenvectors of $B$ as a vector in $\mathbb{C}^{nd}$ are
completely delocalized.

######  Theorem 3.3 (Eigenvector delocalization for $B$).

Let $d\geq 3$ be fixed, and $B$ be the non-backtracking matrix of a random
$d$-regular graph. Let $\mu_{i},i\in[2n]$ be the eigenvalues of $\tilde{B}$.
Then there exist absolute constants $C_{1},C_{2}>0$ such that with probability
at least $1-n^{-C_{1}}$, all the unit eigenvectors $w_{i}$ of $B$ associated
with $\mu_{i}\not=1$ satisfy

(7) $\displaystyle\|w_{i}\|_{\infty}\leq\frac{\log^{C_{2}}(n)}{\sqrt{nd}}.$

Figure 4. The infinite norms of eigenvectors for $A$, $\tilde{B}$, $B$ of the
graph shown in Figure 2.

The $\ell_{\infty}$ norms of eigenvectors of $B$ and $\tilde{B}$ are
illustrated in Figure 4. Similar results for non-backtracking eigenvectors of
random $d$-regular graph for growing $d$ can be proved in the same way by
using the eigenvector delocalization bounds for the adjacency matrix in [15,
32, 7, 6].

######  Remark 3.4 (Eigenvectors corresponding to eigenvalues $\pm 1$ in
Equation (5)).

In Equation (5), we see there are trivial eigenvalues $\pm 1$ of $B$ with an
extra multiplicity $\frac{nd}{2}-n$ which are not given by eigenvalues of
$\tilde{B}$. The structure of the eigenspaces of $\pm 1$ was discussed in [24,
Proof of Proposition 3.1].

##  4\. Proof of Theorem 3.1

We first derive an algebraic relation between eigenvalues and eigenvectors of
$A$ and $B$ when $G$ is a $d$-regular graph.

######  Lemma 4.1 (Spectral relation between $A$ and $\tilde{B}$).

Let $G$ be a $d$-regular graph with adjacency matrix $A$ and non-backtracking
matrix $B$ and $d\geq 2$. Let $v_{i}$ be an eigenvector of $A$ with respect to
an eigenvalue $\lambda_{i},i\in[n]$. Then

(8)
$\mu_{i},\mu_{i}^{\prime}=\frac{\lambda_{i}\pm\sqrt{{\lambda_{i}}^{2}-4(d-1)}}{2}$

are two eigenvalues of $\tilde{B}$, with two corresponding eigenvectors of the
form

(9) $u_{i}=\begin{bmatrix}v_{i}\\\ \frac{\mu_{i}}{d-1}v_{i}\end{bmatrix},\quad
u_{i}^{\prime}=\begin{bmatrix}v_{i}\\\
\frac{\mu_{i}^{\prime}}{d-1}v_{i}\end{bmatrix}.$

###### Proof.

Let $\begin{bmatrix}x\\\ y\end{bmatrix},x,y\in\mathbb{C}^{n}$ be an
eigenvector of $\tilde{B}$ with respect to an eigenvalue $\mu$, we have

$\mu\begin{bmatrix}x\\\ y\end{bmatrix}=\begin{bmatrix}0&(d-1)I\\\
-I&A\end{bmatrix}\begin{bmatrix}x\\\ y\end{bmatrix}.$

This equation implies

(10) $\displaystyle Ax=\left(\mu+\frac{d-1}{\mu}\right)x,\quad
y=\frac{\mu}{d-1}x.$

For each eigenvalue $\lambda_{i}$ of $A$ with a corresponding eigenvector
$v_{i}\in\mathbb{R}^{n}$, there are two corresponding eigenvalue
$\mu_{i},\mu^{\prime}_{i}$ of $\tilde{B}$ satisfying the quadratic equation

(11) $x^{2}-\lambda_{i}x+d-1=0.$

In particular, $\lambda_{1}=d$ gives two real eigenvalues $\mu_{1}=d-1$ and
$\mu_{1}^{\prime}=1$.

By counting multiplicity, Equation (11) gives all the $2n$ eigenvalues of
$\tilde{B}$. Thus, by Equations (10) and (11), there are $n$ eigenvectors of
$\tilde{B}$ of the form

(12) $u_{i}=\begin{bmatrix}v_{i}\\\
\frac{\mu_{i}}{d-1}v_{i}\end{bmatrix}\text{, where
}\mu_{i}=\frac{\lambda_{i}+\sqrt{{\lambda_{i}}^{2}-4(d-1)}}{2},$

and other $n$ eigenvectors of the form,

(13) $u_{i}^{\prime}=\begin{bmatrix}v_{i}\\\
\frac{\mu_{i}^{\prime}}{d-1}v_{i}\end{bmatrix}\text{, where
}\mu_{i}^{\prime}=\frac{\lambda_{i}-\sqrt{{\lambda_{i}}^{2}-4(d-1)}}{2}.$

∎

With Lemma 4.1, we are ready to prove Theorem 3.1 for random $d$-regular
graphs.

###### Proof of Theorem 3.1.

From Equation (11), when $|\lambda_{i}|^{2}\leq 4(d-1)$, we have

$\mu_{i}+\mu_{i}^{\prime}=\lambda_{i},\quad|\mu_{i}|=|\mu_{i}^{\prime}|=\sqrt{d-1},$

and the real parts of $\mu_{i},\mu_{i}^{\prime}$ satisfy

(14) $\displaystyle x_{i}=x_{i}^{\prime}=\frac{\lambda_{i}}{2}.$

The Kesten-McKay law in [25] shows that with high probability, the limiting
eigenvalue distribution of the adjacency matrix for a random $d$-regular graph
converges weakly to a probability measure with density

$f(x)=\frac{d\sqrt{4(d-1)-x^{2}}}{2\pi(d^{2}-x^{2})}\mathbf{1}\left\\{|x|\leq
2\sqrt{d-1}\right\\}.$

From the Kesten-McKay law above, with high probability, there are $n-o(n)$
many eigenvalues $\lambda_{i}$ satisfies $|\lambda_{i}|^{2}\leq 4(d-1)$.
Therefore, with high probability, there are $2n-o(n)$ eigenvalues of
$\tilde{B}$ on the circle of radius $\sqrt{d-1}$, and $n$ pairs of
$x_{i},x_{i}^{\prime}$ satisfies (14). Therefore, the empirical measure of the
real parts of $\mu_{i}$ converges weakly in probability to a rescaled Kesten-
McKay law given by

(15)
$\mu_{\mathrm{KM}}(x)=\frac{2d\sqrt{(d-1)-x^{2}}}{\pi(d^{2}-4x^{2})}\mathbf{1}\left\\{|x|\leq\sqrt{d-1}\right\\}.$

When $d\to\infty$ as $n\to\infty$, it is shown in [32, 15] that the empirical
spectral distribution of $\frac{A}{\sqrt{d-1}}$ converges weakly in
probability to a semicircle law. With Equation (14), following the same
argument as in the case for fixed $d$, the empirical measure
$\frac{1}{2n}\sum_{i=1}^{2n}\delta_{2x_{i}/\sqrt{d-1}}$ converges weakly in
probability to the semicircle law. ∎

##  5\. Proof of Theorem 3.2

Our proof is based on the following eigenvector delocalization bound from
[21].

######  Lemma 5.1 (Theorem 1.4 in [21]).

For any fixed $d\geq 3$, there exists absolute constants $C_{1},C_{2}>0$ such
that with probability at least $1-n^{-C_{1}}$, all unit eigenvectors $v_{i}$
of $A$ for a random $d$-regular graph satisfies

(16) $\displaystyle\|v_{i}\|_{\infty}\leq\frac{\log^{C_{2}}(n)}{\sqrt{n}}.$

With Lemma 5.1, we are able to prove Theorem 3.2.

###### Proof of Theorem 3.2.

We consider real and complex eigenvalues of $\tilde{B}$ separately.

Case 1: Real eigenvalues. There are two types of real eigenvalues for
$\tilde{B}$.

  1. (1)

$\lambda_{1}=d$ as a deterministic eigenvalue for $A$ yields two eigenvalues
for $\tilde{B}$: $\mu_{1}=d-1$ and $\mu_{1}^{\prime}=1$. The corresponding
eigenvector of $\lambda_{1}=1$ is $v_{1}=(1,\dots,1)^{\top}$, which gives two
corresponding eigenvectors for $\tilde{B}$ as

(17)
$\displaystyle\tilde{u}_{1}=[1,\dots,1]^{\top},\quad\tilde{u}_{1}^{\prime}=\left[1,\dots,1,\frac{1}{d-1},\dots,\frac{1}{d-1}\right]^{\top}\in\mathbb{R}^{2n}.$

The corresponding unit eigenvector then satisfies

(18)
$\displaystyle\|\tilde{u}_{1}\|_{\infty}=\frac{1}{\sqrt{2n}},\quad\|\tilde{u}_{1}\|_{\infty}=\frac{1}{\sqrt{\left(1+(d-1)^{-2}\right)n}}\leq\frac{1}{\sqrt{n}}.$

  2. (2)

Another type of eigenvalue is $\lambda_{i}$ of $A$ such that
$|\lambda_{i}|>2\sqrt{d-1}$. From Lemma 4.1, their corresponding eigenvalues
of $\tilde{B}$ are real. By the spectral gap bound from [9, Theorem 1], for
any eigenvalue $\mu$ of $\tilde{B}$ with $\mu\not=d-1$, we have with
probability at least $1-n^{-C}$,

(19) $\displaystyle|\mu|\leq\sqrt{d-1}+C\left(\frac{\log\log n}{\log
n}\right)^{2}.$

From (9), the corresponding eigenvectors $u_{i},u_{i}^{\prime}$ satisfies

$\|u_{i}\|_{\infty},\|u_{i}^{\prime}\|_{\infty}=\max\left\\{1,\frac{|\mu_{i}|}{d-1}\right\\}\|v_{i}\|_{\infty}.$

Therefore with the spectral gap bound (19) and eigenvector delocalization for
$v_{i}$ in Lemma 5.1, we obtain that with probability at least $1-n^{-C}$ for
some absolute constant $C>0$, the corresponding normalized eigenvectors
$\tilde{u}_{i},\tilde{u}_{i}^{\prime}$ satisfy

(20)
$\displaystyle\|u_{i}\|_{\infty},\|u_{i}^{\prime}\|_{\infty}\leq\frac{\log^{C_{2}}(n)}{\sqrt{n}}.$

Case 2: Complex eigenvalues. If $|\lambda_{i}|^{2}<4(d-1)$, the corresponding
complex eigenvalues $\mu_{i},\mu_{i}^{\prime}$ satisfies
$|\mu_{i}|=|\mu_{i}^{\prime}|=\sqrt{d-1}$. Then, with Lemma 5.1 and Equation
(9), the same bound (20) holds.

From the two cases above, Theorem 3.2 holds for all eigenvectors of
$\tilde{B}$. ∎

##  6\. Proof of Theorem 3.3

We first introduce a lemma contained in [24] for a deterministic relation
between eigenvectors of $B$ and eigenvectors of $A$ and provide a proof for
completeness.

######  Lemma 6.1 (Remark 3.4 in [24]).

Let $\vec{E}$ be oriented edge set for a $d$-regular graph $G=(V,E)$ defined
in (1). Let $\mu_{i},\mu_{i}^{\prime}$ denote each eigenvalue pair of
$\tilde{B}$ corresponding to eigenvalue $\lambda_{i}$ of A, where
$\lambda_{i}\neq 1$, through the quadratic equation (11). Let $v_{i}$ be a
unit eigenvector of $A$ associated $\lambda_{i}$. Define
$\tilde{w}_{i},\tilde{w}_{i}^{\prime}\in\mathbb{C}^{nd}$ such that

(21)
$\tilde{w}_{i}(x,y):=\mu_{i}v_{i}(y)-v_{i}(x),\quad\tilde{w}_{i}^{\prime}(x,y):=\mu_{i}^{\prime}v_{i}(y)-v_{i}(x).$

Then, each $\tilde{w}_{i},\tilde{w}_{i}^{\prime}$ is an eigenvector of B.

###### Proof.

It suffices to check $\tilde{w}_{i},i\in[n]$, and the same argument works for
$\tilde{w}_{i}^{\prime}$. We have for any $(x,y)\in\vec{E}$,
$(B\tilde{w}_{i})(x,y)$ presents an unique entry in
$B\tilde{w}_{i}\in\mathbb{C}^{nd}$. With (21), we are able to calculate

$\displaystyle(B\tilde{w}_{i})(x,y)$
$\displaystyle=\sum_{\begin{subarray}{c}z:(y,z)\in\vec{E}\\\ z\neq
x\end{subarray}}(\mu_{i}v_{i}(z)-v_{i}(y))$
$\displaystyle=\mu_{i}[(Av_{i})(y)-v_{i}(x)]-(d-1)v_{i}(y)$
$\displaystyle=[\mu_{i}\lambda_{i}-(d-1)]v_{i}(y)-\mu_{i}v_{i}(x)$
$\displaystyle=[\mu_{i}^{2}+d-1-(d-1)]v_{i}(y)-\mu_{i}v_{i}(x)$
$\displaystyle=\mu_{i}\tilde{w}_{i}(x,y),$

where in the fourth line $\mu_{i}\lambda_{i}$ is replaced by $\mu_{i}^{2}+d-1$
due to the quadratic equation (11). Then, we have
$B\tilde{w}_{i}=\mu_{i}\tilde{w}_{i}$, which implies $w_{i}$ is an eigenvalue
of $B$ with the corresponding eigenvalue $\mu_{i}$. ∎

We can derive a deterministic $\ell_{\infty}$ norm bound for eigenvectors of
$B$ constructed from Lemma 6.1 as follows.

######  Lemma 6.2.

Let $d>3$ be fixed, and let $\tilde{w}_{i},\tilde{w}_{i}^{\prime}$ be
eigenvectors of $B$ associated with $\mu_{i},\mu_{i}^{\prime}$ constructed in
Lemma 6.1 where $\mu_{i},\mu_{i}^{\prime}\not=1$. Then

(22)
$\displaystyle\frac{\|\tilde{w}_{i}\|_{\infty}}{\|\tilde{w}_{i}\|_{2}}\leq\frac{\|v_{i}\|_{\infty}(|\mu_{i}|+1)}{\sqrt{d}||\mu_{i}|-1|}.$

The same bound holds for $\tilde{w}_{i}^{\prime}$.

###### Proof.

It suffices to prove the estimate for $\tilde{w}_{i}$. From (21),

(23)
$\|\tilde{w}_{i}\|_{\infty}\leq\sup_{(x,y)\in\overrightarrow{E}}|\mu_{i}v_{i}(y)-v_{i}(x)|\leq\|v_{i}\|_{\infty}(|\mu_{i}|+1).$

On the other hand, since $\mu_{i},\mu_{i}^{\prime}$ is a conjugate pair,

(24) $\displaystyle\|\tilde{w}_{i}\|_{2}^{2}$
$\displaystyle=\sum_{(x,y)\in\overrightarrow{E}}(\mu_{i}v_{i}(y)-v_{i}(x))(\mu_{i}^{\prime}v_{i}(y)-v_{i}(x))$
(25)
$\displaystyle=|\mu_{i}|^{2}\sum_{(x,y)\in\overrightarrow{E}}(v_{i}(y))^{2}+\sum_{(x,y)\in\overrightarrow{E}}(v_{i}(x))^{2}-(\mu_{i}+\mu_{i}^{\prime})\sum_{(x,y)\in\overrightarrow{E}}v_{i}(y)v_{i}(x)$
(26)
$\displaystyle=d|\mu_{i}|^{2}\|v_{i}\|_{2}+d\|v_{i}\|_{2}-(\mu_{i}+\mu_{i}^{\prime})\sum_{(x,y)\in\overrightarrow{E}}v_{i}(y)v_{i}(x)$
(27)
$\displaystyle=d|\mu_{i}|^{2}+d-(\mu_{i}+\mu_{i}^{\prime})\sum_{(x,y)\in\overrightarrow{E}}v_{i}(y)v_{i}(x).$

where we use the fact that $v_{i}$ is a unit vector, and for each fixed $y$,
there are $d$ many oriented edges in $G$ ending with $y$. By Cauchy’s
inequality,

(28)
$\displaystyle(\mu_{i}+\mu_{i}^{\prime})\sum_{(x,y)\in\overrightarrow{E}}v_{i}(y)v_{i}(x)$
$\displaystyle\leq
2|\mu_{i}|\sqrt{\sum_{(x,y)\in\overrightarrow{E}}v_{i}(y)^{2}\sum_{(x,y)\in\overrightarrow{E}}v_{i}(x)^{2}}=2d|\mu_{i}|.$

This gives us a lower bound

(29) $\displaystyle\|\tilde{w}_{i}\|_{2}^{2}\geq
d|\mu_{i}|^{2}+d-2d|\mu_{i}|=d(|\mu_{i}|-1)^{2}.$

Inequalities (23) and (29) imply (22). ∎

Now, we are ready to prove Theorem 3.3.

###### Proof of Theorem 3.3.

We consider eigenvectors associated with real and complex eigenvalues of
$\tilde{B}$ separately.

Case 1: Real eigenvalues.

  1. (1)

The deterministic eigenvalue $\lambda_{1}=d$ of $A$ with eigenvector
$v_{1}=[1,\dots,1]^{\top}\in\mathbb{R}^{n}$ yields two eigenvalues
$\mu_{1}=d-1,\mu_{1}^{\prime}=1$. Let $\tilde{w}_{1}$ be the eigenvector given
by (21) associated with $\mu_{1}$. One can directly check that

(30)
$\displaystyle\frac{\|\tilde{w}_{1}\|_{\infty}}{\|\tilde{w}_{1}\|_{2}}=\frac{1}{\sqrt{nd}}.$

  2. (2)

For eigenvalues of $A$ such that $|\lambda_{i}|>2\sqrt{d-1}$, the
corresponding eigenvalues $\mu_{i},\mu_{i}^{\prime}$ of $B$ are real and
$\mu_{i}\mu_{i}^{\prime}=(d-1)$ from (14). By (19), with probability at least
$1-n^{-C}$,

(31) $\displaystyle\frac{d-1}{\sqrt{d-1}+C\left(\frac{\log\log n}{\log
n}\right)^{2}}\leq|\mu_{i}|\leq\sqrt{d-1}+C\left(\frac{\log\log n}{\log
n}\right)^{2}.$

Therefore with (22), with probability at least $1-n^{-C}$,

(32)
$\displaystyle\frac{\|\tilde{w}_{i}\|_{\infty}}{\|\tilde{w}_{i}\|_{2}}\leq\frac{C^{\prime}\|v_{i}\|_{\infty}}{\sqrt{d}}$

for an absolute constant $C^{\prime}>0$. The delocalization bound of $v_{i}$
in Lemma 5.1 implies that with probability at least $1-n^{-C_{1}}$,

(33)
$\displaystyle\frac{\|\tilde{w}_{i}\|_{\infty}}{\|\tilde{w}_{i}\|_{2}}\leq\frac{\log^{C_{2}}(n)}{\sqrt{nd}}.$

Case 2: Complex eigenvalues. If $|\lambda_{i}|^{2}<4(d-1)$, the corresponding
$\mu_{i},\mu_{i}^{\prime}$ satisfies
$|\mu_{i}|=|\mu_{i}^{\prime}|=\sqrt{d-1}$. From (22),

(34)
$\displaystyle\frac{\|\tilde{w}_{i}\|_{\infty}}{\|\tilde{w}_{i}\|_{2}}\leq\frac{\|v_{i}\|_{\infty}(\sqrt{d-1}+1)}{\sqrt{d}(\sqrt{d-1}-1)}\leq\frac{\sqrt{2}+1}{\sqrt{2}-1}\frac{\|v_{i}\|_{\infty}}{\sqrt{d}}.$

Therefore, (33) also holds due to Lemma 5.1.

With the two cases discussed above, taking $\ell_{2}$-normalized eigenvector
$w_{i}=\frac{\tilde{w}_{i}}{\|\tilde{w}_{i}\|_{2}}$ for $i\in[n]$ completes
the proof. ∎

##  7\. Erdős-Rényi graphs

Studying bulk non-backtracking eigenvectors for Erdős-Rényi graphs $G(n,p)$ is
an interesting open problem. Previously, using the fact that when
$np=\omega(\log n)$ the degrees in $G(n,p)$ are concentrated around $np$, [33,
14] studied the eigenvalues of $\tilde{B}$ by approximating the degree matrix
$D$ with $(np)I_{n}$, and analyzed a perturbed quadratic equation from the
Ihara-Bass formula. However, it is not clear to us how to generalize this
partial de-randomization argument introduced in [33] to study eigenvector
delocalization. We conjecture that non-backtracking eigenvector delocalization
holds when $np=\omega(\log n)$ and does not hold for $np=o(\log n)$, similar
to the results established for the adjacency matrix [16, 20, 2, 1, 3]. See
Figures 5 for simulation results.

### Acknowledgments

Y.Z. thanks Ludovic Stephan and Ke Wang for helpful discussions. X.Z. is
supported by an AMS Undergraduate Travel Grant to present this work at Joint
Mathematics Meetings 2024. Y.Z. was partially supported by the NSF-Simons
Research Collaborations on the Mathematical and Scientific Foundations of Deep
Learning. Part of this work was completed during Y.Z.’s visit to Université
Paris Cité, for which Y.Z. thanks Simon Coste for his warm hospitality.

Figure 5. Simulation for Erdős-Rényi graphs. The upper graph is a simulation
for $np=\log^{2}n$ with $n=300$; the lower graph is a simulation for $np=3$
with $n=2000$.

## References

  * [1] Johannes Alt, Raphael Ducatez, and Antti Knowles.  Delocalization transition for critical Erdős–Rényi graphs.  Communications in Mathematical Physics, 388(1):507–579, 2021.
  * [2] Johannes Alt, Raphaël Ducatez, and Antti Knowles.  The completely delocalized region of the Erdős-Rényi graph.  Electronic Communications in Probability, 27:1–9, 2022.
  * [3] Johannes Alt, Raphael Ducatez, and Antti Knowles.  Poisson statistics and localization at the spectral edge of sparse erdős–rényi graphs.  The Annals of Probability, 51(1):277–358, 2023.
  * [4] Omer Angel, Joel Friedman, and Shlomo Hoory.  The non-backtracking spectrum of the universal cover of a graph.  Transactions of the American Mathematical Society, 367(6):4287–4318, 2015.
  * [5] Hyman Bass.  The Ihara-Selberg Zeta function of a tree lattice.  International Journal of Mathematics, 03(06):717–797, December 1992\.
  * [6] Roland Bauerschmidt, Jiaoyang Huang, and Horng-Tzer Yau.  Local Kesten–McKay law for random regular graphs.  Communications in Mathematical Physics, 369:523–636, 2019.
  * [7] Roland Bauerschmidt, Antti Knowles, and Horng-Tzer Yau.  Local semicircle law for random regular graphs.  Communications on Pure and Applied Mathematics, 70(10):1898–1960, 2017.
  * [8] Florent Benaych-Georges, Charles Bordenave, and Antti Knowles.  Spectral radii of sparse random matrices.  Annales de l’Institut Henri Poincaré, Probabilités et Statistiques, 56(3):2141–2161, 2020.
  * [9] Charles Bordenave.  A new proof of Friedman’s second eigenvalue theorem and its extension to random lifts.  In Annales Scientifiques de l’École Normale Supérieure, volume 4, pages 1393–1439, 2020.
  * [10] Charles Bordenave and Djalil Chafai.  Around the circular law.  Probability Surveys, 9:1–89, 2012.
  * [11] Charles Bordenave, Simon Coste, and Raj Rao Nadakuditi.  Detection thresholds in very sparse matrix completion.  Foundations of Computational Mathematics, 23(5):1619–1743, 2023\.
  * [12] Charles Bordenave, Marc Lelarge, and Laurent Massoulié.  Nonbacktracking spectrum of random graphs: Community detection and nonregular ramanujan graphs.  Annals of probability, 46(1):1–71, 2018.
  * [13] Gerandy Brito, Ioana Dumitriu, and Kameron Decker Harris.  Spectral gap in random bipartite biregular graphs and applications.  Combinatorics, Probability and Computing, 31(2):229–267, 2022.
  * [14] Simon Coste and Yizhe Zhu.  Eigenvalues of the non-backtracking operator detached from the bulk.  Random Matrices: Theory and Applications, 10(03):2150028, 2021.
  * [15] Ioana Dumitriu and Soumik Pal.  Sparse regular random graphs: Spectral density and eigenvectors.  Annals of Probability, 40(5):2197–2235, 2012.
  * [16] Ioana Dumitriu and Yizhe Zhu.  Sparse general Wigner-type matrices: Local law and eigenvector delocalization.  Journal of Mathematical Physics, 60(2), 2019.
  * [17] Ioana Dumitriu and Yizhe Zhu.  Spectra of random regular hypergraphs.  Electronic Journal of Combinatorics, 28(3):P3.36, 2021.
  * [18] Ioana Dumitriu and Yizhe Zhu.  Extreme singular values of inhomogeneous sparse random rectangular matrices.  arXiv preprint arXiv:2209.12271, 2022.
  * [19] Yukun He.  Edge universality of sparse Erdős-Rényi digraphs.  arXiv preprint arXiv:2304.04723, 2023.
  * [20] Yukun He, Antti Knowles, and Matteo Marcozzi.  Local law and complete eigenvector delocalization for supercritical Erdős–Rényi graphs.  2019\.
  * [21] Jiaoyang Huang and Horng-Tzer Yau.  Spectrum of random $d$-regular graphs up to the edge.  arXiv preprint arXiv:2102.00963, 2021.
  * [22] Florent Krzakala, Cristopher Moore, Elchanan Mossel, Joe Neeman, Allan Sly, Lenka Zdeborová, and Pan Zhang.  Spectral redemption in clustering sparse networks.  Proceedings of the National Academy of Sciences, 110(52):20935–20940, 2013.
  * [23] Alexander Litvak, Anna Lytova, Konstantin Tikhomirov, Nicole Tomczak-Jaegermann, and Pierre Youssef.  Structure of eigenvectors of random regular digraphs.  Transactions of the American Mathematical Society, 371(11):8097–8172, 2019.
  * [24] Eyal Lubetzky and Yuval Peres.  Cutoff on all Ramanujan graphs.  Geometric and Functional Analysis, 26(4):1190–1216, 2016.
  * [25] Brendan D. McKay.  The expected eigenvalue distribution of a large regular graph.  Linear Algebra and its Applications, 40:203–216, 1981.
  * [26] Raffaella Mulas, Dong Zhang, and Giulio Zucal.  There is no going back: Properties of the non-backtracking laplacian.  Linear Algebra and its Applications, 680:341–370, 2024.
  * [27] Ludovic Stephan and Laurent Massoulié.  Robustness of spectral methods for community detection.  In Conference on Learning Theory, pages 2831–2860. PMLR, 2019.
  * [28] Ludovic Stephan and Laurent Massoulié.  Non-backtracking spectra of weighted inhomogeneous random graphs.  Mathematical Statistics and Learning, 5(3):201–271, 2022.
  * [29] Ludovic Stephan and Yizhe Zhu.  Sparse random hypergraphs: Non-backtracking spectra and community detection.  In 2022 IEEE 63rd Annual Symposium on Foundations of Computer Science (FOCS), pages 567–575. IEEE, 2022.
  * [30] Ludovic Stephan and Yizhe Zhu.  A non-backtracking method for long matrix and tensor completion.  arXiv preprint arXiv:2304.02077, 2023.
  * [31] Audrey Terras.  Zeta functions of graphs: A stroll through the garden, volume 128\.  Cambridge University Press, 2010.
  * [32] Linh V Tran, Van H Vu, and Ke Wang.  Sparse random graphs: Eigenvalues and eigenvectors.  Random Structures & Algorithms, 42(1):110–134, 2013.
  * [33] Ke Wang and Philip Matchett Wood.  Limiting empirical spectral distribution for the non-backtracking matrix of an erdős-rényi random graph.  Combinatorics, Probability and Computing, 32(6):956–973, 2023.


|  An exact expression of three-body system for the complex shear modulus of
frictional granular materials
---|---
|  Michio Otsukia and Hisao Hayakawab
|  We propose a simple model comprising three particles to study the nonlinear
mechanical response of jammed frictional granular materials under oscillatory
shear. Owing to the introduction of the simple model, we obtain an exact
analytical expression of the complex shear modulus for a system including many
monodispersed disks, which satisfies a scaling law in the vicinity of the
jamming point. These expressions perfectly reproduce the shear modulus of the
many-body system with low strain amplitudes and friction coefficients. Even
for disordered many-body systems, the model reproduces results by introducing
a single fitting parameter.

##

††footnotetext: a Graduate School of Engineering Science, Osaka University,
Toyonaka, Osaka 560-8531, Japan. E-mail:
<EMAIL_ADDRESS>b Yukawa Institute for Theoretical
Physics, Kyoto University, Kitashirakawaoiwake-cho, Sakyo-ku, Kyoto 606-8502,
Japan.

##  1 Introduction

The rheological property of densely dispersed grains, e.g., granular
materials, colloidal suspensions, and emulsions, plays an important role in
physics and engineering. This rheological property mainly depends on the
packing fraction $\phi$ of the grains. The materials behave like fluids for
$\phi<\phi_{\rm J}$ with jamming fraction $\phi_{\rm J}$ and exhibit a solid-
like elastic response above $\phi_{\rm J}$1, 2. In the linear response regime
(i.e., for small strains), the shear modulus is characterized by the density
of states 3, 4, 5 and satisfies scaling laws 6, 7, 8, 9. However, the linear
response region becomes narrower as $\phi$ approaches $\phi_{\rm J}$ 10, 11,
and the nonlinear response becomes relevant due to the plastic deformation
associated with the yielding 12, 13, 14, 15, 16, 17, 18, 19, 20.

If we are interested in a nonlinear response to an applied oscillatory shear
strain, it exhibits a complicated stress-strain curve. Although the storage
and loss moduli $G^{\prime}$ and $G^{\prime\prime}$ were originally introduced
to characterize the linear viscoelasticity of materials, they can use to
characterize nonlinear viscoelasticity or visco-elastoplastic responses to
applied strains 21. In this case, G’ and G" are no longer constants but
strongly depend on the strain amplitude $\gamma_{0}$. In particular, we have
recognized that $G^{\prime}$ decreases with $\gamma_{0}$ 11, 22, 23, 24, 25
and $G"$ remains non-zero in the low frequency limit24, 25 for densely
dispersed grains.

The theoretical analysis of densely dispersed grains is challenging as a
typical many-body problem in non-equilibrium systems. To date, a few
theoretical approaches have been proposed for systems related to frictionless
particles. The scaling laws for the linear elastic response were derived in
terms of the vibrational density of states 7, 8. The Fourier analysis of
particle trajectories helps to generate semi-analytical expressions for
$G^{\prime}$ and $G^{\prime\prime}$ 25. Unfortunately, these theories can not
apply to frictional particles because of the history-dependent contact force
9, 24.

It is helpful to analyze a simple model with small degrees of freedom to
understand the behavior of many-body systems, including densely dispersed
grains. This approach has been used in statistical mechanics. The mean-field
approximation of the Ising model is a typical example in which the system
contains only one Ising spin under the influence of a self-consistently
determined mean field 26. For atomic liquids, a cell model, in which a single
atom exists in a cage, was used to derive the equation of state 27, 28. The
coherent potential approximation for disordered solids has been used to
understand electronic band structures 29. The effective medium theory reveals
the elastic response of random spring networks 30. In addition, a simple model
consisting of two particles was proposed to reproduce the liquid-solid phase
transition 31. The advantage of such few-body models is that we can obtain
exact solutions. The qualitative behavior of the corresponding many-body
systems can be determined based on the solutions of the few-body models. Thus,
we adopt this approach to determine the nonlinear responses of the frictional
dispersed grains.

This study proposes a model consisting of three identical particles to
describe the mechanical response of jammed frictional granular materials under
oscillatory shear. In Section 2, we introduce the three-body system (TBS).
This model can be analytically solved for low-strain amplitudes and friction
coefficients near the jamming point in Section 3. In Section 4, we demonstrate
that the analytical solution reproduces the storage and loss moduli of many-
body systems (MBSs) without any fitting parameter if there is no disorder in
the particle configuration. Even if disorder exists, a scaling law for the
complex shear modulus for the TBS semi-quantitatively agrees with the
numerical simulations of the MBS by introducing a fitting parameter. We
discuss and conclude our results in Section 5. In Appendix A, we show the
details of the MBS when the particles are initially placed on a triangular
lattice. The effect of particle rotation is described in Appendix B. In
Appendix C, we derive the analytical expressions for the shear stress and
pressure in the TBS. In Appendix D, we relate the complex shear modulus with
the hysteresis loop of the stress–strain curve. The details of the disordered
MBS are presented in Appendix E. We present the numerical shear modulus for
the TBS in Appendix F.

##  2 Three-Body System

We consider two-dimensional granular materials consisting of many grains under
oscillatory shear (Fig. 1). Here, the grains constituting granular materials
are modeled as frictional spherical particles. Moreover, we introduce a system
of three identical particles to simply describe the MBS (Fig. 2). The MBS can
contain polydisperse particles, while we assume that the TBS is a monodisperse
system. In the TBS, the position $\bm{r}_{i}(t)=(x_{i}(t),y_{i}(t))$ of
particle $i$ with diameter $d$ at time $t$ is given by

$\displaystyle\bm{r}_{1}(t)$ $\displaystyle=$
$\displaystyle\left(\frac{\sqrt{3}\gamma(\theta(t))\ell}{4},\frac{\sqrt{3}\ell}{4}\right),$
(1) $\displaystyle\bm{r}_{2}(t)$ $\displaystyle=$
$\displaystyle\left(-\frac{\sqrt{3}\gamma(\theta(t))\ell}{4}-\frac{\ell}{2},-\frac{\sqrt{3}\ell}{4}\right),$
(2) $\displaystyle\bm{r}_{3}(t)$ $\displaystyle=$
$\displaystyle\left(-\frac{\sqrt{3}\gamma(\theta(t))\ell}{4}+\frac{\ell}{2},-\frac{\sqrt{3}\ell}{4}\right),$
(3)

where $\ell$ is the initial distance between particles. We also introduce
$\varepsilon:=1-\ell/d$ as the compressive strain. The compressive strain
$\epsilon$ in the TBS corresponds to $\phi-\phi_{\rm J}$ in the MBS as shown
in Appendix A. We apply shear strain as

$\gamma(\theta)=\gamma_{0}\sin\theta$ (4)

with strain amplitude $\gamma_{0}$, phase $\theta=\omega t$, and angular
frequency $\omega$. Note that we need at least three particles to realize a
stable interlocking state.

Fig. 1:  Schematics of the ordered MBS (a) and the disordered MBS (b).  Fig.
2:  A schematic of the TBS.

We adopt the interaction force $\bm{f}_{ij}$ between particles $i$ and $j$
given by

$\bm{f}_{ij}=\left(f_{ij}^{\rm(n)}\bm{n}_{ij}+f_{ij}^{\rm(t)}\bm{t}_{ij}\right)H(r_{ij}-d),$
(5)

where $f_{ij}^{\rm(n)}$ and $f_{ij}^{\rm(t)}$ denote the normal and tangential
forces between the particles $i$ and $j$ 32. The distance between the
particles $i$ and $j$ is $r_{ij}=|\bm{r}_{ij}|$ with
$\bm{r}_{ij}:=\bm{r}_{i}-\bm{r}_{j}=(x_{ij},y_{ij})$. Here, $H(x)$ is
Heviside’s step function satisfying $H(x)=1$ for $x>0$ and $H(x)=0$ otherwise.
The normal and tangential unit vectors are denoted by
$\bm{n}_{ij}:=\bm{r}_{ij}/r_{ij}=(n_{ij,x},n_{ij,y})$ and
$\bm{t}_{ij}:=(-n_{ij,y},n_{ij,x})$, respectively. For simplicity, we do not
consider the torque balance and, thus, the rotation of the particles. See
Appendix B for the effect of the rotation.

The normal force is assumed to be

$f_{ij}^{\rm(n)}=-k_{\rm n}u^{\rm(n)}_{ij}$ (6)

with the normal elastic constant $k_{\rm n}$ and normal relative displacement
$u^{\rm(n)}_{ij}:=r_{ij}-d$. Moreover, the tangential force is assumed to be

$f_{ij}^{\rm(t)}={\rm min}\left(|\tilde{f}_{ij}^{\rm(t)}|,\mu
f_{ij}^{\rm(n)}\right){\rm sgn}(\tilde{f}_{ij}^{\rm(t)}),$ (7)

where $\tilde{f}_{ij}^{\rm(t)}=-k_{\rm t}u_{ij}^{\rm(t)}$; $k_{\rm t}$ denotes
the tangential elastic constant, and $\mu$ denotes the friction coefficient.
Here, ${\rm min}(a,b)$ selects the smaller value between $a$ and $b$, ${\rm
sgn}(x)=1$ for $x\geq 0$, and ${\rm sgn}(x)=-1$ for $x<0$. The tangential
displacement $u_{ij}^{\rm(t)}$ satisfies
$\frac{d}{dt}u_{ij}^{\rm(t)}=v_{ij}^{\rm(t)}$ for
$|\tilde{f}_{ij}^{\rm(t)}|<\mu f_{ij}^{\rm(n)}$ with the tangential velocity
$v_{ij}^{\rm(t)}=(\frac{d}{dt}{\bm{r}}_{i}-\frac{d}{dt}{\bm{r}}_{j})\cdot\bm{t}_{ij}$,
whereas $u_{ij}^{\rm(t)}$ remains unchanged for
$|\tilde{f}_{ij}^{\rm(t)}|\geq\mu f_{ij}^{\rm(n)}$. We refer to the contact
with $|\tilde{f}_{ij}^{\rm(t)}|<\mu f_{ij}^{\rm(n)}$ as the stick contact and
the contact with $|\tilde{f}_{ij}^{\rm(t)}|\geq\mu f_{ij}^{\rm(n)}$ as the
slip contact. The tangential displacement, $u_{ij}^{\rm(t)}$, is initially set
to zero.

The (symmetric contact) shear stress is given by

$\sigma(\theta;\gamma_{0},\mu)=\sigma^{\rm(n)}(\theta;\gamma_{0},\mu)+\sigma^{\rm(t)}(\theta;\gamma_{0},\mu)$
(8)

with the normal component of $\sigma$

$\sigma^{\rm(n)}(\theta;\gamma_{0},\mu)=-\frac{1}{A}\sum_{i}\sum_{j>i}\frac{x_{ij}y_{ij}}{r_{ij}}f^{\rm(n)}_{ij}$
(9)

and tangential component of $\sigma$

$\sigma^{\rm(t)}(\theta;\gamma_{0},\mu)=-\frac{1}{2A}\sum_{i}\sum_{j>i}\frac{x_{ij}^{2}-y_{ij}^{2}}{r_{ij}}f^{\rm(t)}_{ij}.$
(10)

Here, $A$ corresponds to the area of the system, and we choose
$A=\sqrt{3}\ell^{2}/2$ as shown in Appendix A. The pressure is given by

$P(\theta;\gamma_{0},\mu)=\frac{1}{2A}\sum_{i}\sum_{j>i}(x_{ij}f_{ij,x}+y_{ij}f_{ij,y}).$
(11)

In the right-hand sides of eqns. (9)-(11), we have omitted the arguments
$\theta$, $\gamma_{0}$, and $\mu$. Similar abbreviations are used below. As we
are interested in quasistatic processes, we do not consider the kinetic parts
of $\sigma$ and $P$ and the dependence on $\omega$. After several cycles of
oscillatory shear, $\sigma(\theta)$ becomes periodic. The storage and loss
moduli are given by 33

$\displaystyle G^{\prime}(\gamma_{0},\mu)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}d\theta\
\sigma(\theta;\gamma_{0},\mu)\sin\theta/\gamma_{0},$ (12) $\displaystyle
G^{\prime\prime}(\gamma_{0},\mu)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\int_{0}^{2\pi}d\theta\
\sigma(\theta;\gamma_{0},\mu)\cos\theta/\gamma_{0}.$ (13)

##  3 Theoretical analysis

Assuming $\gamma_{0}\ll\epsilon\ll 1$, we analytically obtain $G^{\prime}$ and
$G^{\prime\prime}$ for the TBS. The derivation of the analytical results can
be found in Appendix C.

First, the normal component of the shear stress is given by

$\displaystyle\sigma^{\rm(n)}(\theta)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}k_{\rm n}\gamma(\theta)}{4}.$ (14)

The tangential component of the shear stress is given by

$\displaystyle\sigma^{\rm(t)}(\theta)$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}k_{\rm t}\gamma(\theta)}{4}$ (15)

for $\gamma_{0}<\gamma_{c}(\mu)$ with a critical amplitude

$\gamma_{c}(\mu)=\frac{4\mu k_{\rm n}\epsilon}{3k_{\rm t}},$ (16)

which characterizes the transition from stick to slip states in the contact
between the particles. For $\gamma_{0}\geq\gamma_{c}(\mu)$, the tangential
component of the shear stress is given by

$\displaystyle\sigma^{\rm(t)}(\theta)=\left\\{\begin{array}[]{ll}\dfrac{\mu
k_{\rm n}\epsilon}{\sqrt{3}},&0\leq\theta<\dfrac{\pi}{2}\\\ \dfrac{\mu k_{\rm
n}\epsilon}{\sqrt{3}}+\dfrac{\sqrt{3}k_{\rm
t}(\gamma(\theta)-\gamma_{0})}{4},&\dfrac{\pi}{2}\leq\theta<\dfrac{\pi}{2}+\Theta\\\
-\dfrac{\mu k_{\rm
n}\epsilon}{\sqrt{3}},&\dfrac{\pi}{2}+\Theta\leq\theta<\dfrac{3\pi}{2}\\\
-\dfrac{\mu k_{\rm n}\epsilon}{\sqrt{3}}+\dfrac{\sqrt{3}k_{\rm
t}(\gamma(\theta)+\gamma_{0})}{4},&\dfrac{3\pi}{2}\leq\theta<\dfrac{3\pi}{2}+\Theta\\\
\dfrac{\mu k_{\rm
n}\epsilon}{\sqrt{3}},&\dfrac{3\pi}{2}+\Theta\leq\theta<2\pi,\end{array}\right.$
(22)

where $\Theta=\cos^{-1}\left(1-2\gamma_{c}(\mu)/\gamma_{0}\right)$. Regions
with $\frac{\pi}{2}\leq\theta<\frac{\pi}{2}+\Theta$ and
$\frac{3\pi}{2}\leq\theta<\frac{3\pi}{2}+\Theta$ correspond to the stick
state, and the other regions correspond to the slip state. Owing to this
transition in the contact, the stress–strain curve given by eqns. (14)–(22)
exhibits a hysteresis loop. Equation (22) does not exhibit a viscoelastic
response but a typical elastoplastic response without viscous effect.

Figure 3 shows the scaled shear stress $\sigma/\gamma_{0}$ against the scaled
strain $\gamma/\gamma_{0}$ using eqns. (4), (8), and (14)-(22) for various
values of $\gamma_{0}$ with $k_{\rm t}/k_{\rm n}=1.0$ and $\mu=0.01$. The
shape of the scaled stress–strain curve is characterized by a parallelogram as
a typical elastoplastic response. As $\gamma_{0}$ increases, the maximum value
$\tilde{\sigma}_{\rm max}=(\sigma/\gamma_{0})|_{\gamma/\gamma_{0}=1}$
decreases from a larger value $\sqrt{3}(k_{\rm n}+k_{\rm t})/4$ to a smaller
value $\sqrt{3}k_{\rm n}/4$. As shown in Appendix D, the storage modulus
$G^{\prime}$ is approximately given by $\tilde{\sigma}_{\rm max}$. Hence, the
decrease of $\tilde{\sigma}_{\rm max}$ in Fig. 3 indicates the decrease of
$G^{\prime}$. For $\gamma_{0}=0.00003$ and $0.0001$, the hysteresis loop
exists, but the area of the loop is negligible for $\gamma_{0}=0.00001$ and
$0.001$. The loss modulus $G^{\prime\prime}$ is proportional to the area of
the loop as shown in Appendix D. Hence, the dependence of the area on
$\gamma_{0}$ indicates that there is a peak in $G^{\prime\prime}$ as
$\gamma_{0}$ increases.

Fig. 3:  Scaled shear stress $\sigma/\gamma_{0}$ against $\gamma/\gamma_{0}$
using eqns. (4), (8), and (14)-(22) for various values of $\gamma_{0}$ with
$k_{\rm t}/k_{\rm n}=1.0$, $\epsilon=0.001$, and $\mu=0.01$.

Substituting eqns. (8) and (14)-(22) into eqn. (12), we obtain the storage
modulus as

$\displaystyle
G^{\prime}=\left\\{\begin{array}[]{ll}\dfrac{\sqrt{3}\left(k_{\rm n}+k_{\rm
t}\right)}{4},&\gamma_{0}\leq\gamma_{c}(\mu)\\\
\dfrac{\sqrt{3}}{4}\left\\{k_{\rm n}+\dfrac{k_{\rm
t}}{\pi}\left(\Theta-\sin\Theta\cos\Theta\right)\right\\},&\gamma_{0}>\gamma_{c}(\mu).\end{array}\right.$
(25)

As $\gamma_{0}$ increases beyond $\gamma_{c}(\mu)$, $G^{\prime}$ decreases
from a higher value to a lower value. The corresponding behavior has been
observed in the MBS in previous studies 9, 24.

Substituting eqns. (8) and (14)-(22) into eqn. (13), the loss modulus is given
by

$\displaystyle
G^{\prime\prime}=\left\\{\begin{array}[]{ll}0,&\gamma_{0}\leq\gamma_{c}(\mu)\\\
\dfrac{\sqrt{3}k_{\rm
t}}{4\pi}\left(1-\cos^{2}\Theta\right),&\gamma_{0}>\gamma_{c}(\mu).\end{array}\right.$
(28)

The loss modulus $G^{\prime\prime}$ is zero for $\gamma_{0}<\gamma_{c}(\mu)$,
whereas $G^{\prime\prime}$ sharply increases with $\gamma_{0}$ when
$\gamma_{0}$ exceeds $\gamma_{c}(\mu)$ and decreases to $0$ after a peak. The
behavior of $G^{\prime\prime}$ for the TBS qualitatively reproduces that of
the MBS in previous studies 24.

We adopt the abbreviation for the pressure at $\gamma=0$ as

$P_{0}:=P(\theta=0;\gamma_{0},\mu),$ (29)

which is also obtained as

$P_{0}=\sqrt{3}k_{\rm n}\epsilon.$ (30)

From eqns. (16), (25), (28), and (30), we derive scaling laws for a given
$\epsilon$ as

$\displaystyle G^{\prime}(\mu,\gamma_{0})=G^{\prime}_{\rm
M}(\mu)\mathcal{G}^{\prime}\left(\frac{k_{\rm t}\gamma_{0}}{\mu
P_{0}(\gamma_{0},\mu)}\right),$ (31) $\displaystyle
G^{\prime\prime}(\mu,\gamma_{0})=G^{\prime\prime}_{\rm
M}(\mu)\mathcal{G}^{\prime\prime}\left(\frac{k_{\rm t}\gamma_{0}}{\mu
P_{0}(\gamma_{0},\mu)}\right),$ (32)

where $\mathcal{G}^{\prime}(x)$ and $\mathcal{G}^{\prime\prime}(x)$ denote
scaling functions. The maximum values of $G^{\prime}$ and $G^{\prime\prime}$
are denoted as $G^{\prime}_{\rm M}$ and $G^{\prime\prime}_{\rm M}$,
respectively. In the TBS, they are given as

$\displaystyle G^{\prime}_{\rm M}=\sqrt{3}\left(k_{\rm n}+k_{\rm t}\right)/4,\
\ G^{\prime\prime}_{\rm M}=\sqrt{3}k_{\rm t}/(4\pi),$ (33)
$\displaystyle\mathcal{G}^{\prime}(x)=\left\\{\begin{array}[]{ll}1,&x\leq
x_{c},\\\ \left(1+\dfrac{k_{\rm t}}{k_{\rm
n}}\dfrac{T(x)-S(x)}{\pi}\right)/\left(1+\dfrac{k_{\rm t}}{k_{\rm
n}}\right),&x>x_{c},\end{array}\right.$ (36)
$\displaystyle\mathcal{G}^{\prime\prime}(x)=\left\\{\begin{array}[]{ll}0,&x\leq
x_{c},\\\ 1-\cos^{2}T(x),&x>x_{c}\end{array}\right.$ (39)

with $T(x)=\cos^{-1}(1-2x_{c}/x)$, $S(x)=\sin(2T(x))/2$, and
$x_{c}=4/(3\sqrt{3})$.

##  4 Comparison with the MBS

We demonstrate the relevance of the TBS analysis based on the simulation of a
two-dimensional MBS consisting of $N$ frictional grains. First, we consider a
system corresponding to the TBS, where all the particles are identical and
initially placed on the triangular lattice with a unit length $\ell$ (Fig.
1(a)). The details are shown in Appendix A. Next, we consider a bidisperse
system where the number of particles with diameter $d$ is equal to that of
particles with diameter $d/1.4$, and the particles are randomly placed with
packing fraction $\phi$ (Fig. 1(b)). The mass densities of the particles are
identical. The details of the disordered MBS are shown in Appendix E. In both
systems, the shear strain given by eqn. (4) is applied for $N_{\rm c}$ cycles
using the SLLOD equation under the Lees–Edwards boundary condition 34. In the
MBS, we replace the normal force as

$f_{ij}^{\rm(n)}\to-\left(k_{\rm n}u^{\rm(n)}_{ij}+\eta_{\rm
n}v^{\rm(n)}_{ij}\right)$ (40)

with the normal viscous constant $\eta_{\rm n}$ and the normal velocity
$v_{ij}^{\rm(n)}=(\frac{d}{dt}{\bm{r}}_{i}-\frac{d}{dt}{\bm{r}}_{j})\cdot\bm{n}_{ij}$
to include the viscous force depending on the relative velocity. The
tangential force is replaced by

$f_{ij}^{\rm(t)}\to{\rm min}\left(|\tilde{f}_{ij}^{\rm(t)}|,\mu
f_{ij}^{\rm(n,el)}\right){\rm sgn}(\tilde{f}_{ij}^{\rm(t)}),$ (41)

with

$\tilde{f}_{ij}^{\rm(t)}\to-\left(k_{\rm t}u_{ij}^{\rm(t)}+\eta_{\rm
t}v_{ij}^{\rm(t)}\right),$ (42)

where $f_{ij}^{\rm(n,el)}=-k_{\rm n}u^{\rm(n)}_{ij}$ denotes the elastic part
of the normal force with a tangential viscous constant $\eta_{\rm t}$. We
measure $G^{\prime}$, $G^{\prime\prime}$, and $P_{0}$ in the last cycle using
eqns. (11)-(13). For the ordered MBS, we use $N=64$, $k_{\rm t}/k_{\rm
n}=1.0$, and $\epsilon=0.001$, whereas $N=1000$, $k_{\rm t}/k_{\rm n}=0.2$,
and $\phi=0.87$ are used for the disordered MBS. In both systems, the other
parameters are identical: $N_{\rm c}=20$, $\eta_{\rm t}=\eta_{\rm
n}=\sqrt{mk_{\rm n}}$, and $\omega=0.0001\sqrt{m/k_{\rm n}}$ with a mass $m$
of larger particles.

Fig. 4:  Storage modulus $G^{\prime}$ against $\gamma_{0}$ with $k_{\rm
t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various values of $\mu$. The points
represent the results of the ordered MBS. The thin solid lines represent the
analytical result given by eqn. (25). The vertical dashed lines represent the
critical amplitude $\gamma_{c}(\mu)$ given by eqn. (16) for
$\mu=10^{-4},10^{-3},10^{-2},10^{-1}$, and $1$ from left to right.

As shown in Fig. 4, we plot $G^{\prime}$ for the ordered MBS against
$\gamma_{0}$ with $k_{\rm t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various
values of $\mu$ as points. Moreover, we plot the analytical results of the TBS
obtained using eqn. (25) as thin solid lines. Surprisingly, the results of the
TBS agree with those of the MBS for $\gamma_{0}<0.003$ without any fitting
parameters. As $\gamma_{0}$ increases beyond $\gamma_{c}(\mu)$ shown by the
vertical dashed lines, $G^{\prime}$ for $\mu>0$ decreases and converges to a
constant, which is equal to $G^{\prime}$ for $\mu=0$. For larger $\gamma_{0}$,
$G^{\prime}$ for the MBS decreases again, whereas the theoretical $G^{\prime}$
for the TBS is constant. This discrepancy results from the violation of
condition $\gamma_{0}\ll\epsilon$ for the analytical calculation. If we
numerically solve the TBS to obtain $G^{\prime}$ without the assumption
$\gamma_{0}\ll\epsilon$, $G\textquoteright$ decreases after a plateau again as
in the case of MBS, although its value in the TBS for $\gamma_{0}\to 0.1$
slightly deviates from that of the MBS, as shown in Appendix F.

As shown in Fig. 5, we plot $G^{\prime\prime}$ for the MBS on the triangular
lattice against $\gamma_{0}$ with $k_{\rm t}/k_{\rm n}=1.0$ and
$\epsilon=0.001$ for various values of $\mu$ as points. Moreover, we plot the
analytical results of the TBS obtained using eqn. (28) as thin solid lines.
The analytical result agrees perfectly with the MBS for $\gamma_{0}<0.003$
without any fitting parameters. As $\gamma_{0}$ increases beyond
$\gamma_{c}(\mu)$ shown by the vertical dashed lines, $G^{\prime\prime}$ for
$\mu>0$ increases from $0$ and decreases after reaching a peak. The peak
position of $G^{\prime\prime}$ against $\gamma_{0}$ increases with $\mu$.
Thus, our analytical results fail to capture the behavior of
$G^{\prime\prime}$ for $\mu=1$.

Fig. 5:  Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ with $k_{\rm
t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various values of $\mu$. The points
represent the results of the ordered MBS. The thin solid lines represent the
analytical results obtained using eqn. (28). The vertical dashed lines
represent the critical amplitude $\gamma_{c}(\mu)$ given by eqn. (16) for
$\mu=10^{-4},10^{-3},10^{-2},10^{-1}$, and $1$ from left to right.

Consider the disordered MBS shown in Fig. 1(b). Figure 6 shows the scaled
shear stress $\sigma/\gamma_{0}$ against the scaled strain $\gamma/\gamma_{0}$
in the disordered MBS with $\mu=0.0001$. The maximum value
$\tilde{\sigma}_{\rm max}$ decreases as $\gamma_{0}$ increases. The area $S$
of the curve is the largest for $\gamma_{0}=0.00003$. It is interesting that
stress-strain curves are not characterized by parallelograms in this case in
contrast to Fig. 3. This means that the disordered configuration of particles
creates an effective viscosity, and thus, the response to an applied strain
becomes visco-elastoplastic.

Fig. 6:  Scaled shear stress $\sigma/\gamma_{0}$ against $\gamma/\gamma_{0}$
in the disordered MBS for various values of $\gamma_{0}$ with $\mu=0.01$ and
$\phi=0.870$.

The behaviors of $G^{\prime}$ and $G^{\prime\prime}$ in the disordered MBS are
similar to those of the TBS as shown in Appendix E. Therefore, it is expected
that the scaling laws in eqns. (31) and (32) for a given $\epsilon$ in the TBS
can be used even in this system with corresponding $\phi$. This expectation is
verified by Fig. 7, in which we plot the scaled moduli
$G^{\prime}/G^{\prime}_{\rm M}$ and $G^{\prime\prime}/G^{\prime\prime}_{\rm
M}$ against the scaled strain $k_{\rm t}\gamma_{0}/(\mu
P_{0}(\gamma_{0},\mu))$ for various values of $\mu$ in the disordered MBS.
Moreover, we plot the analytical results for the TBS obtained using eqns. (36)
and (39) as solid lines, which qualitatively reproduce the MBS results for
small scaled strain, while the scaling is apparently violated for large scaled
strain. Here, we choose $k_{\rm t}/k_{\rm n}=1.5$ for the TBS to fit the
second plateau to that of the MBS. At present, we do not know the relationship
between $\phi$ and the fitting parameter.

Fig. 7:  (a) Scaled storage modulus $G^{\prime}/G^{\prime}_{\rm M}$ against
the scaled strain $k_{\rm t}\gamma_{0}/(\mu P_{0}(\gamma_{0},\mu))$ with
$\phi=0.87$ and $k_{\rm t}/k_{\rm n}=0.2$ for various values of $\mu$ in the
disordered MBS. The solid line represents the analytical result of the TBS
given by eqn. (36) with $k_{\rm t}/k_{\rm n}=1.5$. (b) Scaled loss modulus
$G^{\prime\prime}/G^{\prime\prime}_{\rm M}$ against the scaled strain $k_{\rm
t}\gamma_{0}/(\mu P_{0}(\gamma_{0},\mu))$ with $\phi=0.87$ and $k_{\rm
t}/k_{\rm n}=0.2$ for various values of $\mu$ in the disordered MBS. The solid
line represents the analytical result of the TBS given by eqn. (39) with
$k_{\rm t}/k_{\rm n}=1.5$.

##  5 Conclusions

We demonstrated the relevancy of a model of the TBS to describe the complex
modulus of jammed frictional granular materials under oscillatory shear. We
obtained the analytical expressions for the $\gamma_{0}$-dependence of
$G^{\prime}$ and $G^{\prime\prime}$ as shown in eqns. (16), (25), and (28),
which predict the $\mu$-dependence of the critical amplitude $\gamma_{c}$, the
decrease of $G^{\prime}$, and the peak of $G^{\prime\prime}$ above
$\gamma_{c}$ for crystalline solids. The analytical expressions lead to the
scaling laws given by eqns. (31) and (32). Although we have ignored the non-
affine motion for crystalline solids, these analytical results quantitatively
agree with those of the ordered MBS. Surprisingly, some characteristic
features of disordered solids for low strain (or high pressure) can be
captured. These results indicate that the analysis of the toy model gives a
basis for understanding the nonlinear rheology of frictional granular
materials under small strain.

Although the values of the plateaus in $G^{\prime}$ for disordered MBS
depended on $\phi-\phi_{\rm J}$ 6, 9, the corresponding values of the TBS are
independent of $\phi-\phi_{\rm J}$, as expressed in eqn. (25). In addition,
our analytical expressions cannot reproduce the second decrease of
$G^{\prime}$ and increase of $G^{\prime\prime}$ near $\gamma_{0}=10^{-2}$ in
the MBS. The discrepancy should result from the disorder because it leads to
the $\phi$-dependence of $G^{\prime}$ 7. To include the disorder effect, we
regarded $k_{\rm t}/k_{\rm n}$ as a fitting parameter. In previous studies on
models with small degrees of freedom, e.g., the coherent potential
approximation 26, 29, 30, the corresponding fitting parameters were self-
consistently determined. In future studies, we will discuss the self-
consistent determination of the parameter for the TBS.

Some researchers are interested in contributions from higher harmonics
characterizing the nonlinear response to oscillatory shear 21, but the
nonlinear viscoelastic moduli characterizing the higher harmonics are
negligibly small for jammed frictionless particles 25. However, the higher
harmonics in the frictional granular materials require further careful
investigation.

## Conflicts of interest

There are no conflicts to declare.

## Appendix A: Details of Ordered MBS

This section explains the details of the ordered MBS consisting of
monodispersed particles initially placed on a triangular lattice. We consider
a two-dimensional assembly of $N$ frictional particles in a periodic box with
sizes along the $x$ and $y$ directions $L_{x}$ and $L_{y}$, respectively.
Here, we initially place $N=2N_{x}N_{y}$ particles of diameter $d$ with
integers $n_{x}$ and $n_{y}$ at $\bm{r}_{i}$ as

$\bm{r}_{i}=\left(n_{x}\ell-L_{x}/2,\sqrt{3}n_{y}\ell-L_{y}/2\right)$ (43)

for $0\leq i<N_{x}N_{y}$ with integers $n_{x}$, $n_{y}$, and
$i=n_{x}+N_{x}n_{y}$. For $N_{x}N_{y}\leq i<2N_{x}N_{y}$, $\bm{r}_{i}$ is
defined as

$\bm{r}_{i}=\left((n_{x}+1/2)\ell-L_{x}/2,\sqrt{3}(n_{y}+1/2)\ell-
L_{y}/2\right)$ (44)

with $i=n_{x}+N_{x}n_{y}+N_{x}N_{y}$. The initial configuration is illustrated
in Fig. 8. We choose $L_{x}=N_{x}\ell$ and $L_{y}=\sqrt{3}N_{y}\ell$ with
$\ell=d(1-\epsilon)$.

Fig. 8:  Initial configuration of mono-dispersed particles on a triangular
lattice. The red rectangle, including interactions represented by the blue
lines, corresponds to the TBS.

The position $\bm{r}_{i}$ and peculiar momentum ${\bm{p}}_{i}$ of particle $i$
with mass $m_{i}$ and diameter $d_{i}$ are driven by the SLLOD equation under
the Lees-Edwards boundary condition as 34

$\displaystyle\frac{d}{dt}{\bm{r}}_{i}$ $\displaystyle=$
$\displaystyle\dot{\gamma}(t)y_{i}\bm{e}_{x}+\frac{\bm{p}_{i}}{m_{i}},$ (45)
$\displaystyle\frac{d}{dt}{\bm{p}}_{i}$ $\displaystyle=$
$\displaystyle-\dot{\gamma}(t)p_{i,y}\bm{e}_{x}+\bm{f}_{i},$ (46)

where $\dot{\gamma}(t)=\gamma_{0}\omega\cos\omega t$ and $\bm{e}_{x}=(1,0)$ is
the unit vector along the $x$ direction. The interaction force $\bm{f}_{i}$ is
defined as

$\bm{f}_{i}=\sum_{j\neq
i}\left(f_{ij}^{\rm(n)}\bm{n}_{ij}+f_{ij}^{\rm(t)}\bm{t}_{ij}\right)H(d_{ij}-r_{ij})$
(47)

with $d_{ij}=(d_{i}+d_{j})/2$, $\bm{n}_{ij}=\bm{r}_{ij}/r_{ij}$,
$\bm{t}_{ij}=(-n_{ij,y},n_{ij,x})$, and
$\bm{r}_{ij}=\bm{r}_{i}-\bm{r}_{j}=(x_{ij},y_{ij})$. The normal force is given
by

$f_{ij}^{\rm(n)}=-\left(k_{\rm
n}u^{\rm(n)}_{ij}+\eta_{n}v^{\rm(n)}_{ij}\right)$ (48)

with a normal viscous constant $\eta_{\rm n}$ and

$\displaystyle
v^{\rm(n)}_{ij}=\left(\bm{v}_{i}-\bm{v}_{j}\right)\cdot\bm{n}_{ij},$ (49)

where the velocity of particle $i$ is given by
$\bm{v}_{i}=\frac{d}{dt}\bm{r}_{i}$. The following model is adopted for the
tangential force:

$f_{ij}^{\rm(t)}={\rm min}\left(|\tilde{f}_{ij}^{\rm(t)}|,\mu
f_{ij}^{\rm(n,el)}\right){\rm sgn}(\tilde{f}_{ij}^{\rm(t)}),$ (50)

where $f_{ij}^{\rm(n,el)}=-k_{\rm n}u^{\rm(n)}_{ij}$ denotes the elastic part
of the normal force. Here, $\tilde{f}_{ij}^{\rm(t)}$ is given by

$\tilde{f}_{ij}^{\rm(t)}=-\left(k_{\rm t}u_{ij}^{\rm(t)}+\eta_{\rm
t}v_{ij}^{\rm(t)}\right)$ (51)

with a tangential viscous constant $\eta_{\rm t}$. The tangential velocity
$v_{ij}^{\rm(t)}$ is given by

$v_{ij}^{\rm(t)}=(\bm{v}_{i}-\bm{v}_{j})\cdot\bm{t}_{ij}.$ (52)

The tangential displacement $u_{ij}^{\rm(t)}$ satisfies
$\frac{d}{dt}u_{ij}^{\rm(t)}=v_{ij}^{\rm(t)}$ for
$|\tilde{f}_{ij}^{\rm(t)}|<\mu f_{ij}^{\rm(n,el)}$, whereas $u_{ij}^{\rm(t)}$
remains unchanged for $|\tilde{f}_{ij}^{\rm(t)}|\geq\mu f_{ij}^{\rm(n,el)}$.
The tangential displacement $u_{ij}^{\rm(t)}$ is set to zero if $i$ and $j$
are detached.

If all the particles are separated, the packing fraction $\phi$ for the
ordered MBS is defined as

$\displaystyle\phi=\frac{\sum_{i}\pi d_{i}^{2}}{4L_{x}L_{y}}.$ (53)

Even if contacts exist between the particles, we use eqn. (53) by assuming
that the contact length $d_{ij}-r_{ij}$ is sufficiently lower than $d_{ij}$.
Using eqn. (53), $\phi$ is defined as

$\displaystyle\phi=\frac{\pi}{2\sqrt{3}(1-\epsilon)^{2}}.$ (54)

The jamming point of this system is

$\displaystyle\phi_{\rm J}=\frac{\pi}{2\sqrt{3}}$ (55)

with $\epsilon=0$. The distance from the jamming point is proportional to
$\epsilon$.

$\displaystyle\phi-\phi_{\rm J}\simeq\frac{\pi}{\sqrt{3}}\epsilon$ (56)

for $\epsilon\ll 1$.

The shear stress $\sigma$ is defined by eqn. (8) in the main article with the
normal component

$\sigma^{\rm(n)}=-\frac{1}{L_{x}L_{y}}\sum_{i}\sum_{j>i}\frac{x_{ij}y_{ij}}{r_{ij}}f^{\rm(n)}_{ij}$
(57)

and tangential component

$\sigma^{\rm(t)}=-\frac{1}{2L_{x}L_{y}}\sum_{i}\sum_{j>i}\frac{x_{ij}^{2}-y_{ij}^{2}}{r_{ij}}f^{\rm(t)}_{ij}.$
(58)

The pressure is defined as

$P=\frac{1}{2L_{x}L_{y}}\sum_{i}\sum_{j>i}(x_{ij}f_{ij,x}+y_{ij}f_{ij,y}).$
(59)

We use $N_{x}=8$, $N_{y}=4$, $N_{\rm c}=20$, $k_{\rm t}=k_{\rm n}$, and
$\eta_{\rm n}=\eta_{\rm t}=k_{\rm n}\sqrt{m/k_{\rm n}}$, where $m$ denotes the
mass of a particle with diameter $d$. This model corresponds to a restitution
coefficient $e=0.043$. We adopt the leapfrog algorithm considering a time step
of $\Delta t=0.05t_{0}$. We choose $\omega=1.0\times 10^{-4}\sqrt{k_{\rm
n}/m}$ as the quasistatic shear deformation because $G^{\prime}$ and
$G^{\prime\prime}$ are almost independent of $\omega$ for $\omega\leq
1.0\times 10^{-3}\sqrt{k_{\rm n}/m}$.

As shown in Figs. 4 and 5, the behaviors of $G^{\prime}$ and
$G^{\prime\prime}$ of the TBS agree with that of the MBS. We explain the
theoretical background of the TBS. The initial configuration is shown in Fig.
8; it contains the unit cell represented by the red rectangle with length
$\ell$ and height $\sqrt{3}\ell/2$. It contains interactions between the three
particles represented by blue lines. Here, we assume that the particles move
affinely as

$\displaystyle\bm{r}_{i}(t)=\bm{r}_{i}(0)+\gamma(\theta(t))y_{i}(0)\bm{e}_{x}.$
(60)

In this case, the corresponding relative distances between the particles in
any unit cell are identical.

In particular, in a unit cell containing particles $i=i_{1},i_{2}$, and
$i_{3}$ with $i_{1}=N_{x}N_{y}$, $i_{2}=0$, and $i_{3}=1$, the positions of
the particles are given by

$\displaystyle\bm{r}_{i_{1}}(t)=\left(\gamma(\theta(t))\left(\frac{\sqrt{3}\ell-
L_{y}}{2}\right)+\frac{\ell-L_{x}}{2},\frac{\sqrt{3}\ell-L_{y}}{2}\right),$
(61)
$\displaystyle\bm{r}_{i_{2}}(t)=\left(-\gamma(\theta(t))\frac{L_{y}}{2}-\frac{L_{x}}{2},-\frac{L_{y}}{2}\right),$
(62)
$\displaystyle\bm{r}_{i_{3}}(t)=\left(-\gamma(\theta(t))\frac{L_{y}}{2}+\ell-\frac{L_{x}}{2},-\frac{L_{y}}{2}\right).$
(63)

The relative distances between these particles are identical to those of the
TBS, given by eqns. (1)-(3), which indicates that the TBS provides the
interaction forces among the three particles. This system includes
$2N_{x}N_{y}$ unit cells with identical interaction forces. Hence, the normal
and tangential components of $\sigma$ are given by

$\displaystyle\sigma^{\rm(n)}=-\frac{2N_{x}N_{y}}{L_{x}L_{y}}\sum_{i=i_{1},i_{2},i_{3}}\left\\{\sum_{\begin{subarray}{c}j=i_{1},i_{2},i_{2}\\\
(j>i)\end{subarray}}\frac{x_{ij}y_{ij}}{r_{ij}}f^{\rm(n)}_{ij}\right\\},$ (64)
$\displaystyle\sigma^{\rm(t)}=-\frac{N_{x}N_{y}}{L_{x}L_{y}}\sum_{i=i_{1},i_{2},i_{3}}\left\\{\sum_{\begin{subarray}{c}j=i_{1},i_{2},i_{2}\\\
(j>i)\end{subarray}}\frac{x_{ij}^{2}-y_{ij}^{2}}{r_{ij}}f^{\rm(t)}_{ij}\right\\}.$
(65)

The pressure is also given by:

$P=\frac{N_{x}N_{y}}{L_{x}L_{y}}\sum_{i=i_{1},i_{2},i_{3}}\left\\{\sum_{\begin{subarray}{c}j=i_{1},i_{2},i_{2}\\\
(j>i)\end{subarray}}(x_{ij}f_{ij,x}+y_{ij}f_{ij,y})\right\\}.$ (66)

Using the relation $L_{x}L_{y}/(2N_{x}N_{y})=\sqrt{3}\ell^{2}/2$ corresponding
to $A=\sqrt{3}\ell^{2}/2$, $\sigma^{\rm(n)}$, $\sigma^{\rm(t)}$, and $P$
coincide with eqns. (8)-(11). Hence, if the assumptions of the affine motion,
i.e., eqns. (61)–(63), are satisfied, $G^{\prime}$ and $G^{\prime\prime}$ in
the ordered MBS coincide with those in the TBS.

## Appendix B: Effect of particle rotation

In this section, we illustrate the effect of particle rotation, which was not
described in Appendix A. In the model with rotation, the tangential velocity
$v_{ij}^{\rm(t)}$ is given by

$v_{ij}^{\rm(t)}=(\bm{v}_{i}-\bm{v}_{j})\cdot\bm{t}_{ij}-(d_{i}\omega_{i}+d_{j}\omega_{j})/2$
(67)

instead of eqn. (52), where $\omega_{i}$ denotes the angular velocity of
particle $i$. The time evolution of $\omega_{i}$ is given by

$I_{i}\frac{d}{dt}\omega_{i}=T_{i}$ (68)

with the moment of inertia $I_{i}=m_{i}d_{i}^{2}/8$ and torque
$T_{i}=-\sum_{j}\frac{d_{i}}{2}\bm{F}_{ij}^{\rm(t)}\cdot\bm{t}_{ij}$.

As shown in Fig. 9, we plot $G^{\prime}$ in the ordered MBS with and without
rotation with $k_{\rm t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various
values of $\mu=0.01$. The values of other parameters are the same as those in
Appendix A. The effect of particle rotation is negligible, except for the
region near $\gamma_{c}$.

Fig. 9:  Storage modulus $G^{\prime}$ against $\gamma_{0}$ for the ordered MBS
with $\mu=0.01$ and $\epsilon=0.001$. The open and filled symbols represent
the results of the particles with and without rotation, respectively.

As shown in Fig. 10, we plot $G^{\prime\prime}$ in the ordered MBS with and
without rotation with $k_{\rm t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for
$\mu=10^{-2},10^{-3},10^{-4},10^{-5}$ and $0.00001$. The values of other
parameters are the same as those in Appendix A. There are slight deviations in
the peak position near $\gamma_{c}$ between particles with and without
rotation.

Fig. 10:  Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ for the ordered
MBS with $\mu=10^{-2},10^{-3},10^{-4},10^{-5}$ and $\epsilon=0.001$. The open
and filled symbols represent the results of the particles with and without
rotation, respectively.

In Fig. 11, we plot $G^{\prime\prime}$ in the ordered MBS with and without
rotation with $k_{\rm t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for $\mu=1.0$ and
$0.1$. There are slight deviations in the peak position near $\gamma_{c}$
between particles with and without rotation even for these higher $\mu$. In
addition, the second increase in $G^{\prime\prime}$ around $\gamma_{0}=0.1$
for particles without rotation disappears for those with rotation.

Fig. 11:  Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ for the ordered
MBS with $\mu=1.0,0.1$ and $\epsilon=0.001$. The open and filled symbols
represent the results of the particles with and without rotation,
respectively.

## Appendix C: Analytical calculation of shear stress and pressure

This section briefly explains the derivation of the normal and tangential
components of shear stress and pressure for a small value of $\gamma_{0}$.
From eqns. (1)–(3), the relative distance $\bm{r}_{ij}(t)$ is given by

$\displaystyle\bm{r}_{12}(\theta(t))=\left(\frac{\sqrt{3}\gamma(\theta(t))+1}{2}\ell,\frac{\sqrt{3}\ell}{2}\right),$
(69)
$\displaystyle\bm{r}_{13}(\theta(t))=\left(\frac{\sqrt{3}\gamma(\theta(t))-1}{2}\ell,\frac{\sqrt{3}\ell}{2}\right),$
(70) $\displaystyle\bm{r}_{23}(\theta(t))=\left(-\ell,0\right).$ (71)

Substituting these equations into $u^{\rm(n)}_{ij}=r_{ij}-d$, the normal
displacements are given by

$\displaystyle u^{\rm(n)}_{12}(t)=-\epsilon
d+\frac{\sqrt{3}}{4}\ell\gamma(\theta(t))+O(\gamma_{0}^{2}),$ (72)
$\displaystyle u^{\rm(n)}_{13}(t)=-\epsilon
d-\frac{\sqrt{3}}{4}\ell\gamma(\theta(t))+O(\gamma_{0}^{2}),$ (73)
$\displaystyle u^{\rm(n)}_{23}(t)=-\epsilon d.$ (74)

Substituting these equations into eqn. (6), we obtain the normal force as

$\displaystyle f_{12}^{\rm(n)}=k_{\rm n}\left(\epsilon
d-\frac{\sqrt{3}}{4}\gamma(\theta)\ell\right),$ (75) $\displaystyle
f_{13}^{\rm(n)}=k_{\rm n}\left(\epsilon
d-\frac{\sqrt{3}}{4}\gamma(\theta)\ell\right),$ (76) $\displaystyle
f_{23}^{\rm(n)}=k_{\rm n}\epsilon d$ (77)

up to $O(\gamma_{0})$.

By differentiating eqns. (69)–(71) with time $t$, we obtain the relative
velocity as

$\displaystyle\bm{v}_{12}(t)=\left(\frac{\sqrt{3}\dot{\gamma}(\theta(t))\ell}{2},0\right),$
(78)
$\displaystyle\bm{v}_{13}(t)=\left(\frac{\sqrt{3}\dot{\gamma}(\theta(t))\ell}{2},0\right),$
(79) $\displaystyle\bm{v}_{23}(t)=\left(0,0\right)$ (80)

with the strain rate $\dot{\gamma}(\theta(t))=\frac{d}{dt}\gamma(\theta(t))$.
The tangential unit vector is given by

$\displaystyle\bm{t}_{12}(t)=\left(-\frac{\sqrt{3}\ell}{2},\frac{\sqrt{3}\gamma(\theta(t))+1}{2}\ell\right)/\left|\bm{r}_{12}\right|,$
(81)
$\displaystyle\bm{t}_{13}(t)=\left(-\frac{\sqrt{3}\ell}{2},\frac{\sqrt{3}\gamma(\theta(t))-1}{2}\ell\right)/\left|\bm{r}_{13}\right|,$
(82) $\displaystyle\bm{t}_{23}(t)=\left(0,-1\right).$ (83)

By considering the inner product of $\bm{v}_{ij}$ and $\bm{t}_{ij}$, the
tangential velocity is given by

$\displaystyle
v^{\rm(t)}_{12}(t)=-\frac{3}{4}\ell\dot{\gamma}(\theta(t))+O(\gamma_{0}^{2}),$
(84) $\displaystyle
v^{\rm(t)}_{13}(t)=-\frac{3}{4}\ell\dot{\gamma}(\theta(t))+O(\gamma_{0}^{2}),$
(85) $\displaystyle v^{\rm(t)}_{23}(t)=0.$ (86)

If the transition from the stick state to the slip state does not occur under
oscillatory shear, the tangential displacement is obtained by integrating
$v^{\rm(t)}_{ij}(t)$ as

$\displaystyle
u^{\rm(t)}_{12}(t)=u^{\rm(t)}_{13}(t)=-\frac{3}{4}\ell\gamma(\theta(t))+O(\gamma_{0}^{2}),$
(87) $\displaystyle u^{\rm(t)}_{23}(t)=0.$ (88)

Substituting these equations into $f_{ij}^{\rm(t)}=-k_{\rm t}u_{ij}^{\rm(t)}$
yields

$\displaystyle f_{12}^{\rm(t)}=f_{13}^{\rm(t)}=3k_{\rm
t}\gamma(\theta(t))\ell/4,$ (89) $\displaystyle f_{23}^{\rm(t)}=0$ (90)

up to $O(\gamma_{0})$. The condition that the transition does not occur is
satisfied when $f_{12}^{\rm(t)}<\mu f_{12}^{\rm(n)}$ for $\gamma=\gamma_{0}$.
Using eqns. (75) and (89) with the assumption $\gamma_{0}\ll\epsilon$, the
condition is replaced by $\gamma_{0}<\gamma_{c}$ with $\gamma_{c}$ given by
eqn. (16).

For $\gamma_{0}>\gamma_{c}$, there exist regions where $u_{ij}^{\rm(t)}$ is
unchanged in the slip state as

$\displaystyle u_{12}^{\rm(t)}=\left\\{\begin{array}[]{ll}-\dfrac{\mu k_{\rm
n}\epsilon d}{k_{\rm t}},&0\leq\theta(\theta)<\dfrac{\pi}{2}\\\ -\dfrac{\mu
k_{\rm n}\epsilon d}{k_{\rm
t}}-\dfrac{3d(\gamma(\theta)-\gamma_{0})}{4},&\dfrac{\pi}{2}\leq\theta<\dfrac{\pi}{2}+\Theta\\\
\dfrac{\mu k_{\rm n}\epsilon d}{k_{\rm
t}},&\dfrac{\pi}{2}+\Theta\leq\theta<\dfrac{3\pi}{2}\\\ \dfrac{\mu k_{\rm
n}\epsilon d}{k_{\rm
t}}-\dfrac{3d(\gamma(\theta)+\gamma_{0})}{4},&\dfrac{3\pi}{2}\leq\theta<\dfrac{3\pi}{2}+\Theta\\\
-\dfrac{\mu k_{\rm n}\epsilon d}{k_{\rm
t}},&\dfrac{3\pi}{2}+\Theta\leq\theta<2\pi,\end{array}\right.$ (96)
$\displaystyle u_{13}^{\rm(t)}=u_{12}^{\rm(t)},$ (97) $\displaystyle
u_{23}^{\rm(t)}=0,$ (98)

where $\Theta$ satisfies

$\displaystyle-\frac{\mu k_{\rm n}\epsilon d}{k_{\rm
t}}-3d\frac{\gamma\left(\frac{\pi}{2}+\Theta\right)-\gamma_{0}}{4}=\frac{\mu
k_{\rm n}\epsilon d}{k_{\rm t}}.$ (99)

This equation provides
$\Theta=\cos^{-1}\left(1-2\gamma_{c}/\gamma_{0}\right)$. Substituting these
equations into $f_{ij}^{\rm(t)}=-k_{\rm t}u_{ij}^{\rm(t)}$ yields

$\displaystyle f_{12}^{\rm(t)}=\left\\{\begin{array}[]{ll}-\mu k_{\rm
n}\epsilon d,&0\leq\theta(\theta)<\dfrac{\pi}{2}\\\ -\mu k_{\rm n}\epsilon
d-\dfrac{3k_{\rm
t}d(\gamma(\theta)-\gamma_{0})}{4},&\dfrac{\pi}{2}\leq\theta<\dfrac{\pi}{2}+\Theta\\\
\mu k_{\rm n}\epsilon d,&\dfrac{\pi}{2}+\Theta\leq\theta<\dfrac{3\pi}{2}\\\
\mu k_{\rm n}\epsilon d-\dfrac{3k_{\rm
t}d(\gamma(\theta)+\gamma_{0})}{4},&\dfrac{3\pi}{2}\leq\theta<\dfrac{3\pi}{2}+\Theta\\\
-\mu k_{\rm n}\epsilon
d,&\dfrac{3\pi}{2}+\Theta\leq\theta<2\pi,\end{array}\right.$ (105)
$\displaystyle f_{13}^{\rm(t)}=f_{12}^{\rm(t)},$ (106) $\displaystyle
f_{23}^{\rm(t)}=0.$ (107)

The normal component of $\sigma$ in eqn. (9) is given by

$\sigma^{\rm(n)}=\sigma^{\rm(n)}_{12}+\sigma^{\rm(n)}_{13}$ (108)

with

$\displaystyle\sigma^{\rm(n)}_{12}=-\frac{1}{A}\frac{x_{12}y_{12}}{r_{12}}f_{12}^{\rm(n)}$
(109)
$\displaystyle\sigma^{\rm(n)}_{13}=-\frac{1}{A}\frac{x_{13}y_{13}}{r_{13}}f_{12}^{\rm(n)}.$
(110)

Substituting eqns. (69) and (70) with eqns. (75) and (76) into eqns. (109) and
(110) and using eqn. (108), we obtain $\sigma^{\rm(n)}$ as eqn. (14).

The tangential component of $\sigma$ in eqn. (10) is given by

$\sigma^{\rm(t)}=\sigma^{\rm(t)}_{12}+\sigma^{\rm(t)}_{13}$ (111)

with

$\displaystyle\sigma^{\rm(t)}_{(12)}=-\frac{1}{2A}\frac{x_{12}^{2}-y_{12}^{2}}{r_{12}}f_{12}^{\rm(t)}$
(112)
$\displaystyle\sigma^{\rm(t)}_{(13)}=-\frac{1}{2A}\frac{x_{13}^{2}-y_{13}^{2}}{r_{13}}f_{12}^{\rm(t)}.$
(113)

Substituting eqns. (69) and (70) with eqn. (89) into eqns. (111), (112), and
(113), we obtain $\sigma^{\rm(t)}$ as eqn. (15) for $\gamma_{0}<\gamma_{c}$.
Using eqns. (105) and (106) instead of eqn. (89), we obtain $\sigma^{\rm(t)}$
as eqn. (22) for $\gamma_{0}\geq\gamma_{c}$.

The pressure, i.e., $P$, in eqn. (11) is defined as

$P=P_{12}+P_{13}+P_{23}$ (114)

with

$\displaystyle P_{ij}=\frac{1}{2A}r_{ij}f_{ij}^{\rm(n)}.$ (115)

Substituting eqns. (69)-(71) with eqns. (75)-(77) into eqns. (114) and (115)
with $\gamma=0$, we obtain $P_{0}(\gamma_{0},\mu)$ as eqn. (30).

## Appendix D: Relation between shear modulus and stress–strain curve

In this section, we relate the shape of the stress–strain curve to the complex
shear modulus. The shear stress $\sigma(\theta)$ is expanded using the Fourier
series as

$\displaystyle\sigma(\theta)=\gamma_{0}\sum_{n=1}^{\infty}G^{\prime}_{n}\sin(n\theta)+\gamma_{0}\sum_{n=1}^{\infty}G^{\prime\prime}_{n}\cos(n\theta),$
(116)

where $G^{\prime}_{n}$ and $G^{\prime\prime}_{n}$ with $n>1$ denote the higher
harmonics, $G^{\prime}=G^{\prime}_{1}$, and
$G^{\prime\prime}=G^{\prime\prime}_{1}$. By neglecting $G^{\prime}_{n}$ and
$G^{\prime\prime}_{n}$ for $n>1$,

$\displaystyle
G^{\prime}\simeq\frac{\sigma\left(\theta=\pi/2\right)}{\gamma_{0}}=\left.\frac{\sigma}{\gamma_{0}}\right|_{\gamma/\gamma_{0}=1}=\tilde{\sigma}_{\rm
max},$ (117)

which is the maximum value of the scaled stress–strain curve illustrated in
Fig. 3(b). This expression and the scaled stress–strain curve in Fig. 3(b)
explain the decrease of $G^{\prime}$ defined by eqn. (25).

The area $S$ of the curve for $\sigma(\theta)/\gamma_{0}$ against
$\gamma(\theta)/\gamma_{0}$ is given by

$\displaystyle
S=\int_{0}^{2\pi}d\theta\frac{1}{\gamma_{0}}\frac{d\gamma(\theta)}{d\theta}\frac{\sigma(\theta)}{\gamma_{0}}.$
(118)

Substituting eqn. (4) into eqn. (118) with eqn. (12), we obtain

$\displaystyle S=\int_{0}^{2\pi}d\theta\sigma(\theta)\cos\theta/\gamma_{0}=\pi
G^{\prime\prime},$ (119)

which results in $G^{\prime\prime}=S/\pi$. As $\gamma_{0}$ increases, the area
$S$ of the scaled stress–strain curve in Fig. 3(b) increases first and
decreases later, which explains the $\gamma_{0}$-dependence of
$G^{\prime\prime}$ provided by eqn. (28).

## Appendix E: Details of Disordered MBS

In this section, we present the details of the disordered MBS. This model is
an extension of the monodisperse model used in Appendix A, including the
dispersion of the particles and disordered initial configuration.

The system is bidisperse and includes an equal number of particles with
diameters $d$ and $d/1.4$. To simulate the disordered MBS, we randomly place
the particles in a rectangular box with an initial packing fraction of
$\phi_{\rm I}=0.75$. The system is slowly compressed until the packing
fraction reaches $\phi$ 24. In each compression step, the packing fraction is
increased by $\Delta\phi=1.0\times 10^{-4}$ with an affine transformation.
Thereafter, the particles are relaxed to a mechanical equilibrium state with
the kinetic temperature $T_{\rm K}=\sum_{i}p_{i}^{2}/(mN)<T_{\rm th}$. Here,
we choose $T_{\rm th}=1.0\times 10^{-8}k_{\rm n}d^{2}$. After compression, the
oscillatory shear strain given by eqn. (4) is applied for $N_{\rm c}$ cycles.
In the last cycle, we measure $G^{\prime}$ and $G^{\prime\prime}$ using eqns.
(12) and (13) with eqns. (8)–(10). The pressure, $P_{0}(\gamma_{0},\mu)$ is
obtained using eqn. (11) after the last cycle. We use $\phi=0.87$, $N=1000$,
$N_{\rm c}=20$, $L_{y}/L_{x}=1$, $k_{\rm t}=0.2k_{\rm n}$, and $\eta_{\rm
n}=\eta_{\rm t}=k_{\rm n}\sqrt{m/k_{\rm n}}$.

Figure 12 shows the storage modulus $G^{\prime}$ against $\gamma_{0}$ in the
disordered MBS for various values of $\mu$. The storage modulus $G^{\prime}$
is almost independent of $\gamma_{0}$ for a small $\gamma_{0}$ and decreases
as $\gamma_{0}$ increases. The endpoint of the first plateau increases with
$\mu$ except for $\mu=0$. A second plateau of $G^{\prime}$ exists for
$\mu=10^{-4}$ and $10^{-5}$. The behavior of $G^{\prime}$ for relatively small
$\gamma_{0}$ is similar to that of crystalline solids as depicted in Fig. 4.
On the other hand, the decrease of $G^{\prime}$ for larger $\gamma_{0}$ cannot
be captured by the analytical results of the TBS. Note that $G^{\prime}$ for
$\mu=0.1$ in the limit $\gamma_{0}\to 0$ is different from that for $\mu\leq
0.01$, which results from the $\mu$ dependence of the jamming point $\phi_{\rm
J}$ 9.

Fig. 12:  Storage modulus $G^{\prime}$ in the disordered MBS against
$\gamma_{0}$ with $\phi=0.870$ for various values of $\mu$.

Figure 13 shows the loss modulus $G^{\prime\prime}$ in the disordered MBS
against $\gamma_{0}$ for various values of $\mu$. For sufficiently small
$\gamma_{0}$, $G^{\prime\prime}$ is zero, while $G^{\prime\prime}$ becomes
non-zero as $\gamma_{0}$ increases. The loss modulus $G^{\prime\prime}$ starts
to increase for smaller $\gamma_{0}$ as $\mu$ decreases. Similar to the case
of $G^{\prime}$, TBS captures only the behavior of relatively small
$\gamma_{0}$ (see Figs. 5 and 13).

Fig. 13:  Loss modulus $G^{\prime\prime}$ in the disordered MBS against
$\gamma_{0}$ with $\phi=0.870$ for various values of $\mu$.

## Appendix F: Numerical shear modulus for TBS

In this section, we show the behaviors of $G^{\prime}$ and $G^{\prime\prime}$
in the TBS without the assumption used to obtain the analytical solution.
Here, we numerically obtain $G^{\prime}$ and $G^{\prime\prime}$ under
quasistatic oscillation using eqns. (12) and (13) based on the left Riemann
sum, where the integration of $\Psi(\theta)$, i.e.,

$\displaystyle\int_{0}^{2\pi}d\theta\ \Psi(\theta),$ (120)

is approximated as

$\displaystyle\int_{0}^{2\pi}d\theta\
\Psi(\theta)\simeq\sum_{n=1}^{M}\Psi(\theta_{n})\Delta\theta$ (121)

with $\Delta\theta=2\pi/M$ and $\theta_{n}=(n-1)\Delta\theta$. We use
$\epsilon=0.001$ and $\Delta\theta=5.0\times 10^{-5}$ in our simulation.

As shown in Fig. 14, we plot the storage modulus $G^{\prime}$ numerically
obtained from the TBS against $\gamma_{0}$ with $k_{\rm t}/k_{\rm n}=1.0$ for
various values of $\mu$ as points. Moreover, we plot the analytical results
derived from eqn. (25) as thin solid lines. The numerical results agree with
the analytical results for $\gamma_{0}<0.003$ and reproduce the second plateau
of the MBS shown in Fig. 4.

Fig. 14:  Storage modulus $G^{\prime}$ against $\gamma_{0}$ with $k_{\rm
t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various values of $\mu$. The points
represent the numerical results of the TBS, while the thin solid lines
represent the analytical result given by eqn. (25). The vertical dashed lines
represent the critical amplitude $\gamma_{c}(\mu)$ given by eqn. (16) for
$\mu=10^{-4},10^{-3},10^{-2},10^{-1},10^{0}$ from left to right.

Figure 15 shows the loss modulus $G^{\prime\prime}$ numerically obtained from
the TBS against $\gamma_{0}$ with $k_{\rm t}/k_{\rm n}=1.0$ for various values
of $\mu$ as points. We also plot the analytical results given by eqn. (28) as
thin solid lines. The numerical results agree with the analytical results for
$\gamma_{0}<0.003$.

Fig. 15:  Loss modulus $G^{\prime\prime}$ against $\gamma_{0}$ with $k_{\rm
t}/k_{\rm n}=1.0$ and $\epsilon=0.001$ for various values of $\mu$. The points
represent the numerical results of the TBS, while the thin solid lines
represent the analytical results obtained from eqn. (28). The vertical dashed
lines represent the critical amplitude $\gamma_{c}(\mu)$ obtained from eqn.
(16) for $\mu=10^{-4},10^{-3},10^{-2},10^{-1},10^{0}$ from left to right.

## Acknowledgements

The authors thank K. Saitoh, D. Ishima, and S. Takada for fruitful
discussions. This study was supported by JSPS KAKENHI under Grant Nos.
JP19K03670 and JP21H01006.

## Notes and references

  * van Hecke 2010 M. van Hecke, _J. Phys.: Condens. Matter_ , 2010, 22, 033101.
  * Behringer and Chakraborty 2019 R. P. Behringer and B. Chakraborty, _Rep. Prog. Phys._ , 2019, 82, 012601\.
  * Maloney and LemaÃ®tre 2004 C. Maloney and A. LemaÃ®tre, _Phys. Rev. E_ , 2004, 93, 195501\.
  * Maloney and LemaÃ®tre 2006 C. E. Maloney and A. LemaÃ®tre, _Phys. Rev. E_ , 2006, 74, 016118\.
  * 5 D. Ishima, K. Saitoh, M. Otsuki and H. Hayakawa, 2022, arXiv:2207.06632.
  * O’Hern _et al._ 2002 C. S. O’Hern, S. A. Langer, A. J. Liu and S. R. Nagel, _Phys. Rev. Lett._ , 2002, 88, 075507.
  * Wyart 2005 M. Wyart, _Ann. Phys._ , 2005, 30, 1.
  * Tighe 2011 B. P. Tighe, _Phys. Rev. Lett._ , 2011, 107, 158303.
  * Otsuki and Hayakawa 2017 M. Otsuki and H. Hayakawa, _Phys. Rev. E_ , 2017, 95, 062902.
  * Coulais _et al._ 2014 C. Coulais, A. Seguin and O. Dauchot, _Phys. Rev. Lett._ , 2014, 113, 198001.
  * Otsuki and Hayakawa 2014 M. Otsuki and H. Hayakawa, _Phys. Rev. E_ , 2014, 90, 042202.
  * Nagamanasa _et al._ 2014 K. H. Nagamanasa, S. Gokhale, A. K. Sood and R. Ganapathy, _Phys. Rev. E_ , 2014, 89, 062308.
  * Knowlton _et al._ 2014 E. D. Knowlton, D. J. Pine and L. Cipelletti, _Soft Matter_ , 2014, 10, 6931–6940.
  * Kawasaki and Berthier 2016 T. Kawasaki and L. Berthier, _Phys. Rev. E_ , 2016, 94, 022615\.
  * Leishangthem _et al._ 2017 P. Leishangthem, A. D. S. Parmar and S. Sastry, _Nat. Commun._ , 2017, 8, 14653.
  * Clark _et al._ 2018 A. H. Clark, J. D. Thompson, M. D. Shattuck, N. T. Ouellette and C. S. O’Hern, _Phys. Rev. E_ , 2018, 97, 062901.
  * Boschan _et al._ 2019 J. Boschan, S. Luding and B. P. Tighe, _Granular Matter_ , 2019, 21, 58.
  * Boschan _et al._ 2016 J. Boschan, D. VÃ¥gberg, E. Somfai and B. P. Tighe, _Soft Matter_ , 2016, 12, 5450–5460.
  * Nakayama _et al._ 2016 D. Nakayama, H. Yoshino and F. Zamponi, _J. Stat. Mech._ , 2016, 2016, 104001.
  * 20 T. Kawasaki and K. Miyazaki, 2020, arXiv:2003.10716.
  * Hyun _et al._ 2011 K. Hyun, M. Wilhelm, C. O. Klein, K. S. Cho, J. G. Nam, K. H. Ahn, S. J. Lee, R. H. Ewoldt and G. H. McKinley, _Prog. Polym. Sci._ , 2011, 36, 1697\.
  * Dagois-Bohy _et al._ 2017 S. Dagois-Bohy, E. Somfai, B. P. Tighe and M. van Hecke, _Soft Matter_ , 2017, 13, 9036–9045.
  * Ishima and Hayakawa 2020 D. Ishima and H. Hayakawa, _Phys. Rev. E_ , 2020, 101, 042902\.
  * Otsuki and Hayakawa 2021 M. Otsuki and H. Hayakawa, _Eur. Phys. J. E_ , 2021, 44, 70.
  * Otsuki and Hayakawa 2022 M. Otsuki and H. Hayakawa, _Phys. Rev. Lett._ , 2022, 128, 208002\.
  * Goldenfeld 1992 N. Goldenfeld, _Lectures on Phase Transitions and the Renormalization Group_ , CRC Press, 1992.
  * Lennard-Jones and Devonshire 1937 J. E. Lennard-Jones and A. F. Devonshire, _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , 1937, 163, 53–70.
  * Lennard-Jones and Devonshire 1938 J. E. Lennard-Jones and A. F. Devonshire, _Proc. R. Soc. A: Math. Phys. Eng. Sci._ , 1938, 165, 1–11.
  * Yonezawa and Morigaki 1973 F. Yonezawa and K. Morigaki, _Prog. Theor. Phys. Suppl._ , 1973, 53, 1–76.
  * Feng _et al._ 1985 S. Feng, M. F. Thorpe and E. Garboczi, _Phys. Rev. B_ , 1985, 31, 276–280.
  * Awazu 2001 A. Awazu, _Phys. Rev. E_ , 2001, 63, 032102.
  * Cundall and Strack 1979 P. A. Cundall and O. D. L. Strack, _GÃ©otechnique_ , 1979, 29, 47–65.
  * Doi and Edwards 1986 M. Doi and S. F. Edwards, _The theory of polymer dynamics_ , Oxford University Press, 1986.
  * Evans and Morriss 2008 D. J. Evans and G. Morriss, _Statistical Mechanics of Nonequilibrium Liquids_ , Cambridge University Press, 2008.