# Creep response of athermal amorphous solids under imposed shear stress

Suman Dutta  International Centre for Theoretical Sciences - Tata Institute of
Fundamental Research
Survey No 51, Hessaraghatta Hobli, Shivakote, Bangalore 560089, India  Kirsten
Martens  Université Grenoble Alpes, CNRS, LIPhy, 38000 Grenoble, France
Pinaki Chaudhuri  The Institute of Mathematical Sciences, Taramani, Chennai
600113, India

( August 27, 2024)

###### Abstract

Yield stress materials fail when the imposed stress crosses a critical
threshold. A well-known dynamical response to the applied stress is the
phenomenon of creep where the cumulative deformation grows sublinearly with
time, prior to failure or arrest. Using extensive molecular dynamics
simulations, we study such response for a model amorphous system, in the
athermal limit, and probe how the annealing history of the initial state
determines the observed behaviour to an applied shear stress. Further, we
analyze the microscopic dynamics in the vicinity of the yield threshold, using
large systems, and characterize the spatiotemporal signatures towards arrest
or flow, at different scales.

##  I Introduction

The mechanical properties of amorphous materials determine the occurrence of
diverse natural phenomena (e.g. landslides, avalanches, etc.) as well a large
number of products and applications (e.g in the form of glasses, gels, foams,
emulsions, colloids etc.). Therefore, developing a physical understanding,
from a microscopic perspective, of the processes that lead to the observed
mechanical response is necessary and hence draws a significant amount of
research efforts utilizing a combination of experiments, numerical modelling
at micro as well as meso scales, and analytical frameworks (e.g. see reviews
Schuh  _et al._ (2007); Rodney  _et al._ (2011); Bonn  _et al._ (2017);
Nicolas  _et al._ (2018)). The characteristic feature of amorphous solids,
soft or hard, is the existence of a yield stress, exceeding which the material
fails or enters steady flow. The approach to yielding occurs via the
proliferation of local plastic events which self-organise to form extensive
avalanches. Even after yielding, such spatially organised events persist, e.g.
in the form of shear bands Fielding (2014). It is now well-established that
the transient response displayed by these materials depends on the preparation
history of the amorphous state, and recent theoretical work has focused on
probing this in terms of ductile/brittle response depending upon the degree of
annealing Ozawa  _et al._ (2018); Popović  _et al._ (2018); Barlow  _et al._
(2020); Richard  _et al._ (2021a, b); Rossi  _et al._ (2022); Pollard and
Fielding (2022).

One of the typical mechanical protocols involves the application of an
externally imposed step stress on the material and in response, amorphous
solids display an interesting dynamical feature known as creep, wherein the
time evolution of the emergent strain, resulting from the applied loading,
increases sub-linearly Betten (2008). In the case of thermal glasses or gels,
the creep process can in principle proceed indefinitely with diverse sub-
linear exponents reported in literature (e.g. see Bauer  _et al._ (2006);
Divoux  _et al._ (2011a); Leocmach  _et al._ (2014); Ballesta and Petekidis
(2016); Chaudhuri and Horbach (2013); Sentjabrskaja  _et al._ (2015); Landrum
_et al._ (2016); Cipelletti  _et al._ (2020); Merabia and Detcheverry (2016);
Joshi and Petekidis (2018); Cabriolu  _et al._ (2019); Dijksman and Mullin
(2022); Bhowmik  _et al._ (2022)) and a logarithmic growth has been proposed
in the asymptotic regime Siebenbürger  _et al._ (2012). However, in the case
of athermal or non-Brownian systems, the sub-linear growth can suddenly
terminate in the form of a dynamical arrest, if the amount of applied stress
is below the static yield stress threshold Paredes  _et al._ (2013); Bonn  _et
al._ (2017). Above the threshold, the material eventually reaches a steady
flow, as in the case of soft glassy materials, or can fail catastrophically,
as in the case of brittle solids. Although there are several experimental
studies of creep in non-Brownian materials, numerical or analytical
investigations have been limited. Only recently, using mean field or elasto-
plastic models, athermal creep has been analysed in some details Liu  _et al._
(2018a, b, 2021); Popović  _et al._ (2022). However, most of these studies
provide the overall macroscopic analysis and the spatio-temporal aspects need
to be explored in a systematic way to analyse the processes at play.

In this work, we study how a two-dimensional dense athermal amorphous solid
responds to applied shear stress, using extensive molecular dynamics
simulations. Our study is two-fold. The initial part focuses on demonstrating
how the annealing history for the preparation of the athermal solid determines
the macroscopic response, specially the transient regime, including the
determination of the static yield threshold. In the second half, we analyse
the mechanism of dynamical arrest or flow, around the yield threshold, and
reveal the spatio-temporal plasticity at work.

This manuscript is organised as follows. After the introductory discussion in
Section I, we provide details of the model and the numerical methods in
Section II. This is followed an extensive discussion of our findings in
Section III. Finally, we have a concluding discussion in Section IV.

##  II Model and Methods

In this study, we consider a model two-dimensional glass-forming binary
mixture having two different sizes Lançon and Billard (1988); Lançon  _et al._
(1986), whose mechanical response has been extensively studied Falk and Langer
(1998); Patinet  _et al._ (2016); Barbot  _et al._ (2018). The ratio of such
large (L) and small (S) is given by $\frac{N_{L}}{N_{S}}=\frac{1+\sqrt{5}}{4}$
where $N_{L}$ and $N_{S}$ are the number of L and S particles respectively.
The interaction between any pair of particles is the following:

$\displaystyle V_{\alpha\beta}(r)$ $\displaystyle=$ $\displaystyle
4\epsilon_{\alpha\beta}\left[\left(\sigma_{\alpha\beta}/r\right)^{12}-\left(\sigma_{\alpha\beta}/r\right)^{6}\right]\,.$
(1)

where $r=|{\bf r}_{i}-{\bf r}_{j}|$, and $\alpha,\beta$ correspond to the
identities $S$ or $L$. The values of the interaction parameters are set to
$\epsilon_{\textrm{SL}}=1.0$,
$\epsilon_{\textrm{SS}}=\epsilon_{\textrm{LL}}=0.5\epsilon_{\textrm{LS}}$,
$\sigma_{\textrm{LS}}=1$, $\sigma_{\textrm{LL}}=2\sin{\pi/5}$,
$\sigma_{\textrm{ss}}=2\sin{\pi/10}$. In the following, we use
$\epsilon_{\textrm{LS}}$ and $\sigma_{\textrm{LS}}$ as the unit for energy and
length, respectively. The cutoff radius in Eq. (1) is chosen as
$R_{c}=2.5\sigma_{\textrm{LS}}$ and the potential is smoothened out near the
cutoff Patinet  _et al._ (2019). As the time unit, we use
$\sqrt{{m\sigma_{\textrm{LS}}^{2}}/\epsilon_{\textrm{LS}}}$, where $m$ is the
mass of a particle that is considered to be equal for both type of particles,
i.e. $m=m_{\textrm{L}}=m_{\textrm{S}}=1.0$. More details about the model can
be found in Ref.Falk and Langer (1998); Patinet  _et al._ (2016); Barbot  _et
al._ (2018).

The equation of motion for any particle within the athermal assembly is given
by:

${m_{i}}\ddot{{\bf r}_{i}}=\sum^{N}_{j\neq i}{\bf f}^{int}_{(i,j)}+{\bf
f}^{D}_{(i,j)}$ (2)

Here, the total force acting on the particle is a sum over the interaction
forces coming from the potential discussed above, viz. ${\bf
f}^{int}_{(i,j)}=-\vec{\nabla}V_{(i,j)}$, and a dissipative force,
$\bf{f}^{D}_{(i,j)}=-\zeta w^{2}(r_{(i,j)})\leavevmode\nobreak\ (\hat{{\bf
r}}_{(i,j)}\leavevmode\nobreak\ .\leavevmode\nobreak\ {\bf
v}_{(i,j)})\hat{{\bf r}}_{(i,j)}$ which stabilizes the system against runaway
heating effects resulting from the continuously imposed external drive in the
form of shear. Apart from the equations of motion corresponding to the
particle co-ordinates, we also have an equation of motion of the macroscopic
shear-rate ($\dot{\gamma}$), resulting from the constant shear stress
($\sigma_{0}$) imposed along the $xy$-plane. The constant stress protocol is
implemented via a feedback control scheme Vezirov  _et al._ (2015); Cabriolu
_et al._ (2019), which leads to the following form for the time evolution of
shear rate due to the imposed stress:

${\ddot{\gamma}(t)}=B[\sigma_{0}-\sigma_{xy}(t)]$ (3)

where
$\sigma_{xy}=[m_{i}(v^{x}_{i}(t)v^{y}_{j}(t))+\sum^{N}_{j>i}(r^{x}_{(i,j)}\leavevmode\nobreak\
f^{y}_{(i,j)}(t))]/V$, is the Irving-Kirkwood expression for the shear stress,
using the velocity components, $v^{x}$ and $v^{y}$ at single particle level.
Here, we choose a suitable value of damping parameter, $B=1$ such that it can
capture the correct long time dynamics and the fluidization of the material.
Previously, as mentioned earlier, the rheological properties of this binay
mixture model has been studied and the dynamical yield threshold has been
reported to be $\sigma_{d}=0.5107$ Liu  _et al._ (2021) All imposed shear
stress values ($\sigma_{0}$) in this work will be expressed in units of
$\sigma_{d}$.

The molecular dynamics simulations were done using LAMMPS Plimpton (1995). The
cosnidered system sizes are $N=102400,409600$ which lead to respective box-
lengths of $L=316.174,632.348$ for the density at which we probe the
properties of the model system. The time step for the numerical integration
was chosen to be $\tau=0.001$. We prepare equilibrium thermal states at
$T=2.98,0.351$ and subsequently, we generate inherent structures states,
labelled $HTL$ and $ESL$, via energy minimization of states sample from the
respective equilibrated ensembles. The third ensemble, labelled $GQ$, was
prepared by slowly quenching the thermal sample at $T=0.351$ to $T=0.0299$ in
$10^{6}\tau$.

Figure 1: Response of states having different preparation histories.
$N=102400$. The three columns correspond to histories labelled as HTL, ESL, GQ
(see text for details). In each column, we show time evolution of ensemble-
averaged (a) macroscopic strain, $\langle\gamma(t)\rangle$ (b) macroscopic
shear-rate, strain-rate $\dot{\gamma}$ (labelled as $\langle
d\gamma(t)/dt\rangle$) (c) strain fluctuations, $\Delta^{2}\gamma(t)$ and (d)
potential energy fluctuations, $\Delta^{2}E_{p}(t)$. For HTL states, data is
shown for $\sigma_{0}/\sigma_{d}=0.39,0.489,0.587,0.685,0.783,0.881,1.077$.
For ESL states, data is shown for
$\sigma_{0}/\sigma_{d}=0.587,0.783,1.038,1.077,1.126,1.175,1.273$. For GQ
states, data is shown for
$\sigma_{0}/\sigma_{d}=1.468,1.547,1.566,1.615,1.664,1.713,1.762$.

Figure 2: Comparison across different preparation histories labelled as HTL,
ESL, GQ : (A) For $\sigma_{0}/\sigma_{d}=0.783$, where eventually flow is
arrested in all cases; and (B) for $\sigma_{0}/\sigma_{d}=1.664$ where there
is a long-time steady flow in all cases. In each case, time evolution of (a)
strain $\langle\gamma(t)\rangle$, (b) strain rate strain-rate $\dot{\gamma}$
(labelled as $\langle{d\gamma/dt}\rangle$), (c) strain fluctuations
$\Delta^{2}\gamma(t)$, and (d) MSD along gradient direction
$\langle{r^{2}_{y^{2}}}\rangle$, are shown. The data is shown for a system
size of $N=102400$, for the three different preparation protocols, in each
case, as marked.

##  III Results

###  III.1 Creep data: averages and fluctuations

We start our discussion of results by first analysing the macroscopic response
to applied shear stress. Since we are imposing a shear stress on the system,
our observable of interest is the cumulative deformation that the system
undergoes with time which is quantified via the shear strain, $\gamma(t)$, and
the corresponding rate of deformation quantified via the shear-rate
$\dot{\gamma}(t)$. When the rate of deformation, i.e. $\dot{\gamma}(t)$,
becomes constant, then $\gamma(t)\sim t$; this would correspond to a non-
equilibrium steady state. If the system does not eventually flow but gets
stuck, then the amount of deformation will remain fixed after that, i.e.
$\gamma(t)$ reaches a plateau and $\dot{\gamma}(t)$ goes to zero. An
intermediate regime is where the $\gamma(t)$ grows sublinearly with time,
which is termed as creep. To remind, below an applied stress threshold, we
expect the dynamics to get arrested, viz. $\dot{\gamma}(t)$ will vanish, and
above the threshold, steady flow, i.e. a constant $\dot{\gamma}(t)$, is
expected in compliance with the steady state macroscopic rheological flow
curve.

In the first and second rows of Fig.1, respectively, we show the data for
$\gamma(t)$ and $\dot{\gamma}(t)$, averaged over the ensemble of states
considered, with the three columns corresponding to the data for the three
different preparation histories, viz. HTL, ESL, GQ. The first important
finding is that the yield threshold, under applied stress, is dependent on
preparation history. For states labelled HTL and ESL, within the resolution of
our study, we observe flow for $\sigma_{0}/\sigma_{d}\geq 1.077$, and for
smaller values of $\sigma_{0}$, the system gets stuck at long times. However,
for GQ, the threshold changes dramatically; there is no flow for
$\sigma_{0}/\sigma_{d}\leq 1.468$. Thus, with improved annealing, i.e. for
more low lying inherent structures states on the energy landscape, the yield
threshold shifts to higher values. This is consistent with findings in
athermal quasistatic shear, where the peak stress during yield is observed to
increase with increase annealing Ozawa  _et al._ (2018). Since the yielding
stress threshold increases with preparation history, the flow rate at which
steady state is observed also changes; lower lying energy states have a larger
onset flow rate, to be consistent with the underlying macroscopic flow curve.

Note that the shape of the transient $\dot{\gamma}(t)$ have different trends
for the three different histories. For HTL, $\dot{\gamma}(t)$ decreases with
time, and then eventually vanishes or reaches a steady state. For the other
two histories, $\dot{\gamma}(t)$ decreases in the case when the system gets
arrested at long times. On the other hand, when there is eventual steady flow,
$\dot{\gamma}(t)$ has a non-monotonic shape, with $\dot{\gamma}(t)$ initially
increasing, goes through a minimum and then increases to the final steady
state. Broadly, these features in $\dot{\gamma}(t)$ observed in our numerical
simulations are consistent with experimental data Divoux  _et al._ (2011b);
Siebenbürger  _et al._ (2012).

Next we study the fluctuations in these macroscopic quantities, as measured
across the different trajectories within the ensemble starting from
independent amorphpus states. In this context, while studying creep in thermal
glasses, we have shown that the onset of large scale plasticity that leads to
steady flow is demarcated by a peak in strain fluctuations across
trajectories, and the timescale of the peak depends upon applied stress
Cabriolu  _et al._ (2019). These peaks are clearly visible for the GQ, ESL
states; see third row in Fig.1 . Further, we now see that the fluctuations in
the potential energy across trajectories also exhibit a similar behaviour,
with the peak approximately occurring when the peak in strain fluctatuations
occur; see fourth row in Fig.1. This allows for using energy fluctuations as a
simpler way to identify onset of plasticity. Note that the peak location is at
a time later than the location of the minimum in $\dot{\gamma}(t)$.

In Fig.2, we further analyse the nature of the transient behaviour, for the
three different preparation histories, by focusing on the response to the same
applied stress, viz. two different magnitudes, viz.
$\sigma_{0}/\sigma_{d}=0.783,1.664$. In the case of
$\sigma_{0}/\sigma_{d}=0.783$, the deformation gets arrested independent of
history (see left panel of Fig.2), whereas for $\sigma_{0}/\sigma_{d}=1.664$,
there is a steady flow at long times in all the three cases (see right panel
of Fig.2). In each case, we show data for the time evolution of (a) strain
$\langle\gamma(t)\rangle$, (b) strain rate strain-rate $\dot{\gamma}$, (c)
strain fluctuations $\Delta^{2}\gamma(t)$, and (d) single particle mean
squared displacement (MSD) measured along gradient direction
$\langle{r^{2}_{y^{2}}}\rangle$ which is a measure of non-affine motion
generated due to plasticity originating from the imposed shear. For the
smaller stress value, the system gets arrested very quickly for the ESL and GQ
states, whereas for the HTL states, there is a prolonged transient regime
before the system gets arrested; this is evident from both the macroscopic
strain data as well as the MSD data which shows a long time plateau indicating
that the motion of the particles is frozen in time. The situation is very
different for the larger applied stress, where there is a long time steady
state with a well-defined constant $\dot{\gamma}$ for all the three cases and
the MSD data also show identical long-time diffusive regime. We note that the
more low lying inherent structure states require longer time to reach the
steady-state, which is also known for the case of imposed shear rate (e.g.
Lamp  _et al._ (2022)). Interestingly, only in this case, there is pronounced
signal in the fluctuations data, with the heterogeneity across trajectories
larger for the low-lying states. Such non-monotonic behaviour in fluctuations
is not present in the case of the smaller applied stress. Thus, this exercise
also validates the use the fluctuations signal as a diagonistic for onset of
large scale steady flow.

###  III.2 Comparing paths to dynamical arrest and fluidization: a microscopic
analysis

After exploring the macroscopic response across the different ensembles
corresponding to different preparation histories, we now focus on probing at a
microscopic scale how fluidization proceeds in the vicinity of the estimated
yield threshold. For this study, we consider only the low-lying GQ states,
where the onset of plasticity is dramatic; see the data of shear-rate in the
third column of Fig.1 where we observe that typically $\dot{\gamma}(t)$
initially decreases, goes through a minimum, before a sudden increase to reach
the steady state.

In particular, we discuss the following thought experiment where we start with
the same initial state and follow the trajectories corresponding to two close-
by applied stress values of $\sigma_{0}/\sigma_{d}=1.547,1.557$; see
Fig.3(a)-(c) for data on strain, strain-rate and MSD. We can clearly see that
for this initial state and applied shear stress of
$\sigma_{0}/\sigma_{d}=1.547$, the single particle dynamics (quantified via
MSD) gets arrested at long times and thereby the macroscopic deformation also
stops; labelled NF. However, for $\sigma_{0}/\sigma_{d}=1.557$, the system
flows with a constant shear-rate at long times; labelled F. Thus, with a very
small increase in the magnitude of the applied shear stress, the fate of the
system, starting from the same initial state, changes from being arrested at
long times to that of a steady flow. However, the interesting and significant
feature is that the evolution of strain, shear-rate, MSD with time are very
similar for both cases, till around $10^{2}$, after which the trajectories
diverge. If one follows $\dot{\gamma}(t)$, then beyond this point of
divergence, for $\sigma_{0}/\sigma_{d}=1.557$, there is suddenly a rapid
increase and the system enters a steady flow regime whereas for
$\sigma_{0}/\sigma_{d}=1.547$, it continues to decrease and eventually
vanishes to become arrested, i.e. the system gets absorbed in a local minima
of its underlying potential energy landscape. The MSD data also show similar
signatures – for the smaller stress, it evolves to the asymptotic plateau,
whereas for the larger stress, there is a jump after $t\approx{10^{2}}$
indicating largescale breaking of cages and then proceeds towards diffusion.
Thus, the arrested state obtained for $\sigma_{0}/\sigma_{d}=1.547$ could be
labelled as a marginally stable absorbed state, which gets destabilized and
goes into eventual steady flow at the slightest increase of imposed stress.

To analyse at a spatial scale the divergence of trajectories for the two
applied shear stresses, we use two-dimensional local MSD maps, using a coarse-
graining scale of $\ell=2.47$ which gives a mesh of $256\times 256$; see the
panels in Fig.3(d), where the maps (A)-(E) in each row corresponds to five
different times as marked in Fig.3(a)-(c). The top row corresponds to the case
of $\sigma_{0}/\sigma_{d}=1.547$, where there is no flow at the end and the
system gets arrested, and the bottom row corresponds to the case of
$\sigma_{0}/\sigma_{d}=1.557$ where we have a steady flow at the end. All
displacements are measured relative to the initial state, which is the same
for both cases. We clearly observe that even at the local scale, the MSD maps
are very similar for (A)-(C), i.e. for $t=0.06,87.6,167.7$. From (D) onward,
i.e. $t=377.3$, the maps differ extensively. This indicates that for both the
applied stresses the system evolves through similar trajectories in
configuration space till the point of divergence. For the eventually flowing
state, we clearly saw strong heterogeneity in the flow pattern in (D), with a
vertical band of large mobility around a possible slip line.

To characterize the local plasticity, we compute maps of local vorticity
$\omega=\nabla\times\tilde{\bf v}$, where $\tilde{\bf v}$ is the coarse-
grained velocity field, using the same mesh size detailed above; these maps
are shown in the panels in Fig.3(e), for the both the stress magnitudies and
for the same time points for which the MSD maps have been constructed. At
earlier times (e.g. see (C)), there are only some sporadic vortices scattered
here and there, in both cases. However, the main map of interest is (D), which
for the case of $\sigma_{0}/\sigma_{d}=1.547$, has no spatial signal, whereas
for the case of $\sigma_{0}/\sigma_{d}=1.557$, we observe local vortices
aligned along the percolating slip line discussed above. Thus, large scale
flow for the larger stress proceeds via the formation of these spatially
spanning structure consisting of local vortices. Other such vortical structure
can also been seen elsewhere across the system in this case, indicating that
plastic events have started to occur across the system as it transits to
steady flow.

Figure 3: Comparative study of trajectories for $\sigma_{0}/\sigma_{d}=1.547$
and $\sigma_{0}/\sigma_{d}=1.557$, starting from same initial state.
$N=409600$. Time evolution of (a) strain, $\gamma(t)$, (b) strain-rate
$\dot{\gamma}$ (labelled as $d\gamma/dt$), (c) mean-squared displacement,
$r^{2}_{y}(t)$: trajectory for $\sigma_{0}/\sigma_{d}=1.557$ gets fluidized
(F) whereas it gets arrested, i.e. non-fluidized (NF) for
$\sigma_{0}/\sigma_{d}=1.547$. (d) Local MSD maps (see text for details),
using a coarse-graining scale of $\ell^{2}$, with $\ell=2.47$, for
$\sigma_{0}/\sigma_{d}=1.547$ (first row; labelled NF) and
$\sigma_{0}/\sigma_{d}=1.557$ (second row; labelled F). The maps are shown at
the time points $t=0.06$ (A), $87.6$ (B), $167.7$ (C), $377.3$ (D) and
$3369.7$ (E), as marked in (a), (b) (c). (e) Local vorticity maps (see text
for details), using coarse-graining scale used in (d), computed at the same
time points for both NF (first row) and F (second row) trajectories.

###  III.3 Probing response at different coarse-graining scales

Finally, we probe and characterise the response to stress at different scales,
to obtain a multi-scale flavour of the ensuing plasticity.

As a first step in this direction, we compute spatially coarse-grained self-
overlap function, $Q_{\ell}(X,Y;t)$, measured at different times $t$. For a
given coarse-grained mesh-block of size $\ell$, centered at position $(X,Y)$,
the self-overlap function is defined as
$Q_{\ell}(X,Y;t)=\sum_{i=1}^{N_{XY}}w(\delta{r}(t))$ where
$\delta{r}(t)=|\vec{r_{i}}(t)-\vec{r_{i}}(0)|$ is the displacement from the
initial state at $t=0$, $N_{XY}$ is the number of particles located in the
mesh at $(X,Y)$ at time $t=0$, and the window-function, $w(\delta{r}(t))=1$ if
$\delta{r}(t)\leq a$ or otherwise $w(\delta{r}(t))=0$. A threshold of
$a\approx 0.2$ allows to pick up positional changes due to plasticity.
Further, we define that a mesh has fluidized when $Q_{\ell}(X,Y;t)\leq 0.1$,
which provides some kind of first passage time-scale for initiation of plastic
activity. Via this, we obtain the timescale maps, $\tau_{F}(X,Y)$ where
$\tau_{F}$ is the fluidization timescale of each individual mesh-block.. This
exercise is done for different coarse-graining scales, viz.
$\ell=19.76,9.88,4.94,2.47$ which gives number of mesh-blocks to be,
respectively, $n_{L}=32,64,128,256$.

For the F and NF trajectories discussed above, the corresponding timescale
maps are shown in Fig.4(a); see Appendix (Fig. A1 and Fig. A2) for examples of
corresponding time evolution of $Q_{\ell}(X,Y;t)$. In the case of
$\sigma_{0}/\sigma_{d}=1.547$, we observe an interesting feature, viz. that
there exists a percolating backbone on the right hand side which does not
yield over the observation time window $(t\approx 3.39X10^{3})$. This implies
that this region is resistant to yield at this applied stress, and therefore
we infer that this rigid backbone leads to the final dynamical arrest as
discussed in the context of Fig.3. In contrast, this backbone fluidizes for
$\sigma_{0}/\sigma_{d}=1.557$ within the observation window, which makes this
a flowing state. Thus, it is this resistant backbone which makes the
difference between arrest and flow.

A pertinent question to ask is what triggers the breakdown of this backbone
for a small increase in stress. Our analysis of the difference in activity
(see Fig. A3 in Appendix) between the two trajectories seem to indicate that
increased plastic activity, at early times, in the area around where the
eventual slip-line emerges at later times (see sub-panel (D) in Fig.3(d)),
seems to contribute to the toppling of the region which is rigid for the
marginally stable absorbed state at $\sigma_{0}/\sigma_{d}=1.547$. Whether
this is the generic mechanism that triggers yielding of a marginally absorbed
state needs to be examined in more details, which will be explored in future
work.

The diversity of local fluidization timescales is highlighted for smaller
coarse-graining scales, e.g. in the case of the F trajectory, but gets
smoothened out over larger scales, as expected. On the other hand, coarse-
graining over larger $\ell$ brings into focus the robust backbone visible for
the NF trajectory. Further, for $\sigma_{0}/\sigma_{d}=1.557$, when we
visualize the vorticity map at $t=377.3$, using these different coarse-
graining scales, complex spatio-temporal patterns are visible at nearly
granular scales, which get smoothened out with increasing coarse-graining,
with the slip line remaining the dominant feature; see Fig.4(b).

It is possible to obtain a measure of the dynamical heterogeneity, $\chi(t)$,
from the instantaneous spatial fluctuations in $Q_{\ell}(X,Y;t)$. The data for
$\chi(t)$ is shown in in Fig.4(c) (i)-(ii); for the different coarse graining
scales $\ell$ (labelled I-IV) the data is shown in the top two panels for the
F and NF trajectories discussed above. We observe that $\chi(t)$ is non-
monotonic in both cases. For the trajectory where complete fluidization
occurs, $\chi(t)$ decays to zero, while for the the case where the system gets
arrested, $\chi(t)$ reaches a plateau capturing the extent of fluidization
undergone till arrest. The peak scale of fluctuations is of course more for
the least coarse-graining scale and for the choice of $a=0.2$, the peak in
$\chi(t)$ occurs around $t\approx{30}$. These fluctuations are indicators of
the plastic events that lead the system towards some local minima (see Fig.3,
but not the large-scale plasticity that leads to the system to steady flow;
this would correspond to larger $a$.) The peak shifts to larger times for
increasing choice of $a$, as can be seen in Fig.4(c) (iii)-(iv) where we show
the fluctuations, labelled as $\chi_{a}(t)$, for different $a$ values.

Figure 4: Probing response at various lengthscales. (a) Maps of local
fluidization timescale $\tau_{F}$ computed at different coarse-grained scales
$\ell=19.76$ (I), $9.88$ (II), $4.94$ (III), $2.47$ (IV), for (left) the non-
fluidized (NF) state and (middle) fluidized (F) state, as indicated in Fig.3.
(b) Vorticity maps for the fluidized trajectory, measured at $t=377.3$ using
the same $\ell$ scales as listed in (a). (c) Spatial dynamical fluctuations,
$\chi(t)$, computed via the overlap function $Q_{\ell}(X,Y;t)$, for the (i)
non-fluidized and (ii) fluidized trajectories, using the $\ell$ scales
labelled I-IV in (a). (iii)-(iv) Variation in fluctuations, $\chi_{a}(t)$,
with choice of cutoff ($a$) in defining $Q_{\ell}(X,Y;t)$ for F and NF
trajectories for different choices of $a$ as marked.

##  IV Summary and perspectives

To summarize, we have studied the response of athermal amorphous states, which
correspond to inherent structures of the model system, to applied shear
stress.

At first, we have studied the macroscopic response as well as corresponding
fluctuations for different preparation histories of the inherent structures,
which populate lower energy levels with increased annealing. The main
observation, in this context, is that the yield threshold ($\sigma_{y}$)
increases with decrease in the energy level of the initial undeformed state.
For imposed stress values above the respective thresholds, the system reaches
a long-time steady state flowing with a finite a shear-rate, which is of
course independent of the history and compliant with the bulk flow curve.
Below the threshold, the system reaches arrested states after some transient
activity. Interestingly, the transient lifetime to eventual arrest or flow
also depends upon history of preparation; more low lying states get arrested
quickly if $\sigma_{0}<\sigma_{y}$ and take longer to reach steady state when
$\sigma_{0}>\sigma_{y}$. The onset of fludization is also more sudden for such
states. We show that the onset time for the large-scale plastic flow can be
obtained by measuring the energy or strain fluctuations, across the ensemble
of initial states: maximal fluctuation mark the onset of large-scale
plasticity and this increases as the imposed stress decreases towards
$\sigma_{y}$.

Subsequently, we have done a spatio-temporal analysis of the response of the
system for imposed stresses around the yield threshold, for the case of the
lowest lying energy levels of the model that we could access via conventional
annealing dynamics. For an imposed stress just below $\sigma_{y}$, we
identified that the occurrence of percolating clusters resist yielding and
this leads to dynamical arrest. These arrested states are marginally stable,
since with a small increase in the applied stress, this rigid cluster breaks
down due to the imposed stress and the system eventually fluidizes. Our
preliminary analysis indicates that the increases activity emanating from the
soft spots leads to the local yielding and eventual flow. The onset of
largescale flow occurs via spatio-temporal patterns that can be visualised in
non-affine velocities as well as local vorticity, forming slip-plane like
structures spanning the system. Furthermore, there is a wide variation of
local yielding timescales during the onset of flow, which we characterise by
quantifying the spatiotemporal fluctuations.

Finally, we note that our work can also be contextualised to the recent
explorations of the relaxation dynamics within the energy landscapes of
amorphous solids; e.g. see Nishikawa  _et al._ (2022); Chacko  _et al._
(2019); Vasisht  _et al._ (2022); Mandal and Sollich (2020). The distinction
in our case is that the imposed stress alters the energy landscape and thus
one studies the search dynamics in this altered space. However, similar to
observations in some of these earlier studies, there is a power-law approach
in the dynamics to arrest, as is manifested in the observed macroscopic shear-
rate in the protocol that we study. We also note that a recent mean field
study of fatigue behaviour in cyclic shear of amorphous systems reported yield
rate curves which are very similar to shear-rate curves that we observe and
similar bifurcation in the vicinity of an yield threshold Parley  _et al._
(2022), suggesting some generic scenario of marginal stability around yield
that needs to be explored further.

##  V Acknowledgements

KM and PC acknowledge funding for this project via the CEFIPRA Project 5604-1.
SD acknowledges support of the Department of Atomic Energy, Government of
India, under project no RTI4001. We thank the HPC facility at IMSc Chennai for
providing computational resources. We also thank P. Sollich, M. Ozawa for
useful discussions.

## References

  * Schuh  _et al._ (2007) C. A. Schuh, T. C. Hufnagel, and U. Ramamurty, Acta Materialia 55, 4067 (2007).
  * Rodney  _et al._ (2011) D. Rodney, A. Tanguy, and D. Vandembroucq, Modelling and Simulation in Materials Science and Engineering 19, 083001 (2011).
  * Bonn  _et al._ (2017) D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Reviews of Modern Physics 89, 035005 (2017).
  * Nicolas  _et al._ (2018) A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Reviews of Modern Physics 90, 045006 (2018).
  * Fielding (2014) S. M. Fielding, Reports on Progress in Physics 77, 102601 (2014).
  * Ozawa  _et al._ (2018) M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Proceedings of the National Academy of Sciences 115, 6656 (2018).
  * Popović  _et al._ (2018) M. Popović, T. W. de Geus, and M. Wyart, Physical Review E 98, 040901 (2018).
  * Barlow  _et al._ (2020) H. J. Barlow, J. O. Cochran, and S. M. Fielding, Physical Review Letters 125, 168003 (2020).
  * Richard  _et al._ (2021a) D. Richard, C. Rainone, and E. Lerner, The Journal of Chemical Physics 155, 056101 (2021a), https://doi.org/10.1063/5.0053303 .
  * Richard  _et al._ (2021b) D. Richard, E. Lerner, and E. Bouchbinder, MRS Bulletin 46, 902 (2021b).
  * Rossi  _et al._ (2022) S. Rossi, G. Biroli, M. Ozawa, G. Tarjus, and F. Zamponi, Physical Review Letters 129, 228002 (2022).
  * Pollard and Fielding (2022) J. Pollard and S. M. Fielding, Physical Review Research 4, 043037 (2022).
  * Betten (2008) J. Betten,  _Creep mechanics_ (Springer Science & Business Media, 2008).
  * Bauer  _et al._ (2006) T. Bauer, J. Oberdisse, and L. Ramos, Physical review letters 97, 258303 (2006).
  * Divoux  _et al._ (2011a) T. Divoux, C. Barentin, and S. Manneville, Soft Matter 7, 9335 (2011a).
  * Leocmach  _et al._ (2014) M. Leocmach, C. Perge, T. Divoux, and S. Manneville, Physical Review Letters 113, 1 (2014), arXiv:1401.8234 .
  * Ballesta and Petekidis (2016) P. Ballesta and G. Petekidis, Physical Review E 93, 042613 (2016).
  * Chaudhuri and Horbach (2013) P. Chaudhuri and J. Horbach, Physical Review E 88, 040301 (2013).
  * Sentjabrskaja  _et al._ (2015) T. Sentjabrskaja, P. Chaudhuri, M. Hermes, W. Poon, J. Horbach, S. Egelhaaf, and M. Laurati, Scientific reports 5, 11884 (2015).
  * Landrum  _et al._ (2016) B. J. Landrum, W. B. Russel, and R. N. Zia, Journal of Rheology 60, 783 (2016).
  * Cipelletti  _et al._ (2020) L. Cipelletti, K. Martens, and L. Ramos, Soft matter 16, 82 (2020).
  * Merabia and Detcheverry (2016) S. Merabia and F. Detcheverry, EPL (Europhysics Letters) 116, 46003 (2016).
  * Joshi and Petekidis (2018) Y. M. Joshi and G. Petekidis, Rheologica Acta 57, 521 (2018).
  * Cabriolu  _et al._ (2019) R. Cabriolu, J. Horbach, P. Chaudhuri, and K. Martens, Soft matter 15, 415 (2019).
  * Dijksman and Mullin (2022) J. A. Dijksman and T. Mullin, Physical Review Letters 128, 238002 (2022).
  * Bhowmik  _et al._ (2022) B. P. Bhowmik, H. Hentschel, and I. Procaccia, Physical Review E 106, 034906 (2022).
  * Siebenbürger  _et al._ (2012) M. Siebenbürger, M. Ballauff, and T. Voigtmann, Physical review letters 108, 255701 (2012).
  * Paredes  _et al._ (2013) J. Paredes, M. A. Michels, and D. Bonn, Physical review letters 111, 015701 (2013).
  * Liu  _et al._ (2018a) C. Liu, K. Martens, and J.-L. Barrat, Phys. Rev. Lett. 120, 028004 (2018a).
  * Liu  _et al._ (2018b) C. Liu, E. E. Ferrero, K. Martens, and J.-L. Barrat, Soft matter 14, 8306 (2018b).
  * Liu  _et al._ (2021) C. Liu, S. Dutta, P. Chaudhuri, and K. Martens, Physical Review Letters 126, 138005 (2021).
  * Popović  _et al._ (2022) M. Popović, T. W. de Geus, W. Ji, A. Rosso, and M. Wyart, Physical Review Letters 129, 208001 (2022).
  * Lançon and Billard (1988) F. Lançon and L. Billard, Journal de Physique 49, 249 (1988).
  * Lançon  _et al._ (1986) F. Lançon, L. Billard, and P. Chaudhari, EPL (Europhysics Letters) 2, 625 (1986).
  * Falk and Langer (1998) M. L. Falk and J. S. Langer, Physical Review E 57, 7192 (1998).
  * Patinet  _et al._ (2016) S. Patinet, D. Vandembroucq, and M. L. Falk, Physical review letters 117, 045501 (2016).
  * Barbot  _et al._ (2018) A. Barbot, M. Lerbinger, A. Hernandez-Garcia, R. García-García, M. L. Falk, D. Vandembroucq, and S. Patinet, Physical Review E 97, 033001 (2018).
  * Patinet  _et al._ (2019) S. Patinet, A. Barbot, M. Lerbinger, D. Vandembroucq, and A. Lemaître, arXiv preprint arXiv:1906.09818 (2019).
  * Vezirov  _et al._ (2015) T. A. Vezirov, S. Gerloff, and S. H. Klapp, Soft Matter 11, 406 (2015).
  * Plimpton (1995) S. Plimpton, Journal of computational physics 117, 1 (1995).
  * Divoux  _et al._ (2011b) T. Divoux, C. Barentin, and S. Manneville, Soft Matter 7, 9335 (2011b).
  * Lamp  _et al._ (2022) K. Lamp, N. Küchler, and J. Horbach, The Journal of Chemical Physics 157, 034501 (2022).
  * Nishikawa  _et al._ (2022) Y. Nishikawa, M. Ozawa, A. Ikeda, P. Chaudhuri, and L. Berthier, Physical Review X 12, 021001 (2022).
  * Chacko  _et al._ (2019) R. N. Chacko, P. Sollich, and S. M. Fielding, Phys. Rev. Lett. 123, 108001 (2019).
  * Vasisht  _et al._ (2022) V. V. Vasisht, P. Chaudhuri, and K. Martens, Soft Matter 18, 6426 (2022).
  * Mandal and Sollich (2020) R. Mandal and P. Sollich, Phys. Rev. Lett. 125, 218001 (2020).
  * Parley  _et al._ (2022) J. T. Parley, S. Sastry, and P. Sollich, Phys. Rev. Lett. 128, 198001 (2022).

## Appendix

Figure A1: (Top) Maps of coarse-grained overlap function, $Q_{\ell}(X,Y;t)$
for $\sigma_{0}=0.79$, computed at $t=0.06$ (A), $87.6$ (B), $167.7$ (C),
$377.3$ (D) and $3369.7$ (E), using $\ell=4.94$. (Bottom) Distribution of
coarse-grained overlap at corresponding times. Figure A2: Maps of coarse-
grained overlap function, $Q_{\ell}(X,Y;t)$ for $\sigma_{0}=0.795$ computed at
$t=0.06$ (A), $87.6$ (B), $167.7$ (C), $377.3$ (D) and $3369.7$ (E), using
$\ell=4.94$. (Bottom) Distribution of coarse-grained overlap at corresponding
times.

Figure A3: Spatial maps of difference in activity (quantified via coarse-
grained local MSDs), computed for responses to $\sigma=0.79$ and
$\sigma=0.795$, using $\ell=4.94$, starting from same initial condition.