# Blow-up for the 1D compressible Euler equations with separated time and
space dependent damping coefficient

Yuusuke<EMAIL_ADDRESS>The University of Shiga
Prefecture

###### Abstract

In this paper, we consider the 1D Euler equation with time and space dependent
damping term $-a(t,x)v$. We suppose that the coefficient $a$ satisfies the
following timewise and spatial separated decaying condition

$|a(t,x)|\leq a_{1}(t)+a_{2}(x),$

where $a_{1}$ and $a_{2}$ are integrable functions with respect to $t$ and
$x$. Under this condition, we show solutions blow up with small perturbed
initial data. The upper estimate of the life-span is also given.

##  1 Introduction

###  1.1 Problem and main theorem

In this paper, we consider the following Cauchy problem of the compressible
Euler equation with time and space dependent damping

$\displaystyle\left\\{\begin{array}[]{ll}u_{t}-v_{x}=0,\\\
v_{t}+p(u)_{x}=-a(t,x)v,\\\
(u(0,x),v(0,x))=(1+\varepsilon\varphi(x),\varepsilon\psi(x)).\end{array}\right.$
(1.4)

Here $x\in\mathbb{R}$ is the Lagrangian spatial variable and
$t\in\mathbb{R}_{+}$ is time. $u=u(t,x)$ and $v=v(t,x)$ are the real valued
unknown functions, which stand for the specific volume and the fluid velocity.
$\varepsilon$ is a positive and small parameter. We impose that the damping
coefficient is integrable $t$ and $x$ separably. Namely, for the damping
coefficient $a(t,x)\in C^{1}_{b}([0,\infty)\times\mathbb{R})$, we suppose that
there exist two non-negative functions $a_{1}\in C^{1}_{b}([0,\infty))$ and
$a_{2}\in C^{1}_{b}(\mathbb{R})$ such that

$\displaystyle|a(t,x)|+|a_{x}(t,x)|+|a_{t}(t,x)|\leq a_{1}(t)+a_{2}(x),$ (1.5)
$\displaystyle\int_{0}^{\infty}a_{1}(t)dt+\int_{-\infty}^{\infty}a_{2}(x)dx=C_{a}<\infty,$
(1.6) $\displaystyle x\dfrac{d}{dx}a_{2}(x)\leq 0.$ (1.7)

The typical example of $a(t,x)$ is

$\displaystyle a(t,x)=(1+t)^{-\lambda_{1}}+(1+|x|)^{-\lambda_{2}},$

where $\lambda_{1},\lambda_{2}>1$. From the Young inequality, the above
condition also covers the following example

$\displaystyle a(t,x)=(1+t)^{-\lambda_{1}}(1+|x|)^{-\lambda_{2}},$

where $\lambda_{1}$ and $\lambda_{2}$ are non-negatives such that
$\lambda_{1}+\lambda_{2}>1$. In the case with $a\equiv 0$, the equations in
(1.4) are the compressible Euler, which is a fundamental model for the
compressible inviscid fluids. In the case with $a\not\equiv 0$, this system
describes the flow of fluids in porous media. We assume that the flow is
barotropic ideal gases. Namely the pressure $p$ satisfies that

$\displaystyle p(u)=\frac{u^{-\gamma}}{\gamma}\ \ \mbox{for}\ \gamma>1.$ (1.8)

For initial data, in order to avoid the singularity of $p^{\prime}$, we assume
that there exists constant $\delta_{0}>0$ such that for all $x\in\mathbb{R}$

$\displaystyle u(0,x)\geq\delta_{0}.$ (1.9)

The local existence and the uniqueness of solutions of (1.4) hold with
$C^{1}_{b}$ initial data (e.g. Friedrichs [6], Lax [11] and Li and Wu [14]).
If $u(0,\cdot),v(0,\cdot)\in C^{1}_{b}(\mathbb{R})$ and $a\in
C^{1}_{b}([0,\infty)\times\mathbb{R})$ and (1.9) is assumed, then (1.4) has a
local and unique solution $(u,v)\in C^{1}_{b}([0,T]\times\mathbb{R})$ such
that $u(t,x)>0$.

To state the main theorem, we define the life-span $T^{*}$ by

$\displaystyle T^{*}=$ $\displaystyle\sup\\{T>0\ |\
\sup_{t\in[0,T)}\\{\|(u_{t},v_{t})(t)\|_{L^{\infty}}$ (1.10)
$\displaystyle+\|(u_{x},v_{x})(t)\|_{L^{\infty}}+\|p^{\prime}(u(t))\|_{L^{\infty}}\\}<\infty\\}.$

For $t\geq 0$, we set

$\displaystyle\Omega(t)=\\{x\in\mathbb{R}\ |\ x\geq x_{+}(t;0,0)\ \mbox{or}\
x\leq x_{-}(t;0,0)\\},$

and also set for $\pm x_{0}>0$ respectively

$\displaystyle\Omega_{+,x_{0}}(t)=\\{x\in\mathbb{R}\ |\ x\geq
x_{+}(t;0,x_{0})\\},$ $\displaystyle\Omega_{-,x_{0}}(t)=\\{x\in\mathbb{R}\ |\
x\leq x_{-}(t;0,x_{0})\\},$

where $x_{\pm}(t;t_{0},x_{0})$ is the plus and minus characteristics through
$(t_{0},x_{0})$ (see (2.4) for the definition). Now we state the main theorem
of this paper.

######  Theorem 1.1.

Let $\gamma>1$ and $(\varphi,\psi)\in C^{1}_{b}(\mathbb{R})$. Suppose that

$\displaystyle\psi_{x}(x_{0})\leq-K$ (1.11)

for some point $x_{0}\in\mathbb{R}$ and a constant
$K=K(\|\varphi\|_{C^{1}(\mathbb{R})},\|\psi\|_{L^{\infty}(\mathbb{R})})>0$.
Then there exists a number $\varepsilon_{0}>0$ such that if
$\varepsilon\leq\varepsilon_{0}$ then

$\displaystyle T^{*}\leq C\varepsilon^{-1}.$

Furthermore, it holds that if $x_{0}=0$

$\displaystyle\limsup_{t\rightarrow
T^{*}}\|(u_{t},v_{t})(t)\|_{L^{\infty}(\Omega(t))}+\|(u_{x},v_{x})(t)\|_{L^{\infty}(\Omega(t))}=\infty.$

If $\pm x_{0}>0$, it holds that respectively

$\displaystyle\limsup_{t\rightarrow
T^{*}}\|(u_{t},v_{t})(t)\|_{L^{\infty}(\Omega_{\pm,x_{0}}(t))}+\|(u_{x},v_{x})(t)\|_{L^{\infty}(\Omega_{\pm,x_{0}}(t))}=\infty.$

######  Remark 1.2.

Assumption on the smallness of initial data One can show the same statement as
in Theorem 1.1 without the smallness of derivatives of initial data. To holds
theorem 1.1, it is enough to assume $\|r(0,\cdot)\|_{L^{\infty}}$ and
$\|s(0,\cdot)\|_{L^{\infty}}$ are small, where $r$ and $s$ are the Riemann
invariants defined in (2.3). In stead of (1.11), we assume that with
$\varepsilon=1$

$\displaystyle r_{x}(0,0)>-K^{*}\ \mbox{or}\ s_{x}(0,0)>-K^{*},$

where $K^{*}$ is a positive constant depending on
$\|r(0,\cdot)\|_{L^{\infty}}$, $\|s(0,\cdot)\|_{L^{\infty}}$, $\gamma$ and
$C_{a}$ in (1.5). Then the solution blows up in finite time, and the life-span
is bounded by

$\displaystyle\frac{C}{r_{x}(0,0)-K^{*}}\ \mbox{or}\
\frac{C}{s_{x}(0,0)-K^{*}}.$

This can be shown in same way as in the proof of Theorem 1.1, since the key
estimate in Lemma 2.4 is valid under the assumption of the smallness of
$\|r(0,\cdot)\|_{L^{\infty}}$ and $\|s(0,\cdot)\|_{L^{\infty}}$ only (not
derivative of $r$ and $s$).

######  Remark 1.3.

Possibility of the blow-up of $p^{\prime}$ Here we only consider the case that
$x_{0}=0$ in Theorem 1.1. From the proof of Theorem 1.1, we can see that the
blow-up of $p^{\prime}$ does not occurs in the region $\cup_{t\geq
0}\\{t\\}\times\Omega(t)$. Theorem 1.1 states that the derivative blow-up
occurs in this region. While, in the region $\left(\cup_{t\geq
0}\\{t\\}\times\Omega(t)\right)^{c}$, we do not obtain the boundedness of
$p^{\prime}$, which means that the the blow-up of $p^{\prime}$ may occurs
before the derivative blow-up, where we denote the complement set of $X$ in
$[0,T^{*})\times\mathbb{R}$ by $X^{c}$. From the same reason, we do not have a
lower estimate (optimality) of the estimate of the life-span. However, the
upper estimate of the life-span is optimal, if $a_{1}\equiv 0$ or $a_{2}\equiv
0.$

###  1.2 Known results

If $a\equiv 0$, it is well-known that solutions can blow up in finite time
with even small initial data (e.g. Lax [12] and Zabusky [28]). It is not hard
to show that $C^{\prime}\varepsilon^{-1}\leq T^{*}\leq C\varepsilon^{-1}$ in
this case. Here the lower estimate ”$T^{*}\geq$” means the minimal existence
time of solutions for arbitrary function $\varphi$ and $\psi$. On the other
hand hand, the upper estimate ”$T^{*}\leq$” means maximal existence time for
some $\varphi$ and $\psi$ such that a condition as (1.11) is satisfied. In the
case that $a(t,x)\equiv 1$, Hsiao and Liu in [8] have proved that classical
solutions exist globally in time, if initial data are small perturbations near
constant states and that small solutions asymptotically behave like that to
the following porous media system as $t\rightarrow\infty$:

$\displaystyle\left\\{\begin{array}[]{ll}\bar{u}_{t}=-p(\bar{u})_{xx},\\\
\bar{v}=-p(\bar{u})_{x}.\end{array}\right.$

After Hsiao and Liu’s work, many improvements and generalizations of this work
have been investigated (see Hsiao and Liu [9], Nishihara [18], Hsiao and Serre
[10], Marcati and Nishihara [16] and Mei [17]). Let us review previous results
for the following 1D Euler equation with time dependent damping case:

$\displaystyle\left\\{\begin{array}[]{ll}u_{t}-v_{x}=0,\\\
v_{t}+p(u)_{x}=-\dfrac{\mu v}{(1+t)^{\lambda}},\\\ \end{array}\right.$ (1.15)

where $\lambda\geq 0$ and $\mu\in\mathbb{R}$. In [20, 21, 22], Pan has
established critical phenomena for $\mu$ and $\lambda$. Namely, in the case
with $0\leq\lambda<1$ and $\mu>0$ or $\lambda=1$ and $\mu>2$, Pan [21] has
proved that solutions of (1.15) exist globally in time, while, in the case
with $\lambda>1$ or $\lambda=1$ and $0\leq\mu\leq 2$, he has proved that
solutions of (1.15) can blow up under some conditions on initial data in [20,
22]. In the previous author’s papers [23], the optimal life-span estimates are
given. To summarize, the following estimates are given:

$\displaystyle C^{\prime}\varepsilon^{-1}\leq T^{*}\leq C\varepsilon^{-1}\ \
\mbox{for}\ \lambda>1\ \mbox{and}\ \lambda\geq 0,$ $\displaystyle
C^{\prime}\varepsilon^{-\frac{2}{2-\lambda}}\leq T^{*}\leq
C\varepsilon^{-\frac{2}{2-\lambda}}\ \ \mbox{for}\ \ \lambda=1\ \mbox{and}\
0\leq\mu<2,$ $\displaystyle e^{\frac{C^{\prime}}{\varepsilon}}\leq T^{*}\leq
e^{\frac{C}{\varepsilon}}\ \ \mbox{for}\ \lambda=1\ \mbox{and}\ \mu=2,$
$\displaystyle T^{*}=\infty\ \ \mbox{for}\ \mu>0\ \mbox{and}\ 0\leq\lambda<1\
\mbox{or}\ \mu>2\ \mbox{and}\ \lambda=1$

for small $\varepsilon>0$. We remark that the lower estimates hold under the
assumption (1.11). In [23], the author also has proved that
$C^{\prime}\varepsilon^{-1}\leq T^{*}\leq C\varepsilon^{-1}$ holds when
$a(t,x)$ satisfies that

$\displaystyle|a(t,x)|+|a_{x}(t,x)|+|a_{t}(t,x)|\leq C(1+t)^{-\lambda}$ (1.16)

for $\lambda>1$. It is not difficult to extend the above results to the case
that more general time dependent damping $a(t)$. However, there have not been
many papers dealing with space-time dependent damping coefficients. In Chen et
al. [4], they have considered (1.4) and (1.15) in the large data regime. In
[4] they have also show that the global existence for small perturbed initial
data, when $a(t,x)$ fulfills that with $0<\lambda<1$

$\displaystyle a(t,x)\cong(1+t)^{-\lambda},$
$\displaystyle|a_{x}(t,x)|+|a_{t}(t,x)|\leq C(1+t)^{-2\lambda}.$

The method in Chen et al. [4] to prove uniform estimate for the Riemann
invariant is based on a new maximum principle. They impose the positivity of
the damping coefficient. They also suppose that $u_{0}(x),v_{0}(x)$ and
$a(t,x)$ converge as $|x|\rightarrow\infty$. For $x$ dependent case, Chen,
Young and Zhang’s result in [2] can be applicable (see Remark 2.1). By
applying their method to the case that $a_{1}\equiv 0$, we can estimate the
life-span as $C^{\prime}\varepsilon^{-1}\leq T^{*}\leq C\varepsilon^{-1}$. In
[2], they have considered the following non-isentropic compressible Euler
equations (see also Chen, Pan and Zhu [3], Zheng [29] and Pan and Zhou [19]):

$\displaystyle\left\\{\begin{array}[]{ll}u_{t}-v_{x}=0,\\\
v_{t}+p(S,u)_{x}=0,\\\ S_{t}=0,\end{array}\right.$ (1.20)

where $p(S,u)=e^{S(x)}u^{-\gamma}/\gamma$ with $1<\gamma<3$ and $S=S(x)$ is a
given function such that$\int_{\mathbb{R}}|S^{\prime}(x)|dx<\infty$. They have
established that the formation of singularities (blow-up) occurs in finite
time. In [2], the key for the proof is the uniform boundedness of Riemann
invariant by using non-trivial characteristic method. It seems hard to apply
the method in [2] to the case that $a_{1}\not\equiv 0$ in (1.5), since the use
of $x$ variable characteristics $t_{\pm}(x)$ plays an important role to show
the uniform boundedness of Riemann invariant. We also should mention that
there are many studies on asymptotic behavior results for (1.15) (see [18, 1,
5, 13, 17] and references therein).

######  Remark 1.4.

The linear wave equation with space independent damping The study of the Euler
equation with the time-dependent damping (1.15) is motivated by Wirth’s
papers. In [25], he studies the behavior of solutions to

$u_{tt}-\Delta u=-\dfrac{\lambda u_{t}}{(1+t)^{\mu}}.$

In [26, 27], he gives a threshold of $\mu$ separating solutions to this
equation asymptotically behave like that to the corresponding hear or wave
equation. In [25], the critical case of the threshold is studied.

######  Remark 1.5.

Large data regime If $a\equiv 0$, to show the occurrence of the blow-up,
restrictions on $\|r(0,x)\|_{L^{\infty}}$ and $\|s(0,x)\|_{L^{\infty}}$ are
not necessary. Namely the necessary and sufficient condition for the global
existence has been given by Chen, Pan and Zhu’s paper [3]. They has shown that
the global classical solution exists if and only if

$\displaystyle r_{x}(0,x)\geq 0\ \mbox{and}\ s_{x}(0,x)\geq 0.$

In Chen et al. [4], they give a new sufficient condition of global existence
for the 1D damped Euler equation. They show the solution globally exists with
monotonically decreasing or increasing and large initial data.

###  1.3 Idea of the proof

In the previous results ([2, 23, 4]), the key for the proof is uniform
estimates of the Riemann invariant. More precisely, for time-decaying damping
case ($a_{2}\equiv 0$), the proof is based on the method of characteristic
with $t$ valuable. On the other hand, Characteristic with $x$ valuable is used
for spacial decaying damping case ($a_{1}\equiv 0$) to control the Riemann
invariant Thus, it seems difficult to balance the two method. The idea of the
proof is to estimate the Riemann invariant in the region outside a
characteristic curve. In the proof, we use the following two important
properties of this domain. The first is the property of trapping backward
characteristic curves by the region. Second, it is possible to replace the
spatially decaying damping coefficient with a time-decaying one in the region.

### Notation

For $\Omega\subset\mathbb{R}^{n}$, $C^{1}_{b}(\Omega)$ are the set of bounded
and continuous functions whose first partial derivatives are also bounded on
$\Omega$. The norm of $C^{1}_{b}(\Omega)$ is
$\|f\|_{C^{1}_{b}(\Omega)}=\|f\|_{L^{\infty}(\Omega)}+\|(\partial_{x_{1}}f,\ldots,\partial_{x_{n}}f)\|_{L^{\infty}(\Omega)}$.
When $\Omega=\mathbb{R}$, for abbreviation, we denote
$\|\cdot\|_{L^{\infty}(\Omega)}$ and $\|\cdot\|_{C^{1}_{b}(\Omega)}$ by
$\|\cdot\|_{L^{\infty}}$ and $\|\cdot\|_{C^{1}_{b}}$ respectively.

##  2 Preliminaries

###  2.1 Formulas for Riemann invariants

Now we prepare some identities for Riemann invariants. We set
$c=\sqrt{-p^{\prime}(u)}$ and
$\eta=\int_{u}^{\infty}c(\xi)d\xi=\frac{2}{\gamma-1}u^{-(\gamma-1)/2}$ and
define Riemann invariants as follows:

$\displaystyle\begin{array}[]{ll}r=v-\eta+\frac{2}{\gamma-1},\\\
s=v+\eta-\frac{2}{\gamma-1}.\end{array}$ (2.3)

This type of the Riemann invariant is used in [24]. For
$c=\sqrt{-p^{\prime}(u)}$, the plus and minus characteristic curves are
solutions to the following deferential equations:

$\displaystyle\frac{dx_{\pm}}{dt}(t)=\pm c(t,u(t,x_{\pm}(t))).$ (2.4)

If necessary, we denote the characteristic curves through $(t,x)$ by
$x_{\pm}(\cdot;t,x)$. Riemann invariants $r$ and $s$ are solutions to

$\displaystyle\left\\{\begin{array}[]{ll}\partial_{-}r=-\dfrac{a(t,x)}{2}(r+s),\\\
\partial_{+}s=-\dfrac{a(t,x)}{2}(r+s),\end{array}\right.$ (2.7)

where $\partial_{\pm}=\partial_{t}\pm c\partial_{x}$.We set
$A_{\pm}(\tau,t,x)=\exp\left(\int_{0}^{\tau}a(\bar{\tau},x_{\pm}(\bar{\tau};t,x))/2d\bar{\tau}\right)$
and abbreviate $A_{\pm}(t,t,x)=A_{\pm}(t,x)$ On the characteristic curves
through $(t,x)$, we have

$\displaystyle
A_{-}(t,x)r(t,x)=r(0,x_{-}(0))-\int_{0}^{t}\dfrac{A_{-}(\tau,t,x)a(t,x_{-}(\tau))s(\tau,x_{-}(\tau))}{2}d\tau,$
(2.8) $\displaystyle
A_{+}(t,x)s(t,y)=s(0,x_{+}(0))-\int_{0}^{t}\dfrac{A_{+}(\tau,t,x)a(t,x_{+}(\tau))r(\tau,x_{+}(\tau))}{2}d\tau.$
(2.9)

Differentiating the equations in (2.8) and (2.9) with $x$, from the identity
$s_{x}-r_{x}=2\eta_{x}=-2cu_{x}$, we have

$\displaystyle\left\\{\begin{array}[]{ll}\partial_{-}r_{x}=\dfrac{c^{\prime}}{2c}r_{x}(s_{x}-r_{x})-\dfrac{a(t,x)}{2}(r_{x}+s_{x})-\dfrac{a_{x}(t,x)}{2}(r+s),\\\
\partial_{+}s_{x}=\dfrac{c^{\prime}}{2c}r_{x}(r_{x}-s_{x})-\dfrac{a(t,x)}{2}(r_{x}+s_{x})-\dfrac{a_{x}(t,x)}{2}(r+s).\end{array}\right.$
(2.12)

Now we rewrite (2.12) as integral equalities. We define $\theta_{\gamma}(u)$
as follows:

$\displaystyle\theta_{\gamma}(u)=\left\\{\begin{array}[]{ll}\frac{4}{3-\gamma}u^{\frac{3-\gamma}{4}}-\frac{4}{3-\gamma}\
\ \mbox{for}\ \gamma\not=3,\\\ \log u\ \ \mbox{for}\
\gamma=3.\end{array}\right.$

From (1.4) and the definitions of $c$ and $\eta$, $r_{x}$ and $s_{x}$ can be
written by

$\displaystyle\begin{array}[]{ll}r_{x}=u_{t}+c(u)u_{x},\\\
s_{x}=u_{t}-c(u)u_{x}.\end{array}$ (2.16)

To deal with (2.16) on the characteristic curves through $(0,x)$, we set

$\displaystyle Y(t)=A_{-}(t,0,x)\sqrt{c(u(t,x_{-}(t)))}r_{x}(t,x_{-}(t)),$
$\displaystyle Q(t)=A_{+}(t,0,x)\sqrt{c(u(t,x_{+}(t)))}s_{x}(t,x_{+}(t)).$

Since it holds that

$\displaystyle\sqrt{c}s_{x}(t,x_{-}(t))=\frac{d}{dt}\theta_{\gamma}(u(t,x_{-}(t)))$
(2.17)

from (2.8), we have

$\displaystyle Y(t)=$ $\displaystyle
Y(0)+\int_{0}^{t}\dfrac{d}{d\tau}(A_{-}(\tau,0,x)a(t,x_{-}(\tau)))\theta_{\gamma}(\tau,x_{-}(\tau))d\tau$
$\displaystyle-\frac{1}{2}(A_{-}(t,0,x)a(t,x_{-}(t))\theta_{\gamma}(t,x_{-}(t))-a(0,x_{-}(0))\theta_{\gamma}(0,x_{-}(0)))$
$\displaystyle-\int_{0}^{t}A_{-}(\tau,0,x)\dfrac{a_{x}(\tau,x_{-}(\tau))}{2}\sqrt{c}(r+s)d\tau$
$\displaystyle-\int_{0}^{t}A_{-}(\tau,0,x)^{-1}\dfrac{\gamma+1}{4}u^{\frac{\gamma-3}{4}}Y(\tau)^{2}d\tau.$
(2.18)

Similarly we have

$\displaystyle Q(t)=$ $\displaystyle
Q(0)+\int_{0}^{t}\dfrac{d}{d\tau}(A_{+}(\tau,0,x)a(t,x_{+}(\tau)))\theta_{\gamma}(\tau,x_{+}(\tau))d\tau$
$\displaystyle-\frac{1}{2}(A_{+}(t,x)a(t,x_{+}(t))\theta_{\gamma}(t,x_{+}(t))-a(0,x_{+}(0))\theta_{\gamma}(0,x_{+}(0)))$
$\displaystyle-\int_{0}^{t}A_{+}(\tau,0,x)\dfrac{a_{x}(\tau,x_{+}(\tau))}{2}\sqrt{c}(r+s)d\tau$
$\displaystyle-\int_{0}^{t}A_{+}(\tau,0,x)^{-1}\dfrac{\gamma+1}{4}u^{\frac{\gamma-3}{4}}Q(\tau)^{2}d\tau.$
(2.19)

######  Remark 2.1.

Non-isentropic Euler equation Using Riemann invariants, non-isentropic Euler
system (1.20) can be written by

$\displaystyle\left\\{\begin{array}[]{ll}\partial_{-}r=-\dfrac{m^{\prime}}{2m\gamma}(r-s),\\\
\partial_{+}s=\dfrac{m^{\prime}}{2m\gamma}(s-r),\end{array}\right.$ (2.22)

where $m=e^{S/c_{v}}$ and the constant $c_{v}$ is called specific heat. In
[2], to obtain a uniform estimate for $r$ and $s$, they use the characteristic
curve $t=t_{\pm}(x)$ such that

$\displaystyle\frac{d}{dx}t_{\pm}(x)=\frac{\pm 1}{c(u(t_{\pm}(x),x))}.$

We also remark that the integrability of $S^{\prime}$ implies that

$\displaystyle\int_{\mathbb{R}}\frac{|m^{\prime}(x)|}{m(x)}dx=\frac{1}{c_{v}}\int_{\mathbb{R}}|S^{\prime}(x)|dx.$

The method in [2] to control the $L^{\infty}$-norm of the Riemann invariant is
applicable to the case that $m^{\prime}/m$ is time-dependent such that

$\displaystyle\int_{\mathbb{R}}\left|\frac{m_{x}}{m}(t,x)\right|dx\leq C.$

###  2.2 Estimates for Riemann invariants

In this subsection, we always suppose that it holds that

$\displaystyle c_{1}/4\leq c(u(t,x))\leq 4c_{1},$ (2.23)

where we set $c_{1}=c(1)$. Henceforth, we often use the following fact for
backward characteristics $x_{\pm}(\tau;t,x)$ with $\tau\leq t$ (see Figure 1)

$\displaystyle x_{\pm}(\tau;t,x)\in\Omega_{+,x_{0}}(\tau),\ \ \mbox{if}\
x\in\Omega_{+,x_{0}}(t).$ (2.24)

This property also holds for $\Omega_{-,x_{0}}(t)$.

Figure 1 : $\Omega_{+,x_{0}}(t)$ and characteristic curves
$\Omega_{+,x_{0}}(t)$$(x,t)$$x_{+}(\cdot;t,x)$$x_{-}(\cdot;t,x)$$x$$x_{0}$$\Omega_{+,x_{0}}(\tau)$

######  Lemma 2.2.

For $C^{1}_{b}$ solution of (1.4) satisfying (2.23), the following uniform
estimate holds for $0\leq\tau\leq t$ and $x\in\Omega(t)$,
$x\in\Omega_{+,x_{0}}(t)$ or $x\in\Omega_{-,x_{0}}(t)$:

$\displaystyle C^{-1}\leq|A_{\pm}(\tau,t,x)|\leq C,$ (2.25)

where the positive constant $C$ is independent of $x_{0}$.

###### Proof.

We only show this lemma in the case that $x\in\Omega_{+,x_{0}}(t)$. From (1.5)
and the assumption $\tau\leq t$, we estimate $A_{\pm}(\tau,t,x)$ as

$\displaystyle A_{\pm}(\tau,t,x)\leq$
$\displaystyle\exp\left(\int_{0}^{t}|a(t,x_{\pm}(\tau;t,x))|/2d\tau\right)$
$\displaystyle\leq$
$\displaystyle\exp\left(C\int_{0}^{t}a_{1}(\tau)+a_{2}(x_{\pm}(\tau;t,x))d\tau\right)$
$\displaystyle\leq$ $\displaystyle
C\exp\left(C\int_{0}^{t}a_{2}(x_{\pm}(\tau;t,x))d\tau\right),$

where we used the integrability of $a_{1}$. Noting (2.23) and (2.24), by the
change of variable $y=x_{\pm}(\tau)$, we have that

$\displaystyle\int_{0}^{t}a_{2}(x_{\pm}(\tau;t,x))d\tau\leq$
$\displaystyle\frac{C}{c_{1}}\left|\int_{x_{0}}^{x}a_{2}(y)dy\right|$
$\displaystyle\leq$
$\displaystyle\frac{C}{c_{1}}\int_{\mathbb{R}}a_{2}(y)dy<\infty.$

Hence we obtain the upper bounds of $A_{\pm}(\tau,t,x)$. The lower bounds can
be shown in the same as above. ∎

######  Lemma 2.3.

If $f\in C([0,T])$ satisfies that $f(t)\geq 0$ and

$\displaystyle f(t)\leq f(0)+C_{0}\int_{0}^{t}b(\tau)f(\tau)d\tau,$

then we have that

$\displaystyle f(t)\leq e^{C_{b}}f(0),$

where $b$ is a nonnegative integrable function and
$C_{b}=\int_{0}^{\infty}b(t)dt<\infty$.

###### Proof.

Setting

$\displaystyle F(t)=f(0)+\int_{0}^{t}b(\tau)f(\tau)d\tau,$

we have that from the assumption

$\displaystyle F^{\prime}(t)\leq b(t)F(t).$

Thus

$\displaystyle(B^{-1}(t)F(t))^{\prime}\leq 0,$

where $B(t)=\exp\left(\int_{0}^{t}b(\tau)d\tau\right)$. Since $B$ and $B^{-1}$
are bounded from the integrability of $b$, integrating this inequality over
$[0,t]$, we have the desired inequality. ∎

Next we introduce the key inequality to control $r$ and $s$. We define
$\Phi_{\pm,x_{0}}$ and $\Phi$ as follows

$\displaystyle\Phi_{\pm,x_{0}}(t)=\|r(t)\|_{L^{\infty}(\Omega_{\pm,x_{0}}(t))}+\|s(t)\|_{L^{\infty}(\Omega_{\pm,x_{0}}(t))}.$

and

$\displaystyle\Phi(t)=\|r(t)\|_{L^{\infty}(\Omega(t))}+\|s(t)\|_{L^{\infty}(\Omega(t))}.$

######  Lemma 2.4.

Let $x_{0}\in\mathbb{R}$. For $C^{1}_{b}$ solution of (1.4) satisfying (2.23).
If $\pm x_{0}>0$, it holds respectively that

$\displaystyle\Phi_{\pm,x_{0}}(t)\leq
C\Phi_{\pm,x_{0}}(0)=C(\|r(0,\cdot)\|_{L^{\infty}(\pm x\geq
x_{0})}+\|s(0,\cdot)\|_{L^{\infty}(\pm x\geq x_{0})}).$ (2.26)

We also have

$\displaystyle\Phi(t)\leq C\Phi(0)=C(\|r(0)\|_{L^{\infty}(x\geq
0)}+\|s(0)\|_{L^{\infty}(x\geq 0)}).$ (2.27)

###### Proof.

First we show (2.27). We only prove that

$\displaystyle\Phi_{+,0}(t)\leq C\Phi_{+,0}(0).$

The estimate of $\Psi_{-,0}(t)$ can be shown in the similar way. Applying
Lemma 2.2 to (2.8), we have that from (1.5) and (2.23) for $(t,x)$ with
$x\in\Omega_{+,0}(t)$

$\displaystyle|r(t,x)|\leq$ $\displaystyle
C|r(0,x_{-}(0))|+C\int_{0}^{t}|a(\tau,x_{-}(\tau))||s(\tau,x_{-}(\tau))|d\tau$
$\displaystyle\leq$ $\displaystyle C|r(0,x_{-}(0))|$
$\displaystyle+C\int_{0}^{t}(a_{1}(\tau)+a_{2}(x_{-}(\tau)))|s(\tau,x_{-}(\tau))|d\tau,$
(2.28)

where we took $x_{-}(s)=x_{-}(s;t,x)$. Similarly we have that

$\displaystyle|s(t,x)|\leq$ $\displaystyle C|s(0,x_{+}(0))|$
$\displaystyle+C\int_{0}^{t}(a_{1}(\tau)+a_{2}(x_{+}(\tau)))|r(\tau,x_{+}(\tau))|d\tau.$
(2.29)

Thus, taking $L^{\infty}$-norm with $x\in\Omega(t)$ in from (2.24), (2.29) and
(2.28), we have that

$\displaystyle\|r(t)\|_{L^{\infty}(\Omega_{+,0}(t))}\leq$ $\displaystyle
C\|r(0,\cdot)\|_{L^{\infty}(x\geq
0)}+C\int_{0}^{t}a_{1}(\tau)\|s(\tau)\|_{L^{\infty}(\Omega_{+,0}(\tau))}d\tau$
$\displaystyle+C\int_{0}^{t}\sup_{x\in\Omega_{+,0}(t)}a_{2}(x_{+}(\tau;t,x))\|s(\tau)\|_{L^{\infty}(\Omega_{+,0}(\tau))}d\tau$
(2.30)

and

$\displaystyle\|s(t)\|_{L^{\infty}(\Omega_{+,0}(t))}\leq$ $\displaystyle
C\|s(0,\cdot)\|_{L^{\infty}(x\geq
0)}+C\int_{0}^{t}a_{1}(\tau)\|r(\tau)\|_{L^{\infty}(\Omega_{+,0}(\tau))}d\tau$
$\displaystyle+C\int_{0}^{t}\sup_{x\in\Omega_{+,0}(t)}a_{2}(x_{+}(\tau;t,x))\|r(\tau)\|_{L^{\infty}(\Omega_{+,0}(\tau))}d\tau.$
(2.31)

We estimate $\sup_{x\in\Omega_{+}(t)}a_{2}(x_{+}(\tau))$. Putting
$x_{1}=x_{+}(0;t,x)$, we have from the definition of the characteristic curve
(see also Figure 1)

$\displaystyle x_{+}(\tau;t,x)=$ $\displaystyle
x+\int_{t}^{\tau}c(u(s,x_{+}(s;t,x)))ds$ $\displaystyle=$ $\displaystyle
x+\int_{t}^{0}c(u(s,x_{+}(s;t,x)))ds-\int_{\tau}^{0}c(u(s,x_{+}(s;t,x)))ds$
$\displaystyle=$ $\displaystyle x_{1}+\int_{0}^{\tau}c(u(s,x_{+}(s;t,x)))ds.$

Since $x_{1}\in\Omega_{+,0}(0)$ (namely $x_{1}\geq 0$), we obtain from (2.23)

$\displaystyle x_{+}(\tau;t,x)\geq\frac{c_{1}t}{4}.$ (2.32)

From the fact that $x_{+}(\tau;t,x)\leq x_{-}(\tau;t,x)$ with $\tau\leq t$, we
also obtain that

$\displaystyle x_{-}(\tau;t,x)\geq\frac{c_{1}t}{4}.$ (2.33)

Hence, from (1.7), we have that

$\displaystyle\sup_{x\in\Omega_{+}(t)}a_{2}(x_{\pm}(\tau))\leq
a_{2}(\frac{c_{1}t}{4}).$

Figure 2 : $\Omega_{\pm}(t)$, characteristic curves and the line
$x=\frac{c_{1}t}{4}$
$\Omega_{-,0}(t)$$(x,t)$$\Omega_{+,0}(t)$$x_{+}(\cdot;0,0)$$x_{-}(\cdot;0,0)$$x_{1}$$x_{-}(\cdot;t,x)$$x_{+}(\cdot;t,x)$$x$$t$$x=\frac{c_{0}t}{4}$

Applying (2.32) and (2.33) to (2.30) and (2.31) and summing up these
inequalities, we see that

$\displaystyle\Phi(t)\leq C_{1}\Phi(0)+C\int_{0}^{t}b(\tau)\Phi(\tau)d\tau.$

where we set $b(t)=a_{1}(t)+a_{2}(c_{1}t/4).$ From Lemma 2.3 and (1.6) we have
(2.27). Next we show (2.26) in the case that $x_{0}>0$. We remarks that the
property (2.24) holds for $\Omega_{+,x_{0}}$. Thus we can obtain similar
inequalities to in (2.30) and (2.31) hold for $\|r(t)\|_{\Omega_{+,x_{0}}(t)}$
and $\|s(t)\|_{\Omega_{+,x_{0}}(t)}$. For $x\in\Omega_{+,x_{0}}(t)$, one can
easily check that for $\tau\leq t$

$\displaystyle x_{\pm}(\tau;t,x)$ $\displaystyle\geq
x_{\pm}(0;t,x)+\frac{c_{1}\tau}{4}$ $\displaystyle\geq\frac{c_{1}t}{4}.$

Thus, from the same argument as in the case that $x_{0}=0$, we can obtain the
bounds for $\Phi_{+,x_{0}}(t)$. We can also show (2.26) for
$\Phi_{-,x_{0}}(t)$ in the same way as above. ∎

##  3 Proof of Theorem1.1

###  3.1 Continuity argument

Combining Lemma 2.4 and the continuity argument, we can show the following
lemma to control $u,v$ and $c$

######  Lemma 3.1.

If $\pm x_{0}>0$, then the following estimates hold on $[0,T^{*})$ for
sufficient small $\varepsilon>0$ with $t\in[0,T^{*})$ and
$x\in\Omega_{\pm,x_{0}}(t)$ respectively

$\displaystyle|u(t,x)-1|\leq C\varepsilon,$ (3.1) $\displaystyle|v(t,x)|\leq
C\varepsilon,$ (3.2) $\displaystyle|r(t,x)|+|s(t,x)|\leq C\varepsilon.$ (3.3)

Moreover, if $x_{0}=0$, the above three estimates hold for $t\in[0,T^{*})$ and
$x\in\Omega(t)$.

###### Proof.

We show this lemma by the continuity argument. We only consider the case that
$x_{0}>0$. We fix $T\in(0,T^{*})$ arbitrarily. From the continuity of the
solution, (2.23) holds neat $t=0$, if $\varepsilon>0$ is small. Let us denote
by $T_{1}$ the maximal time that (2.23) holds. We suppose that $T_{1}<T$ in
the contradiction argument. From Lemma 2.4, (3.3) holds on $[0,T_{1}]$, which
implies that (3.1) and (3.2) for sufficient small $\varepsilon$. Thus, for
sufficiently small $\varepsilon>0$, we find that it holds on
$(t,x)\in[0,T_{1}]\times\Omega_{+,x_{0}}(t)$

$\displaystyle\frac{c_{1}}{2}\leq c(u)\leq 2c_{1}.$

We note that this smallness of $\varepsilon$ is independent of $T$. Hence,
from the continuity of the solution, the estimate (2.23) can be extend to the
interval $[0,T_{2}]$ for $T_{1}<T_{2}\leq T$, which contradicts to the
definition of $T_{1}$. From the arbitrariness of $T$, we see that (2.23) holds
on $[0,T^{*})$, which implied (3.1), (3.2) and (3.3) hold on $[0,T^{*})$ by
Lemma 2.4. ∎

###  3.2 Conclusion

First we consider the case that $x_{0}=0$. We consider (2.19) on the
characteristic curve $x_{+}(\cdot;0,0)$ . It should be noted that the
following uniform estimate holds

$\displaystyle\int_{0}^{t}|a(\tau,x_{\pm}(\tau;0,0))|+|a_{t}(\tau,x_{\pm}(\tau;0,0))|+|a_{x}(\tau,x_{\pm}(\tau;0,0))|d\tau\leq
C.$

Applying this inequality, (3.1), (3.2), (3.3) and Lemma 2.2 to (2.19), we have
that

$\displaystyle Q(t)\leq Q(0)+K^{*}\varepsilon-C\int_{0}^{t}Q(\tau)^{2}d\tau,$
(3.4)

where $K^{*}$ is a positive constant depending on $M$ abd $\delta_{0}$. From
the assumption (1.11), it follows that

$\displaystyle Q(t)\leq-
C_{1}\varepsilon(K-K^{*})-C_{2}\int_{0}^{t}Q(\tau)^{2}d\tau,$ (3.5)

which implies that $-Q(t)$ blows up in finite time for sufficiently large $K.$
In fact, putting $q(t)$ as a solution of the integral equation:

$\displaystyle
q(t)=C_{1}\varepsilon(K-K^{*})+C_{2}\int_{0}^{t}q(\tau)^{2}d\tau,$

we can easily check that $-Q(t)\geq q(t)$. Solving this integral equation, we
have that

$\displaystyle
q(t)=\frac{C_{1}\varepsilon(K-K^{*})}{1-C_{1}C_{2}\varepsilon(K-K^{*})t}$

Thus we also obtain the upper estimate of the life-span for $Q(t)$

$\displaystyle T^{*}\leq\frac{1}{C_{1}C_{2}\varepsilon(K-K^{*})}.$

We can also show that $-Y(t)$ on the characteristic curve $x_{-}(\cdot;0,0)$
in (2.19) blows up in finite time.

In the case that $x_{0}>0$, we consider consider (2.19) on the characteristic
curve $x_{+}(\cdot;0,x_{0})$. Since
$x_{+}(\tau;0,x_{0})\in\Omega_{+,x_{0}}(\tau)$, we can apply Lemma 3.1 to the
right hand of (2.19) and find that (3.4) holds on $x_{+}(\cdot;0,x_{0})$. In
the case that $x_{0}<0$, we can show that $Y(t)$ on the characteristic curve
$x_{-}(\cdot;0,x_{0})$ in (2.19) blows up in finite time.

## Acknowledgments

The author would like to Thank Prof. Hiroyuki Takamura for his comment on
generalization of the damping coefficient. The research is supported by Grant-
in-Aid for Early-Career Scientists, No. 19K14573.

## References

  * [1] H. Cui, H. Yin, J. Zhang, and C. Zhu, Convergence to nonlinear diffusion waves for solutions of Euler equations with time-depending damping, J. Differential Equations, 264 (2018) 4564-4602.
  * [2] G. Chen, R. Young and Q. Zhang, Shock formation in the compressible Euler equations and related systems, J. Hyperbolic Differ. Equ., 10, (2013) 149-172.
  * [3] G. Chen, R. Pan and S. Zhu, Singularity formation for compressible Euler equation, SIAM J. Math. Anal., 49, (2017) 2591-2614.
  * [4] S. Chen, H. Li, J. Li, M. Mei and K. Zhang, Global and blow-up solutions to compressible Euler equations with time-dependent damping. J. Differential Equations, 268 (2020) 5035-5077.
  * [5] S. Geng and Y. Lin and M. Mei, Asymptotic behavior of solutions to Euler equations with Ttme-dependent damping in critical case, SIAM J. Math. Anal., 33 (2020) 1463-1488.
  * [6] K.O. Friedrichs, Nonlinear hyperbolic differential equations for functions of two independent variables, Amer. J. Math. 70 (1948) 555-589.
  * [7] F. Hou and H. Yin, On the global existence and blowup of smooth solutions to the multi-dimensional compressible Euler equations with time-depending damping, Nonlinearity, 30 (2017) 2485-2517.
  * [8] L. Hsiao and T. P. Liu, Convergence to nonlinear diffusion waves for solutions of a system of hyperbolic conservation laws with damping, Comm. Math. Phys., 143 (1992) 599-605.
  * [9] L. Hsiao and T. P. Liu, Nonlinear diffusive phenomena of nonlinear hyperbolic systems, Chinese Ann. Math. Ser. B 14 (1993) 65-480.
  * [10] L. Hsiao and D. Serre, Global existence of solutions for the system of compressible adiabatic flow through porous media, SIAM J. Math. Anal. 27 (1996) 70-77.
  * [11] P.D. Lax, Nonlinear hyperbolic equations, Comm. Pure Appl. Math. 6 (1953) pp. 231-258.
  * [12] P. D. Lax, Development of Singularities of solutions of nonlinear hyperbolic partial differential equations, J. Math. Phys. 5 (1964) 611-613.
  * [13] H. Li, J. Li, M. Mei and K. Zhang, Convergence to nonlinear diffusion waves for solutions of p-system with time-dependent damping, J. Math. Anal. Appl., 456 (2017) 849-871.
  * [14] T.-T. Li and W.-C. Yu, Boundary Value Problems for Quasilinear Hyperbolic Systems, Duke Univ. Press, Durham (1985).
  * [15] P. Marcati and A. Milani, The one-dimensional Darcy’s law as the limit of a compressible Euler flow, J. Differential Equations, 84 (1990) 129-147.
  * [16] P. Marcati and K. Nishihara, The $L^{p}$-$L^{q}$ estimates of solutions to one-dimensional damped wave equations and their application to the compressible flow through porous media, J. Differential Equations 191 (2003) 445-469.
  * [17] M. Mei, Best asymptotic profile for hyperbolic p-system with damping, SIAM J. Math. Anal., 42 (2010), 1-23.
  * [18] K. Nishihara, Asymptotic behavior of solutions of quasilinear hyperbolic equations with linear damping, J. Differential Equations 137 (1997) 384-395.
  * [19] R. Pan and Y. Zhu, Singularity formation for one dimensional full Euler equations, J. Differential Equations 261 (2016) 7132-7144.
  * [20] X. Pan, Remarks on 1-D Euler equations with time-decayed damping, arXiv:1510.08115.
  * [21] X. Pan, Global existence of solutions to 1-d Euler equations with time-dependent damping, Nonlinear Anal. 132 (2016) 327-336.
  * [22] X. Pan, Blow up of solutions to 1-d Euler equations with time-dependent damping. J. Math. Anal. Appl. 442 (2016) 435-445.
  * [23] Y. Sugiyama, Singularity formation for the 1D compressible Euler equations with variable damping coefficient, Nonlinear Analysis Series A: Theory, Methods & Applications, 170, (2018) 70-87.
  * [24] Y. Sugiyama, Remark on global existence of solutions to the 1D compressible Euler equation with time-dependent damping, Advanced Studies in Pure Mathematics, 85, (2020) 379-389.
  * [25] J. Wirth, Solution representations for a wave equation with weak dissipation, Math. Meth. Appl. Sci. 27 (2004) 101-124.
  * [26] J. Wirth, Wave equations with time-dependent dissipation. I. Non-effective dissipation, J. Differential Equations 222 (2006) 487-514.
  * [27] J. Wirth, Wave equations with time-dependent dissipation. II. Effective dissipation, J. Differential Equations 232 (2007) 74-103.
  * [28] N. J. Zabusky, Exact solution for the vibrations of a nonlinear continuous model string, J. Math. Phys. 3 (1962) 1028-1039.
  * [29] H. Zheng, Singularity formation for the compressible Euler equations with general pressure law, J. Math. Anal. Appl., 438 (2016) 59-72.


|  Residual stress in athermal soft disordered solids: insights from
microscopic and mesoscale models
---|---
|  Vishwas V.Vasisht∗a,b, Pinaki Chaudhuric and Kirsten Martensa
|  In soft amorphous materials, shear cessation after large shear deformation
leads to structures having residual shear stress. The origin of these states
and the distribution of the local shear stresses within the material is not
well understood, despite its importance for the change in material properties
and consequent applications. In this work, we use molecular dynamics
simulations of a model dense non-Brownian soft amorphous material to probe the
non-trivial relaxation process towards a residual stress state. We find that,
similar to thermal glasses, an increase in shear rate prior to the shear
cessation leads to lower residual stress states. We rationalise our findings
using a mesoscopic elasto-plastic description that explicitly includes a long
range elastic response to local shear transformations. We find that after flow
cessation the initial stress relaxation indeed depends on the pre-sheared
stress state, but the final residual stress is majorly determined by newly
activated plastic events occurring during the relaxation process. Our
simplified coarse grained description not only allows to capture the
phenomenology of residual stress states but also to rationalise the altered
material properties that are probed using small and large deformation
protocols applied to the relaxed material.

##

††footnotetext: ∗ E-mail<EMAIL_ADDRESS>a Univ. Grenoble
Alpes, CNRS, LIPhy, 38000 Grenoble, France.††footnotetext: b Dept. of Physics,
Indian Institute of Technology, Palakkad.††footnotetext: b The Institute of
Mathematical Sciences, Taramani, Chennai 600113, India.

##  1 Introduction

Soft glassy materials are ubiquitous in our everyday life and find their
application in various domains, such as food science, pharmaceutical
engineering or medical applications 1, 2. Being easier to handle and visualise
in experimental setups, these materials have become a playground for analysing
diverse mechanical response of amorphous solids in general, since many
features are common to both hard and soft materials.

Despite a large body of empirical knowledge, the fundamental understanding of
the mechanical behaviour of dense disordered materials is still a work in
progress 1, 2, 3, 4, 5. Besides many interesting properties in the deformation
process, like complex yielding phenomena and non-trivial rheology with
possible flow localization, these materials also display fascinating behaviour
once they are let to relax. The most striking observation is that they can
hold internal stresses even long after the cessation of their mechanical
stimulus 6. Such locked-in stresses, known as residual stress, can occur from
diverse procedures, e.g. at the end of thermal quenches (e.g., Rupert’s drop
7, 8), due to chemical processes (e.g., Tempered glass 9) or by the cessation
of some external mechanical drive 10. Understanding how such residual stresses
are created is of fundamental interest, within the broader framework of
identifying specific microscopic processes responsible for the observed
mechanical properties of amorphous systems. Further, knowing how such stresses
build up during the formation of the material, is also of practical
importance, since there can be both detrimental or beneficial aspects
depending upon the specific functional aspects of the material under
consideration 6. For example, in many cases, residual stresses can lead to
materials being more prone to fracture 11. On the other hand, such stresses
also allow some materials to be resistant or control the extent of failure 12,
13. Therefore, the ability to control the extent of residual stresses is
significant for the design of new material properties, and thus there is a
need for better understanding of the processes involved.

In the context of soft amorphous materials, various studies have recently
extensively explored the occurrence of residual stresses, via switching off
the flow driven by an external shear-rate 1, 14, 15, 16, 17, 18, 19, 20, 21.
For the case of non-Brownian systems (e.g. emulsions), upon switch off, the
shear stress, accumulated hitherto during flow, decays rapidly followed by a
slower relaxation to an eventual plateau, i.e. a stuck state with a well-
defined residual stress is obtained 19, 20. Further, a scaling relation
between initial stress, at the point of shear switch off, and the
corresponding locked-in stresses was revealed. For the case of Brownian
systems (e.g. colloidal glasses, gels), the eventual plateau is not observed.
Rather, a power-law decay in stress continues to occur, due to aging processes
characteristic to glasses 18. Mode coupling calculations for model glassy
systems also report similar behaviour 22.

In this work, we probe the origin of residual stresses in model athermal
amorphous systems, via a combination of microscopic and mesoscale simulations.
First, using molecular dynamics simulations of a model dense non-Brownian soft
amorphous material, we demonstrate that dynamically arrested states having
finite residual stresses are obtained upon the cessation of an applied shear,
consistent with experimental observations. As is also demonstrated via
experiments, the final residual stress depends upon the strain-rate at which
the system was being driven, prior to the switch-off of shear; larger shear-
rate for the initially flowing states lead to eventual arrested states having
lower residual stress. We rationalise our findings using a mesoscopic elasto-
plastic description that explicitly includes a long range elastic response to
local shear transformations.

In a serious of works, it has been shown that these type of elasto-plastic
descriptions account for various phenomena in the deformation process of dense
disordered media. Since their introduction for the description of local
plasticity 23, 24, 25, they have been used in a variety of contexts. Notably
they have extensively used to study critical fluctuations in the yielding
transition 26, 27, 28, 29, 30 as well as in the vicinity of a finite shear
rate critical point 31, 32. Further, they have been useful for the
understanding of strain localization 33, 34 and permanent shear banding 35,
36, and also creep phenomena 37, 38. And more recent developments have been
successful to establish a more quantitative link between these coarse grained
elasto-plastic descriptions and the microscopic dynamics measured in particle-
based simulations 39, 40, 41, 42, 43.

Our elastoplastic model, used in this study, exhibits all the phenomenological
observations related to residual stress measurements, observed in experiments
and our microscopic simulations. Using the ability to track local plastic
events in such elastoplastic models, we provide a semi-analytical analysis of
the formation of residual stresses after the flow cessation. Our significant
finding is that the final residual stress is majorly determined by newly
activated plastic events occurring during the relaxation process, after the
shear switch-off. Finally, using both the microscopic and mesoscale model, we
illustrate that states having less residual stresses are more rigid, which
becomes evident via the observed transient response to applied external shear.

The paper is organised as follows. After the introductory discussion in
Section 1, we elaborate the microscopic and mesoscale models that we have
studied in Section 2, along with the the methods involved in the numerical
simulations. In Sections 3 and 4, we report and discuss the measurement of
residual stresses in the microscopic and mesoscale models. In Section 5, we
provide the semi-analytical analysis rationalising the observations regarding
the residual stress measurements. In Section 6, we illustrate some of the
mechanical characteristics of the residual stress states. Finally, in Section
7, we provide a concluding discussion.

##  2 Model and protocol

###  2.1 Molecular dynamics simulations

In our study, the dense amorphous solid is modelled as a non-Brownian
suspension of soft repulsive spheres, at a volume fraction $\phi\approx 70\%$,
consisting of $97556$ particles, with repulsive effective interactions
mimicked via a truncated and shifted Lennard-Jones potential 44 given by
$U(r)=4\epsilon\left[(a_{ij}/r_{ij})^{12}-(a_{ij}/r_{ij})^{6}\right]+\epsilon$,
for $r_{ij}\leq 2^{1/6}a_{ij}$, else $U(r_{ij})=0$. Here $\epsilon$ is the
unit energy in the simulations, $r_{ij}$ being the center to center distance
between the particle $i$ and $j$ and $a_{ij}=0.5(a_{i}+a_{j})$, with $a_{i}$
and $a_{j}$ being the diameter of particles $i$ and $j$ respectively. The
diameters of the particles are drawn from a Gaussian distribution with
variance of $10\%$, whose mean is used as unit length $a$. We prepare the
initial samples by quenching high temperature liquid states
($T=5\epsilon/k_{B}$) to low temperature ($T=0.01\epsilon/k_{B}$) using a NVT
Molecular Dynamics protocol at a fixed cooling rate $\Gamma$
($10^{-4}\epsilon/k_{B}\tau_{0}$). Each sample is subsequently brought to the
closest energy minimum and to $k_{B}T/\epsilon\approx 0$ via energy
minimization. Under athermal conditions, these samples are subjected to a
shear rate $\dot{\gamma}$ using Lees-Edwards boundary conditions (LEBC) and
solving the following dissipative particle dynamics (DPD) based equations of
motion.

$m\frac{d^{2}{\bf r}_{i}}{dt^{2}}=-\zeta_{DPD}\sum_{j(\neq
i)}\omega(r_{ij})(\hat{r}_{ij}.v_{ij})\hat{r}_{ij}-\triangledown_{{\bf
r}_{i}}U$ (1)

where $m$ is the mass of the particle, the first term in the right hand side
(RHS) is the damping force which depends on the damping coefficient
$\zeta_{DPD}$. We have taken $\zeta_{DPD}=1.0$ that guarantees minimum
inertial effects 45. The relative velocity ${\bf v}_{ij}={\bf v}_{j}-{\bf
v}_{i}$ is computed over a cut-off distance $r_{ij}\leq 2.5a_{ij}$, with the
weight factor $\omega(r_{ij})=1$. The second term in the RHS is the force due
to interactions between particles. In all our simulations [x,y,z] dimensions
refers to flow, gradient and vorticity directions respectively. The shear
stress is computed from the complete virial stress tensor as
$\sigma_{xy}\equiv\sigma=\frac{1}{V}\sum_{i}\sum_{j>i}x_{ij}f^{y}_{ij}$, where
V($=lx*ly*lz$) is the volume of the system, $x_{ij}$ represents x-component of
the distance between particle $i$ and $j$ and $f^{y}_{ij}$ is the y-component
of force on the particle $i$ due to $j$.

Also, for the residual stress states, we compute the storage modulus by
performing small strain oscillatory shear simulations via the application of a
shear strain $\gamma(t)=\gamma_{0}sin(\omega t)$, with a strain amplitude of
$\gamma_{0}=1\%$. Computing the stress response for frequencies $\omega$, from
the steady state regime we extract visco-elastic coefficient using

$G^{\prime}(\omega,\gamma_{0})=\frac{\omega}{\gamma_{0}\pi}\int_{t_{0}}^{t_{0}+2\pi/\omega}\sigma_{xy}(t)sin(\omega
t)dt$ (2)

The storage modulus is obtained using $G^{\prime}=G^{\prime}(\omega\rightarrow
0,\gamma_{0})$.

All the microscopic simulations are done using LAMMPS 46, with a modification
to handle polydispersity in size.

###  2.2 Mesoscopic simulations

The mesoscopic simulations we perform in this study are based on an elasto-
plastic model for a yield stress material under steady shear at an applied
shear rate $\dot{\gamma}$. The model consists of a regular square lattice with
each lattice site representing an elasto-plastic element holding a shear-
stress variable evolve according to

$\frac{\partial\sigma_{i}}{\partial
t}=\mu\dot{\gamma}+\mu\sum_{j}G_{ij}\frac{\partial{\gamma}^{pl}_{j}}{\partial
t}$ (3)

where $\sigma_{i}$ is the scalar shear stress component at a site i, $\mu$ is
the material dependent shear modulus, and hence $\mu\dot{\gamma}$ is the
elastic contribution to the evolution of stress. The second term on the right
hand side (RHS) accounts for the change in stress due to local plastic
yielding. Here $\partial{\gamma}^{pl}_{j}/\partial t$ is the rate of plastic
strain deformation and $G_{ij}\propto\cos{(4\theta)}/r^{2}$ is the Eshelby
kernel which accounts for the stress redistribution due to a plastic event.
The sum over $j$ in the second term on the RHS represent the rate of local
elastic deformation associated with the response to a plastic deformation on a
distant site j. At the plastically deforming site, the relaxation dynamics is
modeled as Maxwellian visco-elastic relaxation
$\partial{\gamma}^{pl}_{j}/\partial t=(1/\mu\tau)n_{j}\sigma_{j}$, where
$\tau$ is the characteristic time for stress release in a plastic phase, a
system property which depends on the volume fraction $\phi$, $n_{j}$ is the
local state variable which refer to the local activity. $n_{j}=1$, if the site
is in the plastic phase, otherwise $n_{j}=0$. Similar to Ref. 47, a stochastic
dynamics is followed for the evolution of state variable $n$. If at a given
site $i$ the system has $\sigma_{i}<\sigma_{y}$, with site variable $n_{i}=0$,
according to equation 1, the stress changes due to elastic contribution
($\mu\gamma$). If the local stress $\sigma_{i}$ is greater than the local
yield stress $\sigma^{i}_{y}$, the $n_{i}$ transform from $0$ to $1$ within a
time period of $\tau_{pl}$. Once the system yields, the stress relaxation has
a local contribution as well in the second term in the equation 1. Within the
plastically active element the stress relaxes as $e^{-gt/\tau}$ during a
typical time $\tau_{res}$, where g is numerical value of the propagator at
site i, $G_{ii}$, that insures mechanical equilibrium. In two dimensions in
the large system limit the value of g tends to $g\approx 0.57$.

In our model the active state ($n_{j}=1$) relaxes back to $0$ with a rate
$\tau_{res}^{-1}$, where $\tau_{res}$ is the typical local restructuring time
of the material to regain its local elastic properties after a local yield
event. In our protocol the time scales $\tau_{pl}$ and $\tau_{res}$ are
typically smaller than the driving time scale $\tau_{a}=\dot{\gamma}^{-1}$.
The distribution of local yield stress threshold values follows from an
argument wherein yielding is considered as a jump in the potential energy
landscape (PEL) with energy barrier $E_{y}=\sigma^{2}_{y}/4\mu$. The
distribution is given by $P(E_{y})=\Theta(E_{y}-E^{min}_{y})\lambda
exp(-\lambda(E_{y}-E^{min}_{y}$)), where $E^{min}_{y}$ term considers only
larger jumps ignoring the small jumps within the metabasins of PEL, the
parameters $\lambda$ and $E_{min}$ determine the average yield strain
$\gamma_{y}$. We choose $E_{min}=0.001$ and $\lambda=2015$ such that we obtain
strain overshoot in the shear start-up regime while performing re-shear
simulations, reaching a residual stress state.

We emphasise here that our objective is not develop a mesoscale model of the
microscopic model that we study. Rather, the objective is to demonstrate that
the mesoscale model qualitatively reproduces the phenomenology observed in the
microscale model and thereafter use the simplicity of the mesoscale model to
gain insights into the physical processes involved in the creation of residual
stresses in amorphous solids.

Fig. 1: Microscale particulate model. (a) Relaxation of shear stress
($\sigma$) when the external drive (imposed shear-rate) is switched off, while
in steady state flow. The arrow shows the direction of increasing shear-rate
($\dot{\gamma}$). (inset) Load curve ($\sigma$ vs. strain $\gamma$) showing
the start-up to steady state regime, for different shear rates, prior to
switch-off. (b) The stress $\sigma_{I}$ at the shear switch-off (red square),
the residual stress $\sigma_{R}$ reached at the end of stress relaxation upon
switch-off (black circle) and $\Delta\sigma=\sigma_{I}-\sigma_{R}$ (blue
diamond), shown for two different system sizes (see labels) as a function of
the imposed shear rate.

##  3 Approach to residual stress states in microscopic simulations

In this section, we discuss our microscopic simulation results related
residual stress states and it’s dependence on the shear rate.

As discussed above, we initially impose an external shear deformation at a
chosen shear rate $\dot{\gamma}$ till the system reaches a steady state
regime. The corresponding variation of shear stress ($\sigma$) with shear
strain ($\gamma$), for different imposed $\dot{\gamma}$ is shown in the inset
of Fig.1(a). Once in steady state, we switch off the external shear and allow
the accumulated stress to relax while keeping the strain fixed. The stress
decreases with time and is observed to reach a finite saturation value, which
is referred to as the residual shear stress $\sigma_{R}$ of the system. To
obtain good statistics for $\sigma_{R}$, we not only perform three different
independent runs, but also relax the system from 5 different flow cessation
points for each run. Further, we check for finite size effects by probing the
behaviour for two different system size, viz. N=91556 and N=10976.

The relaxation of stress, for the different imposed shear-rates, is shown in
the main panel of Fig.1 (a). At a given shear rate, upon flow cessation, the
stress relaxes steadily at small times, followed by a dramatic decrease and
eventual arrest to a fixed value. Also, the timescale for the stress
relaxation is fastest for the largest $\dot{\gamma}$ and slows down with
increasing shear-rate. These features are similar to previous reported work on
emulsions 19 and attractive gels 21, and unlike thermal glasses where aging
continues after the initial stress relaxation 18. The dependence of
$\sigma_{R}$ upon the initially imposed $\dot{\gamma}$ is shown in the Fig. 1
(b). Consistent with previous reports, $\sigma_{R}$ decreases with increasing
$\dot{\gamma}$. Also shown in Fig. 1 (b), are the values of the shear stress
at flow cessation $\sigma_{I}$ and the difference between the stress in flow
and and the residual stress at arrest, $\Delta\sigma=\sigma_{I}-\sigma_{R}$.
The shear rate dependence of $\sigma_{I}$ is well described by the Herschel-
Bulkley function $\sigma_{y}+\kappa\dot{\gamma}^{n}$, as expected, where the
dynamic yield stress $\sigma_{y}=2.7\epsilon/a^{3}$ and the exponent $n=0.7$.
In the limit of vanishing $\dot{\gamma}$, as the stress in the system reaches
the dynamic yield stress, $\Delta\sigma$ approaches to zero, i.e. the residual
stress value coincides with the dynamical yield stress. On the other hand,
with the increase in shear rate where the $\sigma_{R}$ decrease,
$\Delta\sigma$ increases following a power-law profile (solid lines in the
Fig. 1 (b)). Finally, note that, we find that there is very little variation
observed with changing system size (as seen from the open symbols with the
opaque symbols in Fig. 1 (b)).

In the case of thermal glasses, Ballauff et al. 18 associate the residual
stresses with the long-lived memory effects of the pre-sheared glasses which
introduce supra-caging lengthscale in the system. Long time stress relaxation
are associated with ageing effects. In the case of non-Brownian colloidal gel
like systems, Mohan et al.19 find microstructural signature of the stress
relaxation which associate with different structural rearrangements. An
initial rapid relaxation is associated with unjamming process where in a
ballistic motion is observed and long time slow relaxation is associated with
cage breaking process. In this work we address the questions related to
approach towards residual stress state using a elasto-plastic model based
approach using mesoscopic simulations.

Fig. 2: Mesoscale model. (Top) Relaxation of shear stress ($\sigma$) when the
external drive (imposed shear-rate) is switched off, while in steady state
flow. The arrow shows the direction of increasing shear-rate ($\dot{\gamma}$).
Overall behaviour is similar to the response observed in microscale
simulations as shown in Fig.1. (inset) Load curve ($\sigma$ vs. strain
$\gamma$) showing the start-up to steady state regime, for different shear
rates, prior to switch-off. (Bottom) The stress $\sigma_{I}$ at the shear
switch-off (red square), the residual stress $\sigma_{R}$ reached at the end
of stress relaxation upon switch-off (black circle) and
$\Delta\sigma=\sigma_{I}-\sigma_{R}$ (blue diamond), shown for $N=1024^{2}$ as
a function of the imposed shear rate (system size dependence is shown in the
inset). Again, overall behaviour similar to what is observed in particulate
model; see Figs.2.

##  4 Approach to residual stress states in mesoscale simulations

Similar to particulate simulations, in the mesoscale simulations, upon flow
cessation, the sheared state relaxes to a state having finite residual stress
(Fig. 2(a)). The start-up to steady state evolution of the stress is shown in
the inset of top panel of Fig.2 (a), and the relaxation of the shear stress,
after the flow cessation, is shown in main panel Fig.2(a). In the mesoscale
model, the flow cessation is mimicked by turning off the first term
($\mu\dot{\gamma}$) in the RHS of Equation 1. The variation of the
$\sigma_{R}$ with the shear-rate during flow, $\dot{\gamma}$, is shown in
Fig.2 (b), along with the stress ($\sigma_{I}$) at which the switch off
occurs, and the corresponding stress difference
$\Delta{\sigma}=\sigma_{I}-\sigma_{R}$. The qualitative behaviour of all these
quantities is very similar to what has been observed for the particulate
system, as shown in Fig.1. The data for $\sigma_{I}$ can also be fitted with a
Herschel Bulkley law, obtaining $\sigma_{y}=0.06$ and $n=0.8$. Again, the
system size dependence for $\Delta{\sigma}$ does not show much variation
(Fig.2 (b) inset). Thus, overall, the mesoscale model can successfully
reproduce all the qualitative features related to the observations on residual
stress that we obtain in the microscopic model.

By construction, the mesoscale simulations provides easy access to local
yielding events or plastic events and its spatio-temporal evolution and hence
we utilize these simulations to analyse the process of approaching the locked-
in or residual stress state and it’s dependence on flow rate prior to shear
switchoff.

Fig. 3: Mesoscale model. (Top) Maps of state of activity at the point of shear
switch-off, for different driving rates (a) $10^{-2}$, (b) $10^{-3}$, (c)
$10^{-2}$. Black color correspond to active or plastic site ($n_{i}$=1) and
Yellow correponds to $n_{i}$=0. (Bottom) Corresponding maps of local stress
difference between initial and final stress, $\delta\sigma_{i}$ for the
corresponding shear rates in top. In all cases, only half of the simulation
domain is shown for clarity purposes.

###  4.1 Mesoscale stress relaxation maps

To begin with we study, in the mesoscale simulations, how well the the
configurations at the flow cessation regulate the final residual state.

In Fig. 3 (a-c) we show the local activity map of the initial configurations
for three different shear rates. In the maps the yielded sites ($n_{i}=1$)
(hereby termed active) is represented by black color and the yellow depicts
the elastic sites ($n_{i}=0$) (henceforth termed inactive). As expected, with
the decrease in shear rates (see Fig. 3 (a-c)), the number and also the
spatial density of active sites decreases. Naively, one expects that when the
shear is switched off, these active mesoblocks release their stress and become
inactive, leading to the residual stress states. However spatial maps of local
stress difference between initial flowing and final arrested states, viz.
$\Delta{\sigma}_{i}=(\sigma_{I}-\sigma_{R})_{i}$ (see Fig. 3 (d-f)) tell a
different story. We observe that $\Delta{\sigma}_{i}$ is finite in sites
beyond the initial active sites, which is very evident even at the low shear
limit, see Fig.3(f). This clearly indicates that the stress relaxation
process, upon flow cessation, involves more sites than the ones which were
active at the time of shear switch-off, and spatial scale of this become more
and more extensive with increasing pre-sheared $\dot{\gamma}$.

###  4.2 Activity statistics in mesoscale model

Following the above analysis vis-a-vis the stress difference maps, we now
proceed to quantify our observations.

We start by identifying the active sites at the time of shear switch-off and
compute the fraction of such sites in the system, $a_{0}$. During the stress
relaxation, we monitor these identified sites and check how many remain active
at any time instant, and thereby arrive at the time evolution of $a_{0}(t)$.
At any time instant, we also count the overall fraction of observed active
sites, $a_{s}(t)$. The difference, $a_{s}(t)-a_{0}(t)$, gives the fraction of
sites that have been newly activated. In the inset of Fig. 4 (a), for system
size of $1024^{2}$, we show how $a_{0}$, $a_{s}$ and the difference
$a_{s}-a_{0}$ evolve in time, for the case of relaxation from a state that was
being sheared at $\dot{\gamma}=10^{-3}\tau^{-1}$. We find that from very early
times, $a_{s}(t)>a_{0}(t)$, indicating the formation of new active sites. This
is evident from the behaviour of $a_{s}(t)-a_{0}(t)$, which grows with time as
relaxation proceeds, goes through a maximum and then eventually decays.
Further, we also note that $a_{0}(t)$, i.e. the original active sites, deplete
very soon, whereas $a_{s}(t)$ takes longer time to die out, implying that the
late time relaxation is due to the decay of the newly activated sites.

In the main panel of Fig. 4 (a), we show the time evolution of
$a_{s}(t)-a_{0}(t)$ for different values of shear-rate at which the system was
flowing prior to shear switch-off. In all cases, we observe the non-monotonic
growth and decay of the newly activated sites. The peak in $a_{s}(t)-a_{0}(t)$
is highest for the case of switch-off from the largest shear-rate and
decreases with decreasing shear-rate. The presence of larger fraction of
active sites at the point of switch off from a flow having large shear-rate,
is likely to lead to generate more active sites and therefore larger peak.
Next, as is the case with residual stress, the eventual decay in the number of
active sites is fastest for the case of shear switch-off from the flowing
state having the highest shear-rate.

In Fig. 4 (b), we show the rate of production of the new active sites,
directly measured from the mesoscale simulations. As is evident, the rate is
higher for the case of the flowing state having largest shear-rate and the
decay is also fastest, as pointed out above.

Both Fig. 4 (a) and (b) show that beyond a certain time the system no longer
creates new active sites (or no new local yielding events are observed) and
the existing active sites continue to relax stress.

Even though identifying plastic sites in molecular simulation is a challenging
task, using the knowledge of mean square displacement and the measure of non-
affine displacements $d^{2}_{min}$ , we confirm that a similar observations
are made in the particulate system.

Hence the emerging scenario is that the residual stress obtained at the end of
stress relaxation in these amorphous states, depends not only on the state of
activity in the initial sheared state but also on the dynamics during the
stress relaxation. We utilize these observations to provide a semi-analytical
rational for the magnitude of residual stresses obtained from the mesoscale
simulations.

Fig. 4: Mesoscale model. (Top) For different imposed shear-rates during the
steady flow, the variation in the number of newly activated sites (see text
for discussion) as a function of time, during stress relaxation after shear
switch-off. Arrow shows the direction of increasing shear-rate. (inset)
Variation of total number of activated sites (in black) with time, during the
relaxation, as well as the decay of pre-existing activated sites within the
configuration at switch-off (in red), apart from the newly activated sites (in
blue) (for $\dot{\gamma}=10^{-3}$). (Bottom) Rate of activation of new plastic
sites as a function of the time during the stress relaxation following switch-
off. Arrow shows the direction of increasing shear-rate.

##  5 Semi-analytical analysis of residual stress

We propose a simple analytical calculation which utilizes the observations we
discussed in the previous section to predict the residual stress values, with
inputs from mesoscale simulations.

In the calculation, we define the residual stress $\sigma_{R}$ as the stress
retained in the system after subtracting from the initial stress $\sigma_{I}$
the stress relaxed through the process of local yielding. In our mesoscopic
model, every time a site yields, the stress relaxes exponentially over a
certain time window as $\exp{(-gt/\tau)}$. Given that the local yielding
process is stochastic in nature, we will have a distribution the duration of
events $P(\tau_{res})$. For an initial sample (at the flow cessation), it is
plausible that fraction of the active sites might have yielded at earlier
times and hence would have already relaxed to an extent. Assuming a uniform
distribution of active sites, we can simply write residual stress as

$\sigma_{R}=\sigma_{I}-\int_{0}^{1}dxP_{1}(x)\int_{0}^{\infty}d\tau_{res}P(\tau_{res})\int_{0}^{x\tau_{res}}dt\frac{n^{ini}_{a}}{N^{2}}\left\langle\sigma^{ini}_{a}\right\rangle\exp(-gt/\tau)$
(4)

where $(n^{ini}_{a}/N^{2})$ gives the average fraction of active sites in the
system when the shear is switched off and $\left<\sigma^{ini}_{a}\right>$ is
the average stress per active site at that time. The product of the two gives
the total plastic stress which will eventually relax. We input this quantity
from the mesoscale simulations. With the simplified expression for fraction of
relaxation before shear cessation of the already active sites with
$P_{1}(x)=1$ being the uniform distribution on the interval zero to one and
$P(\tau_{res})=(1/\tau_{res})\exp{(-\tau_{res})}$ describing the stochastic
elastic recovery, we obtain the residual stress as

$\sigma_{R}=\sigma_{I}-\frac{n^{ini}_{a}}{N^{2}}\left\langle\sigma^{ini}_{a}\right\rangle\tau\tau_{res}\left(\frac{\tau_{res}\log(\tau_{res})}{g^{2}\tau^{2}}-\frac{\tau_{res}\log(\tau_{res}+g\tau)-g\tau}{g^{2}\tau^{2}}\right)$
(5)

Taking $\tau=1$ and $\tau_{res}=1$ as in the simulations we would like to
compare to, this expression simplifies to

$\sigma_{R}=\sigma_{I}-\frac{n^{ini}_{a}}{N^{2}}\left\langle\sigma^{ini}_{a}\right\rangle\left(\frac{g-log(1+g)}{g^{2}}\right)$
(6)

Fig. 5: Mesoscale model. (Left) (top) Fraction of active sites,
$n^{ini}_{a}/N^{2}$ in the flowing state, at the point of shear switch-off,
and (bottom) corresponding average stress per active site,
$\langle\sigma^{ini}_{a}\rangle$, as a function of shear rate, shown for
different system sizes. (Right) (top) The cumulative fraction of active sites,
$N^{new}_{a}/N^{2}$, computed from the flow cessation to the eventual relaxed
arrested state, and (bottom) associated stress per active site,
$\langle\sigma^{rel}_{a}\rangle$, averaged over the whole range of relaxation,
as a function of shear rates, shown for different system sizes. Fig. 6:
Mesoscale model. Comparison of residual stress values measured from
simulations and estimated from semi-analytical calculations (see main text).

In the Fig. 5 (a) and (b) we show $n^{ini}_{a}/N^{2}$ and
$\langle\sigma^{ini}_{a}\rangle$ respectively as a function of $\dot{\gamma}$
for various system sizes. The system size effects are visible at lower shear
rate limits but are within the sample to sample error bars. Using these values
in equation 6, with $g=0.57$, we estimate the residual stress as a function of
the shear-rate of the initial flow; see Fig. 6 (using blue symbols). Clearly,
the residual stress predicted solely from the decay of active site at the flow
cessation does not match with the numerically measured quantity. As discussed
before, we observe production of new plastic events; hence, one has to take
into account the occurrence and decay of the new active sites as well. We can
quite easily extend the argument presented in the form of equation 4, except
that for these newly activated sites we can drop the integral which took care
fraction of relaxed active sites. The area under the curves presented in
Fig.4(a) would give the total fraction newly activated sites
$N^{new}_{a}/N^{2}$, which is shown in Fig. 5 (c). The associated stress per
active sites $\langle\sigma^{rel}_{a}\rangle$ (averaged over the whole
relaxation time) is shown in Fig. 5 (d). Hence the total stress associated
with the newly created active sites is given by

$\sigma^{new}_{tot}=\frac{N^{new}_{a}}{N^{2}}\left\langle\sigma^{new}_{a}\right\rangle\frac{\tau\tau_{res}}{g\tau+\tau_{res}}$
(7)

Taking $\tau=1$ and $\tau_{res}=1$, the residual stress, including the
contribution from the initially active sites and newly activated sites, can be
written as

$\sigma_{R}=\sigma_{I}-\frac{n^{ini}_{a}}{N^{2}}\left\langle\sigma^{ini}_{a}\right\rangle\left(\frac{g-\log(1+g)}{g^{2}}\right)-\frac{N^{new}_{a}}{N^{2}}\left\langle\sigma^{new}_{a}\right\rangle\left(\frac{1}{1+g}\right)$
(8)

The residual stress computed from the above equation is shown in the Fig. 6
(dark green triangles). We find that our simple semi-analytical estimation of
residual stresses match quite well with the mesoscale simulation results.

##  6 Characterising residual stress states

After analysing the approach towards residual stress states, we now probe the
mechanical properties of these states. Understanding the rheological
properties of these states with frozen-in stress is quite relevant to
applications and designing materials .

###  6.1 Response to re-shearing

We probe the response of the residual stress state by imposing a sequence of
shear switch-on and switch-off, which we term as shear $\&$ relax cycles. The
observed evolution of shear stress due to this protocol is shown in Fig. 7 (a)
and (b), for microscopic and mesoscale simulations respectively.

Fig. 7: Reproducible residual stress states, is demonstrated via the variation
of stress with time, during a sequence of shear switch-off and switch-on using
the same imposed shear-rate. Data in the top panel corresponds to microscopic
model and data in the bottom panel corresponds to mesoscale model.

Firstly, we note that if we impose the shear on the residual stress state
using the same shear-rate as was done prior to the first switch-off and then
again switch-off the applied shear, and continue with such shear $\&$ relax
cycles, the intervening residual stress states obtained are always at the same
level; see Fig. 7 (a) and (b). Some minor variations, observed in the
microscopic simulations, are due to the fluctuations in the sampled states
from the steadily flowing states.

Secondly, we observe consistent stress overshoots during the shear switch-on
part of the shear $\&$ relax cycles, in both microscopic as well as mesoscale
simulations. This observation suggest that there is an unique transient shear
response characteristic to a residual stress state. We note that the magnitude
of the stress overshoot can be different from the one obtained from the shear
of a amorphous state prepared from a thermal quench, where the history of the
quench is relevant. Here, as we discuss later, it depends upon the residual
stress level.

Overall, the above discussion highlights that residual stress states can be an
ideal avenue for preparing amorphous states with reproducible response. Such
states can be used to study in greater details the start-up response of
disordered systems, both via elastoplastic as well as microscopic simulations.

Fig. 8: Rigidity of residual of stress states. (Top) Variation of storage
modulus, $G^{\prime}$ of the arrested states, with the imposed shear-rate of
the flowing state prior to switch-off, measured from the microscopic
simulations. Also marked are the values of the obtained residual stresses,
$\sigma_{R}$. (Bottom) Comparing local yield stress distribution for sheared
state and corresponding residual stress state, for two different shear rates
as marked. (a), (b) correspond to data from microscopic simulations and (c),
(d) correspond to data from mesoscale simulations.  Fig. 9: (a) Response
(stress vs strain curves) to imposed shear-rate
($\dot{\gamma}=10^{-3}\tau_{0}^{-1}1$) for states having different residual
stress value, as indicated by the arrow. (b) The stress-strain curves, with
the initial residual stress value deducted by $\sigma_{R}^{0}$, the residual
stress value associated with $\dot{\gamma}=10^{-3}$ demonstrating that states
having least residual stress demonstrates the largest stress overshoot.

###  6.2 Probing rigidity

Next, we analyse the rigidity aspects of the residual stress states by
computing the complex shear modulus using small amplitude oscillatory shear
protocol. The storage modulus $G^{\prime}$ obtained at low frequency limits
are show in the Fig. 8 (a) as a function of $\dot{\gamma}$, the shear-rate of
the flowing state prior to shear switch-off. We find that $G^{\prime}$
decreases with increase in residual stress $\sigma_{R}$ (marked alongside the
data points) or decrease in $\sigma_{I}-\sigma_{R}$. The loss modulus
$G^{\prime\prime}$, which is around two orders of magnitude smaller, also show
a similar dependence on the $\sigma_{R}$. For comparison, we also note that
the $G^{\prime}$ of the initial configuration, prepared via the thermal
quenching at a particular cooling rate, is around $75\epsilon/a^{3}$, which is
higher than that of the residual stress state having the lowest $\sigma_{R}$
in our study.

We further analyse the rigidity by computing the distribution of local yield
thresholds $(X=\sigma_{lcl}^{y}-\sigma_{lcl}^{0})$. In the microscopic
simulations we use the frozen matrix method  to compute X. The distribution
$P(X)$ computed (in both microscopic and mesoscale simulations) at the pre-
sheared state as well as at residual state is shown in Fig.8 (b-e) for two
different $\dot{\gamma}$. Note that we have scaled the local yield thresholds
by the macroscopic dynamic yield stress $\sigma_{y}$. In both microscopic and
mesoscale models, we observe that the mean of P(X) is lower at higher
$\dot{\gamma}$ of pre-sheared samples. This suggests that the residual states
have higher shear rigidity than the flowing state, which is expected. The P(X)
of residual states shows that mean is higher at higher pre-sheared
$\dot{\gamma}$, suggesting that lower $\sigma_{R}$ generate higher mean yield
threshold. This is a microscopic backing of what is observed in macroscopic
storage modulus $G^{\prime}$. Further analysis of local yield stress
distribution both in understanding the local yielding phenomenon  as well as
in the context of improving the mesoscale models  has to carried out
systematically.

Having demonstrated that the low lying residual stress states have more
rigidity, we now rationalise the transient response to the re-shear for a
range of residual states. In Fig. 9 (a) and (c) we show three different
residual states (obtained from three different pre-shear rates) subjected to
re-shear at the same imposed shear rate. The three residual states show
different transient response but eventually end up in the same steady state
stress value. The stress overshoot obtained during the transient response show
an initial residual stress state dependence. From the Fig. 9 (a) and (c) it
seems like lower residual stress states show lower overshoot and is counter
intuitive since the lower $\sigma_{R}$ corresponds to higher $G^{\prime}$ and
one expects from the previous studies  a higher stress overshoot. Once we
subtract the initial residual stress ($\sigma-\sigma_{R}$) such that all the
curves are compared at the same initial stress values (see Fig. 9 (b) and (c))
we find that the overshoot stress values increases with increase in rigidity
state (or decrease in $\sigma_{R}$).

##  7 Conclusion

In this work, we have investigated the origin of residual stresses in
athermally driven amorphous solids after shear cessation, using a combination
of microscopic and mesoscale elastoplastic models.

We first demonstrated that the mesoscopic model qualitatively reproduces the
phenomenology observed in microscopic simulations. In both cases, we obtain
arrested states exhibiting residual stresses after the forced shear flow is
switched off. And, the monotonic dependence of the residual stress on the
shear stress of the flowing state under the applied shear, is also
consistently observed. Thus, this provides the groundwork for using the
mesoscale model to gain further insight into the formation of residual
stresses, since these models have local plasticity as an intrinsic variable
which allows for the identification as well as tracking of active sites, be it
during shear or in the absence of it.

Thereafter, using the information provided by the mesoscale model, we develop
a semi-analytical argument to provide a reasonable estimate of the eventually
measured residual stress. The naive expectation is that when the applied shear
is switched off, the active sites available in the system at that instant
would release their stress and thereby the system would reach the arrested
state. We demonstrate, via the semi-analytical calculations, that such a
scenario does not explain the level of residual stresses observed in the
simulations. On the other hand, the spatial mesoscale maps, constructed from
the cumulative stress relaxation following the shear switch-off, reveal that
stress relaxation not only happens at these active sites but also in the
surroundings even long after the forcing is switched off. This implies that
new plastic events occur during the stress relaxation, which contribute to the
observation of new active sites. The number of these newly activates sites
eventually relaxes after a given time, and thereby the system reaches the
final state with a finite residual stress. The scale of the cascade of the new
events depends upon the initial stress in the steady flow, the larger the
initial stress, the more activity is generated during the relaxation following
shear switch off which leads to a lower residual stress. When the occurrence
of the newly activated events are appropriately included within the semi-
analytical calculations, we can reasonably reproduce the scale of residual
stresses observed in the simulations.

Finally, we have studied the mechanical response of the residual stress
states. We first observe, by imposing a sequence of steady deformation and
relaxation cycles, that the level of residual stress is reproducible,
expectedly, for a state flowing under a particular shear-rate prior to shear
switch-off. More importantly this sequence also enabled us to demonstrate that
the transient shear response of these residual stress states is also
reproducible, with the same level of stress overshoots in each cycle. This is
very significant for preparation of amorphous states with reproducible
rheological response. We also probed the rigidity of the residual stress
states in two ways, viz. by measuring the storage modulus via oscillatory
shear, and then by measuring local yield stress distributions. Both protocols
revealed that the states having lower residual stresses are more rigid and
have a higher yield threshold. Consequently, when we impose the same shear-
rate to states having different residual stresses, the state having the least
residual stress exhibits the largest stress overshoot. All these observations
are consistent across microscopic and mesoscale models.

##  8 Acknowledgements

All authors acknowledge financial support from CEFIPRA Grant No. 5604-1
(AMORPHOUS-MULTISCALE). K. M. acknowledges financial support of the French
Agence Nationale de la Recherche (ANR), under grant ANR-14-CE32-0005 (FAPRES).
We thank Jean-Louis Barrat, Emanuela Del Gado and Alberto Rosso for useful
discussions.

## References

  * Bonn _et al._ 2017 D. Bonn, M. M. Denn, L. Berthier, T. Divoux and S. Manneville, _Rev. Mod. Phys._ , 2017, 89, 035005.
  * Joshi and Petekidis 2018 Y. M. Joshi and G. Petekidis, _Rheologica Acta_ , 2018, 57, 521–549.
  * Coussot 2014 P. Coussot, _Journal of Non-Newtonian Fluid Mechanics_ , 2014, 211, 31–49.
  * Nicolas _et al._ 2018 A. Nicolas, E. E. Ferrero, K. Martens and J.-L. Barrat, _Reviews of Modern Physics_ , 2018, 90, 045006.
  * Rodney _et al._ 2011 D. Rodney, A. Tanguy and D. Vandembroucq, _Modelling and Simulation in Materials Science and Engineering_ , 2011, 19, 083001.
  * Withers 2007 P. Withers, _Reports on progress in physics_ , 2007, 70, 2211\.
  * Brodsley _et al._ 1986 L. Brodsley, F. C. Frank and J. W. Steeds, _Notes and Records of the Royal Society of London_ , 1986, 41, 1–26.
  * Kooij _et al._ 2021 S. Kooij, G. van Dalen, J.-F. Molinari and D. Bonn, _Nature communications_ , 2021, 12, 1–7.
  * Garfinkel and King 1970 H. Garfinkel and C. King, _Journal of the American Ceramic Society_ , 1970, 53, 686–691.
  * Withers and Bhadeshia 2001 P. J. Withers and H. Bhadeshia, _Materials science and technology_ , 2001, 17, 366–375.
  * Reiter _et al._ 2005 G. Reiter, M. Hamieh, P. Damman, S. Sclavons, S. Gabriele, T. Vilmin and E. Raphaël, _Nature materials_ , 2005, 4, 754–758.
  * Green _et al._ 1999 D. Green, R. Tandon and V. Sglavo, _Science_ , 1999, 283, 1295–1297.
  * Zhang _et al._ 2006 Y. Zhang, W. Wang and A. Greer, _Nature materials_ , 2006, 5, 857–860.
  * Chung _et al._ 2006 B. Chung, S. Ramakrishnan, R. Bandyopadhyay, D. Liang, C. Zukoski, J. Harden and R. Leheny, _Physical review letters_ , 2006, 96, 228301\.
  * Bandyopadhyay _et al._ 2010 R. Bandyopadhyay, P. H. Mohan and Y. M. Joshi, _Soft Matter_ , 2010, 6, 1462–1466.
  * Negi and Osuji 2009 A. S. Negi and C. O. Osuji, _Physical Review E_ , 2009, 80, 010404\.
  * Negi and Osuji 2010 A. S. Negi and C. O. Osuji, _Journal of Rheology_ , 2010, 54, 943–958.
  * Ballauff _et al._ 2013 M. Ballauff, J. M. Brader, S. U. Egelhaaf, M. Fuchs, J. Horbach, N. Koumakis, M. Krüger, M. Laurati, K. J. Mutch, G. Petekidis _et al._ , _Physical review letters_ , 2013, 110, 215701.
  * Mohan _et al._ 2013 L. Mohan, R. T. Bonnecaze and M. Cloitre, _Physical Review Letters_ , 2013, 111, 268301.
  * Mohan _et al._ 2015 L. Mohan, M. Cloitre and R. T. Bonnecaze, _Journal of Rheology_ , 2015, 59, 63–84.
  * Moghimi _et al._ 2017 E. Moghimi, A. R. Jacob and G. Petekidis, _Soft matter_ , 2017, 13, 7824–7833.
  * Fritschi _et al._ 2014 S. Fritschi, M. Fuchs and T. Voigtmann, _Soft matter_ , 2014, 10, 4822–4832.
  * Argon 1979 A. Argon, _Acta metallurgica_ , 1979, 27, 47–58.
  * Baret _et al._ 2002 J.-C. Baret, D. Vandembroucq and S. Roux, _Physical review letters_ , 2002, 89, 195506.
  * Picard _et al._ 2004 G. Picard, A. Ajdari, F. Lequeux and L. Bocquet, _The European Physical Journal E_ , 2004, 15, 371–381.
  * Lin _et al._ 2014 J. Lin, E. Lerner, A. Rosso and M. Wyart, _Proceedings of the National Academy of Sciences_ , 2014, 111, 14382–14387.
  * Liu _et al._ 2016 C. Liu, E. E. Ferrero, F. Puosi, J.-L. Barrat and K. Martens, _Physical review letters_ , 2016, 116, 065501.
  * Budrikis _et al._ 2017 Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser and S. Zapperi, _Nature Communications_ , 2017, 8, 15928.
  * Karimi _et al._ 2017 K. Karimi, E. E. Ferrero and J.-L. Barrat, _Physical Review E_ , 2017, 95, 013003.
  * Ferrero and Jagla 2019 E. E. Ferrero and E. A. Jagla, _Soft matter_ , 2019, 15, 9041–9055.
  * Le Goff _et al._ 2019 M. Le Goff, E. Bertin and K. Martens, _Physical review letters_ , 2019, 123, 108003.
  * Le Goff _et al._ 2020 M. Le Goff, E. Bertin and K. Martens, _Journal of Physics: Materials_ , 2020, 3, 025010.
  * Talamali _et al._ 2012 M. Talamali, V. Petäjä, D. Vandembroucq and S. Roux, _Comptes Rendus Mecanique_ , 2012, 340, 275–288.
  * Jagla 2007 E. Jagla, _Physical Review E_ , 2007, 76, 046119.
  * Martens _et al._ 2012 K. Martens, L. Bocquet and J.-L. Barrat, _Soft Matter_ , 2012, 8, 4197–4205.
  * Barbot _et al._ 2020 A. Barbot, M. Lerbinger, A. Lemaitre, D. Vandembroucq and S. Patinet, _Physical Review E_ , 2020, 101, 033001.
  * Merabia and Detcheverry 2016 S. Merabia and F. Detcheverry, _EPL (Europhysics Letters)_ , 2016, 116, 46003.
  * Liu _et al._ 2018 C. Liu, E. E. Ferrero, K. Martens and J.-L. Barrat, _Soft matter_ , 2018, 14, 8306–8316.
  * Puosi _et al._ 2014 F. Puosi, J. Rottler and J.-L. Barrat, _Physical Review E_ , 2014, 89, 042302.
  * Puosi _et al._ 2015 F. Puosi, J. Olivier and K. Martens, _Soft matter_ , 2015, 11, 7639–7647.
  * Albaret _et al._ 2016 T. Albaret, A. Tanguy, F. Boioli and D. Rodney, _Physical Review E_ , 2016, 93, 053002.
  * Patinet _et al._ 2016 S. Patinet, D. Vandembroucq and M. L. Falk, _Physical review letters_ , 2016, 117, 045501.
  * Liu _et al._ 2021 C. Liu, S. Dutta, P. Chaudhuri and K. Martens, _Physical Review Letters_ , 2021, 126, 138005.
  * Weeks _et al._ 1971 J. D. Weeks, D. Chandler and H. C. Andersen, _The Journal of Chemical Physics_ , 1971, 54, 5237–5247.
  * Vasisht and Del Gado 2020 V. V. Vasisht and E. Del Gado, _Physical Review E_ , 2020, 102, 012603\.
  * Plimpton 1995 S. Plimpton, _Journal of computational physics_ , 1995, 117, 1–19.
  * Martens _et al._ 2012 K. Martens, L. Bocquet and J.-L. Barrat, _Soft Matter_ , 2012, 8, 4197–4205.