\section{Introduction}
Solving partial differential equations (PDEs) through deep learning approach has attracted extensive attention recently. Thanks to the universal approximation theorem of neural networks, it is natural to approximate solutions of PDEs using neural network. Many popular neural network based methods have been proposed recently, such as Deep Ritz Method \cite{yu2018deep}, Deep Galerkin Method \cite{sirignano2018dgm}, Physics Informed Neural Networks (PINNs) \cite{raissi2019physics} and the Weak Adversarial Networks \cite{zang2020weak}. In spite of the great success of these methods in solving various PDEs, the neural networks need to be re-trained if one seeks solutions corresponding to different initial conditions (ICs), boundary conditions (BCs) or parameters for the same PDEs. Instead, the recently proposed parametric operator learning methods, such as DeepONet \cite{lu2019deeponet} and FNO \cite{li2020fourier} enable learning of PDEs corresponding to varying BCs or ICs without re-training the networks.
However, there is one important caveat in the aforementioned operator neural networks. Namely they are essentially supervised learning and often require solving large number of PDEs to form the training data, which can be extremely expansive, especially when PDEs of interest lie in high dimensional spaces. To overcome this issue, Wang et al. \cite{wang2021learning,wang2021long} proposed the physics-informed DeepONet, which uses only the physical information (for instance the governing law of the PDEs) to construct loss function and thus making DeepONet self-supervised. Nevertheless, in practice the physics-informed DeepONets are more difficult to train compared to its vanilla version since the exact differential operators act on the networks and make the convergence behavior highly depends on the underlying physics problem.

Recently DeepONet has also been applied to learning the propagators of evolution equations; see e.g. \cite{liu2022deeppropnet,wang2021long}. The basic idea is to employ DeepONets to learn the solution operator of a PDE within a short time interval subject to a collection of (random) initial conditions. The solution of the PDE at later times can be computed as recursive actions of the trained network operator on solutions obtained at the prior steps. However, the approximation accuracy of solutions can deteriorate in the long-run for at least two reasons. First, due to the approximation error, the trained DeepONet, as a surrogate propagator, may be expansive even if the exact propagator is non-expansive, which leads to the accumulation of approximation error in time and hence makes it difficult to predict the solution in the long-run. Second, during the time-evolution of PDEs, the functions that a propagator inputs and outputs can vary in time, even though the form of the propagator within a fixed time-slot may remain unchanged (e.g. when the dynamics is autonomous).  Taking diffusion equation as an example, one observes that the functions in the range space of the propagator or the semigroup are much smoother  than those in the domain, and for this reason, the solutions in later times become increasingly more regular than those in earlier times. Furthermore, some evolution equations may develop various complexities in a long time-horizon, such as turbulence and scale separations. For those equations,  iterating a DeepONet surrogate that is usually only trained in a single (short) time-frame using a finite collection of initial functions may fail to capture the correct regularity or complexity of the solutions in the long time.

Transfer learning \cite{bozinovski1976influence,do2005transfer} is an important class of  machine learning techniques that use one neural network trained for one task for a new neural network trained for a different related task. The idea is that the knowledge or important features of one problem gained by training the former neural nets  can be transferred to other problems. Transfer learning has been widely used in image recognition \cite{yin2019feature,jin2020deep}, natural language processing \cite{ruder2019transfer} and recently in PINNs \cite{goswami2019transfer,obiols2021surfnet,song2022transfer,desai2021one}. To the best of our knowledge, the present work is the first work to employ transfer learning for learning solution operators of evolutionary PDEs.

\subsection{Our contributions}
We propose a novel physics-informed DeepONet approach based on the transfer learning for predicting time-dependent PDEs. Different from the existing usage of DeepONets in learning the propagators of PDEs where the learned propagators are treated constant in time, we use transfer learning to sequentially update the learned propagators as time evolves. The resulting time-changing DeepONets offer several advantages compared to the vanilla counterparts: (1) the evolving DeepONets can better adapt to the varying complexities associated to the evolution equations; (2) the DeepONets are updated in a computationally efficient way that the hidden layers are frozen once trained and only the parameters in the last layer are re-trained.

We hereby highlight the major contributions of the proposed method:\begin{itemize}
	\item Time marching with the transfer-learning tuned DeepONet gives more accurate and robust long-time prediction of solutions of PDEs while still maintaining low computational cost.

	\item The proposed method is applied to various types of evolutionary PDEs, including the reaction diffusion equations, Allen-Cahn and Cahn-Hilliard equations, the Navie-Stokes equation and multiscale linear radiative transfer equations.

	\item Through extensive numerical results, we show that  our method can significantly reduce the training sample size needed by DeepONet to achieve the same (or even higher) accuracy.

\end{itemize}

\subsection{Related works}\label{sec:relatedwork}
Transfer-learning has been previously combined with physics informed neural networks for solving PDEs problems arising from diverse fields, including the phase-field modeling of fracture \cite{goswami2019transfer}, super-resolution of turbulent flows \cite{obiols2021surfnet}, training of CNNs on multi-fidelity data (e.g. multi-resolution images of PDE solutions on fine and coarse meshes) \cite{song2022transfer}, etc. In \cite{chakraborty2022domain}, transfer-learning was also applied
as a domain adaption method for learning solutions of PDEs defined on complex geometries. The recent paper \cite{desai2021one} proposed a one-shot transfer learning strategy that freezes the hidden layers of a pre-trained PINN and reduces the training neural networks for solving new differential equations to optimizing only the last (linear) layer. This approach  eliminates the need of re-training the whole network parameters while still produces  high-quality solutions by tuning a small fraction of parameters in the last layer.  The present paper marry this transfer learning idea with  DeepONet for learning the propagators of evolution equations in order to predict the long time evolution.

While we are finalizing the current paper, we are aware of a recent preprint  \cite{goswami2022deep} where transfer learning was exploited together with DeepONet for learning  PDEs under conditional shift. The purpose there is to train a source PDE model with sufficient labeled data from one source domain and transfer the learned parameter to a target domain with limited labeled data. The technology developed there is mainly applied for transferring the knowledge of a solution operator trained on a system of PDEs from one domain to another. Different from \cite{goswami2022deep}, we leverage transfer learning to successively tuning the surrogate models of  propagators learned via physics-informed DeepONet so that the tuned operator networks can adaptively track the evolving propagators that carry evolving inputs and outputs.
The proposed approach is proven to be more accurate and robust for learning the long-time evolution of PDEs.

\section{Numerical method}\label{sec:numerical_method}

\paragraph{Problem set-up}
Consider the initial boundary value problem for a general evolution equation:
\begin{equation} \label{111}
	\left\{\begin{array}{l}
		\partial_t f(t, \boldsymbol{x}) = \mathcal{L}(f(t, \boldsymbol{x}))  , \\
		f(t, \boldsymbol{x}) = \phi(\boldsymbol{x}), ~ \boldsymbol{x} \in \partial \Omega_x, \\
		f(0, \boldsymbol{x}) = f_0(\boldsymbol{x})\,, ~ \boldsymbol{x} \in  \Omega_x\,.
	\end{array}\right.
\end{equation}
Throughout the paper, we assume that the equations are dissipative in the sense that $\int_{\Omega_x} f \mathcal{L} f d\boldsymbol{x} \leq 0$.
Given a time step size $\Delta t$, we consider the semi-discrete approximation $f^n(\boldsymbol{x})$ of the solution $f(n \Delta t,\boldsymbol{x})$ to \eqref{111} defined by the backward Euler discretization:
\begin{equation}\label{eqn:mP}
	f^{n+1}(\boldsymbol{x}) = (I - \Delta t \mathcal{L})^{-1}f^n(\boldsymbol{x})
	:=\mathcal{P}^{\Delta t} f^n(\boldsymbol{x})\,.
\end{equation}
Our goal is to approximate the propagator
\[
\mathcal{P}^{\Delta t}: f^n(\boldsymbol{x}) \mapsto f^{n+1}(\boldsymbol{x})
\]
by an operator neural network $\mathcal{P}_{NN}$ so that only one forward pass of the neural network achieves time-marching solutions from one step to the next, and that the evolution dynamics can be captured in a long time-horizon.

It is important to point out that the backward Euler scheme is not the only choice for time-marching. One can extend it to high order time discretization schemes such as Runge-Kutta methods, as long as $\Delta t$ is chosen such that $\mathcal{P}^{\Delta t}$ is a non-expanding operator. We will make this point more clear in Section 3 and Appendix~\ref{2nd}. To ease the notation, the superscript $n$, $n+1$ and $\Delta t$ will be omitted in the following context if it does not cause any confusion.

\subsection{Physics-informed DeepONet}

Let $\Omega_x$ be a compact set in $\mathbb R^d$ and let $\mathcal{X}$ be a compact subspace of the space $C(\Omega_x)$ of continuous function defined on $\Omega_x$. Then according to the universal approximation theorem \cite{CC95}, an operator $\mathcal{P}: \mathcal X \rightarrow \mathcal X$ can be approximated by a parametric operator $\mathcal{P}_{NN}$ with arbitrary accuracy. That is, for any $\varepsilon>0$, there exists a sufficiently large parametric neural network $\mathcal{P}_{NN}$, such that
\[
\int_{\mathcal{X}} \int_{\Omega_x} | \mathcal{P}(f)(\boldsymbol{x}) - \mathcal{P}_{NN}(f)(\boldsymbol{x}) |^2 d \boldsymbol{x} d \mu(f) < \varepsilon\,.
\]
Here $\mu$ denotes a probability measure on $\mathcal{X}$.
In practice,  $\mu$ is chosen as a Gaussian measure induced by the law of a Gaussian random field. Several operator networks have been proposed recently, including DeepONets \cite{lu2019deeponet,wang2021learning} and various neural operators \cite{bhattacharya2020model,li2020fourier,kovachki2021neural}. In this paper, we adopt DeepONet as the basic architecture and refine it with transfer learning. The vanilla DeepONet  takes the following form:
\begin{align}\label{eqn:approx2}
	\mathcal{P}_{NN} (f)(&\boldsymbol{x};\theta, \xi) \nonumber \\
	&= \sum_{k=1}^p b^{NN}_k(f(\boldsymbol{y}_1), \cdots, f(\boldsymbol{y}_N); \theta) t^{NN}_k(\boldsymbol{x}; \xi) \nonumber
	\\ & =:\sum_{k=1}^p b^{NN}_k(f; \theta) t^{NN}_k(\boldsymbol{x}; \xi).
\end{align}
The operator network $\mathcal{P}_{NN}$ consists of two sub-networks: the branch net $b^{NN}$ is paramterized by $\theta$ and maps an encoded input function $\{f(\boldsymbol{y}_i)\}_{i=1}^N$ to $p$ scalars $b_k^{NN}$, and the trunk net $t^{NN} = \{t_k^{NN}\}_{k=1}^p$  is parameterized by $\xi$ and forms a directionary of functions in the output space. Both networks can be modified in practice to fit with various set-ups of PDEs, such as the boundary condition. A diagram of the DeepONets architecture we use in this paper is shown in Figure~\ref{fig:tl_pidon_arch}. The vanilla DeepONets is often trained in a supervised fashion and  requires pairs of input-output functions. To be more specific, given $N_s$ randomly sampled functions $\{ f_s(\boldsymbol{x}) \}_{s=1}^{N_s}$ one needs to prepare reference solutions $\{ \mathcal{P} (f_s)(\boldsymbol{x}) \}_{s=1}^{N_s}$ either analytically or using conventional high-fidelity numerical solvers. Then one trains the $\mathcal{P}_{NN}$ by minimizing the  loss function
\[
\begin{aligned}
	& \frac{1}{2N_s} \sum_{s=1}^{N_s}\Big( \| \mathcal{P}_{NN}(f_s)(\cdot; \theta, \xi) - \mathcal{P} (f_s)(\cdot) \|_{L^2(\Omega_x)}^2 \\
	\qquad & + \|\mathcal{P}_{NN}(f_s)(\cdot; \theta, \xi)  - \phi_s(\cdot)\|_{L^2(\partial \Omega_x)}^2\Big).
\end{aligned}
\]
However, in reality it can be extremely expensive to obtain the outputs $\mathcal{P}(f_s)$, especially when the underlying physical principles are complicated and the dimension of the problem is high. To this end, \cite{wang2021learning} proposed a physics-informed DeepONet which makes the  learning procedure above  self-supervised. More precisely, we turn to minimizing the new loss function
\begin{equation}
	\begin{aligned}
		& \frac{1}{2N_s} \sum_{s=1}^{N_s}\Big(  \| \mathcal{P}^{-1}\big(\mathcal{P}_{NN}(f_s)(\cdot; \theta, \xi)\big) -  f_s(\cdot) \|_{L^2(\Omega_x)}^2 \\
		\qquad & + \|\mathcal{P}_{NN}(f_s)(\cdot; \theta,\xi)  - \phi_s(\cdot)\|_{L^2(\partial \Omega_x)}^2\Big).  \label{eqn:emp_loss}
	\end{aligned}
\end{equation}
Note that the introduction of $\mathcal{P}^{-1}$ in \eqref{eqn:emp_loss} completely avoids the evaluations of $\mathcal{P}(f_s)$. The boundary term in \eqref{eqn:emp_loss} can be further eliminated in practice because the networks can be modified to satify the  boundary conditions (see e.g. \cite{lu2022comprehensive}). We also observe through numerical experiments that eliminating the  the boundary loss can substantially improves the training efficiency.
The physics-informed DeepONet has  been applied to learning evolution equations \cite{wang2021long}. For equation \eqref{111}, instead of first discretizing it in time, they consider time as an additional input variable, and try to learn an operator $\mathcal{P}^I$ that maps the initial condition to the solutions over an time-interval $[0,t_0]$:
\begin{align*}
	\mathcal{P}^I: f(0,\boldsymbol{x}) \mapsto f(t,\boldsymbol{x}), \quad \text{for} \quad t \in [0,t_0]\,.
\end{align*}
The corresponding loss function to be minimized is
\begin{align}
	L(\theta,  w,  \xi) & =  \frac{1}{2N_s} \sum_{s=1}^{N_s}   ( \| \partial_t (\mathcal{P}_{NN}^I (f_s))(t,\boldsymbol{x})  \nonumber \\ & - \mathcal L (\mathcal{P}_{NN}^I (f_s))(t,\boldsymbol{x}) \|_{L^2(\Omega_x \times [0,t_0])}^2 \nonumber \\
	& + \| \mathcal{P}_{NN}^I f_s(\boldsymbol{x}) - g_s(\boldsymbol{x}) \|_{L^2(\partial\Omega_x \times [0,t_0])}^2 ),  \label{eqn:bc_loss}
\end{align}
where $\mathcal{P}_{NN}^I$ is the neural network approximator to $\mathcal{P}^I$. Once trained, $\mathcal{P}_{NN}^I$ can be applied to $f(t_0,\boldsymbol{x})$ to get the solution $f(t,\boldsymbol{x})$ over $[t_0, 2t_0]$. Repeating this process enables one to obtain approximation solutions in any finite time. However, this methodology may suffer from long-time instability. In fact, let us illustrate this using the Allen-Cahn equation \eqref{eqn:ac} in one dimension.  Figure~\ref{fig:failure} shows that the average $L^2$ errors
of approximated solutions learned by  DeepONets using both the single-shot loss \eqref{eqn:emp_loss} (labeled as DeepONet) and the time-integrated loss \eqref{eqn:bc_loss} (labeled as CONT DeepONet). Both errors accumulate rapidly as time increases, indicating the instability of the learned DeepONets in prediction of long-time solution. In contrast, our transfer-learning assisted DeepOnet dramatically reduces the error and stabilizes the prediction. For completeness, we also compare them with the Fourier Neural Operator(FNO) \cite{li2020fourier}, which is another  state-of-the-art operator learning method.
\begin{figure}[th!]
	\centering
	{\includegraphics[width=0.3\textwidth, height=0.15\textheight]{ac_failure_new-eps-converted-to}}
	\caption{The relative $L^2$ error in time (defined in $\eqref{eqn:avg_l2_loss_t}$) for 1D Allen Cahn $\eqref{eqn:ac}$. The networks with ``CONT" refer to networks trained by using \eqref{eqn:bc_loss} whereas the others are trained by using  \eqref{eqn:emp_loss}. The prefix ``TL" means tuned by transfer learning. See implementation details in Appendix \ref{sec:detail_ac_ch}.}

\label{fig:failure}
\end{figure}

\subsection{DeepONet with transfer learning}

The main idea of transfer learning is to train a neural network on a large data set and then partially freeze and apply it to a related but unseen task. Inspired by \cite{desai2021one}, we employ the transfer learning technique to successively correct the trained DeepONet at the prediction steps:  we freeze the majority of the well-trained DeepONet and merely re-train the weights in the last hidden layer of the branch net by fitting the same physics-informed loss \eqref{eqn:emp_loss} defined by the underlying PDEs.
To be more precise, by separating the parameters $\theta$ in the hidden layers and the parameter $w$ in the last layer of the branch net, we rewrite the  branch net as
\begin{align*}
b_k^{NN}(f;\theta, w) = \sum_{j=1}^q w_j h_{k,j}(f;\theta)\,,
\end{align*}
where $h = \{h_{k,j}\}$ are the outputs of the last hidden layer of the branch net and $w = \{w_j\}$ are the  weights in the last layer. Inserting this into \eqref{eqn:approx2} gives
\begin{equation}\label{eqn:approx_trans}
\mathcal{P}_{NN} (f)(\boldsymbol{x};\theta, w, \xi) = \sum_{k=1}^p \sum_{j=1}^q w_j h^{NN}_{k, j}(f; \theta) t^{NN}_k(\boldsymbol{x}; \xi)\,.
\end{equation}
The architecture of the new operator network is illustrated in Figure~\ref{fig:tl_pidon_arch}.
\begin{figure*}[h!]
\centering
{\includegraphics[width=0.78\textwidth, height=0.22\textheight]{new_arc2.png}}
\caption{The architecture of transfer learning aided physics-informed DeepONet. Here $P$ and $D$ are optional layers that  enforce periodic and Dirichlet boundary conditions, respectively. The block named Modified FC is a modified fully connected neural networks architecture introduced in \cite{wang2021learning}. The parameter $w$ (in the red box) denotes the tunable weights in the last hidden layer of the branch net. In the transfer learning step, only $w$ will be re-trained while the $\theta, \xi$ are frozen. }
\label{fig:tl_pidon_arch}
\end{figure*}
In the training step, the optimal  parameters $(\theta^*, w^*,\xi^*)$ of the DeepONet \eqref{eqn:approx_trans} can be obtained by minimizing the empirical loss \eqref{eqn:emp_loss}. Later in each prediction step, we freeze the value of $\theta^*$ and $\xi^*$, but update $w^\ast$ by re-training the loss \eqref{eqn:emp_loss} with newly-predicted solution as the initial condition. Namely with the predicted solution $f_n$ at step $n$, we seek $w^\ast_{n+1}$ defined by
\begin{align}\label{eqn:opt_w}
&  w^*_{n+1}  \in \arg \min_w  \nonumber \\ & \frac{1}{2N_s} \sum_{s=1}^{N_s}\Big(  \| \mathcal{P}^{-1}  \mathcal{P}_{NN} (f)(\boldsymbol{x};\theta^\ast, w, \xi^\ast)  - f_n(\boldsymbol{x}) \|_{L^2(\Omega_x)}^2 \nonumber \\
+ & \|  \mathcal{P}_{NN} (f)(\boldsymbol{x};\theta^\ast, w, \xi^\ast)  - \phi(\boldsymbol{x})\|_{L^2(\partial \Omega_x)}^2\Big), n=1,2,\cdots.
\end{align}
Note that $w^\ast_1 = w^\ast$. The optimal sequence of weights $w_n^\ast$ defines a sequence of operator networks $\mathcal{P}^n_{NN} := \mathcal{P}^{n}_{NN}(\theta^\ast, w^\ast_n, \xi^\ast)$, which can be used to approximate the solution at $t=n\Delta t$ by
$$
f(n\Delta t) \approx \mathcal{P}^{n}_{NN}\circ \mathcal{P}^{n-1}_{NN}\circ \mathcal{P}^{1}_{NN} (f_0).
$$
It is interesting to note that the proposed method shares some similarities with the classical Galerkin approximation. In fact, the operator network \eqref{eqn:approx_trans} can be further rewritten as
\begin{align*}
\mathcal{P}_{NN} (f)(\boldsymbol{x}) & = \sum_{j=1}^q w_j \left( \sum_{k=1}^p  h^{NN}_{k, j}(f; \theta) t^{NN}_k(\boldsymbol{x}; \xi) \right) \\
&=: \sum_{j=1}^q w_j \phi_j(\boldsymbol{x}; f)\,,
\end{align*}
Observe that $\phi_j$ playing the role of basis functions in Galerkin methods, and $w_j$ being the corresponding weight.
However, unlike most Galerkin methods which often use  handcraft bases, such as piecewise polynomials and trigonometric functions, here the bases are learned from the problem itself, and vary with the function they approximate. This seemly minor change reduces substantially the number of bases needed in the output space, as shown by extensive numerical tests in Section \ref{sec:numerical_exp}. To minimize \eqref{eqn:opt_w}, it amounts to solving a system of $N$ equations with $q$ unknowns, where $N$ is the total number of fixed sensors in the branch net. This is achieved by least square minimization. Since $q \ll N$, the computational complexity of finding the least square solution is only $\mathcal{O}(q^2 N)$.  In practice, we further reduce the computational complexity by sub-sampling $N_c$ grid points out of $N$ in the transfer learning step.

\section{Theoretical result}\label{sec:thm}
In this section, we analyze the long time stability of the learned operator $\mathcal{P}_{NN}$. First let $\mathcal{X} $ be a Banach space and assume that the original propagator $\mathcal{P}$ (i.e., $\mathcal{P}^{\Delta t}$ in \eqref{eqn:mP}): $\mathcal{X} \rightarrow \mathcal{X}$ is non-expansive such that
\begin{align} \label{assump}
\norm{\mathcal{P}}_{\mathcal{X}}:= \sup_{f\in \mathcal{X}, \norm{f}_{\mathcal{X}} =1} \norm{\mathcal{P} f}_\mathcal{X}  \leq 1.
\end{align}
In the case that $\mathcal{X} = L^2(\Omega_x)$, assumption\eqref{assump} follows from the dissipative assumption of $\mathcal{L}$; see Appendix \ref{sec:assump} for more details. Let $\mathcal{U} \subseteq \mathcal{X}$ be a linear subspace, we also assume that
\begin{align}\label{assmp2}
\mathcal{P} f \in \mathcal{U}, \quad \forall f \in \mathcal{U}\,.
\end{align}
The theorem below shows that the long-time prediction error of the operator network can be bounded by the loss function.
\begin{theorem}\label{thm:main}
Assume \eqref{assump} and \eqref{assmp2} hold. If the neural network approximator $\mathcal{P}^{NN}$ satisfies that the maximum  loss over the set $\mathcal U$ is less than $\delta$, i.e.
\begin{equation} \label{loss0}
\sup_{f\in \mathcal U, \|f\|_{\mathcal{X}}=1 } \norm{\mathcal{P}^{-1} \mathcal{P}_{NN} f - f}_\mathcal{X} \leq \delta\,,
\end{equation}
and that
\begin{equation} \label{eq:Pnnf}
\mathcal{P}_{NN} f \in \mathcal{U}, \quad \forall f \in \mathcal{U},
\end{equation}
then the following long-time stability holds
\begin{equation} \label{sta1}
\sup_{f\in\mathcal{U}, \|f\|_{\mathcal{X}}=1}	\norm{(\mathcal{P})^K f - (\mathcal{P}_{NN})^K f}_\mathcal{X} \leq \delta K(1+\delta )^K.
\end{equation}
Moreover, if we further assume that
\begin{equation} \label{assump2}
\norm{\mathcal{P}}_{\mathcal{X}} \leq \eta < 1
\end{equation}
and  \eqref{loss0} holds with $\delta \leq \frac{1}{2} (1-\eta)$, then we have
\begin{equation} \label{est2}
\sup_{f\in\mathcal{U}, \|f\|_{\mathcal{X}}=1} 	\norm{(\mathcal{P})^K f - (\mathcal{P}_{NN})^K f}_\mathcal{X} \leq \delta K \Big(\frac{ 1+\eta}{2}\Big)^{K-1}\,.
\end{equation}
\end{theorem}

The proof of Theorem \ref{thm:main} is provided in Appendix~\ref{pf}.

\begin{remark}

\begin{itemize}
\item[1.] If the error tolerance $\delta = \Delta t^2$ with $\Delta t$ being the time-discretization stepsize and the number of iterations $K = \frac{T}{\Delta t}$, then  \eqref{sta1} becomes
$$
\sup_{f\in\mathcal{U}, \|f\|_{\mathcal{X}}=1}	\norm{(\mathcal{P})^K f - (\mathcal{P}_{NN})^K f}_\mathcal{X} \leq e^{T\Delta t} T \Delta t \,.
$$
When  assumption \eqref{assump2} holds, the estimate above improves to
$$\begin{aligned}
	\sup_{f\in\mathcal{U}, \|f\|_{\mathcal{X}}=1} 	\norm{(\mathcal{P})^K f - (\mathcal{P}_{NN})^K f}_\mathcal{X} & \leq \big(\frac{1+\eta}{2}\big)^{\frac{T}{\Delta t}} T \Delta t\\
	& \leq  C\Delta t,
\end{aligned}
$$
where $C$ is a constant independent of $T$ and $\Delta t$, suggesting that the prediction error is of the order $\mathcal{O}(\Delta t)$ uniformly in time.

\item[2.] We comment on the assumptions made in Theorem \ref{thm:main}. In practice, physics-informed loss \eqref{eqn:opt_w} is trained so that condition \eqref{loss0}  is fulfilled for some subspace space $\mathcal{U}$, e.g. $\mathcal{U} = \{e^{i \mathbf{k}\cdot \mathbf{x}}\}_{|\mathbf{k}| \leq K_0}$. Assumption \eqref{assmp2}  holds for such a choice of $\mathcal{U}$ when the evolution equation involves diffusion. Assumption \eqref{eq:Pnnf} holds in particular when sine or cosine activation function is used in the operator network $\mathcal{P}_{NN}$.
\end{itemize}

\end{remark}

\section{Numerical experiments}\label{sec:numerical_exp}
In this section, we demonstrate the effectiveness of transfer learning enhanced DeepONet and show its advantages over the vanilla DeepONet through several evolutionay PDEs, including reaction diffusion equation, Allen-Cahn and Cahn-Hilliard equations, Navier-Stokes equation and multiscale linear radiative transfer equation.
The equations of consideration are equipped with either Dirichlet or periodic boundary conditions. In all the test problems, our goal is to predict the long time evolution of the equations obtained by successive actions of the propagators learned via DeepONets. More concretely, we first build the first-step neural operator approximation $\mathcal{P}^1_{NN}$ to the propagator $\mathcal{P} = \mathcal{P}^{\Delta t}$ by minimizing the physics-informed loss \eqref{eqn:emp_loss} with $M$ training initial data. The operator network $\mathcal{P}^1_{NN}$ is then gradually tuned to $\mathcal{P}^j_{NN},j=2,\cdots K$ via updating the weights $w$ in the last-layer of its trunk nets. With the learned (and adjusted) operators $\mathcal{P}^j_{NN},j=1,\cdots K$, the solution of a PDE at  time $t=K\Delta t$ with an initial condition $f_0$ can then be obtained approximately by $\mathcal{P}^K_{NN}\circ \cdots \circ \mathcal{P}^1_{NN} f_0$.  We remark that  the $M$ training data is constructed as a subset of a larger  training set of size $N_s\times N_p$, which consists of pointwise evaluations of $N_s$ randomly sampled functions at $N_p$ physical locations (sensors). We refer to Appendix \ref{sec:data} for detailed discussions on the data generating process and treatment of boundary conditions in various test problems. Choices of parameters for the operator networks and the training process are discussed in the end of Appendix \ref{sec:data}.   In all the numerical results to follow, we quantify the performance of the proposed method by measuring the aggregated relative prediction error over a time horizon $[0,T]$; see the precise definition of the relative error in Appendix \ref{sec:data}.  The codes used for the numerical experiments will be published on the website  \url{https://github.com/woodssss/TL-PI-DeepONet}.

\subsection{Reaction diffusion equation}
Consider the reaction diffusion equation
\begin{equation}\label{eqn:nrd}
\left\{ \begin{array}{l}
\partial_t f = d  \Delta f + k f^2, x\in \Omega_x:= [0,1]
\\ f(t,\boldsymbol{x}) = 0, ~~\boldsymbol{x} \in \partial \Omega_x = \{0,1\}\,,
\\ f(0,\boldsymbol{x}) = f_0(\boldsymbol{x}) \,,
\end{array}\right.
\end{equation}
where $d=k=0.001$. In this example, we train DeepOnets and our transfer learning enhanced DeepONets using two different loss functions, with one based on the physics-informed loss within a single time-step $\Delta t = 0.05$ (c.f. \eqref{eqn:emp_loss}), and the other based on the aggregated physics-informed loss \eqref{eqn:bc_loss} in the time window $t\in [0,1]$. We refer to the DeepONets trained using the latter loss as {\em CONT DeepONets} and reserve {\em DeepONets} for the one trained by the former loss.  The numerical results of different DeepONets with varying training sample sizes $M$ are shown in Table~\ref{table:nrd}. Our proposed method provides more accurate and more robust prediction of the solutions with little extra computational cost. In particular, the prediction error of vanilla DeepONets and CONT DeepONets increase dramatically as time increases from $0.2$ to $50$, while transfer learning can significantly reduces the error and stabilizes the prediction in the long time.
In addition, our method also substantially reduces the size of training data to achieve the same order of prediction accuracy.
Note that, since trained only within a single time step, DeepONets take far less training time than the corresponding CONT DeepONets while maintain comparable accuracy. The similar trade-off of accuracy and training cost applies to other experiments. For this reason, in subsequent examples we will only report results on our proposed method and the vanilla DeepONet, and exclude the results from CONT DeepONets. Note also that the results obtained in Table ~\ref{table:nrd} are for propagators defined by the backward Euler scheme. One can also consider propagators defined by higher order time-discretization schemes and their neural network approximation. We refer the numerical results obtained using the Crank-Nicolson method to Table ~\ref{table:nrd_2nd} in Appendix \ref{sec:reaction}.

\begin{table}[!ht]
\centering
\begin{tabular}{l|lll|ll}
\toprule
\cmidrule(r){1-2}
Neural network    & $M$  & $t_1$ & $t_2$ &  $T=0.2$ & $T=50$ \\
\midrule
CONT  & 1000  & 43030  & 0.26 & 7.95e-2 & 4.17e-1     \\
DeepONet & 3000   & 79375 & 0.28 & 1.01e-2 & 4.15e-2    \\
& 10000  & 83465  & 0.26 & 3.40e-3 & 1.34e-2      \\
CONT    & 1000  &43030 & 0.93 & \textbf{5.38e-3} & \textbf{1.95e-3}    \\
TL-DeepONet & 3000   & 79375 & 0.98 & \textbf{1.87e-3} & \textbf{6.88e-4}   \\
& 10000  &83465 & 0.93 & \textbf{1.84e-3} & \textbf{4.87e-4}     \\
DeepONet  & 1000  & 2575  & 5.02 & 2.75e-1 & 1.69e0     \\
& 3000   & 4313 & 4.9 & 1.05e-1 & 1.60e0    \\
& 10000  & 5854 & 4.8 & 9.19e-2 & 1.87e0      \\
TL-DeepONet & 1000  & 2575 &  9.86 & \textbf{1.03e-3} & \textbf{1.52e-3}    \\
& 3000   & 4313 & 7.1 & \textbf{8.34e-4}  & \textbf{1.19e-3}    \\
& 10000  & 5854 & 8.2  &  \textbf{8.19e-4} &  \textbf{9.05e-4}     \\
\bottomrule
\end{tabular}
\caption{Results on reaction diffusion equation. Here $t_1$ is the training time and $t_2$ is the averaged  time of predicting the solution trajectories among the time interval $[0, 50]$ based on $30$ test initial conditions. The last two columns to the right are the averaged relative $L^2$ error within $[0,T]$.}
\label{table:nrd}
\end{table}

\subsection{Allen-Cahn and Cahn-Hilliard equations}
In the second example, we consider Allen-Cahn equation
\begin{equation}\label{eqn:ac}
\left\{ \begin{array}{l}
\partial_t f =  d_1 \Delta f +   d_2 f(1-f^2), \\
f(0, \boldsymbol{x}) = f_0(\boldsymbol{x}) ,
\end{array}\right.
\end{equation}
and Cahn-Hilliard equation
\begin{equation}\label{eqn:ch}
\left\{ \begin{array}{l}
\partial_t f =   \Delta g, \\
g = - d_1 \Delta f +   d_2 (f^3 - f), \\
f(0, \boldsymbol{x}) = f_0(\boldsymbol{x}),
\end{array}\right.
\end{equation}
both equipped with periodic boundary conditions. They are prototype models for the motion of anti-phase boundaries in crystalline solids. The computational domain is $ \Omega := [0, 1]^d$ with $d=1,2$. We are interested in learning the propagator $\mathcal{P} = \mathcal{P}^{\Delta t}$ with $\Delta t = 0.05$ and used it to predict the solutions $f(t, \boldsymbol{x})$ for every $t\leq T=50$. The results on Allen-Cahn equation are shown in Table~\ref{table:ac} (1D) and Table~\ref{table:ac_2d} (2D). See also Figure~\ref{fig:failure} for a plot of evolving relative errors on 1D Allen-Cahn equation. Similar results for 1D Cahn-Hilliard equation are presented in
Table~\ref{table:ch} and Figure~\ref{fig:ch2d} compares the snapshots of predicted solutions to the 2D Cahn-Hilliard equation. In all the results, for a fixed $d_2$, the relative errors increase as $d_1$ decreases because the transition layers of solutions are increasingly sharper and hence make the numerical resolution more challenging. Similar to the previous example, our proposed method provides more accurate prediction of solutions than the vanilla DeepONets among all the configurations of parameters. We note that the average trajectory prediction times of TL-DeepONets increase for about 3 times compared to those of the vanilla DeepONets while the prediction errors of the former decrease by at least two orders of magnitude.
\begin{table}[h]
\centering
\begin{tabular}{l|llll}
\toprule
\cmidrule(r){1-2}
Neural network    & $M$    &  $d_1$=1e-3 &  $d_1$=5e-4 &  $d_1$=1e-4  \\
\midrule
DeepONet  & 1000   &  1.13e0 & 1.33e0 &  1.29e0 \\
& 3000   & 1.33e0 & 1.18e0   &  1.23e0 \\
& 10000  & 1.01e0 &  8.95e-1  & 1.43e0     \\
TL-DeepONet  & 1000  & \textbf{9.25e-4} & \textbf{9.64e-4}  & \textbf{2.16e-2}   \\
& 3000 & \textbf{7.78e-4}  & \textbf{8.83e-4} &   \textbf{1.81e-2}   \\
& 10000 & \textbf{5.81e-4} & \textbf{7.94e-4} & \textbf{1.16e-2}       \\
\bottomrule
\end{tabular}
\caption{Results on 1D Allen-Cahn equation: the time-average of relative prediction errors within $[0,50]$. The average trajectory prediction time is 5.1s for DeepONet and 16.1s for TL-DeepONet.}
\label{table:ac}
\end{table}

\begin{table}[h]
\centering
\begin{tabular}{l|llll}
\toprule
\cmidrule(r){1-2}
Neural network    & $M$    &  $d_1$=4e-3 &  $d_1$=2e-3 &  $d_1$=1e-3  \\
\midrule
DeepONet   & 1000   &  9.96e-1 & 1.01e0 &  1.02e0 \\
& 10000   & 9.96e-1 & 1.00e0   &  1.00e0 \\
TL-DeepONet   & 1000  & \textbf{6.54e-3} & \textbf{8.43e-3}  & \textbf{1.01e-2}   \\
& 10000 & \textbf{4.96e-3}  & \textbf{6.46e-3} &   \textbf{9.01e-3}   \\
\bottomrule
\end{tabular}
\caption{Results on 2D Allen-Cahn equation: the time-average of the relative prediction errors within $[0,10]$. The average trajectory prediction time is 6.1s for DeepONet and 29.5s for TL-DeepONet.}
\label{table:ac_2d}
\end{table}

\begin{table}[h]
\centering
\begin{tabular}{l|llll}
\toprule
\cmidrule(r){1-2}
Neural network    & $M$    &  $d_1$=4e-6 &  $d_1$=2e-6 &  $d_1$=1e-6  \\
\midrule
DeepONet  & 1000 & 9.43e-1  & 9.47e-1 & 9.44e-1 \\
& 3000   & 9.54e-1  & 9.58e-1  & 9.39e-1 \\
& 10000  & 9.65e-1 & 9.58e-1  & 9.38e-1    \\
TL-DeepONet  & 1000  & \textbf{1.04e-2}  & \textbf{1.36e-2}  & \textbf{4.25e-2}  \\
& 3000 & \textbf{8.11e-3} & \textbf{9.04e-3} & \textbf{3.86e-2}  \\
& 10000 &\textbf{2.29e-3}  &\textbf{7.89e-3}  &   \textbf{3.03e-2}   \\
\bottomrule
\end{tabular}
\caption{Results on 1D Cahn-Hilliard equation: the time-average of the relative prediction errors within $[0,50]$. the average trajectory prediction time is 4.4s for DeepONet and 12.3s for TL-DeepONet.}
\label{table:ch}
\end{table}

\begin{figure}
\centering
\setlength{\tabcolsep}{-2pt}
\renewcommand{\arraystretch}{-1}
\begin{tabular}{ccc}
\text{Initical condition} \quad & t=4 & t=10\\
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_ref_tz-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_ref_t4-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_ref_t1z-eps-converted-to}}\\
\text{DeepONet}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_DON_t4-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_DON_t1z-eps-converted-to}}\\
\text{TL-DeepONet}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_TL_t4-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{CH2D_TL_t1z-eps-converted-to}}\\
\end{tabular}
\caption{Results on 2D Cahn-Hilliard equation: snapshots of reference solutions (top), and of approximate solutions predicted by   DeepONet (middle) and TL-DeepONet (bottom).}
\label{fig:ch2d}
\end{figure}

\subsection{Navier-Stokes equation}
Consider the 2D Navier-Stokes equation in the vorticity form:
\begin{equation} \label{eqn:ns}
\left\{ \begin{array}{l}
\partial_{t} w(\boldsymbol{x}, t)+u(\boldsymbol{x}, t) \cdot \nabla w(\boldsymbol{x}, t)=\nu \Delta w(\boldsymbol{x}, t)+  f(\boldsymbol{x}) , \\
w(\boldsymbol{x}, 0) = w_0(\boldsymbol{x})
\end{array}\right.
\end{equation}
with periodic boundary condition and source $f(\boldsymbol{x}) = 0.1 (\sin(2 \pi (x + y)) + \cos(2 \pi (x + y)))$. We would like to learn the propagator $\mathcal{P}^{\Delta t}$ with $\Delta t = 0.01$ and apply it to predict the solution $w|_{(0,1)^2 \times (\Delta t, T)}$. Table~\ref{table:ns} shows the results with varying values of viscosity $\nu$. Note that the prediction error increases as $\nu$ decreases. TL-DeepONets reduces the errors of DeepONets by two orders of magnitudes although the prediction time of the former increases for less than 4 times. Figure~\ref{fig:ns_zpzz1} shows the snapshots of solutions to \eqref{eqn:ns} with $\nu=0.001$ at two different times.

\begin{table}[h]
\centering
\begin{tabular}{l|llll}
\toprule
\cmidrule(r){1-2}
Neural network     &  $\nu$=1e-1 &  $\nu$=1e-2 &  $\nu$=1e-3 &  $\nu$=1e-4 \\
\midrule
DeepONet   & 9.95e-1    & 1.02e0 &  9.96e-1 &  1.04e0   \\
TL-DeepONet  & \textbf{1.41e-2} &\textbf{1.07e-2} & \textbf{3.35e-2}  & \textbf{9.42e-2}   \\
\bottomrule
\end{tabular}
\caption{Results on 2D Navier-Stokes equation: the relative prediction errors within $[0, 10]$. The average trajectory prediction time is 5.3s for DeepONet and 24.8s for TL-DeepONet.}
\label{table:ns}
\end{table}

\begin{figure}[ht]
\centering
\begin{center}
\setlength{\tabcolsep}{-1pt}
\renewcommand{\arraystretch}{-1}
\begin{tabular}{ccc}
\text{Initical condition} \quad & $t=4$ &  $t=10$\\
\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_ref_tz-eps-converted-to}&
\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_ref_t4-eps-converted-to}&
\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_ref_t_end-eps-converted-to}\\
\text{DeepONet}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_DON_t4-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_DON_t_end-eps-converted-to}}\\
\text{TL-DeepONet}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_TPIDON_t4-eps-converted-to}}&
\raisebox{-.5\height}{\includegraphics[width=0.16\textwidth, height=0.1\textheight]{NS_nu_zpzz1_TPIDON_t_end-eps-converted-to}}\\
\end{tabular}
\end{center}
\caption{Results on Navier-Stokes equation with $\nu=0.001$: snapshots of reference solutions (top), and of approximate solutions predicted by DeepONet (middle) and TL-DeepONet (bottom).}
\label{fig:ns_zpzz1}
\end{figure}

\subsection{Multiscale linear radiative transfer equation}
Consider the  linear multiscale radiative transfer equation:
\begin{equation}\label{eq:rte}
\left\{ \begin{array}{l}
\varepsilon \partial_t f + \boldsymbol{v} \cdot  \nabla f = \frac{1}{\varepsilon} \mathcal{L} f, ~ t \in [0, T], ~ (\boldsymbol{x}, \boldsymbol{v}) \in \Omega_x \times \mathcal{S}^{d-1}\,, \\
f(t, \boldsymbol{x}, \boldsymbol{v}) = \phi(\boldsymbol{x}) , ~ (\boldsymbol{x}, \boldsymbol{v}) \in  \Gamma_- \,,\\
f(0, \boldsymbol{x}, \boldsymbol{v}) = f_0(\boldsymbol{x}, \boldsymbol{v}).
\end{array}\right.
\end{equation}
Here  $\varepsilon>0$ is the Knudsen number which is a dimensionless parameter that determines the physical regime of the equation, $\displaystyle \mathcal{L}(f) = \frac{1}{|\mathcal{S}^{d-1}|} \int_{\mathcal{S}^{d-1}}  f d \boldsymbol{v} - f =: \average{f} - f$, and $\Gamma_-=\{(x,v): x\in \partial \Omega_x,~ v\cdot n_x <0\}$ is the inflow part of the boundary. In this example, we aim to learn the propagator $\mathcal{P}^{\Delta t}$ with $\Delta t =0.01$ and employ it to predict the solution $f(t, \boldsymbol{x}, \boldsymbol{v})$ for  $t\in [0,10]$. We mainly consider \eqref{eq:rte} in one and two physical dimensions and refer to Appendix \ref{sec:detailsrte} for a detailed discussion on the experiment set-up and the numerical method. Table~\ref{tab:rte} displays the results corresponding to different Knudsen numbers. The transfer learning enhanced DeepONets reduces the relative error by one or two orders of magnitude although increase the prediction time by around 4 times. Figure~\ref{fig:rte_eps_1}  and Figure~\ref{fig:rte_eps_zpzzz1} in Appendix \ref{sec:detailsrte} show several snapshots of solutions to \eqref{eq:rte} with $\varepsilon=1$ and $\varepsilon=10^{-4}$ respectively.

\begin{table}[ht!]
\centering
\begin{tabular}{l|lll}
\toprule
\cmidrule(r){1-2}
RTE & Neural network & $t_2$ & relative error  \\
\midrule
1D $\varepsilon =$ 1 & DeepONet  & 4.6  & 3.06e-1    \\
& TL-DeepONet  & 22.5 & $\textbf{1.52e-2}$      \\
1D $\varepsilon =$ 1e-4 & DeepONet  &5.1  & 3.74e-1    \\
& TL-DeepONet  & 21.9 & $\textbf{5.52e-3}$      \\
2D $\varepsilon =$ 1 & DeepONet   & 79.8 & 3.58e-1   \\
& TL-DeepONet  &   431.3 & $\textbf{2.19e-2}$   \\
2D $\varepsilon =$ 1e-4 & DeepONet   &83.1 & 2.37e0   \\
& TL-DeepONet  & 379.3 &  $\textbf{8.93e-3}$   \\
\bottomrule
\end{tabular}
\caption{Results on the radiative transfer equation: the relative prediction errors over the time-horizon $[0,10]$. }
\label{tab:rte}
\end{table}

\section{Conclusion}
In this paper, we proposed a new physics-informed DeepONet based on transfer learning for learning evolutionary PDEs. This is achieved in two steps: first learn the propagators and then predict the solutions by  successive actions of propagators on the initial condition. The experimental results demonstrated that the proposed method improves substantially upon the vanilla DeepONet in terms of long-time accuracy and stability while maintains low computational cost. The proposed method also reduced the training sample size needed to achieve the same order of prediction accuracy of the vanilla DeepONets.

\section*{Acknowledgement}
L. Wang is partially supported by
NSF grant DMS-1846854. Y. Lu thanks NSF for the support via the award DMS-2107934.


\section{Introduction\label{intro}}

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

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

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

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

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

The remainder of this paper is organized as follows.
In Sec.~\ref{Sec2}, we introduce the model to be analyzed in this study.
In Sec.~\ref{Sec3}, we summarize the theoretical framework for determining the rigidity of an amorphous solid consisting of frictional grains without considering Coulomb's slip process.
In Sec.~\ref{results}, we present the results of the stress-strain relation obtained using the theory formulated in Sec.~\ref{Sec3}.
We also compare the theoretical results with the simulation results to demonstrate the relevancy of our theoretical analysis.
In Sec.~\ref{concluding_remarks}, we summarize the obtained results and address future tasks to be solved.
In Appendix~\ref{App:Absence_of_2nd}, we numerically demonstrate the absence of history-dependent tangential deformations in our system.
In Appendix~\ref{app:gamma_Th}, we explain the detailed behavior of the eigenvalue near the stress-drop points.
In Appendix~\ref{AppHessian}, we present some detailed properties of the Hessian matrix in a harmonic system.
In Appendix~\ref{AppHarmonic}, we present some properties of the Jacobian matrix in a harmonic system and demonstrate its equivalency to the Hessian matrix.
Appendix~\ref{app:rigidity} presents the detailed properties of rigidity.

\section{Model \label{Sec2}}

Our system contains $N$ frictional circular disks embedded in a two-dimensional space.
To prevent the system from crystallization, it contains an equal number of grains with diameters $d$ and $d/1.4$ \cite{Luding01}.
We assume that the mass of grain $i$ is proportional to $d_i^2$, where $d_i$ is the diameter of $i$-th grain.
For later convenience, we introduce $m$ as the mass of grain having a diameter $d$.
In this study, $x_{i}, y_{i}$, and $\theta_{i}$ denote $x$ and $y$ components of the position of $i$-th grain, and the rotational angle of the $i$-th grain, respectively.
We introduce the generalized coordinates of the $i$-th grain as follow:
\begin{equation}\label{def_q}
\bm{q}_{i}:=(\bm{r}_{i}^{\textrm{T}},\ell_{i})^{\textrm{T}},
\end{equation}
where $\bm{r}_{i}:=(x_{i},y_{i})^{\textrm{T}}$, $\ell_{i}:=d_{i}\theta_{i}/2$, and the superscript $\textrm{T}$ denotes the transposition.

Let the force, and $z$-component of the torque acting on the $i$-th grain be $\bm{F}_{i}:=(F_{i}^{x},F_{i}^{y})^{\textrm{T}}$ and $T_{i}$, respectively.
Then, the equations of motion of $i$-th grain are expressed as
\begin{align}
m_{i}\frac{d^{2}\bm{r}_{i}}{dt^{2}} &=  \bm{F}_{i}, \label{tra}\\
I_{i}\frac{d^{2}\theta_{i}}{dt^{2}} &= T_{i}, \label{rot}
\end{align}
with mass $m_i$ and momentum of inertia $I_i:=m_i d_i^2/8$ of $i$-th grain.
In a system without volume forces, such as gravity, we can write
\begin{align}
\bm{F}_{i}&=\sum_{j\neq i} \bm{f}_{ij} - m_{i}\eta_{D}\dot{\bm{r}}_{i}, \label{Fi}\\
T_{i}&=\sum_{j\neq i}T_{ij} - I_{i}\eta_{D}\dot{\theta}_{i}, \label{Ti}
\end{align}
where we have adopted the notation $\dot{A}:=dA/dt$ for arbitrary variable $A$ such as $A=\bm{r}_i$ and ${\theta}_i$.
Here, $\bm{f}_{ij}$ and $T_{ij}$ are the force and $z$-component of the torque acting on the $i$-th grain from the $j$-th grain, respectively.
As a simplified description of the drag terms from the background fluid,
$\eta_D$ is a damping constant uniformly acting on grains, and
$T_{ij}$ is given by
\begin{align}
T_{ij}=-\frac{d_{i}}{2}( {n}_{ij}^{x}{f}_{ij}^{y} - {n}_{ij}^{y}{f}_{ij}^{x} ),
\label{Torque}
\end{align}
where we have introduced the normal unit vector between $i$-th and $j$-th grains as $\bm{n}_{ij}:=\bm{r}_{ij}/|\bm{r}_{ij}|:=(\bm{r}_{i}-\bm{r}_{j})/|\bm{r}_{i}-\bm{r}_{j}|$.
The force $\bm{f}_{ij}$ can be divided into the normal part $\bm{f}_{N,ij}$ and tangential part $\bm{f}_{T,ij}$ as
\begin{align}
\bm{f}_{ij}=(\bm{f}_{N,ij}+\bm{f}_{T,ij} )H(d_{ij}/2-|\bm{r}_{ij}|)
\label{fijSum}
,
\end{align}
where $d_{ij}:=d_i+d_j$ and $H(x)$ is Heaviside's step function, taking $H(x)=1$ for $x>0,$ and $H(x)=0$ otherwise.
We assume that the contact force is expressed as
\begin{align}
\bm{f}_{N,ij}:&=k_{N}\xi_{N,ij}\bm{n}_{ij}
\label{fijn}
,\\
\bm{f}_{T,ij}:&=k_{T}\xi_{T,ij}\bm{t}_{ij}
\label{fijt}
,
\end{align}
where $k_N$ and $k_T$ are the stiffness parameters of normal and tangential contacts, respectively.
The contact force can be derived from the harmonic potential.
In Eqs.~\eqref{fijn} and ~\eqref{fijt} we have introduced
$\xi_{N,ij}:=d_{ij}/2-|\bm{r}_{ij}|$ and
\begin{align}
\bm{\xi}_{T,ij}(t):&=\int_{C_{ij}(t')}dt'\bm{v}_{T,ij}(t')
\notag\\
&\quad
 -\left[\left( \int_{C_{ij}(t')}dt'\bm{v}_{T,ij}(t') \right)\cdot\bm{n}_{ij}(t)\right]\bm{n}_{ij}(t) ,
\label{xiT}
\end{align}
where we have used
\begin{align}
\bm{v}_{T,ij}:=\dot{\bm{r}}_{ij}-\dot{\bm{\xi}}_{N,ij}+\bm{u}_{ij}(d_{i}\dot{\theta}_{i}+d_{j}\dot{\theta}_{j})/2
\end{align}
with $\bm{u}_{ij}:=(n_{ij}^{y}, -n_{ij}^{x})^{\textrm{T}}$, $\bm{t}_{ij}:=-\bm{\xi}_{T,ij}/|\xi_{T,ij}|$,  $\xi_{T,ij}:=|\bm{\xi}_{T,ij}|$,
and the integration over the duration time of contact between $i$-th and $j$-th grains $\int_{C_{ij}(t')}dt'$ with the trajectory $C_{ij}(t')$ of the contact point between $i$-th and $j$-th grains at $t'$.
As shown in Appendix~\ref{App:Absence_of_2nd}, we have confirmed that the second term on the right-hand side (RHS) of Eq.~\eqref{xiT} is zero for harmonic systems, although we do not have any mathematical proof for this statement thus far.
For simplicity, we consider neither the effects of Coulomb's slip in the tangential equation of motion nor the dissipative contact force,
where the tangential forces $\bm{f}_{T,ij}^{C}$ including the effect of Coulomb's slip process, satisfy
\begin{align}
\bm{f}_{T,ij}^{C}
:=
\left\{
\begin{matrix}
\bm{f}_{T,ij} &(\bm{f}_{T,ij} < \mu|\bm{f}_{N,ij}|)\\
\mu|\bm{f}_{N,ij}|\bm{t}_{ij} &(\textrm{otherwise})
\end{matrix}
\right.
\end{align}
where $\mu$ is the friction constant.

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

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

We introduce the rate of nonaffine displacements for $r_i^{{\rm FB},\zeta}(\gamma)$ with $\zeta=x$ or $y$ and $\ell_i^{\rm FB}(\gamma)$ as:
\begin{align}
\frac{d\mathring r_{i}^{\zeta}(\gamma)}{d\gamma}:&=
\lim_{\Delta\gamma\to0}\frac{r_{i}^{\textrm{FB},\zeta}(\gamma+\Delta\gamma)-r_{i}^{\textrm{FB},\zeta}(\gamma)}{\Delta\gamma}-\delta_{\zeta x} y_{i}^{\textrm{FB}}(\gamma), \label{drdG}
\\
\frac{d\mathring \ell_{i}(\gamma)}{d\gamma}:&=
\lim_{\Delta\gamma\to0}\frac{\ell_{i}^{\textrm{FB}}(\gamma+\Delta\gamma)-\ell_{i}^{\textrm{FB}}(\gamma)}{\Delta\gamma}
. \label{dldG}
\end{align}

Our system is characterized by the generalized coordinate $\bm{q}(\gamma):=(\bm{q}_{1}^{\textrm{T}}(\gamma), \bm{q}_{2}^{\textrm{T}}(\gamma), \cdots, \bm{q}_{N}^{\textrm{T}}(\gamma))^{\textrm{T}}$.
The configuration in the FB state at strain $\gamma$ is denoted by $\bm{q}^{\textrm{FB}}(\gamma):=((\bm{q}_{1}^{\textrm{FB}}(\gamma))^{\textrm{T}}, (\bm{q}_{2}^{\textrm{FB}}(\gamma))^{\textrm{T}},\cdots,(\bm{q}_{N}^{\textrm{FB}}(\gamma))^{\textrm{T}})^{\textrm{T}}$.
The shear stress $\sigma(\gamma)$ at $\bm{q}^{\textrm{FB}}(\gamma)$ for one sample is given by:
\begin{align}
\sigma(\bm{q}^{\textrm{FB}}(\gamma))=
-\frac{1}{2L^{2}} \sum_{i}\sum_{j>i}
&\left[ f_{ij}^{x}(\bm{q}^{\textrm{FB}}(\gamma)) r_{ij}^{y}(\bm{q}^{\textrm{FB}}(\gamma))
\right.
\nonumber \\
&\ +\left. f_{ij}^{y}(\bm{q}^{\textrm{FB}}(\gamma)) r_{ij}^{x}(\bm{q}^{\textrm{FB}}(\gamma))
\right]
\label{Sigma}
.
\end{align}
The rigidity $g$ for one sample is defined as
\begin{align}
g&:= \left. \frac{d\sigma(\bm{q}(\gamma))}{d\gamma} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}
\label{g}
,
\end{align}
where the differentiation on the RHS of Eq. \eqref{g} is defined as follows:
\begin{align}
&\left. \frac{d\sigma(\bm{q}(\gamma))}{d\gamma} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}
\nonumber \\
&\quad\quad\quad:=\lim_{\Delta \gamma\to 0}
\frac{\sigma(\bm{q}^{\textrm{FB}}(\gamma+\Delta\gamma)) - \sigma(\bm{q}^{\textrm{FB}}(\gamma))}{\Delta\gamma}
\label{Gnm}
.
\end{align}
In the numerical calculation, we use a non-zero but sufficiently small $\Delta \gamma$ for the evaluation of $g$.
Then, the averaged rigidity $G$ is defined as
\begin{align}
G&:=\Braket{ g }
\label{Gdef}
,
\end{align}
where $\langle \cdot \rangle$ is the ensemble average.

\section{Theoretical Analysis \label{Sec3}}
In this section, we introduce Hessian matrix in Sec.~\ref{Sec3A} and theoretical expressions of rigidity in Sec.~\ref{Sec3B}.

\subsection{Hessian matrix for frictional grains}\label{Sec3A}

Because the Hessian matrix is equivalent to the Jacobian matrix for harmonic grains, as shown in  Appendices~\ref{Hessian=Jacobian} and \ref{AppHkn2},
in this study, we adopt the Hessian matrix ($\mathcal{H}$), where the element is given by~\cite{Saitoh19}:
\begin{align}
\mathcal{H}_{ij}^{\alpha\beta}:=\left. \frac{\partial^{2} \delta e_{ij} (\bm{q}(\gamma))}{\partial q_{i}^{\alpha}\partial q_{j}^{\beta}} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)},
\label{Jacobian}
\end{align}
where $\alpha$ and $\beta$ are any of $x,y$ and $\ell$, while $i$ and $j$ express the grain indices.
Here, we have introduced the effective potential energy $\delta e_{ij}$ between the contacted $i$-th and $j$-th grains as:
\begin{align}
\delta e_{ij}:=\frac{k_{N}}{2}(\delta\bm{r}_{ij}\cdot\bm{n}_{ij})^{2} + \frac{k_{T}}{2}\delta\bm{r}_{ij,\perp}^{2} ,
\label{potential}
\end{align}
where $\delta\bm{r}_{ij,\perp}$ is defined as
\begin{align}
\delta\bm{r}_{ij,\perp}:=\delta\bm{r}_{ij}-(\delta\bm{r}_{ij}\cdot\bm{n}_{ij})\bm{n}_{ij}-\delta\bm{\ell}_{ij}\times\bm{n}_{ij}
\label{tangential}
\end{align}
with
\begin{align}
\delta\bm{r}_{ij}:&=\delta\bm{r}_{i}-\delta\bm{r}_{j},\\
\delta\bm{\ell}_{ij}:&=(\delta\ell_{i}+\delta\ell_{j})\bm{e}_{z}
\end{align}
under the virtual displacements $\delta\bm{r}_i$ and $\delta\ell_i$ from the FB state at $\bm{r}_i^{\rm FB}$ and $\ell_i^{\rm FB}$, respectively.

The Hessian matrix introduced in Eq.~\eqref{Jacobian} can be written as
\begin{align}
\mathcal{H}=
\left[
\begin{matrix}
\mathcal{H}_{11} & \cdots & \mathcal{H}_{1i} & \cdots & \mathcal{H}_{1j} & \cdots & \mathcal{H}_{1N}\\
\vdots & \ddots & \vdots &   & \vdots &  &\vdots \\
\mathcal{H}_{i1} & \cdots & \mathcal{H}_{ii} & \cdots & \mathcal{H}_{ij} & \cdots & \mathcal{H}_{iN}\\
\vdots &  & \vdots & \ddots & \vdots &  &\vdots \\
\mathcal{H}_{j1} & \cdots & \mathcal{H}_{ji} & \cdots & \mathcal{H}_{jj} & \cdots & \mathcal{H}_{jN}\\
\vdots &  & \vdots &   & \vdots & \ddots &\vdots \\
\mathcal{H}_{N1} & \cdots & \mathcal{H}_{Ni} & \cdots & \mathcal{H}_{Nj} & \cdots & \mathcal{H}_{NN}
\end{matrix}
\right]
,
\label{defJ}
\end{align}
where $\mathcal{H}_{ij}$ is a $3\times3$ submatrix of the Hessian $\mathcal{H}$ for a pair of grains $i$ and $j$ satisfying:
\begin{align}
\mathcal{H}_{ij}=
\left[
\begin{matrix}
\mathcal{H}_{ij}^{xx} & \mathcal{H}_{ij}^{xy} & \mathcal{H}_{ij}^{x\ell} \\
\mathcal{H}_{ij}^{yx} & \mathcal{H}_{ij}^{yy} & \mathcal{H}_{ij}^{y\ell} \\
\mathcal{H}_{ij}^{\ell x} & \mathcal{H}_{ij}^{\ell y} & \mathcal{H}_{ij}^{\ell\ell}
\end{matrix}
\right]
\label{33J}
.
\end{align}
See  Appendix \ref{Hessian=Jacobian} for an explicit expression of each component of the Hessian matrix.
Note that $\mathcal{H}_{ij}^{\alpha\beta}=0$ if the $i$-th and $j$-th grains are not in contact with each other.

Because $\mathcal{H}$ is a real symmetric matrix, its eigenvalues and eigenvectors are also real.
Using the decomposition of the potential, the Hessian matrix can be divided into
\begin{align}
\mathcal{H}_{ij}^{\alpha\beta}=\mathcal{H}_{N,ij}^{\alpha\beta}+\mathcal{H}_{T,ij}^{\alpha\beta}
\label{sepNT}
,
\end{align}
where
\begin{align}
\mathcal{H}_{N,ij}^{\alpha\beta}:&=\left. \frac{\partial^{2} \delta e_{N,ij}^{\alpha}(\bm{q}(\gamma))}{\partial q_{i}^{\alpha}\partial q_{j}^{\beta}} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)},\\
\mathcal{H}_{T,ij}^{\alpha\beta}:&=\left. \frac{\partial^{2} \delta e_{T,ij}^{\alpha}(\bm{q}(\gamma))}{\partial q_{i}^{\alpha}\partial q_{j}^{\beta}} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}
\end{align}
for $i\neq j$ and
\begin{align}
\mathcal{H}_{N,ij}^{\alpha\beta}:&
\left. \frac{\partial^{2} \delta e_{N,ik}^{\alpha}(\bm{q}(\gamma))}{\partial q_{i}^{\alpha}\partial q_{i}^{\beta}} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)},\\
\mathcal{H}_{T,ij}^{\alpha\beta}:&
\left. \frac{\partial^{2} \delta e_{T,ik}^{\alpha}(\bm{q}(\gamma))}{\partial q_{i}^{\alpha}\partial q_{i}^{\beta}} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}
\end{align}
for $i=j$.
Here, we have introduced
\begin{align}
\delta e_{N,ij}:&=\frac{k_{N}}{2}(\delta\bm{r}_{ij}\cdot\bm{n}_{ij})^{2},    \\
\delta e_{T,ij}:&=\frac{k_{T}}{2}\delta\bm{r}_{ij,\perp}^{2} .
\end{align}
To determine the explicit expression of each component of the Hessian matrix, refer to Appendix~\ref{Hessian=Jacobian}.

The eigenequation of the Hessian matrix $\mathcal{H}$ is given by
\begin{align}
\mathcal{H}\ket{\Phi_{n}}&= \lambda_{n}\ket{\Phi_{n}}, \label{REieq}
\end{align}
where $\ket{\Phi_{n}}$ is the right eigenvector corresponding to the $n$-th eigenvalue $\lambda_{n}$ of $\mathcal{H}$.
Because the Hessian matrix is a real symmetric matrix, its left eigenequation is equivalent to its right eigenequation.
Such properties remain unchanged even under the Lees-Edwards boundary conditions (see Appendix \ref{app:eff_boundary}).
If all eigenstates are non-degenerate, $\ket{\Phi_n}$ satisfies the orthonormal relation $\braket{\Phi_m|\Phi_n}=\delta_{mn}$ with normalization $\braket{\Phi_n|\Phi_n}=1$,
where $\braket{\Phi_{n}|\Phi_{n}}:=\sum_{i=1}^{N}\sum_{\alpha=x,y,\ell}(\Phi_{n,i}^{\alpha})^{2}$.

\subsection{Expressions of the rigidity via eigenmodes \label{Sec3B}}

In this subsection, we consider the rigidity $g$ introduced in Eq.~\eqref{Gdef}.
See Appendix~\ref{app:rigidity} for the detailed properties of $g$.

Let us introduce $\tilde{\bm{F}}_i := (\tilde{F}^x_i , \tilde{F}^y_i , \tilde{F}^\ell_i)^{\rm T} := (F^x_i, F^y_i,2T_i/d_i)^{\rm T}$ and $\ket{\tilde F(\bm{q}(\gamma))}$ as
\begin{align}
\ket{\tilde F(\bm{q}(\gamma))}:=[ \tilde{\bm{F}}_{1}^{\textrm{T}}(\bm{q}(\gamma)), \tilde{\bm{F}}_{2}^{\textrm{T}}(\bm{q}(\gamma)), \cdots, \tilde{\bm{F}}_{N}^{\textrm{T}}(\bm{q}(\gamma)) ]^{\textrm{T}}.
\end{align}
Because the FB state is the minimum state of the potential energy, as shown in Appendix~\ref{Hessian=Jacobian},
$\ket{\tilde{F}(\bm{q}(\gamma))}|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}$ satisfies
\begin{align}
\left. \ket{\tilde F(\bm{q}(\gamma))} \right|_{ \bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma) }=\left. \frac{d\ket{\tilde F(\bm{q}(\gamma))}}{d\gamma} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)}=\ket{0},
\end{align}
where $\ket{0}$ is the ket vector containing $0$ for all components.

Introducing
\begin{align}
\Ket{\frac{d\mathring q}{d\gamma}}
:=
\left[
\frac{d\mathring r_{1}^{x}}{d\gamma},
\frac{d\mathring r_{1}^{y}}{d\gamma},
\frac{d\mathring \ell_{1}}{d\gamma},
\cdots,
\frac{d\mathring r_{N}^{x}}{d\gamma},
\frac{d\mathring r_{N}^{y}}{d\gamma},
\frac{d\mathring \ell_{N}}{d\gamma}
\right]^{\textrm{T}}
,
\label{NonAff}
\end{align}
one can write
 \begin{align}
\left. \frac{d\ket{\tilde F(\bm{q}(\gamma))}}{d\gamma} \right|_{\bm{q}(\gamma)=\bm{q}^{\textrm{FB}}(\gamma)} =-\ket{\Xi} +\tilde\mathcal{H}\Ket{\frac{d\mathring q}{d\gamma}}
\label{dFdG2}
,
\end{align}
where we have used Eqs. \eqref{drdG} and \eqref{dldG}.
The first and second terms on the RHS of Eq. \eqref{dFdG2} represent
the strain derivatives of the forces for the contributions from the affine and nonaffine displacements, respectively.
In Eq.~\eqref{dFdG2} we have introduced $\ket{\Xi}$, which is defined as:
\begin{align}
\ket{\Xi}:=\sum_{j}
\left[
\begin{matrix}
\mathcal{H}_{N,j1}^{xx}r_{1j}^{y}\\
\mathcal{H}_{N,j1}^{xy}r_{1j}^{y}\\
\mathcal{H}_{N,j1}^{x\ell}r_{1j}^{y}\\
\vdots\\
\mathcal{H}_{N,jN}^{xx}r_{Nj}^{y}\\
\mathcal{H}_{N,jN}^{xy}r_{Nj}^{y}\\
\mathcal{H}_{N,jN}^{x\ell}r_{Nj}^{y}
\end{matrix}
\right]
.
\end{align}
We have used $\tilde\mathcal{H}$ in Eq. \eqref{dFdG2}, which is defined as
\begin{align}
\tilde\mathcal{H}_{ii}^{\alpha\beta}
&:=
\left\{
\begin{matrix}
-\mathcal{H}_{ii}^{\ell x} & (\alpha=\ell,\ \beta=x) \\
-\mathcal{H}_{ii}^{\ell y} & (\alpha=\ell,\ \beta=y) \\
\mathcal{H}_{ii}^{\alpha\beta} & (\textrm{otherwise})
\end{matrix}
\right.
\label{tildeJ1}
\end{align}
and
\begin{align}
\tilde\mathcal{H}_{ij}^{\alpha\beta}
&:=
\left\{
\begin{matrix}
-\mathcal{H}_{ij}^{x\ell} & (\alpha=x,\ \beta=\ell) \\
-\mathcal{H}_{ij}^{y\ell} & (\alpha=y,\ \beta=\ell) \\
\mathcal{H}_{ij}^{\alpha\beta} & (\textrm{otherwise})
\end{matrix}
\right.
\label{tildeJ2}
\end{align}
for $i\neq j$.
Note that $\mathcal{H}_{T,ij}^{\alpha\beta}$ or $\tilde\mathcal{H}_{T,ij}^{\alpha\beta}$ does not affect $\ket{\Xi}$,
because the affine displacements are instantaneously applied to the system as a step strain.
Thus, the integral interval of the tangential displacement during the affine deformation is zero.

Expanding the nonaffine displacements by the eigenvectors of $\tilde\mathcal{H}$ and using the fact that the left-hand side of Eq. \eqref{dFdG2} is zero, we obtain
\begin{align}
\Ket{\frac{d\mathring q}{d\gamma}}=\sideset{}{^{'}}\sum_{n} \frac{\braket{\tilde{\Phi}_{n}|\Xi}}{\tilde{\lambda}_{n}}\ket{\tilde{\Phi}_{n}},
\label{eqEiExp}
\end{align}
where $\tilde{\lambda}_{n}$ and $\ket{\tilde{\Phi}_{n}}$ are the $n$-th eigenvalues of $\tilde\mathcal{H}$, and the eigenvector corresponding to $\tilde{\lambda}_n$, respectively.
Here, $\sum_{n}\nolimits'$ on the RHS of Eq. \eqref{eqEiExp} excludes low-frequency modes for $\tilde\lambda_{n}t_{0}^{2}/m\leq10^{-12}$ to maintain the numerical accuracy.
Note that $\ket{\tilde{\Phi}_n}$ satisfies the orthonormal relation $\braket{\tilde{\Phi}_m|\tilde{\Phi}_n}=\delta_{mn}$, if all eigenstates are non-degenerate.
The expression for $d\mathring q/d\gamma$ in Eq.~\eqref{eqEiExp} leads to a discontinuous change in $d\mathring q/d\gamma$ at a critical strain $\gamma_c$ for a plastic event
because the eigenvectors and eigenvalues are discontinuously changed at this point.

The rigidity is decomposed into two parts:
\begin{align}
g:=g_{\textrm{A}}+g_{\textrm{NA}} ,
\label{expG}
\end{align}
where $g_{\textrm{A}}$ and $g_{\textrm{NA}}$ are the rigidities corresponding to the affine and nonaffine displacements, respectively, for one sample.
With the aid of Eqs. \eqref{Sigma}, \eqref{Gdef}, and \eqref{tildeJ2}, the expressions for $g_{\textrm{A}}$ and $g_{\textrm{NA}}$ can be obtained as:
\begin{align}
g_{\textrm{A}}
:&=
\frac{1}{4 L^{2}}
\sum_{i,j(i\neq j)}
r_{ij}^{y}
\left[
r_{ij}^{y}\mathcal{H}_{N,ji}^{xx}
+
r_{ij}^{x}\mathcal{H}_{N,ji}^{yx}
\right]
\label{GAff}
, \\
g_{\textrm{NA}}
:&=
\frac{1}{4 L^{2}}
\sum_{i,j(i\neq j)}
\Biggl[
\sum_{\zeta=x,y}
 \left( r_{ij}^{y} \tilde\mathcal{H}_{ij}^{x\zeta}  + r_{ij}^{x} \tilde\mathcal{H}_{ij}^{y\zeta}  \right) \frac{d\mathring r_{ij}^{\zeta}}{d\gamma} \nonumber \\
&\quad\quad\quad\quad\quad\quad\quad\quad\quad
-\left( r_{ij}^{y} \tilde\mathcal{H}_{ij}^{x\ell}  + r_{ij}^{x} \tilde\mathcal{H}_{ij}^{y\ell} \right) \frac{d\mathring \ell_{ij}}{d\gamma}
\Biggl]
,
\label{GNAff}
\end{align}
where we have introduced
\begin{align}
\frac{d\mathring r_{ij}^{\zeta}}{d\gamma}:&=\frac{d\mathring r_{i}^{\zeta}}{d\gamma} - \frac{d\mathring r_{j}^{\zeta}}{d\gamma}, \\
\frac{d\mathring \ell_{ij}}{d\gamma}:&=\frac{d\mathring \ell_{i}}{d\gamma} + \frac{d\mathring \ell_{j}}{d\gamma}.
\end{align}

Substituting Eq. \eqref{eqEiExp} into Eq. \eqref{GNAff}, $g_{\textrm{NA}}$ can be rewritten as
\begin{align}
g_{\textrm{NA}}
&=
-
\frac{1}{L^{2}}
\sideset{}{^{'}}\sum_{n}\frac{\braket{\tilde{\Phi}_{n}|\Xi}\braket{\Theta|\tilde{\Phi}_{n}}}{\tilde\lambda_{n}}
\label{GNAffv2}
,
\end{align}
where we have introduced
\begin{widetext}
\begin{align}
\bra{\Theta}:=\frac{1}{2}\sum_{j}
\left[
\left( r_{1j}^{y}\tilde{\mathcal{H}}_{j1}^{xx}+r_{1j}^{x}\tilde{\mathcal{H}}_{j1}^{yx} \right),
\left( r_{1j}^{y}\tilde{\mathcal{H}}_{j1}^{xy}+r_{1j}^{x}\tilde{\mathcal{H}}_{j1}^{yy} \right),
\hdots,
\left( r_{Nj}^{y}\tilde{\mathcal{H}}_{jN}^{xy}+r_{Nj}^{x}\tilde{\mathcal{H}}_{jN}^{yy} \right),
\left( r_{Nj}^{y}\tilde{\mathcal{H}}_{jN}^{x\ell}+r_{Nj}^{x}\tilde{\mathcal{H}}_{jN}^{y\ell} \right)
\right] .
\end{align}
The affine rigidity can be also expressed as
\begin{align}
g_{\textrm{A}}
= \frac{1}{L^2} \braket{Y|\Xi} ,
\end{align}
where
\begin{align}
\bra{Y}:=\left[ y_{1}, 0, 0, y_{2},0, 0,  \cdots, y_{N}, 0, 0 \right]
.
\end{align}

\end{widetext}

To verify the validity of the theoretical treatment,
we introduce the theoretical stress $\sigma^{\textrm{th}}(\gamma)$ with the aid of Eq.\eqref{expG} as
\begin{align}\label{theoretical_sigma}
\sigma^{\textrm{th}}(\gamma+\Delta\gamma):=\sigma(\bm{q}^{\textrm{FB}}(\gamma))+g(\gamma)\Delta\gamma .
\end{align}

\section{Results and discussion\label{results}}

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

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

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{Stress.eps}
\caption{
A stress-strain curve for $0\leq\gamma\leq0.5$ for one sample of the collection of grains ($N=128$), which includes the theoretical results (line) and simulation results (filled symbols) under the condition$\Delta\gamma_{\textrm{Th}}=\Delta\gamma_{\textrm{in}}=10^{-4}$.
The inset is a close-up of the stress-strain curve in the vicinity of a stress-drop event.
}
\label{Stress}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{StressCompProtocol.eps}
\caption{
Stress-strain relations for
$\Delta\gamma_{\textrm{Th}}=1.0\times10^{-4}$ (blue circles) and $\Delta\gamma_{\textrm{Th}}=1.0\times10^{-8}$ (red triangles) with $N=128$.}
\label{StressCmpGammaTh}
\end{figure}

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

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

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

\begin{figure}[htbp]
\centering
\includegraphics[width=7cm]{EiDrop.eps}
\includegraphics[width=7cm]{GStressDrop.eps}
\caption{
Plots of (a) the smallest eigenvalue except for the zero modes and  (b) the rigidity $g$ based on the eigenvalue analysis (open symbols) and numerical shear stress (filled symbols) against $\gamma$ for $0.3408\leq\gamma\leq0.3413$ and $N=128$.
}
\label{Drop}
\end{figure}

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

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{EiFieldDrop-.eps}
\includegraphics[width=8cm]{EiFieldDrop+.eps}
\caption{
Plots of eigenvectors at (a) $\gamma_{c-}$  and at (b) $\gamma_{c+}$  corresponding to the smallest eigenvalue.
For visualization, the magnitudes of the vectors are three times larger than their true values ($N=1024$).
}
\label{EiFieldDrop}
\end{figure}

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

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

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{NAFieldDrop+Th.eps}
\includegraphics[width=8cm]{NAFieldDrop+Nm.eps}
\caption{
Plots of $\Ket{d\mathring q/d\gamma}$ at $\gamma=\gamma_{c+}$ for $N=1024$ and $\Delta\gamma_{\textrm{Th}}=1.0\times10^{-8}$,
where (a) and (b) are based on the eigenvalue analysis and simulation, respectively.
For visualization, we magnify  the magnitudes of the vectors with the factor $10.0$.
}
\label{NAFieldComp}
\end{figure}

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{NAFieldScaledLrgN.eps}
\caption{
Plot of $\Delta\mathring{q}|_{c}$ in the simulation for $N=1024$ and $\Delta\gamma_{\textrm{Th}}=1.0\times10^{-8}$.
For visualization, we magnify $\Delta\mathring q|_{c}$ with the factor $1.0\times10^{3}$.
}
\label{NAFieldDropLrgN}
\end{figure}

Because we cannot use the eigenvalue analysis at the critical strain $\gamma_c$ for a plastic event,
let us analyze the nonaffine displacement $\Delta \mathring{\bm{q}}$ between $\gamma_{c-}$ and $\gamma_{c+}$ caused by an avalanche using the simulation:
\begin{align}
\Delta\mathring{\bm{q}}|_{c}
:=
\left[
\begin{matrix}
\left. \Delta\mathring{\bm{q}}_{1}\right|_{c} \\
\left. \Delta\mathring{\bm{q}}_{2} \right|_{c} \\
\vdots \\
\left. \Delta\mathring{\bm{q}}_{N} \right|_{c}
\end{matrix}
\right]
,
\end{align}
where
\begin{align}
\left.\Delta\mathring{\bm{q}}_{i}\right|_{c}:=
\left[
\begin{matrix}
r_{i}^{\textrm{FB},x}(\gamma_{c+}) - r_{i}^{\textrm{FB},x}(\gamma_{c-}) - \Delta\gamma r_{i}^{\textrm{FB},y}(\gamma_{c-}) \\
r_{i}^{\textrm{FB},y}(\gamma_{c+}) - r_{i}^{\textrm{FB},y}(\gamma_{c-}) \\
\ell_{i}^{\textrm{FB}}(\gamma_{c+}) - \ell_{i}^{\textrm{FB}}(\gamma_{c-})
\end{matrix}
\right]
.
\end{align}

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

\begin{figure}[htbp]
\centering
\includegraphics[width=7cm]{GContact.eps}
\includegraphics[width=7cm]{EiMinContact.eps}
\caption{
Plots of the shear modulus with (a) the theoretical evaluation (open symbols) and the simulation results (filled symbols) except for critical strain with the close-up of $g$ near a yielding point (inset),
and (b) the smallest eigenvalue, except for zero modes
against $\gamma$ for $0\leq\gamma\leq0.002$ for $N=128$.
Note that the rigidity is not plotted at the yielding points, because it diverges there.
}
\label{Contact}
\end{figure}

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

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{StressContact.eps}
\caption{
The plot of the stress-strain curve in the region $0\leq\gamma\leq0.002$ for $N=128$.
}
\label{StressContact}
\end{figure}

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

\begin{figure}[htbp]
\centering
\includegraphics[width=8cm]{Gave.eps}
\caption{
Plots of the rigidity $G$ based on the eigenvalue analysis (line) and on the simulation (filled symbols) with $\Delta\gamma_{\textrm{Th}}=1.0\times10^{-4}$.
We have used $30$ samples for $N=128$, where the error bars represent the standard deviations for $\gamma$.
Note that we have omitted the data if stress-drop events occur.
}
\label{GContactEnsemble}
\end{figure}

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

\section{Conclusion}\label{concluding_remarks}

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

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

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