# Yielding transition of a two dimensional glass former under athermal cyclic
shear deformation

Himangsu Bhaumik  Jawaharlal Nehru Center for Advanced Scientific Research,
Jakkur Campus, Bengaluru 560064, India.  Giuseppe Foffi  Université Paris-
Saclay, CNRS, Laboratoire de Physique des Solides, 91405 Orsay, France
Srikanth Sastry<EMAIL_ADDRESS>Jawaharlal Nehru Center for Advanced
Scientific Research, Jakkur Campus, Bengaluru 560064, India.

###### Abstract

We study numerically the yielding transition of a two dimensional model glass
subjected to athermal quasi-static cyclic shear deformation, with the aim of
investigating the effect on the yielding behaviour of the degree of annealing,
which in turn depends on the preparation protocol. We find two distinct
regimes of annealing separated by a threshold energy. Poorly annealed glasses
progressively evolve towards the threshold energy as the strain amplitude is
increased towards the yielding value. Well annealed glasses with initial
energies below the threshold energy exhibit stable behaviour, with negligible
change in energy with increasing strain amplitude, till they yield.
Discontinuities in energy and stress at yielding increase with the degree of
annealing, consistently with recent results found in three dimensions. We
observe significant structural change with strain amplitude that closely
mirrors the changes in energy and stresses. We investigate groups of particles
that are involved in plastic rearrangements. We analyse the distributions of
avalanche sizes, of clusters of connected rearranging particles, and related
quantities, employing finite size scaling analysis. We verify previously
investigated relations between exponents characterising these distributions,
and a newly proposed relation between exponents describing avalanche and
cluster size distributions.

##  I Introduction

The response of amorphous solids to external deformation or applied stress is
of obvious importance in characterising their behaviour, both from a
fundamental and an applied point of view Bonn  _et al._ (2017); Nicolas  _et
al._ (2018). For small deformations, such response is characterized by
elastic-like behaviour, but for large enough deformation the response begins
to display plasticity and irreversible flow-like behaviour. Yielding
behaviour, and the transition to elasto-plastic flow, can be broadly
classified into two different categories. Upon deformation, the stress may
increase gradually and monotonically, till a steady value is reached for large
deformation. FoamsLauridsen  _et al._ (2002), emulsionsBécu  _et al._ (2006)
and colloidal suspensions Schall  _et al._ (2007) have been shown to display
such phenomenology. On the other hand, one may observe a stress overshoot and
a drop, accompanied by strain localization. This behaviour has been observed,
for example, in metallic glasses Greer  _et al._ (2013); Hufnagel  _et al._
(2016) and window glasses Guin and Wiederhorn (2004). The phenomenon of
yielding is central to understanding the mechanical properties in different
materials spanning several scales, ranging from nano-structured materials Zhao
_et al._ (2016) to large scale events like earthquakes, landslides, and
avalanches Chester and Chester (1998); Uhl  _et al._ (2015). In recent years,
an increasing number of studies have explored the fundamental aspects of the
yielding transition in amorphous solids through experimentsKeim and Arratia
(2013); Knowlton  _et al._ (2014); Hima Nagamanasa  _et al._ (2014); Di
Michele  _et al._ (2014); Keim  _et al._ (2020), numerical simulations Maloney
and Lemaître (2006); Karmakar  _et al._ (2010); Jaiswal  _et al._ (2016);
Fiocco  _et al._ (2013); Shrivastav  _et al._ (2016); Regev  _et al._ (2015);
Leishangthem  _et al._ (2017); Parisi  _et al._ (2017); Parmar  _et al._
(2019); Barbot  _et al._ (2020); Schinasi-Lemberg and Regev (2020) and
theoretical models Urbani and Zamponi (2017); Popović  _et al._ (2018); Barlow
_et al._ (2020); Liu  _et al._ (2020); Sastry (2021). Differently from
crystals Sethna  _et al._ (2017), amorphous solid plasticity is not driven by
specific defects and thus the correlation of plastic rearrangements and
specific structural motifs has been a subject of several investigations, as
well as the characterisation of structural changes that arise from plastic
events Richard  _et al._ (2020); Bonfanti  _et al._ (2019); Parmar  _et al._
(2019); Denisov  _et al._ (2015); Vasisht  _et al._ (2020); Mitra  _et al._
(2021). The distribution of sizes of events – avalanches – involved in plastic
rearrangements Dahmen  _et al._ (2009); Lin  _et al._ (2014); Regev  _et al._
(2015); Leishangthem  _et al._ (2017); Oyama  _et al._ (2020) have also been
widely investigated, which display characteristic power law forms whose
origins have been explored through elasto-plastic models and mean-field
theories Lin  _et al._ (2014); Jagla (2015); Liu  _et al._ (2016); Bouchbinder
_et al._ (2007); Dahmen  _et al._ (2011); Franz and Spigler (2017).

Recent computer simulations, performed largely using the athermal quasi-static
shear deformation (AQS) implementing both cyclic and uniform shear, support
the idea that the nature of the yielding transition is highly sensitive on the
initial degree of annealing of the amorphous Ozawa  _et al._ (2018); Bhaumik
_et al._ (2021a); Yeh  _et al._ (2020); Ozawa  _et al._ (2020). For uniform
shear, a threshold degree of annealing that separates the gradual and
sharp/discontinuous yielding behaviours (termed in this context as ductile and
brittle behaviour) has been argued to correspond to a random critical point
Ozawa  _et al._ (2018, 2020). Investigations of yielding under cyclic shear
reveal that, for a well annealed glass (above a threshold degree of annealing
or below a threshold energy), the drop in the maximum stress over a cycle is
discontinuous at the yielding transition, and the magnitude of the stress jump
increases with an increase in the degree of annealing. On the other hand, for
poorly annealed glasses, irrespective of the initial energy, cyclic shear
progressively anneals the glasses as the strain amplitude increases. In all
such cases, independently of the initial (low) degree of annealing, the
glasses approach a common energy – the threshold energy – at the yielding
point Bhaumik  _et al._ (2021a); Yeh  _et al._ (2020). The energy and stress
jumps observed for initially poorly annealed glasses is small Parmar  _et al._
(2019); Bhaumik  _et al._ (2021a). In previous work, we have investigated such
behaviour both in a strong glass (silica) and a fragile glass (the Kob-
Andersen binary Lennard-Jones mixture) in three dimensions.

Here, we investigate the corresponding properties of yielding under cyclic
shear for a model two-dimensional glass. In particular, we explore whether the
yielding transition displays the same kind of dependence on annealing as the
three dimensional glasses. In addition to studying the behaviour of energy and
stress, we also characterise the mechanically induced annealing and yielding
by analysing the structural changes involved. Above yielding, strain has been
observed to become localised in shear bands under cyclic shear Radhakrishnan
and Fielding (2016); Fiocco  _et al._ (2013); Parmar  _et al._ (2019) as well
as for annealed glasses under uniform shear Shi and Falk (2006); Shi  _et al._
(2007a); Ozawa  _et al._ (2018). We investigate strain localisation in the
model we study, focusing attention on the nature of structural change within
the shear bands, compared to regions outside the shear band. Finally, we
investigate the statistics of avalanches Dahmen  _et al._ (2009); Lin  _et
al._ (2014); Regev  _et al._ (2015); Leishangthem  _et al._ (2017), i.e. the
collection of particles undergoing plastic deformation, both below and above
yielding. In previous work, the nature of rearrangements under cyclic shear
has been investigated through the statistics of clusters of connected
particles that are displaced beyond a threshold value Leishangthem  _et al._
(2017), for a three dimensional glass. The statistics of sizes of such
clusters is qualitatively different below and above yielding. In particular,
the characteristic size of such clusters was shown Leishangthem  _et al._
(2017) to remain finite below the yielding transition, whereas above, their
sizes scale with the system size. These observations are in contrast with
simulations of a two-dimensional model Regev  _et al._ (2015) wherein the
sizes of avalanches diverges as the yielding transition is approached from
below. We investigate whether such a difference arises from dimensionality.
Further, we perform a careful investigation of avalanches above yielding, by
distinguishing the number of all particles participating in a plastic
rearrangement, for which we use the term avalanche, vs. spatially connected
clusters of rearranging particles. We further investigate the distribution of
the number of clusters that participate in an avalanche, to understand the
clustering properties of avalanches recently discussed in Le Priol  _et al._
(2021). These analyses confirm earlier conclusions for three dimensional
glasses Leishangthem  _et al._ (2017), and provide additional insights into
the nature of avalanches in the two dimensional model we investigate. In
particular, we verify a relation between exponents characterising the
distribution of avalanches and clusters, proposed in an accompanying paper
Bhaumik  _et al._ (2021b) and verified for three dimensional glasses.

The rest of the paper is organized as follows: In Section II we describe the
model potential (II.A), methods of initial glass preparation (II.B) and
protocol of cyclic shear that we use in this study (II.C). In Sec. III we
present the results starting with the yielding transition (III.A), followed by
results concerning structural change across the yielding transition (III.B)
and strain localization (III.C). Next, we present a detailed statistical
investigation of avalanches, analyzing various quantities related to plastic
rearrangements (III.D). A rigorous finite size scaling analysis will be
provided in support of various exponents and scaling relation we describe
(III.E). Finally, Sec. IV contains a summary and discussion of our results.
Various supporting results are presented in the Appendix.

##  II Simulation Methods

###  II.1 System potentials

We perform molecular dynamics simulations of a binary Lennard-Jones mixture in
two dimensions, introduced in Ref. Lancon  _et al._ (1984), that has been
extensively used to investigate the mechanical properties of amorphous solids
Widom  _et al._ (1987); Shi and Falk (2005); Shi  _et al._ (2007b); Barbot
_et al._ (2018, 2020). In this model, the two types of atoms interact via. the
Lennard-Jones inter-atomic potential, with a cut-off, given by,

$U_{ij}=\left\\{\begin{array}[]{lr}4\epsilon_{ij}[(\frac{\sigma_{ij}}{r_{ij}})^{12}-(\frac{\sigma_{ij}}{r_{ij}})^{6}]+A,&r_{ij}<R_{in}\\\
\sum_{k=0}^{4}C_{k}(r_{ij}-R_{in})^{k},&R_{in}<r_{ij}<R_{\rm{cut}}\\\
0,&r_{ij}>R_{cut}\end{array}\right.$ (1)

where the constants (A, C’s) can be found in Ref. Barbot  _et al._ (2018). A
plot of the interaction potentials is shown in Fig. 1. The density of the
system is kept constant and equals $1.02$. The reduced units, used throughout,
are expressed in terms of the mass $m=m_{S}=m_{L}$, and the energy scale
$\epsilon_{SL}$, length scale $\sigma_{SL}$ defining the inter species
interaction. For this model $\sigma_{SS}/\sigma_{SL}=2\sin(\pi/10)$,
$\sigma_{SS}/\sigma_{SL}=2\sin(\pi/5)$ and
$\epsilon_{SS}/\epsilon_{SL}=\epsilon_{LL}/\epsilon_{SL}=0.5$. Time is
expressed in units of $\tau=\sigma_{SL}\sqrt{m/\epsilon_{SL}}$ and temperature
is expressed in units of $\epsilon_{SL}/k_{B}$ where $k_{B}$ is the Boltzmann
constant. The composition of the mixture is given by the number ratio of large
($L$) and small ($S$) particles, $N_{L}:N_{S}=(1+\sqrt{5}):4$. The potential
is modified so that its value as well as the first two derivatives go
continuously to zero at a cutoff distance $R_{cut}=2.5\sigma$, and be twice
continuously differentiable at $R_{in}=2\sigma$.

Figure 1:  The interaction potential (Eq. (1)) used in this study, for the
different pairs ($SS$,$SL$ and $LL$) of particles.

###  II.2 Initial glass preparation

We perform constant temperature and volume (NVT) molecular dynamics
simulations using the Nosé-Hoover thermostat with an integration time step
$0.001$. Two different temperatures, $T=2.98$ and $T=0.35$, are considered.
The system was equilibrated in the liquid state for $10^{4}\tau_{\alpha}$ for
$T=2.98$ and $100\tau_{\alpha}$ for $T=0.35$, where $\tau_{\alpha}$ is the
structural relaxation time. We consider different system sizes, i.e.
$N=1024,2500,4900,10000,22500$. We employ two different protocols to prepare
glasses in four different states. The first two sets of glasses are obtained
from the well equilibrated high-temperature liquid (HTL, T=2.98) and the
supercooled liquid (ESL, $T=0.35$) configurations, by performing a local
energy minimization, to obtain the corresponding inherent structures. Two sets
of well annealed glasses with lower energies are obtained by the the finite
temperature and shear rate annealing protocol described in Das  _et al._
(2018). Following this procedure, we subject equilibrated configurations at
$T=0.35$ to cyclic shear deformation following the SLLOD equations of motion
Evans and Morriss (1984) with a shear rate $\dot{\gamma}=10^{-3}$, strain
amplitude $\gamma_{max}=0.035$ and simulation temperature $T=0.3$. The
energies of the inherent structures obtained at the end of each cycle of shear
decrease with cycles, as shown in Fig. A1 of the Appendix, which also shows
the inherent structure energy as a function of temperature for $T=2.98$ to
$T=0.35$ obtained from MD simulations. We obtain two different sets of samples
of such glasses whose energies on average are equal to $E_{IS}=-2.41$ and
$-2.45$, which we refer to as WAL1 and WAL2, respectively. These energies,
based on extrapolation of the inherent structure energies from the MD
simulations, corresponds to $T=0.271$ and $T=0.241$ respectively. We analyze
$10-12$ configurations for HTL glass for different system sizes and $3$
configurations for each of ESL, WAL1 and WAL2.

###  II.3 Athermal cyclic shear

The different initial configurations obtained from the protocols described
above, are subjected to the athermal quasi-static (AQS) shearing protocol
which involves two steps: (i) an affine transformation that increases the
strain by a small amount $d\gamma$ in the $xy$-direction (with coordinate
transformations $x^{\prime}\to x+d\gamma~{}y$, $y^{\prime}\to y$) and (ii) an
energy minimization. For a given strain amplitude $\gamma_{max}$, the strain
$\gamma$ is varied cyclically as :
$0\to\gamma_{\rm{max}}\to-\gamma_{\rm{max}}\to 0$ in each cycle. The
accumulated strain $\gamma_{acc}=\sum_{i}|d\gamma|$ serves as an effective
time variable in the analysis. The accumulated strain at $\gamma=0$ at the end
of each cycle is given by $\gamma_{acc}=4\times\gamma_{max}\times N_{cyc}$,
$N_{cyc}$ being the number of performed cycles. Repeating the deformation
cycle for a fixed amplitude $\gamma_{\rm{max}}$, the shear deformed glasses
evolve through a transient into a steady state, in which their properties
either become invariant from cycle to cycle (termed an absorbing stated,
obtained for $\gamma_{\rm{max}}$ below the yield value) or in which their
properties fluctuate around a mean value, but the particles as a function of
$\gamma_{acc}$ exhibit diffusive behaviour (e. g. Leishangthem  _et al._
(2017)). The steady states are obtained and studied for different samples for
strain amplitudes ranging from $\gamma_{\rm{max}}=0.02$ to $0.12$. We employ
the conjugate-gradient algorithm for energy minimization and execute all the
numerical simulations in LAMMPS Plimpton (1995).

The variation of the energy, stress and mean squared displacements (MSD) vs.
accumulated difference in the steady state are shown in the Appendix, in Figs.
A3, A3, A4, and A5. The asymptotic values of the stroboscopic ($\gamma=0$)
energies $U$ are obtained by stretched exponential fits to the dependence on
$\gamma_{acc}$. In most cases, the extracted asymptotic values correspond to
values around which the energies fluctuate in the steady state. The maximum
$\sigma_{xy}^{max}$ is obtained by averaging the stress at $\gamma_{\rm{max}}$
over $100-200$ cycles in the steady state. From the mean squared displacements
of particles with respect to initial configurations in the steady state, we
calculate diffusion coefficients $D$ using the $\gamma_{acc}$ as the time
variable in the usual definition of $D$. Structural quantities, and the
statistics of avalanches are likewise computed in the steady state, averaging
over $100-200$ cycles.

##  III  Results

###  III.1 Yielding transition

Figure 2:  Steady state potential energies (a, d) ($U/N$) of stroboscopic
configurations, (b,e) maximum value of stress $\sigma_{xy}^{max}$ and (c,f)
diffusion coefficient $D$ against strain amplitude ($\gamma_{max}$). Upper
panels are for different types of glasses: HTL, ESL, WAL1, and WAL2 with
system size $N=10000$. Lower panels are for the HTL glass for differnet system
sizes $N=1024,2500,4900,10000$. Dashed vertical line in (c,f) indicates the
yield strain amplitude $\gamma_{y}=0.06$.

The steady state energies $U/N$ are shown in Figs. 2(a) and 2(d) as a function
of strain amplitude $\gamma_{max}$ for the four types of glasses as well as
for different system sizes for the HTL. The HTL glass displays significant
finite size effects. In particular, while the shape of the curves in Fig. 2(d)
are very similar for different system sizes, the yield strain amplitude
decreases with increasing system size. Below yielding, all HTL glasses (also
ESL) anneal to lower energies with increasing strain amplitude, and we find
that lower energies are obtained for larger system sizes for HTL.

As shown in Fig. 2(a), as the strain amplitude $\gamma_{max}$ increases, $U/N$
of HTL and ESL progressively decreases until $\gamma_{max}$ reaches the yield
amplitude $\gamma_{y}=0.06$. At the yield strain the energy of both HTL and
ESL glasses approach to the same energy state. Such an observation is
consistent with the results of cyclic shear explored before Fiocco  _et al._
(2013); Leishangthem  _et al._ (2017); Bhaumik  _et al._ (2021a). Well
annealed glasses (WAL1 and WAL2), on the other hand, display a different
behavior. Unlike poorly annealed glasses (HTL and ESL), well annealed glasses
remain stable for amplitudes of strain below the yield amplitude and retrace
the same sequences of energies over the successive cycles of deformation (see
Fig. A3 in the Appendix), displaying the same behaviour observed in three
dimensional systems Bhaumik  _et al._ (2021a); Yeh  _et al._ (2020). Thus, the
stroboscopic energies remain unaltered as the strain amplitude increases
towards the yield strain amplitude $\gamma_{y}=0.06$. At the yield point, for
well annealed glasses, we observe a finite jump in energy whose size increases
with the degree of annealing. In previous work for a three dimensional glass
Bhaumik  _et al._ (2021a), a finite jump in energy is found even for poorly
annealed glasses. Within the precision of sampling of $\gamma_{max}$ in this
work, we cannot draw the same conclusion. However, very clearly, for well
annealed samples, a finite jump in energy is observed, which appears to be
inconsistent with a recent investigation of a two dimensional glass Schinasi-
Lemberg and Regev (2020). On the other hand, for the three dimensional
Lennard-Jones glass studied in Ref. Bhaumik  _et al._ (2021a) (as also that
investigated in Ref. Yeh  _et al._ (2020)), the yield strain amplitude depends
appreciably on the degree of annealing, which we do not find in the present
case. The yield amplitude in all cases studied here is roughly the same. This
may be a consequence of the limited range of annealing we investigate here,
and requires further study to clarify. As in previous studies, however, all
glasses investigated have the same energies above yielding, regardless of
initial conditions.

Next we study the variation of maximum stress ($\sigma_{xy}^{max}$) with
strain amplitude for the different cases. As shown in Figs. 2(b) (different
glasses) and 2(e) (for different system sizes of HTL), the $\gamma_{max}$
dependence of $\sigma_{xy}^{max}$ below the yield point exhibits significant
dependence of the annealing of the glasses, and a jump in the case of the well
annealed glasses whose size increases with increasing degree of annealing.
Fig. 2(e) displays substantial finite size effects both in the location of the
yielding point and the value of the stress at yielding. As remarked earlier in
the context of energy change, we do not attempt here to resolve whether a
finite jump is present in the case of poorly annealed glasses. Similar to
energy values, $\sigma_{xy}^{max}$ also shows the same $\gamma_{max}$
dependence in post yield regime for all the glasses.

Figure 3:  (a) Fraction of small particles $f_{n}^{S}$ and (b) Fraction of
large particles $f_{n}^{L}$ with $n$ nearest neighbours against $T$ for
inherent structure configurations for $N=4900$.

We next consider the presence or absence of diffusive motion of particles to
identify the yielding transition for the four types of glasses. We find that
for $\gamma_{max}<\gamma_{y}$, the MSD values are constant, indicating the
absence of cumulative movement from one cycle to the next, whereas for
$\gamma_{max}>\gamma_{y}$, the particle motion is diffusive (see Fig. A5 of
the Appendix). We characterize such diffusive motion as a function of
accumulated strain (relative to a reference point in the steady state) with
$MSD=4D\Delta\gamma_{acc}$ where $D$ is the diffusion coefficient and
$\Delta\gamma_{acc}$ is the difference in $\gamma_{acc}$ between the
configurations used to compute the MSD. The values of $D$ as function of
$\gamma_{max}$ is presented in Fig. 2(c) for the the different glasses for
$N=10000$ and for different system sizes for HTL in Fig. 2(f). The diffusion
coefficients are zero below the yielding transition and exhibit a finite value
at the yielding transition, and an increase with $\gamma_{max}$ above. Fig.
2(f) illustrates the finite size effects present for HTL, with the diffusion
coefficients being bigger for larger systems.

Thus, the change in the $\gamma_{max}$ dependence of the energy, maximum
stress, and the diffusion coefficient, consistently mark the discontinuous
yielding transition, with the observed discontinuities of energy and stress
being clearly dependent on the degree of annealing, and the jump in the
diffusion coefficient exhibiting no such dependence.

###  III.2 Structural change across the yielding transition

Although some attention has been paid in the past to the nature of structural
changes associated with shear deformation and yielding in glasses Shi and Falk
(2006); Parmar  _et al._ (2019); Adhikari, Monoj and Sastry, Srikanth (2018);
Bhaumik  _et al._ (2021b); Kawasaki and Berthier (2016); Denisov  _et al._
(2015); Vasisht  _et al._ (2020); Vasisht and Del Gado (2020), such
investigations have been relatively rare. Also, some of the conclusions drawn
by such studies are at odds with each other. For example, it has been argued
that a sharp change in the symmetry associated with anisotropy is observed
when crossing the yielding transition in a colloidal glass subjected to cyclic
shear experimentally Denisov  _et al._ (2015), but no such anisotropy has been
observed in a computer simulation investigation Adhikari, Monoj and Sastry,
Srikanth (2018). It has been claimed that in a system of jammed repulsive non-
Brownian particles, the static structure does not change in the pre- and post-
yield regime Kawasaki and Berthier (2016). However, clear indications of
structural change have been observed in other simulation studies Shi and Falk
(2006); Parmar  _et al._ (2019); Vasisht  _et al._ (2020); Vasisht and Del
Gado (2020). In order to understand the nature of structural change in the
pre-yield and post-yield regimes for the model glass we study, we perform a
detailed analysis of structure.

Figure 4:  Fractions $f_{n}^{S}$ for $n=5,6$ and $f_{n}^{L}$ for $n=6,7$
against strain amplitude $\gamma_{max}$ for different types of glasses with
$N=10000$ in the steady state. Figure 5:  Average orientational order
parameter for (a-c) small particles ($\langle\psi_{n}^{S}\rangle$) and for
(d-f) large particles ($\langle\psi_{n}^{L}\rangle$) against strain amplitude
$\gamma_{max}$ for different types of glasses. Different panels are for
different values of $n=5,6,7$. System size $N=10000$.

Figure 6:  Average orientational order parameter
$\langle\psi_{6}\rangle=N^{-1}\sum_{i=0}^{N}\psi_{6}^{i}$ averaged over all
the particles against strain amplitude $\gamma_{max}$ for (a) different types
of glasses for a given system size $N=10000$ and (b) for the HTL glass with
different system sizes $N=1024,2500,4900,10000$.

We begin by analysing the structural change in the inherent structures
corresponding to equilibrium liquid configurations, to benchmark our
subsequent analysis of the sheared glasses. Since our glasses are composed of
two types of particles (large, ”L”, and small, ”S”), we investigate the short
range order around each type of particle by considering the number of
neighbors around each type of particle. Fig. 3 shows the fraction of large and
small particles with different numbers of neighbors in the first shell
(defined as the geometric neighbors of a particles in a generalised Delaunay
tessellation for configurations of bidisperse discs Sastry  _et al._ (1997);
Maiti  _et al._ (2013)) as a function of temperature, $f^{L}_{n}$,
$f^{S}_{n}$, where $L$ and $S$ refer to the type of particle and $n$ the
number of neighbors. The data clearly reveal that $f^{S}_{n=5}$, $f^{L}_{n=7}$
and $f^{L}_{n=8}$ grow as temperature is lowered, and all other fractions, in
particular $f^{S}_{n=6}$ and $f^{L}_{n=6}$, decrease. This indicates that the
preferred order is for small particles to be surrounded by $5$ neighbours, and
large particles to be surrounded by $7$ or $8$ neighbors. Local structure with
$6$ neighbours, for either large or small particles, corresponds to
unfavorable arrangements at low temperatures.

With these results as reference, we now consider the fractions $f^{L}_{n=7}$,
$f^{S}_{n=5}$, corresponding to favorable motifs, and $f^{L}_{n=6}$,
$f^{S}_{n=6}$, corresponding to unfavorable motifs, for the different glasses
we study, as a function of strain amplitude $\gamma_{max}$. As shown in Fig.
4, $f^{L}_{n=7}$ and $f^{S}_{n=5}$ are large for the well annealed glasses and
small for the poorly annealed glasses, and they remain constant for well
annealed glasses as $\gamma_{max}$ increases (for $\gamma_{max}<\gamma_{y}$,
whereas they increase with increasing $\gamma_{max}$ for the poorly annealed
glasses, indicating an increase in the fraction of favourable structural
motifs. In contrast, $f^{L}_{n=6}$ and $f^{S}_{n=6}$ show exactly the opposite
behaviour. Above the yield point, $f^{L}_{n=7}$ and $f^{S}_{n=5}$ decrease,
whereas $f^{L}_{n=6}$ and $f^{S}_{n=6}$ increase, and the values do not depend
on the initial annealing degree of the glasses.

Figure 7:  (a) Colour map of particle displacements for a zero strain
configuration from the steady state for strain amplitude $\gamma_{max}=0.07$.
Particles are coloured according to their mean square displacement (MSD)
computed over a strain cycle. Particles with $MSD>0.25\sigma^{2}_{SL}$ are
shown in red and the others according to the colour bar shown. (b) Slab-wise
averaged $MSD_{x}$ in the x-direction showing a Gaussian form with width
$w=8.86\sigma_{SL}$ within which the particles are highly mobile. (c) Slab-
wise averaged $\langle\psi_{6}\rangle_{x}$ (averaged over all particles)
(black dots) and slab-wise density (red dots) against the x-coordinate.
Horizontal dashed line indicates values of $\langle\psi_{6}\rangle=0.302$ for
the whole system and overall system density $\rho=1.02$. All the data are for
an HTL glass configuration with $N=10000$.

In order to capture the overall structural changes in a small number of order
parameters, we compute the two dimensional orientational order parameters

$\psi^{i}_{n}=\frac{1}{N_{i}}\sum_{j=1}^{N_{i}}\exp(in\theta_{ij})$ (2)

where $N_{i}$ is the number of nearest neighbours of a particle $i$ obtained
by Delaunay tessellation as mentioned before, and $\theta_{ij}$ is the angle
made by the vector from particle $i$ to its neighbor $j$ with the $x$-axis. We
compute the average values $\langle\psi^{L}_{n}\rangle$ and
$\langle\psi^{S}_{n}\rangle$ for the large and small particles, as averages
over the respective single particle values, for $n=5,6,7$, which are shown in
Fig. 5. We see that $\langle\psi^{L}_{n}\rangle$ and
$\langle\psi^{S}_{n}\rangle$ capture the same trends as the fractions
$f^{L}_{n}$, $f^{S}_{n}$, with the exception that
$\langle\psi^{L}_{n=6}\rangle$ exhibits roughly the same values for the well
annealed glasses as for the yielded glasses. Nevertheless,
$\langle\psi_{n=6}\rangle$, averaged over all the particles, shown in Fig. 6
is a good indicator of overall structural change, and correlates strikingly
well with the changes in the energy for all the different glasses and for
different system sizes. In investigating the spatial variation of structure in
the presence of strain localisation, therefore, we employ
$\langle\psi_{n=6}\rangle$ as the single indicator of structural change. Fig.
A1 (b) shows the temperature dependence of $\langle\psi_{n=6}\rangle$ from MD
simulations, for reference, for different system sizes.

###  III.3 Strain Localisation

The localisation of strain above the yielding point under cyclic shear has
recently been investigated in three dimensional glasses Parmar  _et al._
(2019); Mitra  _et al._ (2021). The shear bands, defined as regions with
greater mobility from one cycle of strain to the next, have been shown to have
comparatively higher energy as well as lower density. We thus consider whether
similar strain localisation is observed for the two dimensional system we
study, and what aspects of structural change it is accompanied by. In Fig.
7(a) we show a snap shot of the stroboscopic configuration for
$\gamma_{max}=0.07$, in the steady state, for the HTL glass. Particles are
coloured according to their mean square displacement(MSD) computed between two
successive stroboscopic configurations and all the particles that have
$MSD>0.25\sigma^{2}_{SL}$ are coloured in red. We see that the most mobile
particles are highly correlated in space and form system spanning band like
structures. In Fig. 7(b) we show the slab wise averaged (along x-direction)
mean square displacement ($MSD_{x}$), which clearly demonstrates the existence
of shear band, within which particles move substantially between successive
cycles. The $MSD_{x}$ profile can be described, to a good approximation, by a
Gaussian with width $w$ roughly equal to $8.86\sigma_{SL}$. We see that most
mobile particles are highly localized within the width of the shear band, with
$MSD_{x}\approx 0.25\sigma_{SL}^{2}$ or greater within the shear band. We next
investigate the structural change with the shear band. As defined in Eq. (2)
we compute $\psi_{n=6}$ for each particle. The slab wise averaged
$\langle\psi_{6}\rangle_{x}$ along with the slab-wise density are presented in
7(c). We see a reduction of structural order accompanied by a reduction of
density inside the shear band. Such an observation is consistent with
observations for the three dimension model glasses investigated earlier Parmar
_et al._ (2019); Mitra  _et al._ (2021).

###  III.4 Avalanches

Figure 8:  (a) Colour map of particle displacements during a plastic
rearrangement at $\gamma\approx 0.042$ during a strain cycle of amplitude
$\gamma_{max}=0.06$. Particles with the lowest displacement $\delta r$ are
shown in blue and those with the largest $\delta r$ are shown in red. Shown in
maroon are the active particles with $\delta r>0.25$ which form a spatially
correlated cluster with a cluster size $s=1359$ at the center. Distribution of
(b) avalanche size $P(S)$ and (c) cluster size $P(s)$ for different
$\gamma_{max}$ for system size $N=10000$. For both the cases two distinct sets
can be observed, one corresponding to $\gamma_{max}\geq\gamma_{y}$ exhibiting
a power law with exponent $\tau_{a}=1.1$ and $\tau_{c}=1.66$, for avalanche
and cluster size distributions respectively. Inset of (c) shows the
distribution of the number of clusters in an avalanche, $P(n_{cl})$, with a
power law with exponent $\kappa\approx 1.18$.

Next we study the statistics of avalanches accompanying plastic rearrangements
of particles during strain deformation of the system. Such events are
accompanied by a discontinuous drop in stress and in energy. In our analysis
here, we consider only the HTL glass. We consider steady state stroboscopic
configurations and consider plastic rearrangement events in the first quadrant
of the cycle, with $\gamma$ varying from $0\to\gamma_{max}$. To obtain the
configurations just before and just after the plastic events precisely, we
vary $\gamma$ with $d\gamma=10^{-6}$. With a given $\gamma_{max}$, for each
sample we consider $20$ cycles to collect information on non-affine
displacements corresponding to each drop event, identified following
procedures in Refs. Leishangthem  _et al._ (2017); Parmar  _et al._ (2019), by
comparing the change in energy with the expected change if the change is
elastic. Even in the presence of a plastic rearrangement, the long range
elastic strain field leads to a continuous range of single particle
displacements that are power law distributed Eshelby (1957); Dyre (1999), and
one needs to identify particles that are displaced as part of the plastic
rearrangement vs. those that are elastically displaced in response. In Refs.
Schröder  _et al._ (2000); Fiocco  _et al._ (2014); Leishangthem  _et al._
(2017), it was observed that the distribution of single particle displacements
exhibits a power law regime, followed by an exponential cutoff. A value of
$\delta r$ that separates the two regimes was chosen as a cutoff, and
particles with displacements larger than the cutoff were identified as being
part of the plastic core of the event. It was, however, also noted in
Leishangthem  _et al._ (2017) that the power law regime follows the expected
scaling ($p(\delta r)\sim(\delta r)^{-(2D-1)/(D-1)}$) only for small
amplitudes of shear. Thus, we follow the procedure of identifying a cutoff
which falls within the exponential tail, but verifying that the distribution
of avalanches and clusters is not sensitive to the choice. We also
investigated a similar prescription in Ref. Salerno and Robbins (2013) that
employs the same reasoning but employs local deviatoric strain instead of
displacements. We find, however, that the resulting avalanche and cluster
distributions are sensitive to the choice of the cutoff and therefore, pending
further investigation, we do not employ the deviatoric strain method in this
work. In Fig. A7 of the Appendix, we show the distributions of particle
displacements $\delta r$, along with different choice of the cutoff. Fig. A7
also shows that the power law part of the cluster size distribution is
insensitive to a choice of cutoff $\delta r_{c}$ from $0.15$ to $0.3$. In what
follows, we employ the choice $\delta r_{c}=0.25$ to identify what we term
active particles.

In Fig. 8(a) we show a colour map of particle displacements for a plastic
rearrangement. Imposing a cut-off value $\delta r_{c}=0.25$ on the
displacement field we identify all the active particles which are coloured in
maroon. It can be observed that such active particles are highly correlated in
space and form a cluster. In general, the set of all active particles (whose
number yields the avalanche size) may be composed of disconnected clusters of
rearranging particles. We identify such clusters, as well as the number of
such clusters in a given avalanche, and investigate their properties. If an
avalanche of size $(S)$ (total number of active particles) consists $n_{cl}$
number of clusters and $i$th cluster has size $s_{i}$, then one has

$S=\sum_{i=0}^{n_{cl}}s_{i}.$ (3)

We compute the cluster statistics of the active particles, considering two
active particles (or type A, B) to be in the same cluster if they are
separated by less than a distance $r_{AB}^{min}$, the first minimum of the
relevant partial pair correlation function ($r_{AB}^{min}=1.03,1.40,1.70$ for
$SS$,$LL$, and $SL$ pairs respectively). The distribution of the investigated
quantities, namely the avalanche size ($S$), cluster size ($s$), and the
number of clusters in an avalanche, $n_{cl}$, are expected to exhibit power
law regimes Le Priol  _et al._ (2021), characterised by exponents identified
below:

$P(S)\sim S^{-\tau_{a}},\ \ \ P(s)\sim s^{-\tau_{c}},\ \ \ P(n_{cl})\sim
n_{cl}^{-\kappa}.$ (4)

We study the size distribution of avalanches and clusters for various strain
amplitudes as shown in Fig. 8(b) and 8(c), respectively. The distribution of
the number of clusters is shown in the inset of Fig. 8(c). For all the three
quantities, one can identify two distinct sets of distributions that
correspond respectively to strain amplitudes $\gamma_{max}$ below, and above,
the yield strain $\gamma_{y}$ as also observed in a previous study for a three
dimensional glass for cluster sizes Leishangthem  _et al._ (2017). For
$\gamma_{max}>\gamma_{y}$ in particular, the distributions exhibit clear
power-law regimes, with exponents $\tau_{a}=1.1\pm 0.1$, $\tau_{c}=1.66\pm
0.06$ and $\kappa=1.18\pm 0.07$. The errors indicated represent the range of
variation of the exponent values obtained by fitting the different data sets
and ranges of sizes. Results are found to be similar over a range of choice of
the cutoff $\delta r_{c}$ as already mentioned (and shown in the Appendix,
Fig. A7). We discuss the relationship between $\tau_{a}$, $\tau_{c}$ and
$\kappa$, following the arguments in Ref. Le Priol  _et al._ (2021) for models
of crack propagation, in the presence of long range interactionsLaurson  _et
al._ (2010), and as discussed for three dimensional glasses in Ref. Bhaumik
_et al._ (2021b).

Considering an avalanche of size $S$, we represent the number of the size of
the clusters present as $n(s|S)$. Since the sum of the sizes of the clusters
must equal the size of the avalanche, we have, by definition,

$\int_{1}^{S}~{}s~{}n(s|S)~{}ds=S.$ (5)

Likewise, by definition, the total number of clusters present is given by

$\int_{1}^{S}~{}n(s|S)~{}ds=n_{cl}(S).$ (6)

The distribution of cluster sizes for a given $S$ exhibits a power law form up
to the maximum size of $S$, as also noted in Le Priol  _et al._ (2021) (The
distributions, however, are noisy within the statistics we have and thus we do
not rely on them directly in the discussion below, but instead on average
quantities as described later.) We thus assume $n(s|S)=As^{\tau}$, and the
condition Eq. 5 requires $A=(2-\tau)S^{\tau-1}$ (for large $S$). Note that we
do not assume $\tau=\tau_{c}$ for reasons that will be discussed further
later. With this choice, and assuming that $\tau>1$, we obtain for the number
of clusters

$n_{cl}(S)\sim S^{\tau-1}.$ (7)

Defining the exponent $\gamma_{ns}$ by $\langle n_{cl}\rangle_{S}\sim
S^{\gamma_{ns}}$, we have $\gamma_{ns}=\tau-1$, or

$\tau=\gamma_{ns}+1$ (8)

We note that the above relations also imply that the mean cluster size

$\langle s\rangle_{S}\sim S^{2-\tau}$ (9)

as along noted in Le Priol  _et al._ (2021). Fig. 9(a) shows that indeed, the
mean cluster size scales for large system sizes with exponent $2-\tau=0.45$,
corresponding to $\tau=1.55$.

Normalising $n(s|S)$ we obtain the distribution $P(s|S)\sim s^{-\tau}$ and
with $P(S)\sim S^{-\tau_{a}}$, we obtain, the full distribution of $s$ to be

$P(s)=\int_{s}^{\infty}P(s|S)P(S)dS\sim s^{-(\tau+(\tau_{a}-1))}.$ (10)

Thus, the exponent describing the cluster sizes is

$\tau_{c}=\tau+(\tau_{a}-1)$ (11)

which, proposed in Bhaumik  _et al._ (2021b) has not, to our knowledge,
previously been discussed. We discuss the validity of this relation in the
next section.

We next consider the distribution of the number of clusters. Since the number
of clusters $n_{cl}$ for a given avalanche size $S$ are distributed with an
average value $\langle n_{cl}\rangle_{S}\sim S^{\gamma_{ns}}$, the conditional
probability $P(n_{cl}|S)$ should have the functional form $P(n_{cl}|S)\sim
S^{-\gamma_{ns}}g(n_{cl}/S^{\gamma_{ns}})$, which results in $\langle
n_{cl}\rangle_{S}\sim S^{\gamma_{ns}}$, through $\langle
n_{cl}\rangle_{S}=\int n_{cl}P(n_{cl}|S)dn_{cl}$. The scaling function $g$,
shown in Fig. 9, is peaked around the mean value $\langle n_{cl}\rangle_{S}$,
as also seen in Le Priol  _et al._ (2021). The full distribution of $n_{cl}$
can be obtained as

$\displaystyle P(n_{cl})$ $\displaystyle=$ $\displaystyle\int
P(n_{cl}|S)P(S)dS$ (12) $\displaystyle\sim$ $\displaystyle
n_{cl}^{-(1+(\tau_{a}-1)/\gamma_{ns})}.$

Using Eq. 4, 8, and 12 we obtain the scaling relation

$\kappa=1+\frac{\tau_{a}-1}{\tau-1}$ (13)

This relation is well satisfied with the measured values of the exponents
described above, namely $\tau_{a}=1.1\pm 0.1$, $\tau=1.55\pm 0.06$ and
$\kappa=1.18\pm 0.07$. To further characterize the scaling behaviour we
perform a system size analysis in the next section.

Figure 9:  (a) The mean cluster size $\langle s\rangle_{S}$ for a given
avalanche size $S$ for different system sizes at $\gamma_{max}=0.09$. For the
small $S$ range, $\langle s\rangle_{S}$ increases as a power law with exponent
converging to $0.45$ for the largest system size. Legends indicate the system
size and the exponent value. (b) Scaled conditional probability distribution
of the number of clusters for given avalanche size $S$ for a system of size
$N=10000$. Legends indicate the range of $S$ values and the mean number of
clusters.

###  III.5 Finite size scaling

Figure 10:  (a) Moments of the cluster size $\langle s^{m}\rangle$ against
system size $N$ for different values of $m$. Lines are the fitted straight
lines obtained from linear regression. (b) Moment exponent $\sigma(m)$ and its
derivative with respect to $m$. (c) Scaled distribution of cluster size
against scaled cluster size for different system sizes for
$\gamma_{max}=0.09$. Inset shows the unscaled data.

Figure 11:  (a) Scaled avalanche size distributions against scaled avalanche
size for different system sizes with $\tau_{a}=1.1$ and
$\beta_{a}/\tau_{a}=1.0$. Inset shows the unscaled data which can be well
described by the power-law $P(S)\sim S^{-1.1}$. (b) Scaled cluster number
distributions for different system sizes. Best data collapse is obtained using
$\beta_{n}/\kappa=0.5$ and $\kappa=1.18$. Inset shows the unscaled data with
increasing cutoff as $N$ increases. The straight line through the data points
indicates the power law $P(n_{cl})\sim n_{cl}^{-1.18}$. (c) Scaled average
number of clusters $\langle n_{cl}\rangle_{S}N^{\gamma_{ns}\delta}$ against
scaled avalanche size $S/N^{\delta}$ for different system sizes. From fitting
of the data for the largest system sizes we obtain $\gamma_{ns}=0.55\pm 0.05$,
as indicated.

Figure 12:  (a) Average cluster size $\langle s\rangle$ against system size
for different $\gamma_{max}$. (b) Average cluster size against strain
amplitude $\gamma_{max}$ for different system sizes. Inset: scaled cluster
size against $\gamma_{max}$ to show the $\langle s\rangle\sim N^{1/3}$
scaling.

To investigate the system size dependence, we collect the statistics of
avalanche sizes, cluster sizes and the number of clusters, for a range of
system sizes $N=1024,2500,4900,10000,22500$ for the HTL glass. In the
following we will describe in detail the finite size scaling analysis of
various quantities. We begin by analysing the distribution of cluster sizes
assuming a finite size scaling form

$P(s,N)=N^{-\beta_{c}}f\left[\frac{s}{N^{\beta_{c}/\tau_{c}}}\right]$ (14)

where the scaling function has the characteristics $f(x)\to x^{-\tau_{c}}$ for
$x\to 0$, so that the power-law behaviour retrieved as $P(s)\sim
s^{-\tau_{c}}$ and $f(x)$ decreases to zero rapidly when $x\to 1$. To extract
the value of $\beta_{c}$ and $\tau_{c}$ we calculate the various moments of
the distribution defined as:

$\displaystyle\langle s^{m}\rangle$ $\displaystyle=$
$\displaystyle\int_{0}^{\infty}s^{m}P(s,N)ds=\int_{0}^{\infty}s^{m}N^{-\beta_{c}}f(s/N^{\beta_{c}/\tau_{c}})ds$
(15) $\displaystyle=$ $\displaystyle
N^{\beta_{c}/\tau_{c}(m+1-\tau_{c})}\int_{0}^{\infty}z^{m}f(z)dz$

where, in the last step we replace the variable $z=s/N^{\beta_{c}/\tau_{c}}$.
Since the integration $\int z^{m}dz$ is a constant the moment $\langle
s^{m}\rangle$ should vary with system size as

$\langle s^{m}\rangle\sim N^{\frac{\beta_{c}}{\tau_{c}}(m+1-\tau_{c})}\sim
N^{\alpha(m)},$ (16)

where the moment exponent $\alpha(m)=\frac{\beta_{c}}{\tau_{c}}(m+1-\tau_{c})$
is evaluated for different values of $m$ ranging from $0$ to $4$ in intervals
of $0.01$. For some integer values of $m$ the data of $\alpha(m)$ against $N$
are presented in Fig.10(a) for different strain amplitudes, larger than the
yield value. For $P(s)$ to obey finite size scaling (FSS) for a given
$\gamma_{max}$, the moment exponent $\alpha(m)$ should have a constant gap
between two successive values of $m$, i.e.,
$\alpha(m+1)-\alpha(m)=\beta_{c}/\tau_{c}$ and a derivative of $\alpha(m)$
with respect to $m$ should converge to $\beta_{c}/\tau_{c}$ for higher values
of $m$. In Fig. 10(b) we present the variation of $\alpha(m)$ and its
derivative $\partial\alpha(m)/\partial m$(obtained from the central
difference) as a function of $m$ for several strain amplitudes. These data
reveal that $\partial\alpha(m)/\partial m$ converges satisfactorily for all
$\gamma_{max}$ considered here in the post yield regime providing the value of
$\beta_{c}/\tau_{c}=1$. By a linear fit of $\sigma(m)$ for large $m$ we
determine $\beta_{c}/\tau_{c}=1.02\pm 0.02$ and $\beta_{c}=1.66\pm 0.08$ which
yield $\tau_{c}=1.69\pm 0.07$ consistent with the direct measurement described
earlier. In Fig.10(c) we present the scaling plot of the cluster size
distribution that demonstrates the excellent quality of the data collapse for
different system sizes.

Similar to cluster sizes $s$, we perform finite size scaling for the avalanche
sizes $S$ and for the number of clusters $n_{cl}$ by writing

$P(S,N)=N^{-\beta_{a}}f_{a}\left[\frac{S}{N^{\beta_{a}/\tau_{a}}}\right]$ (17)

and

$P(n_{cl},N)=N^{-\beta_{n}}f_{n}\left[\frac{n_{cl}}{N^{\beta_{n}/\kappa}}\right]$
(18)

where $f_{a}$ and $f_{n}$ are the scaling function and $\beta_{a}$ and
$\beta_{n}$ are the associated scaling exponents. In Figs. 11(a) and 11(b) we
present the scaled distribution for $S$ and $n_{cl}$ respectively. The best
data collapse obtained for $\beta_{a}/\tau_{a}=1$ and $\tau_{a}=1.1$ for
avalanche size $S$, whereas for $n_{cl}$ we use $\beta_{n}/\kappa=0.5$ and
$\kappa=1.18$. Given the limited range of the power laws, we find that the
moment analysis along the lines performed for cluster sizes is less
satisfactory, and we thus do not use such a procedure in these cases.

We next investigate the system size scaling of average number of clusters
$\langle n_{cl}\rangle_{S}$ for a given avalanche size $S$. Keeping in mind
Eq. 7 and Eq. 8, we assume the finite size scaling form as

$\langle n_{cl}\rangle_{S}=N^{\gamma_{ns}\delta}g_{n}(S/N^{\delta})$ (19)

where the scaling function $g_{n}(x)$ grows as $x^{\gamma_{ns}}$ for small
$x$. In Fig. 11(c) we show the scaled data of average cluster number according
to Eq. (19). A reasonable data collapse obtained for different system sizes
using $\gamma_{ns}=0.55$ and $\delta=1$. The $\gamma_{ns}$ value is consistent
with the estimated exponent $\tau=1.55$. Note that the exponent value we
estimate, of $\tau_{c}\approx 1.66$, with $\gamma_{ns}=0.55$ requires Eq. 11
to be true, and thus, the we conclude that the exponents we estimate

$\tau=1.55,\ \ \gamma_{ns}=0.55,\ \ \tau_{c}=1.66,\ \ \tau_{a}=1.1,\ \
\kappa=1.18$ (20)

are all satisfactorily and consistently explained by the analysis we presented
in the previous section.

Next we consider the variation of average cluster size $\langle s\rangle$ with
system size across the yield strain amplitude. In Fig. 12(a) we show $\langle
s\rangle$ against $N$, which exhibits different scaling behaviour in the post-
and pre-yield regimes: $\langle s\rangle\sim N^{1/3}$ when
$\gamma_{max}>\gamma_{y}$ and the $N$ dependence is negligible below
$\gamma_{y}$. Fig. 12(b) shows the data of $\langle s\rangle$ against strain
amplitude for different $N$. It can be seen that $\langle s\rangle$ does not
depend on system size for $\gamma_{max}<\gamma_{y}=0.06$, whereas it exhibits
a $\langle s\rangle\sim N^{1/3}$ dependence for $\gamma_{max}>\gamma_{y}$ as
described before. The scaled averaged size is plotted against $\gamma_{max}$
in inset of Fig. 12(b) in order to verify the scaling behaviour of $\langle
s\rangle$ in post-yield regime. The observed data collapse confirms the
$N^{1/3}$ scaling, which is consistent with the first moment obtained, as
shown in Fig. 10(b).

Finally, we study the statistics of energy drops during a plastic
rearrangement, which has been investigated in the past as a quantifier of
avalanche size. The plastic component of energy drop ($\Delta U_{pl}$) (see
Appendix Fig. A7) follows power-law scaling $P(\Delta U_{pl})\sim\Delta
U_{pl}^{-1.25}$ in agreement with observations in Ref. Salerno and Robbins
(2013); Leishangthem  _et al._ (2017) for three dimensional glasses. The
variation of $\langle\Delta U_{pl}\rangle$ with system size for
$\gamma_{max}>\gamma_{y}$ displays similar scaling behaviour, $\Delta
U_{pl}\sim N^{1/3}$, as observed for the mean cluster size.

##  IV Summary and Discussion

To summarize, we have studied the yielding behaviour of a two dimensional
glass former under athermal quasi-static cyclic deformation, focusing in
particular on the role of annealing of the glasses on the yielding behaviour.
Upon repeated cycles of shear, the glasses reach steady states, which are
either invariant from one cycle to the other (absorbing states) or whose
properties fluctuate around a mean value, with the particles exhibiting
diffusive behaviour from cycle to cycle. The former corresponds to strain
amplitudes below yielding, whereas the latter to strain amplitudes above the
yield strain amplitude. Our results demonstrate that the nature of yielding
changes with the degree of annealing, consistently with results for three
dimensional glasses, and that well annealed glasses exhibit a strongly
discontinuous yielding transition, with the discontinuous changes in energy,
stress and other properties increasing with the degree of annealing. We
characterise the structural changes during shear deformation, and show that
substantial changes in structure occur both above and below yielding. In
particular, our results show that structural heterogeneities accompany the
formation of shear bands above yielding. We have performed a detailed
investigation of the statistics of avalanches, including an analysis of the
decomposition of avalanches into clusters. We have analysed the relationship
between exponents characterising the distribution of avalanches, cluster sizes
and the number of clusters, and verify predictions from recent investigations
Le Priol  _et al._ (2021) of these exponents in the context of models for
crack propagation with long range interactions. In addition, we propose (in an
accompanying paper Bhaumik  _et al._ (2021b) and here) a new relationship
between the exponents characterising avalanches and clusters and demonstrate
the relationship to be valid. It will be interesting for this relationship to
be investigated in other glasses. These analyses have been performed using
detailed finite size scaling. We further verify for the two dimensional system
studied that avalanche sizes do not grow with system size below yielding,
whereas they do above, mirroring results for three dimensional glasses
Leishangthem  _et al._ (2017).

## Appendix

Figure A1:  (a) Variation of potential energy of the inherent structures
obtained from NVT MD simulation against temperature for different system
sizes. Cross marks indicate the energy of the four glasses considered in this
study. The dashed extrapolated line, obtained from fitting the data for
$N=10000$ with $E_{IS}=a/T^{b}+c$ in the $T<1$ region, maps the $E_{IS}\ \ vs\
\ T$ in the low temperature regime. (b) Average orientational order parameter
$\langle\psi_{6}\rangle$ of inherent structure configurations against
temperature for different system sizes. (c) Potential energy of the inherent
structure of cyclically sheared stroboscopic configurations as a function of
the number of cycles. Dashed horizontal lines represent the energies of the
four glasses which have been investigated in this work. From the extrapolation
of $E_{IS}\ \ vs\ \ T$ data we find the corresponding temperatures of WAL1 and
WAL2 to be $T=0.271$ and $0.241$, respectively.

###  A-1 Glass preparation through finite rate cyclic shear

We explore the yielding transition as a function of the degree of annealing of
glasses subjected to cyclic shear deformation. To generate poorly annealed
glasses (HTL and ESL) we use conventional constant volume and temperature
(NVT) molecular dynamics, while for well annealed glasses (WAL1 and WAL2)
below the energy $-2.40$, we use a finite temperature, finite shear rate,
cyclic shear protocol which we describe in this section. Taking equilibrated
configurations at $T=0.35$ trajectory, the system evolved using the SLLOD
equations of motion under repetitive cyclic shear for $T=0.3$, shear rate
$\dot{\gamma}=10^{-3}$ and amplitude $\gamma_{max}=0.035$. The inherent
structure energy of cyclically sheared stroboscopic configurations are shown
in Fig. A1. With increasing number of cycles, the system explores deeper and
deeper energy configurations. We have collected two sets of configurations
from these well annealed liquid configurations, with energy corresponding to
$-2.41$ (WAL1) and $-2.45$ (WAL2).

###  A-2 Evolution of potential energy under oscillatory deformation

In Fig. A3, we present the evolution of the stroboscopic potential energy
configurations with accumulated strain for different cases (HTL for different
system sizes, along with ESL, WAL1, and WAL2 for $N=10000$). The results
display a trend similar to what have been observed in previous work Fiocco
_et al._ (2013); Bhaumik  _et al._ (2021a) for three dimensional glasses. For
HTL, we have considered $10-12$ samples for each system size for all the
strain amplitudes, whereas for ESL, WAL1, and WAL2 we consider $3$ samples.
Note that the evolution of energy in the post yield regime has a very slow
convergence for WAL2 glass. For each cases we fit the energies with a
stretched exponential function to extract the asymptotic value of the steady
state energy.

In Fig. A3, we present the evolution of the energy during the shear cycles
leading to the steady state. In this case, the energy is presented as a
function of strain, over the full cycles for a set of representative
amplitudes below and above yielding. In the top panel, the evolution for the
ESL is presented. The data illustrate the annealing that occurs with repeated
cycles, with the energy becoming lower with successive cycles. When yielding
occurs and the glass is in a diffusive state, the limit cycles present two
minima in each cycle in agreement with previous observations Leishangthem  _et
al._ (2017). For the WAL1 case, presented in the lower panels of Fig. A3, we
observe no evolution below yielding and, as in the previous case, the
appearance of two minima above yielding. For WAL1 and $\gamma=0.55$, close to
the yielding point, we observe two different final states that can be obtained
depending on the initial configuration.

###  A-3 Stress profile over the strain cycle

In Fig. A4 we present averaged stress-strain curves over the complete
deformation cycle in the steady state for different strain amplitudes. We
average the data over 200-300 cycles in the steady state in each case. For low
strain amplitude we see the stress-strain curves do not enclose a finite area,
whereas for amplitudes larger than the yield strain they clearly do. For each
strain amplitude $\gamma_{max}$, the maximum stress value $\sigma_{xy}^{max}$
are collected and presented in Fig. 2(b)-(c). We see that with increasing
$\gamma_{max}$, $\sigma_{xy}^{max}$ increases and attains its maximum value at
the yield strain amplitude beyond which it drops suddenly. Such behaviour is
more prominent for well annealed glasses.

###  A-4 Mean Square Displacement in the steady state

To characterize the motion of the particles, we calculate the mean square
displacements (MSD) of the particles for stroboscopic configurations, as a
function of the number of elapsed cycles. We calculate this quantity in the
steady state, using different initial configurations, averaging over them, for
each accumulated strain difference $\gamma$:

${\rm MSD}(\gamma)=\left\langle\frac{1}{N}\sum_{i}^{N}|\boldsymbol{\rm
r}_{i}(\gamma_{acc}+\gamma)-\boldsymbol{\rm
r}_{i}(\gamma_{acc})|^{2}\right\rangle$ (A1)

where the averaging is performed over several zero strain configurations as
the initial configurations and over different samples. In Fig. A5, we present
the sample averaged data of MSD for different glasses we study. In Figs. A5
(a)-(c) we show the MSDs for different systems sizes of the HTL system, for
several strain amplitudes. In Figs. A5 (d)-(f) we present the data ESL, WAL1,
and WAL2 glasses for system size $N=10000$. Two regimes are clearly observed.
Above the yield strain amplitude, the MSD grows linearly with $\gamma$ whereas
it remains substantially flat below it. This sharp change corresponds to a
jump in diffusivity presented in Figs. 2(c) and 2(f).

###  A-5 Single particle displacement field and cluster size distribution

In Fig. A7(a) we show the distribution of single particle non affine
displacements $P(\delta r)$ during a plastic event for various strain
amplitudes. Such distributions are expected to have a power-law form with
exponent $-3$ followed by an exponential tail at large $\delta r$, separated
by a cutoff $\delta r_{c}$. The power-law part corresponds to elastic response
whereas the exponential tail is associated with the plastic rearrangements. We
fit the exponential part for the large $\delta r$ range for different
$\gamma_{max}$ and observe that the cutoff value $\delta r_{c}$ varies
approximately from $0.15$ to $0.35$ for the range of strain amplitudes we
consider. Choosing various cut-off $\delta r_{c}=0.15,0.20,0.25$, and $0.30$
we identify the active particles undergoing plastic rearrangements that have
$\delta r>\delta r_{c}$. For each choice of $\delta r_{c}$ we compute the
cluster size distributions, as shown in Figs. A7 (b)-(d), which reveal that
the results are almost independent of the choice of $\delta r_{c}$ in this
range. We use $\delta r_{c}=0.25$ in the analysis presented.

Figure A2:  Evolution of potential energy with accumulated strain
($\gamma_{acc}=4\times\gamma_{max}\times N_{cycle}$) for different system
sizes for (a)-(c) HTL glasses for various strain amplitudes. (d)-(f) Evolution
data of potential energy for ESL, WAL1, and WAL2 for a system of size
$N=10000$ for different strain amplitudes. Solid lines through the data points
are fits to a stretched exponential form.

Figure A3:  Potential energy against strain over the full cycle for ESL (upper
panels) and WAL1 (lower panels) and different strain amplitudes, for system
size $N=10000$. The energies from the early cycles are shown in grey while the
averaged value of the energy in the steady state is shown in red. Dashed blue
lines indicate the initial IS energy. For the case of ESL, the system
gradually anneals with increasing of strain amplitude up to the critical
amplitude $\gamma_{max}$, whereas for the case of WAL1, the system remains at
the same energy below the yield amplitude. Close to the yield strain
amplitude, different samples either remain in the absorbing state or make a
transition to the diffusive state. This bi-stable behaviour is shown for
$\gamma_{max}=0.055$ where two distinct final states are displayed.

Figure A4:  Stress-strain curves averaged over many cycles in the steady state
for different $\gamma_{max}$, (a)-(c) for different system sizes
($N=2500,4900,10000$) of the HTL glass, and (d)-(f) for ESL, WAL1 and WAL2
respectively for $N=10000$.

Figure A5:  Mean squared displacements computed in the steady state against
accumulated strain difference $\gamma$ for different strain amplitudes
$\gamma_{max}$. Different panels are for different types of glasses and system
sizes we study. To highlight the different behaviour across yield strain, we
use the same colour (orange) for the highest $\gamma_{max}$ before yield and
the lowest $\gamma_{max}$ above.

Figure A6:  (a) Distribution of single particle non-affine displacements
$\delta r$ during an avalanche for several strain amplitudes. The fitted
dashed lines for $\delta r>\delta r_{c}$, represent the exponential tail.
Vertical lines indicate the choice of various $\delta r_{c}$ for which we
perform avalanche analysis. (b)-(d) The cluster size distribution of active
particles ($\delta r>\delta r_{c}$) for several strain amplitudes. Different
panels are for different choice of $\delta r_{c}=0.15,0.20,0.30$. Data for
$\delta r_{c}=0.25$ is shown in main text.

Figure A7:  (a) Distribution of the total energy drop $\Delta U_{tot}$ during
a plastic rearrangement, for different strain amplitudes, for $N=10000$. (b)
Average $\langle\Delta U_{tot}\rangle$ against $\gamma_{max}$ for different
system sizes. (c) $\langle\Delta U_{tot}\rangle$ against system size for
different strain amplitudes. (d) Distribution of the plastic component of the
energy drop $\Delta U_{pl}$ for different $\gamma_{max}$. (e) Variation of
$\langle\Delta U_{pl}\rangle$ against $\gamma_{max}$ for different system
sizes. Inset shows the scaling behaviour of $\langle\Delta U_{pl}\rangle\sim
N^{1/3}$ in the post-yield regime. (f) $\langle\Delta U_{pl}\rangle$ against
system size for different strain amplitudes.

###  A-6 Statistics of Energy Drops

The plastic rearrangements lead to drops in energy that are often used as a
quantifier of avalanche size. However, the total energy drop $\Delta U_{tot}$
for the whole system during a plastic rearrangement contains both an elastic
and a plastic component. The energy drop arising from the plastic core,
$\Delta U_{pl}$, is computed considering only active particles that are
identified by the procedure described before. In Figs. A7 (a)-(c) we present
the statistics of $\Delta U_{tot}$ and in A7 (d)-(f) we show the results for
$\Delta U_{pl}$. Both the distributions $P(\Delta U_{tot})$ and $P(\Delta
U_{pl})$ display power-law regimes followed by a $\gamma_{max}$ dependent
cutoff. However, the associated exponent describing the power-law regime is
found to be different, $-1$ for $\Delta U_{tot}$ and $-1.25$ for $\Delta
U_{pl}$. The mean energy drops, $\langle\Delta U_{tot}\rangle$ and
$\langle\Delta U_{pl}\rangle$, show a similar dependence on strain amplitude
$\gamma_{max}$, but with $\langle\Delta U_{tot}\rangle$ displaying a system
size dependence below yielding, which is not present for $\langle\Delta
U_{pl}\rangle$. The system size dependence of $\langle\Delta U_{tot}\rangle$
can be described by power laws ranging from $N^{1/2}$ to $N^{1.25}$, whereas
$\langle\Delta U_{pl}\rangle$ does not display a system size dependence for
$\gamma_{max}<\gamma_{y}$ and a $N^{1/3}$ dependence, similar to the mean
cluster size, above $\gamma_{y}$.

#### Acknowledgements.

We thank J. Horbach, S. Zapper and especially A. Rosso for useful discussions
and comments on the manuscript. We acknowledge Indo-French Center for the
Promotion of Advanced Research (IFCPAR/CEFIPRA) Project 5704-1 for support,
the Thematic Unit of Excellence on Computational Materials Science, and the
National Supercomputing Mission facility (Param Yukti) at the Jawaharlal Nehru
Center for Advanced Scientific Research for computational resources. S.S.
acknowledges support through the J. C. Bose Fellowship (JBR/2020/000015) SERB,
DST (India).

## References

  * 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).
  * Lauridsen  _et al._ (2002) J. Lauridsen, M. Twardos, and M. Dennin, Phys. Rev. Lett. 89, 098303 (2002).
  * Bécu  _et al._ (2006) L. Bécu, S. Manneville, and A. Colin, Phys. Rev. Lett. 96, 138302 (2006).
  * Schall  _et al._ (2007) P. Schall, D. A. Weitz, and F. Spaepen, Science 318, 1895 (2007), https://science.sciencemag.org/content/318/5858/1895.full.pdf .
  * Greer  _et al._ (2013) A. Greer, Y. Cheng, and E. Ma, Materials Science and Engineering: R: Reports 74, 71 (2013).
  * Hufnagel  _et al._ (2016) T. C. Hufnagel, C. A. Schuh, and M. L. Falk, Acta Materialia 109, 375 (2016).
  * Guin and Wiederhorn (2004) J.-P. Guin and S. M. Wiederhorn, Phys. Rev. Lett. 92, 215502 (2004).
  * Zhao  _et al._ (2016) Q. Zhao, C. E. Finlayson, D. R. E. Snoswell, A. Haines, C. SchÃ¤fer, P. Spahn, G. P. Hellmann, A. V. Petukhov, L. Herrmann, P. Burdet, P. A. Midgley, S. Butler, M. Mackley, Q. Guo, and J. J. Baumberg, Nature Communications 7, 11661 (2016).
  * Chester and Chester (1998) F. M. Chester and J. S. Chester, Tectonophysics 295, 199 (1998).
  * Uhl  _et al._ (2015) J. T. Uhl, S. Pathak, D. Schorlemmer, X. Liu, R. Swindeman, B. A. W. Brinkman, M. LeBlanc, G. Tsekenis, N. Friedman, R. Behringer, D. Denisov, P. Schall, X. Gu, W. J. Wright, T. Hufnagel, A. Jennings, J. R. Greer, P. K. Liaw, T. Becker, G. Dresen, and K. A. Dahmen, Scientific Reports 5, 16493 (2015).
  * Keim and Arratia (2013) N. C. Keim and P. E. Arratia, Soft Matter 9, 6222 (2013).
  * Knowlton  _et al._ (2014) E. D. Knowlton, D. J. Pine, and L. Cipelletti, Soft Matter 10, 6931 (2014).
  * Hima Nagamanasa  _et al._ (2014) K. Hima Nagamanasa, S. Gokhale, A. K. Sood, and R. Ganapathy, Phys. Rev. E 89, 062308 (2014).
  * Di Michele  _et al._ (2014) L. Di Michele, D. Fiocco, F. Varrato, S. Sastry, E. Eiser, and G. Foffi, Soft Matter 10, 3633 (2014).
  * Keim  _et al._ (2020) N. C. Keim, J. Hass, B. Kroger, and D. Wieker, Phys. Rev. Research 2, 012004 (2020).
  * Maloney and Lemaître (2006) C. E. Maloney and A. Lemaître, Phys. Rev. E 74, 016118 (2006).
  * Karmakar  _et al._ (2010) S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E 82, 055103 (2010).
  * Jaiswal  _et al._ (2016) P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Phys. Rev. Lett. 116, 085501 (2016).
  * Fiocco  _et al._ (2013) D. Fiocco, G. Foffi, and S. Sastry, Phys. Rev. E 88, 020301 (2013).
  * Shrivastav  _et al._ (2016) G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Phys. Rev. E 94, 042605 (2016).
  * Regev  _et al._ (2015) I. Regev, J. Weber, C. Reichhardt, K. A. Dahmen, and T. Lookman, Nature Communications 6, 8805 (2015).
  * Leishangthem  _et al._ (2017) P. Leishangthem, A. D. S. Parmar, and S. Sastry, Nature Communications 8, 14653 (2017), article.
  * Parisi  _et al._ (2017) G. Parisi, I. Procaccia, C. Rainone, and M. Singh, Proceedings of the National Academy of Sciences 114, 5577 (2017), https://www.pnas.org/content/114/22/5577.full.pdf .
  * Parmar  _et al._ (2019) A. D. S. Parmar, S. Kumar, and S. Sastry, Phys. Rev. X 9, 021018 (2019).
  * Barbot  _et al._ (2020) A. Barbot, M. Lerbinger, A. Lemaître, D. Vandembroucq, and S. Patinet, Phys. Rev. E 101, 033001 (2020).
  * Schinasi-Lemberg and Regev (2020) E. Schinasi-Lemberg and I. Regev, Phys. Rev. E 101, 012603 (2020).
  * Urbani and Zamponi (2017) P. Urbani and F. Zamponi, Physical Review Letters 118, 38001 (2017).
  * 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, Phys. Rev. Lett. 125, 168003 (2020).
  * Liu  _et al._ (2020) C. Liu, E. E. Ferrero, E. A. Jagla, K. Martens, A. Rosso, and L. Talon, (2020), arXiv:2012.15310 [cond-mat.soft] .
  * Sastry (2021) S. Sastry, Phys. Rev. Lett. 126, 255501 (2021).
  * Sethna  _et al._ (2017) J. P. Sethna, M. K. Bierbaum, K. A. Dahmen, C. P. Goodrich, J. R. Greer, L. X. Hayden, J. P. Kent-Dobias, E. D. Lee, D. B. Liarte, X. Ni, _et al._ , Annual Review of Materials Research 47, 217 (2017).
  * Richard  _et al._ (2020) D. Richard, M. Ozawa, S. Patinet, E. Stanifer, B. Shang, S. A. Ridout, B. Xu, G. Zhang, P. K. Morse, J.-L. Barrat, L. Berthier, M. L. Falk, P. Guan, A. J. Liu, K. Martens, S. Sastry, D. Vandembroucq, E. Lerner, and M. L. Manning, Phys. Rev. Materials 4, 113609 (2020).
  * Bonfanti  _et al._ (2019) S. Bonfanti, R. Guerra, C. Mondal, I. Procaccia, and S. Zapperi, Phys. Rev. E 100, 060602 (2019).
  * Denisov  _et al._ (2015) D. V. Denisov, M. T. Dang, B. Struth, A. Zaccone, G. H. Wegdam, and P. Schall, Scientific Reports 5 (2015), 10.1038/srep14359.
  * Vasisht  _et al._ (2020) V. V. Vasisht, G. Roberts, and E. Del Gado, Phys. Rev. E 102, 010604 (2020).
  * Mitra  _et al._ (2021) S. Mitra, A. D. S. Parmar, P. Leishangthem, S. Sastry, and G. Foffi, Journal of Statistical Mechanics: Theory and Experiment 2021, 033203 (2021).
  * Dahmen  _et al._ (2009) K. A. Dahmen, Y. Ben-Zion, and J. T. Uhl, Phys. Rev. Lett. 102, 175501 (2009).
  * Lin  _et al._ (2014) J. Lin, E. Lerner, A. Rosso, and M. Wyart, Proceedings of the National Academy of Sciences 111, 14382 (2014), https://www.pnas.org/content/111/40/14382.full.pdf .
  * Oyama  _et al._ (2020) N. Oyama, H. Mizuno, and A. Ikeda, (2020), arXiv:2009.02635 [cond-mat.soft] .
  * Jagla (2015) E. A. Jagla, Phys. Rev. E 92, 042135 (2015).
  * Liu  _et al._ (2016) C. Liu, E. E. Ferrero, F. Puosi, J.-L. Barrat, and K. Martens, Physical review letters 116, 065501 (2016).
  * Bouchbinder  _et al._ (2007) E. Bouchbinder, J. S. Langer, and I. Procaccia, Phys. Rev. E 75, 036108 (2007).
  * Dahmen  _et al._ (2011) K. A. Dahmen, Y. Ben-Zion, and J. T. Uhl, Nature Physics 7, 554 (2011).
  * Franz and Spigler (2017) S. Franz and S. Spigler, Phys. Rev. E 95, 022139 (2017).
  * 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).
  * Bhaumik  _et al._ (2021a) H. Bhaumik, G. Foffi, and S. Sastry, Proceedings of the National Academy of Sciences 118 (2021a), 10.1073/pnas.2100227118.
  * Yeh  _et al._ (2020) W.-T. Yeh, M. Ozawa, K. Miyazaki, T. Kawasaki, and L. Berthier, Phys. Rev. Lett. 124, 225502 (2020).
  * Ozawa  _et al._ (2020) M. Ozawa, L. Berthier, G. Biroli, and G. Tarjus, Phys. Rev. Research 2, 023203 (2020).
  * Radhakrishnan and Fielding (2016) R. Radhakrishnan and S. M. Fielding, Phys. Rev. Lett. 117, 188001 (2016).
  * Shi and Falk (2006) Y. Shi and M. L. Falk, Phys. Rev. B 73, 214201 (2006).
  * Shi  _et al._ (2007a) Y. Shi, M. B. Katz, H. Li, and M. L. Falk, Physical review letters 98, 185505 (2007a).
  * Le Priol  _et al._ (2021) C. Le Priol, P. Le Doussal, and A. Rosso, Phys. Rev. Lett. 126, 025702 (2021).
  * Bhaumik  _et al._ (2021b) H. Bhaumik, G. Foffi, and S. Sastry, (preprint) (2021b).
  * Lancon  _et al._ (1984) F. Lancon, L. Billard, and A. Chamberod, Journal of Physics F: Metal Physics 14, 579 (1984).
  * Widom  _et al._ (1987) M. Widom, K. J. Strandburg, and R. H. Swendsen, Phys. Rev. Lett. 58, 706 (1987).
  * Shi and Falk (2005) Y. Shi and M. L. Falk, Physical Review Letters 95, 95502 (2005).
  * Shi  _et al._ (2007b) Y. Shi, M. B. Katz, H. Li, and M. L. Falk, Phys. Rev. Lett. 98, 185505 (2007b).
  * Barbot  _et al._ (2018) A. Barbot, M. Lerbinger, A. Hernandez-Garcia, R. García-García, M. L. Falk, D. Vandembroucq, and S. Patinet, Phys. Rev. E 97, 033001 (2018).
  * Das  _et al._ (2018) P. Das, A. D. Parmar, and S. Sastry, arXiv preprint arXiv:1805.12476 (2018).
  * Evans and Morriss (1984) D. J. Evans and G. P. Morriss, Phys. Rev. A 30, 1528 (1984).
  * Plimpton (1995) S. Plimpton, Journal of Computational Physics 117, 1 (1995).
  * Adhikari, Monoj and Sastry, Srikanth (2018) Adhikari, Monoj and Sastry, Srikanth, Eur. Phys. J. E 41, 105 (2018).
  * Kawasaki and Berthier (2016) T. Kawasaki and L. Berthier, Phys. Rev. E 94, 022615 (2016).
  * Vasisht and Del Gado (2020) V. V. Vasisht and E. Del Gado, Phys. Rev. E 102, 012603 (2020).
  * Sastry  _et al._ (1997) S. Sastry, D. S. Corti, P. G. Debenedetti, and F. H. Stillinger, Phys. Rev. E 56, 5524 (1997).
  * Maiti  _et al._ (2013) M. Maiti, A. Lakshminarayanan, and S. Sastry, The European Physical Journal E 36, 5 (2013).
  * Eshelby (1957) J. D. Eshelby, Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 241, 376 (1957).
  * Dyre (1999) J. C. Dyre, Phys. Rev. E 59, 2458 (1999).
  * Schröder  _et al._ (2000) T. B. Schröder, S. Sastry, J. C. Dyre, and S. C. Glotzer, The Journal of Chemical Physics 112, 9834 (2000).
  * Fiocco  _et al._ (2014) D. Fiocco, G. Foffi, and S. Sastry, Phys. Rev. Lett. 112, 025702 (2014).
  * Salerno and Robbins (2013) K. M. Salerno and M. O. Robbins, Phys. Rev. E 88, 062206 (2013).
  * Laurson  _et al._ (2010) L. Laurson, S. Santucci, and S. Zapperi, Phys. Rev. E 81, 046116 (2010).


# An imminent collision monitoring system with safe stopping interventions for
autonomous aerial flights

Jasmine Cheng, Xuning Yang and Nathan Michael The Authors are with the
Robotics Institute at Carnegie Mellon University, Pittsburgh PA, 15213, USA
{xuning<EMAIL_ADDRESS>
###### Abstract

Collision avoidance requires tradeoffs in planning time horizons. Depending on
the planner, safety cannot always be guaranteed in uncertain environments
given map updates. To mitigate situations where the planner leads the vehicle
into a state of collision or the vehicle reaches a point where no trajectories
are feasible, we propose a continuous collision checking algorithm. The
imminent collision checking system continuously monitors vehicle safety, and
plans a safe trajectory that leads the vehicle to a stop within the observed
map. We test our proposed pipeline alongside a teleoperated navigation in
real-life experiments, and in simulated random-forest and warehouse
environments where we show that with our method, we are able to mitigate
collisions with a success rate of at least 90%.

##  I Introduction

In tasks such as navigation and exploration in dense environments, static and
dynamic collision avoidance require trade-offs in planning time horizons.
Planning at a high frequency results in increased demand for computation at
unrealistic rates, while planning at low frequencies may result in inaccurate
plans if map updates between planning iterations. This work seeks to improve
aerial vehicle resilience and robustness in a variety of environments by
introducing a safety monitoring system that can supplement any autonomous or
semi-autonomous aerial system. Prior work has explored collision avoidance and
replanning with respect to goal-seeking, where the vehicle is autonomously
controlled and the goal is known a priori [1]. The proposed method mitigates
possible collisions at map update rate in a lightweight way using the vehicle
state, preventing the vehicle from reaching an unrecoverable state in
situations where the planner fails to mitigate the collision in time.

The proposed system prevents collisions by continuously monitoring collision
safety and generating a stopping trajectory if the vehicle is at risk of
collision. This system is first included in [2]. In this work, we demonstrate
the imminent collision monitoring system in a variety of scenarios and show
its applicability and robustness with respect to a wide variety of vehicle and
environment conditions. We show that, with the inclusion of the proposed
system across a variety of autonomous/semi-autonomous flights, we are able to
detect unsafe states and safely avoid collisions.

##  II Method

The proposed method proceeds as follows: the system continuously monitors
collision possibility by evaluating the closest obstacles to the vehicle’s
current location at a fixed rate. Then, we generate a dynamically feasible and
safe trajectory and bring the vehicle to a stop by first selecting an escape
point along the direction of motion, then generating a trajectory to the
escape point. The method is detailed as follows.

###  II-A Imminent collision checking

We first find the set of possible collision points nearby by evaluating the
normalized vector projection between a set of closest obstacle locations
$\\{\mathbf{p}\\}$ given by a known map and the vehicle velocity:

$\displaystyle\text{proj}(\mathbf{x},\mathbf{p})=\bigg{\langle}\frac{\mathbf{\dot{x}}}{\|\mathbf{\dot{x}}\|},\frac{\mathbf{p}-\mathbf{x}}{\|\mathbf{p}-\mathbf{x}\|}\bigg{\rangle}=\bigg{\langle}\frac{\mathbf{\dot{x}}}{\|\mathbf{\dot{x}}\|},\frac{\mathbf{r}}{\|\mathbf{r}\|}\bigg{\rangle}$

Where $\mathbf{x}$, $\dot{\mathbf{x}}$ is the vehicle position and velocity
respectively, and the vector to an obstacle point $\mathbf{p}$ is
$\mathbf{r}=\mathbf{p}-\mathbf{x}$. The points away from the direction of
motion, e.g., $\text{proj}(\mathbf{x},\mathbf{p}){<}0$ are discarded.

Then, for all points $\\{\mathbf{p}\\}$ in remaining set, we compute a
combined stop criterion based on distance to the obstacle, speed, and the
angle offset:

$\displaystyle
C_{\text{stop}}(\mathbf{x},\mathbf{p})=w_{1}\left\|\mathbf{r}\right\|-w_{2}\left\|\dot{\mathbf{x}}\right\|+w_{3}\arccos({\text{proj}(\mathbf{x},\mathbf{p})})$
$\displaystyle\text{if}\;{\text{proj}(\mathbf{x},\mathbf{p})}$
$\displaystyle>=0$

The third term represents the relative angle difference between the heading of
the vehicle and the direction of the obstacle point, as shown in Figure 3. A
imminent stop trajectory is issued if for any obstacle point $\mathbf{p}$,
$C_{\text{stop}}(\mathbf{x},\mathbf{p})<\beta$.

Figure 1: Trajectory-based teleoperation, current trajectory is collision
free, but imminent collision after the trajectory. Figure 2: Illustration of
imminent collision monitoring and safe stop. (A) Initial assessment: possible
collision points are sampled from nearby obstacles and vectors to collision
are assessed. If high concern vectors exist, stopping trajectory generation
begins. (B) Initial grid of escape points with costs computed and colliding
points discarded (C) A stratified sampling method is used to down sample the
possible escape points, with candidate trajectories generated. Dynamically
infeasible candidate trajectories are discarded. (D) Lowest cost, dynamically
feasible trajectory is selected and sent to the controller.  Figure 3: (A)
Vehicle is heading in the direction of an obstacle. Vehicle is in higher risk
of collision, $C_{\text{stop}}$ is lower. (B) Vehicle is heading between
obstacles/ Vehicle is in lower risk of collision, $C_{\text{stop}}$ higher.

###  II-B Escape points sampling

Escape points are defined as possible goal points where the vehicle could
stop. We first generate escape points and then plan trajectories to them. The
first feasible trajectory incurring the lowest cost is sent to the controller.
Thus, to minimize search time, we aim to check escape points that are more
likely to produce feasible trajectories.

An initial group of escape points $\\{\mathbf{e}\\}$ is generated using a
uniform 3D grid, scaled proportionally to the vehicle’s velocity. Points are
sampled along the direction of motion using stratified sampling for points
near the vehicle’s original heading via the following cost:

$C_{\text{escape}}(\mathbf{x},\mathbf{e})=w_{1}q+w_{2}d$

Where
$q{=}||(\mathbf{x}{-}\mathbf{e}{-}((\mathbf{x}{-}\mathbf{e}){\cdot}\hat{\dot{\mathbf{x}}})\hat{\dot{\mathbf{x}}}||$
is the shortest distance from the escape point $\mathbf{e}$ to the line
$l(s){=}\hat{\dot{\mathbf{x}}}s{+}\mathbf{x}$ along the velocity vector,
$\hat{\dot{\mathbf{x}}}$ is the normalized velocity vector, and $d$ is the
distance from the the escape point to its nearest obstacle.

First, points farther from obstacles are lower cost to encourage the vehicle
towards an area with more free space. Second, points closer to the vehicle’s
velocity vector are given lower costs. Points farther from the vehicle’s path
require the robot to turn severely, which results in undesirable dynamics.

Next, points are downsampled before trajectory generation via stratified
sampling, which samples a fixed number of points from each strata based on
cost. To ensure enough high quality points are chosen, more points are sampled
from higher strata (low cost). In this implementation, 10 points are sampled
from top 1%, 40 points from next 9%, 30 from next 40% and 20 from bottom 50%.
The stratified method guarantees a set of points with both high and low costs
such that in the event low cost trajectories are not dynamically feasible, we
are able to sacrifice cost for dynamic feasibility.

###  II-C Stopping trajectory generation

Given the set of escape points, we generate candidate stop trajectories
starting with the lowest $C_{\text{escape}}$.

A trajectory is a time-parameterized function $\xi(t)$, defined over a time
interval $t\in[0,T]$ that maps a given time $t$ to a position
$\mathbf{x}_{t}$. A stop trajectory is a trajectory that brings a vehicle from
in-motion to at-rest:

s.t. $\displaystyle\xi:[0,T]\rightarrow X$ $\displaystyle\xi^{(j)}(T)=0,j>0$
$\displaystyle\xi\in\Xi_{\text{safe}}$

A single polynomial trajectory is computed by solving an 8th order polynomial
on the current vehicle state and desired vehicle state. The current vehicle
state is defined by the current position, velocity, acceleration, and jerk in
each dimension ($x,y,z,yaw$). The desired state is the position of the escape
point, with all higher-order derivatives of 0. Each dimension’s polynomial is
computed independently.

We solve for the coefficients, $c_{i}$, of the polynomial

$f(x)=\sum_{i=0}^{8}c_{i}x_{i}$

with constraints

$\displaystyle f(0)$
$\displaystyle=x_{0},\dot{f}(0)=v_{0},\ddot{f}(0)=a_{0},f^{(3)}(0)=j_{0}$
$\displaystyle f(t)$ $\displaystyle=x_{f},\dot{f}(t),\ddot{f}(t),f^{(3)}(t)=0$

where $x_{0},v_{0},a_{0},j_{0}$ are the initial state and $x_{f}$ is the final
position given by the escape point.

Next, the trajectory is interpolated and waypoints are checked against a known
map to ensure it is collision-free. The trajectory is also checked for dynamic
feasibility via an acceleration bound along the trajectory, i.e.,
$\ddot{\mathbf{x}}{<}\bar{\alpha}$. For our experiments, we select
$\bar{\alpha}=10\text{m/s}^{2}$.

The first trajectory that is safe and feasible is executed, and the search
stops. If all potential escape points have been exhausted without finding a
valid trajectory, then a stopping trajectory with no goal position is executed
as a last effort.

Input: Vehicle position $\mathbf{x}$ and sample set of points $\mathcal{S}$

1 for  _$s\in\mathcal{S}$_ do

2 Compute $C_{escape}$;

3

4$S_{escape}$:=StratifiedSample($\mathcal{S}$);

5 sort(escapePoints, $C_{escape})$;

6

7for  _$\mathbf{e}_{i}$ in $S_{escape}$_ do

8 $\xi$:=SolvePolynomial($\mathbf{x}$, $\mathbf{e}_{i}$);

9 if  _feasible( $\xi$) and collisionFree($\xi$)_ then

10 executeTrajectory($\xi$);

11

Algorithm 1 Stopping Trajectory Generation

Figure 4: Trial run in the warehouse environment with imminent collision
checking enabled. Left: Vehicle odometry with queried obstacles. Locations A,
B, C indicate where stopping trajectories were issued corresponding with the
left plot. Right: Stopping trajectories are issued when stop cost drops below
a threshold, ensuring that the vehicle is safe. The operator recovers the
vehicle with an in-place yaw before normal flight is resumed.

##  III Experiment and Results

We evaluate our method in a simulated random forest environment and in a dense
realistic warehouse environment in a teleoperated navigation task minimizing
time. We conduct 20 trials in each environment, 10 with the proposed pipeline
enabled, and 10 without. In our experiments, the stop criterion parameters
$w_{1}$, $w_{2}$, $w_{3}$, $\beta$ are chosen to be $0.6$, $0.4$, $1.2$, $0.3$
respectively.

Figure 5: Testing scenarios. (top, left) Simulated warehouse environment.
(bottom, left) Real life flight arena. (right) Simulated random forest
environment.

We discuss task completion with and without enabling collision monitoring in
Table I. With collision avoidance enabled, the operator completes the task
with a success rate of 90% for the forest scenario and 100% for the warehouse
scenario. When collision monitoring is disabled, the vehicle collides under
aggressive flight, leading to 20% task completion rate. The average velocity
of enabled collision monitoring is comparable to disabled collision
monitoring, but the vehicle safety is maintained. “Distance” refers to the
average distance from obstacles throughout the entire trial, while “Min
Distance” is the minimum distance from an obstacle on a given trial.

TABLE I: Task completion with and without collision monitoring  Scenario | Success rate | Velocity (m/s) | Distance (m) | Min Distance (m)
---|---|---|---|---
Forest |  |  |  |
Enabled | $90\%$ | $1.37\pm 0.18$ | $1.61\pm 0.11$ | $0.54\pm 0.13$
Disabled | $20\%$ | $1.69\pm 0.18$ | $1.44\pm 0.22$ | $0.39\pm 0.16$
Warehouse |  |  |  |
Enabled | $100\%$ | $1.73\pm 0.07$ | $1.89\pm 0.11$ | $0.47\pm 0.08$
Disabled | $20\%$ | $1.80\pm 0.15$ | $1.82\pm 0.20$ | $0.32\pm 0.13$

Table II shows the average number of stops issued over the course of
completing task with collision monitoring enabled. We prevent 6 potential
collisions in the forest scenario on average, 4 collisions in the warehouse
scenario and 9 collisions in the hardware experiment. The table shows the
vehicle velocity and obstacle distance conditions that trigger these stops.

TABLE II: Results for trials with collision detecting enabled  Environment | # Stops Issued | Velocity (m/s) | Obstacle Distance (m)
---|---|---|---
Forest | $6\pm 1$ | $1.84\pm 0.12$ | $1.39\pm 0.09$
Warehouse | $4\pm 1$ | $1.97\pm 0.04$ | $1.37\pm 0.08$
Hardware | $9\pm 2$ | $0.38\pm 0.04$ | $0.73\pm 0.17$

In the hardware experiments, we teleoperate the vehicle in a flight arena with
two pillars with imminent collision monitoring enabled. The vehicle flies
between and around the pillars, but successfully stops if it approaches either
pillar.

Figure 6: The stop criterion parameters include distance and relative angle to
obstacle. An approximately linear boundary is formed between stop and no-stop
conditions. Red circles represent conditions for which stops were triggered.

Figure 6 shows the correlation between the obstacle offset angle and the
distance to the obstacle when determining whether to stop from hardware
trials. A rougly linear boundary is formed between stop and no-stop scenarios.
Some no-stop points lie below the boundary due to the third velocity factor,
because the vehicle does not stop when it near obstacles at low velocity.

##  IV Conclusion and Future Work

In this work, we presented an imminent collision avoidance system for safe
teleoperation in dense environments. This system revolves around three main
ideas. First, a continuous monitor on the vehicle state and nearby obstacles,
which will trigger a stop if conditions are met. Second, escape point sampling
to select potential goal locations in free space for the vehicle to stop.
Third, a polynomial trajectoy generator which iterates over escape points to
select a feasibile and collision-free stopping trajectory. The system can be
integrated in an existing teleoperation framework.

Future work can improve upon the stopping condition heuristic. This work used
a KD-tree map representation, but could be expanded to use other map
representations.

Figure 7: Trial run with imminent collision checking enabled in simulated
random forest environment (Left) and real-life flight arena (Right). Vehicle
odometry with queried obstacles, locations in red indicate where stopped
trajectories were issued.

## References

  * Ferreira et al. [2008] André Ferreira, Flávio Garcia Pereira, Raquel Frizera Vassallo, Teodiano Freire Bastos Filho, and Mário Sarcinelli Filho.  An approach to avoid obstacles in mobile robot navigation: the tangential escape.  _Sba: Controle & Automação Sociedade Brasileira de Automatica_, 19:395 – 405, 12 2008.  ISSN 0103-1759.
  * Yang et al. [2021] X. Yang, J. Cheng, and N. Michael.  An intention guided hierarchical framework for trajectory-based teleoperation of mobile robots.  _2021 IEEE/RSJ International Conference on Robotics and Automation (ICRA)_ , 2021.