# Soft Matrix: Extracting Inherent Length Scales in Sheared Amorphous Solids

Monoj Adhikari  Tata Institute of Fundamental Research, 36/P, Gopanpally
Village, Serilingampally Mandal,Ranga Reddy District, Hyderabad, 500046,
Telangana, India  Pinaki Chaudhuri  The Institute of Mathematical Sciences,
CIT Campus, Taramani, Chennai 600113, India  Homi Bhabha National Institute,
Anushaktinagar, Mumbai 400094, India  Smarajit Karmakar  Tata Institute of
Fundamental Research, 36/P, Gopanpally Village, Serilingampally Mandal,Ranga
Reddy District, Hyderabad, 500046, Telangana, India  Vishnu V. Krishnan  Tata
Institute of Fundamental Research, 36/P, Gopanpally Village, Serilingampally
Mandal,Ranga Reddy District, Hyderabad, 500046, Telangana, India  Nandlal
Pingua  Department of Physics Indian Institute of Technology Roorkee 247667
Uttarakhand, India  Shilditya Sengupta  Department of Physics Indian Institute
of Technology Roorkee 247667 Uttarakhand, India  Aparna Sreekumari  Department
of Physics, Indian Institute of Technology Palakkad, Nila Campus, Kanjikode,
Palakkad 678623, Kerala, India  Vishwas V. Vasisht 111The authors are listed
in alphabetical order of their last names. To whom correspondence should be
addressed. E-mail<EMAIL_ADDRESS>Department of Physics, Indian
Institute of Technology Palakkad, Nila Campus, Kanjikode, Palakkad 678623,
Kerala, India

###### Abstract

Amorphous solids yield upon crossing a strain threshold, after an initial
elastic response, when subjected to mechanical deformation. The yielding
process is characterized by local plastic events leading to non-affine
displacements, and their interactions. Despite the lack of long-range
structural order, these disordered materials exhibit long-range spatial
correlations in the non-affine displacement fields, which stems from the
underlying elasticity. Measuring a correlation length scale in deformed
amorphous solids, during the plastic process, is a non-trivial challenge,
often requiring an ad-hoc definition of localized regions. In this paper, we
introduce a novel computational approach called the ”soft matrix” that enables
systematic analysis of mechanical response of local regions within a
disordered solid. In this method, we subject the amorphous solid to a
quasistatic shear and allow a local region of interest to relax freely while
allowing for elastic relaxation of the background. The dependence of the yield
strain upon the size of the probe region naturally reveals the existence of an
intrinsic length scale ($\zeta$) that governs the elasto-plastic properties,
as observed in four distinct model amorphous solids. This finding demonstrates
the universality of this characteristic length scale across a wide range of
materials. We investigate the dependence of this length scale on the
material’s preparation history and find that $\zeta$ increases with better
annealing. Furthermore, the local mechanical properties measured within this
framework provide more accurate estimates compared to existing techniques. Our
study paves the way for a comprehensive understanding of amorphous solids and
facilitates improved characterization and design of these materials.

##  I Introduction

Amorphous solids are disordered systems which are lacking long-range
structural ordering. They encompass a broad range of materials ubiquitous in
our daily lives, including colloids, foams, emulsions, granular matter as well
as metallic and silicate glasses [1]. When subjected to an external load, such
as shear deformation, they display an elastic response at low strain values
and yield beyond a threshold strain [2, 3, 4, 5, 6, 7, 1, 8]. Unlike the
crystalline systems where defects are known to be careers of plasticity, the
physical mechanisms underlying the shear response of an amorphous material
still lacks a complete understanding. The last two decades of research in
theory and experiments have revealed that the plasticity in amorphous solids
occurs through local non-affine displacements or rearrangements resulting in
the redistribution of elastic stresses within the system [9, 10, 11].
Simulation works on athermal systems performed under quasistatic shearing
conditions have shown that the displacement field exhibits a quadrupolar
symmetry resembling the Eshelby inclusion model [12, 11, 13, 14].

Amorphous materials exhibit structural and dynamical heterogeneities [15]. In
recent years, there has been a growing interest in understanding the possible
existence of a characteristic length scale associated with the relaxation
process in driven amorphous solids. This interest stems from its impact on
various aspects, including shear start-up properties [7, 16, 17, 18, 19, 20],
yielding mechanisms [21, 22, 23, 24], and flow heterogeneities [25, 26, 7, 18,
27]. Previously, static correlation lengths have been extracted using non-
affine shear deformation protocols for colloidal glasses and supercooled
liquids [28, 29, 30], which demonstrated that the measured lengthscale is
highly dependent on the applied shear deformation and likely describes the
structural correlations in the initial configuration. Similar protocols have
been employed to extract the correlation length in athermal colloidal
suspensions under steady-state flow with finite-rate shearing [31], showing
dependence on both shear strain and shear rate. Based on ideas of non-local
fluidization, a similar dynamical correlation length has also been predicted
and measured [32, 33, 34]. However, a comprehensive comparison between the
extracted correlation length and other measures of static length scale that do
not involve shear deformation protocols, such as the point-to-set length
scale[35], length scales obtained from finite-size scaling of relaxation times
[36], or minimum eigenvalues [37], is lacking.

Despite the absence of a clear measure of the characteristic length scale,
various approaches have been proposed and utilized to characterize the size of
plastic events or the number of particles involved in non-affine
rearrangements. These methods include analyzing structural properties (e.g.,
free volume [38], local ordering [39]) and linear response measures (e.g.,
elastic moduli [40], soft modes [41]). The frozen matrix method has recently
gained widespread use in computing local mechanical properties, including
yielding thresholds and viscoelastic properties. Initially proposed for
assessing local elastic properties [42], this method involves allowing only
the specific ”target” local region of interest to relax independent of the
background. Under external loading, the entire system undergoes a uniform
affine deformation, except that the local region is permitted to relax freely,
while the rest of the system is restricted from undergoing non-affine
relaxation. By analyzing the stress response in the target region to the
applied strain, the local mechanical properties can be determined.

The frozen matrix method has been widely utilized in various studies on
amorphous solids [43, 44, 45, 46, 47, 48, 49]. However, in all cases, the size
of the target region is chosen in an ad-hoc manner, and the background
material is not allowed to undergo any relaxation and is made rigid - hence
frozen. Although this approach provides reasonable estimates of the local
yield stability distribution, it tends to overestimate the local yield stress
and elastic modulus significantly due to the presence of a rigid background
[43, 46, 44, 48]. This overestimation can have an impact on the predictability
of local plastic events and the correlations between them. Additionally, the
frozen matrix method does not enable the extraction of the length scale
associated with the plastic event.

In this study, we introduce a novel approach called the ”Soft matrix method”
to investigate the local mechanical properties, as described in the following
section. This method builds upon the original proposition of the frozen matrix
approach and it actually provides us a route to extract a characteristic
length scale $\zeta$ in various types of amorphous materials, overcoming the
limitations imposed by rigid or periodic boundary conditions [50, 51]. We
observe that this length scale is influenced by the aging conditions of the
material, and we demonstrate that it exhibits scale-free power law behavior
when studying larger sub-systems. Consequently, this length scale provides a
natural solution to the question of appropriate coarse-graining length,
serving two main purposes: (a) enabling accurate computation of local
mechanical properties crucial for understanding the yielding phenomenon and
providing inputs for mesoscale modeling of amorphous materials and (b)
facilitating comparisons between this length scale and other measures of
static length scale in disordered systems.

Figure 1: Soft Matrix Methodology. (a) Schematic demonstrating the soft matrix
method. The sub-system (central blue box) is subjected to unconstrained
relaxation, whereas the background (surrounding pink region) relaxes with a
constraint imposed by the spring tether (see text for details). (b.) Snapshot
from 2D BMLJ (65:35) simulation illustrating the soft matrix: plastic event is
confined to subsystem, while the background is allowed to elastically relax.
Here X is the shearing direction, and Y is the gradient direction. (bottom)
Variation of stress ($\sigma$) as a function of strain ($\gamma$) for (c)
fixed subsystem size ($L_{s}=10a$) and varying $k$ value and (d) fixed $k$
value ($k=10\epsilon/a^{2}$) and varying subsystem size $L_{s}$.

##  II Soft matrix method

In this approach, within a system having linear dimension $L$, we choose a
local region having size $L_{s}$, whose mechanical properties we plan to
measure. The rest of the system from hereon will be referred to as the
background. The mechanical response of the solid is probed via the athermal
quasi-static shear (AQS) protocol (see Methods). The entire system is first
subjected to a simple shear deformation of magnitude $\delta\gamma\approx
10^{-5}$, and hence an affine displacement is imposed on all particles in the
system. Post deformation, the system is subjected to relaxation or energy
minimization, and here we introduce an essential improvement over the existing
frozen matrix method [43, 44, 45, 46] which is key to extracting the
heterogeneity length scale $\zeta$. Unlike the frozen matrix method, during
the energy minimization, the background can relax by affine rearrangements,
but no plastic events are permitted in the background. On the other hand, the
sub-system is allowed to relax without any constraints. Hence we term the
procedure a soft matrix method. We implement this partial relaxation of the
background by adding an additional spring force of the type $F_{r}=-kr$,
before minimization. Here the spring constant $k$ is a tunable parameter
determining the coefficient of restoring force, and $r$ is the displacement of
a background particle from its reference position, taken to be the position
after applying the step of affine displacement.

We test out the soft matrix method via extensive numerical simulations using
the 2D BMLJ model and later apply our understanding to study three other model
amorphous solids, namely 3D Soft-Rep, 3D SiO2, and 3D Metallic Glass. See
Methods for further details. All the simulations were carried out using LAMMPS
molecular dynamics package [52].

In Fig.1, we show (a) schematic and (b) actual simulation configurations for
the 2D BMLJ model with $L=90a$ and the number of particles $N=10000$. Here $a$
represent the unit length scale. The background region (pink in Fig.1a)) is
subjected to additional restoring force during the relaxation with an
appropriate $\it k$ value such that (almost) no irreversible plastic motion is
observed in the background particles. In Fig.1(b), the colored particles in
the sub-system relax with both affine and non-affine displacement and can
involve in plastic events. The color bar is based on the magnitude of
displacement.

Figure 2: Measuring characteristic length scale. Relative yield strain
$\Delta\gamma_{y}$ as a function of subsystem size $L_{s}$ for four different
models: (a) 2D BMLJ (b) 3D Soft-Rep (c) 3D SiO2 and (d) 3D Metallic glass. The
open symbols are data obtained from the soft matrix method. The opaque symbols
are from the frozen matrix method. Different symbols correspond to different
system sizes. The thick dashed line is the exponential fit, and the thin
dashed line is just a guide to the eye.

##  III Shear startup response - first plastic event

In Fig.1(c), we show the load curves (shear stress $\sigma$ vs. applied strain
$\gamma$) for a range of k values at a fixed sub-system size ($L_{s}=10a$) for
the 2D BMLJ model. For a chosen sample, prepared by the standard quench
protocol (see Methods), we initially perform an AQS simulation of the entire
system and obtain the bulk stress-strain curve (load curve), via which we
recognise the first local yielding event (referred to herein as the first
plastic event) and also identify the spatial location of the event as well as
the particles participating in the yielding event (see Methods). We construct
a sub-box of linear length $L_{s}$ around the first plastic event and perform
the partial relaxation of the surrounding. The bulk load curve corresponds to
$k=0$. In Fig.1(c), we show that with an increase in chosen k values, the
system’s shear modulus increases, as clearly observed via an increase in the
slope of the stress-strain curve, which implies that the system becomes more
and more rigid. Further, with increasing k, the first plastic event is also
delayed and occurs at a higher strain value; thereby the local yield stress
value also increases. As $k\rightarrow\infty$, we should reach the frozen
matrix limit, and that is consistent with our observed trend of larger yield
stress and larger strain thresholds with increasing $k$.

In the current work, we are primarily interested in determining the
heterogeneity length scale $\zeta$. Hence, for further analysis, we fix the
value of $k$ by the following two competing requirements: (i) $k$ should be
small enough so that the probe itself does not significantly alter the local
mechanical properties of the system, and (ii) $k$ should be large enough such
that (almost) no plasticity occurs in the background system during the stress
relaxation. We verify (i) by monitoring the bulk modulus and pressure, which
remain unchanged for a range of $k$ values. For condition (ii), we monitor the
displacement in the sub-system and surrounding and find that for appropriate
$k$ values, the background shows no large displacements. (See SI for further
details).

Having the fixed $k$, we now examine the effect of sub-system size $L_{s}$ on
the load curve, specifically on the yielding strain. In Fig. 1(d) we show load
curves for $\it k=10\epsilon/a^{2}$ for varying $L_{s}$. The solid black line
represents the load curve for an unconstrained system which shows a yield
strain of $\gamma_{y}\approx 0.002$ at the first plastic event. With varying
$L_{s}$, the initial linear response does not show much variation, but we
observe a systematic change in local yield strain $\gamma_{y}$ and local yield
stress $\sigma_{y}$. With the increase in $L_{s}$, $\gamma_{y}$ and
$\sigma_{y}$ decrease. Approaching $L_{s}\rightarrow L$, we recover the
unconstrained load curve. We find similar features in the other three model
systems that we have studied. In each case, the value of $k$ was fixed
according to the above-prescribed protocol. One can rationalise the increasing
yield stress and strain with decreasing $L_{s}$ in terms of interface effects,
with the background curtailing the motion of particles in the sub-system and
hence increasing the barrier height for the first plastic event. However, we
find a more interesting observation when we look at the variation of
$\gamma_{y}$ with $L_{s}$, which we discuss below.

##  IV A characteristic length scale $\zeta$

In Fig. 2, we show the variation of relative yield strain $\Delta\gamma_{y}$
with the sub-system size $L_{s}$ for the four different models of amorphous
solids. The relative yield strain $\Delta\gamma_{y}$ is defined as the
deviation of $\gamma_{y}$ from the bulk value:
$\Delta\gamma_{y}=\gamma_{y}(L_{s})-\gamma_{y}(L)$ where $\gamma_{y}(L)$ is
the yield strain for the bulk system ($k=0$). The data shown in Fig. 2
corresponds to poorly aged samples, and we obtain the average
$\Delta\gamma_{y}$ over multiple initial configurations as well as different
system sizes (whenever possible). All four systems show $\Delta\gamma_{y}$
decreases exponentially as we increase the size of the probing zone $L_{s}$,
suggesting the existence of a characteristic length scale. The exponential fit
function $a_{0}\exp(-L_{s}/\zeta)$ is shown as dashed lines. Here the fit
parameter $a_{0}$ represents the limiting value of the relative yield strain
in the limit of zero sub-system size, and $\zeta$ is the estimated value of
the characteristic heterogeneity length scale intrinsic to each model system.
In 2D BMLJ and 3D Soft-rep models, which have short-range interactions, our
fit finds $\zeta\approx 5a$. Interestingly, this is similar in magnitude to
values often assumed in the literature for local rheological measurements in
similar models [40, 44, 46], where it is interpreted as the length scale below
which Hooke’s law for linear continuum elasticity breaks down. For the first
time, we have provided the rationality of this choice by directly extracting
an intrinsic scale via the proposed soft matrix method. For 3D $SiO_{2}$ and
3D Metallic glass, we find $\zeta$ to be $\approx 12a$ and $\approx 8a$
respectively. The higher values in these models are possibly due to the long-
range nature of inter-particle interactions, unlike the Lennard-Jones-type
models typically studied. In Fig.2, we also show corresponding data obtained
using the frozen method for the purpose of comparison. It is important to note
that at small sub-system lengths, the stress shows a significant increase for
the frozen matrix method, indicating the influence of how the background
relaxation is modeled. This rise is expected to depend on the system size, and
in larger system simulations, the impact of the interface will diminish, which
is what we observe. Thus, we argue that this approach is not suitable for
accurately extracting any characteristic scale which is much less than system
size.

Figure 3: Age dependence. (a) Relative yield strain as a function of subsystem
size for different ages of the initial sample (shown for 2D BMLJ). The inset
has the same data in log-linear scale, for selected ages. Dashed lines
represent fits to exponential functions, from which respective characteristic
length $\zeta$ is extracted. (b) Characteristic length $\zeta$ as a function
of energy of the initial state (i.e. at zero strain), $U_{\rm init}$; note
that $U_{\rm init}$ decreases with an increase in age. The dashed line is a
guide to the eye.

##  V Age dependence of characteristic length scale, $\zeta$

After extracting the characteristic length scale from the exponential
relationship between the local yield strain and the probe window size, we now
investigate the influence of the preparation protocol on this length scale.

It is well known that the mechanical response varies dramatically with the
increasing age of the amorphous solid samples [2, 53, 54, 20, 55]. For
example, a poorly aged or annealed sample often shows more homogeneous
ductile-like yielding behavior, while a well-aged or ultra-stable glassy
sample exhibit heterogeneous brittle-like yielding via shear banding [21, 20,
56]. Thus it is essential to understand how the characteristic length scale
that we obtained using our soft-matrix method depends on the sample age or
annealing history, which is presented here for 2D BMLJ glass. We study samples
quenched by infinite cooling rate (obtained by energy minimization of
configurations generated at different initial temperatures $T_{W}$ to
instantaneous inherent structures) as well as the ones prepared via finite
cooling rate $\Gamma$ 222Around $250$ samples were prepared at infinite to
finite but higher $\Gamma$, $100$ samples at moderate $\Gamma$, and $50$
samples at the lowest $\Gamma$. We then compute the $\Delta\gamma_{y}(L_{s})$
via the soft matrix method protocol - see Fig. 3 (a) main panel; note that we
have scaled the abscissa by the age-dependent fit parameter $a_{0}$ for better
clarity. We find that $\Delta\gamma_{y}$ has a slower decay with $L_{s}$, for
increasing age of the system. The exponential nature of
$\Delta\gamma_{y}(L_{s})$ is clearly shown in the inset of Fig. 3 (a), via
linear-log scale. Therefore, it is evident that with increasing age, the
characteristic length scale increases significantly The extracted $\zeta$ as a
function of the energy of the sample at zero strain is shown in Fig. 3 (b);
lower energy of the initial state indicates higher age of the sample. We see
that the value $\zeta\approx 5a$ occur only for a range of the poorly aged
sample, and $\zeta$ systematically increases with an increase in age. Hence
our analysis clearly suggests that any measure of local mechanical properties
should consider the $\zeta$ dependence with age. Further work needs to be done
to understand the functional form of the age dependence on $\zeta$.

Figure 4: Large length scale behaviour. Absolute yield strain $\gamma_{y}$ as
a function of subsystem size $L_{s}$ for two different system sizes (opaque
circle: L=90 and open circle: L=320) for the 2D BMLJ model (using poorly aged
samples). Fit lines (red dashed line - exponential fit and solid blue line -
power law fit with exponent $1.57$) show a cross-over (at around $15a$) from
an exponential regime at small sub-systems to a power law regime at large sub-
systems.

##  VI Large length scale behaviour - beyond $\zeta$.

The macroscopic response of amorphous solids shows long-range elastic
correlations [11] and power law avalanche statistics in the steady state [50].
Via numerical simulations done using periodic boundary conditions, it has been
shown that, in the steady state, the finite size dependence of the average
strain interval between two stress drops has the following behaviour:
$\langle\Delta\gamma\rangle\sim N^{-\beta}=L^{-\beta d}$, where $d$ is the
dimension, and the exponent $\beta\simeq 0.7$ [50]. One would thus expect that
in $2D$, the average strain interval for the first plastic drop should go as
$\langle\gamma_{y}\rangle\sim L^{-1.4}$. We test whether this behaviour can be
recovered at large sizes of the sub-box, for which we have utilized large
system size soft matrix simulations ($N=125$K; $L=360$) and obtained average
yield strain as a function of $L_{s}$. In Fig.4, we show the yield strain
variation with changing $L_{s}$ for the 2D BMLJ model (in a poorly aged
sample) using two different system sizes ($L=90$ and $L=320$). We observe two
interesting features. Firstly, with the increase in $L_{s}$, the average yield
strain, $\gamma_{y}(L_{s})$, crosses over from an exponential behaviour to a
power law dependence having an exponent around $1.57$. The crossover length is
around $10a$. Secondly, we find that in the limit of $L_{s}\rightarrow L/2$,
$\gamma_{y}$ deviates from power law and seems to saturate, which likely
originates from the effects due to periodic boundary condition. These
observations highlight effectiveness of our the proposed soft-matrix in not
only extracting the inherent length scale associated with sheared amorphous
solids but also reproducing the large length scale behaviour. The power law
exponent from our fit is not very off from the expected value [50], but we
note that the soft elastic confinement employed in our method is different
from the standard periodic boundary condition. This might affect the power-law
exponent, which is an interesting direction for future investigations. Also,
finite size effects cannot be ruled out. One needs to be sufficiently far away
from the cross-over length scale to suppress all the effects of the underlying
length scale. Hence for well annealed samples, one should go to larger system
sizes to extract the same power-law exponent. Further it has been recently
pointed out [51] that the asymptotic power-law behaviour with an exponent
close to $\beta\simeq 0.7$ can only be obtained at system sizes much larger
for well annealed (ultra stable) glasses than currently accessible via
computer simulation. Further studies across various models in different
dimensions might shed more light on this puzzle.

##  VII Local mechanical properties

After extracting the characteristic length scale $\zeta$ using the soft matrix
method, we employ this approach to determine the local yield threshold and
local elastic modulus $\mu_{i}$. The local yield threshold (X) is defined as
$X=\sigma^{i}_{y}-\sigma^{i}_{0}$, where $\sigma^{i}_{y}$ and $\sigma^{i}_{0}$
are the yield stress and zero strain stress at the local region $i$, whose
linear size is of the order of $\zeta$. The local storage modulus $\mu_{i}$ is
obtained from the slope of the local load curve.

In Fig.5(a), we show for the 3D Soft-Rep system (using $N=97556$, $L=42a$ and
$k=100$), the probability distribution function P(X) (inset) and cumulative
distribution function F(X) (main panel) obtained from soft matrix as well as
frozen matrix method 333We find that in 2D systems, the fluctuations dominate
the local load curve and hence noisy for $N=10K$. Such issues are minimized
when we study 3D systems. We find that the frozen matrix not only
overestimates the yield threshold, which would be expected, but also the
functional behaviour is different. The fit lines in Fig.5(a) correspond to
Weibull distribution $F(X)=1-exp[-(X/s_{1})^{s_{2}}]$ with scale parameter
$s_{1}$ being 1.9 (soft matrix) and 1.5 (frozen matrix). The shape parameter
$s_{2}$, which is associated with the instantaneous yield rate [59, 60],
varies by a factor of 3, with $s_{2}=2.3$ (soft matrix) and $s_{2}=6.2$
(frozen matrix). Examining the implications of these differences on the
mesoscale modeling results (where P(X) is a key input) would be an interesting
aspect to explore.

In the Fig.5(b), we show $P(\mu)$, the distribution of local storage modulus
along with a Gaussian fit. Unlike the frozen matrix method, which perturbs the
elasticity of the background drastically, the soft matrix method allows for
having the background elasticity close to the unconstrained bulk value and
hence $P(\mu)$ is distributed around the bulk value of the storage modulus
(marked by green upper triangle in Fig.5(b)), which is measured from the slope
of the load curve obtained from unconstrained simulations.

In conclusion, the soft matrix method introduced in this study not only
enables the extraction of the length scale $\zeta$, but also provides
improvements over the frozen matrix approach in terms of estimating local
mechanical properties.

Figure 5: Local Mechanical Properties. (a.) Distribution of local yield stress
threshold ($X=\sigma_{y}-\sigma_{0}$) for 3D soft-rep model (poorly aged),
comparing the soft matrix (circle) with frozen matrix (square) methods. Inset
shows the probability distribution P(X), and the main figure show the
cumulative distribution F(X). The fit lines are drawn using Weibull
distribution (see main text for details). (b.) Distribution of local shear
modulus ($\mu$) for 3D soft-rep model (poorly aged) comparing the two methods,
shown along with Gaussian fits (using lines). Also marked is the bulk shear
modulus (green triangle).

##  VIII Summary and Discussion

We have devised a numerical technique for measuring local mechanical
properties in amorphous solids. Labelled as the “soft matrix method”, we have
utilized it to successfully extract the characteristic length scale associated
with plastic events or non-affine displacements occurring during the
deformation process. Our approach provides an opportunity to explore the local
mechanical response within a purely elastic background. By applying this
method to four different models of amorphous solids, covering a wide range of
materials, we have demonstrated the universal existence of this intrinsic
length scale. Furthermore, we show that this length scale increases with the
age of the sample, indicating its dependence on the preparation history.
Importantly, our method allows us to accurately estimate the local yield
threshold and local elastic modulus, greatly improving upon the accuracy of
the existing frozen matrix method.

A natural question that arises from our results is the origin of the
characteristic length scale, $\zeta$, in the material. The simplest
interpretation is that $\zeta$ corresponds to the spatial scale over which
local plastic relaxations are devoid of any interfacial effects; i.e. it
captures the intrinsic scale over which the first localized irreversible event
would occur in a bulk system under mechanical loading. Thereby, it also
relates to the spatial scale beyond which interference from other plastic
events would come into play leading to screening of power-law correlations, as
has been proposed recently [61]. Previous studies have extensively explored
the concept of an intrinsic static correlation length scale associated with a
putative amorphous order in supercooled liquids and glasses [15]. In the
context of shear amorphous solids, similar length scales have been extracted
using non-affine shear deformation protocols [28, 29, 30, 31]. However, a
comprehensive comparison with other measures of static length scale that do
not involve shear deformation protocols, such as the point-to-set length scale
[35], length scale derived from finite size scaling of relaxation times [36],
and minimum eigenvalues [37], is currently lacking. It is expected that these
different measures of static length scale in the supercooled liquid regime are
proportional to each other [62].

Our work offers an opportunity to compare and investigate various length
scales, shedding light on the origin of such correlations. The observed
increase in $\zeta$ with better annealing is consistent with the corresponding
increase in the static correlation length obtained through analysis of non-
affine displacement fields or finite size scaling of different observables in
supercooled liquids. It also resembles the increase in the typical distance
between soft spots as inferred from the measured number density of quasi-
localized modes [63, 64]. However, without a detailed comparison, it is not
immediately evident whether the length scale $\zeta$ obtained using the soft
matrix method is identical to these other length scales, although it is highly
likely that $\zeta$ is closely related to the length scale of quasi-localized
modes. The connection to the static length scale of glass transition as
determined by the point-to-set method is not apparent unless we assume that
the point-to-set length scale and the length scale derived from finite size
scaling of minimum eigenvalues are fundamentally the same [62].

Despite the aforementioned limitations, it is noteworthy that the
characteristic length scale $\zeta$ obtained through the soft matrix method
exhibits remarkable similarity to the coarse-graining length scales reported
in various studies on elasto-plastic models. These models often define an
appropriate coarse-graining scale, typically around $5$ in suitable units, to
denote the size of elementary meso-block beyond which the propagation of
plasticity occurs via the Eshelby kernel [40, 44, 47, 45, 46, 49]. It is
interesting to observe that in order to study the effect of ageing or
annealing in these elasto-plastic models, systematically increasing the
coarse-graining length scale is often required to achieve a comparable
correspondence with microscopic simulations [51]. Hence, it appears that the
length scale derived from the soft matrix method may serve as the desired
length scale for the development of a robust mesoscopic elasto-plastic model
capable of incorporating ageing or annealing effects in a more fundamental
manner.

Overall, our work provides valuable insights into the behavior of amorphous
solids at a local level, enhancing our understanding of their macroscopic
properties and paving the way for improved mesoscale models required for
material design.

##  IX Methods

The four different models studied in this work are as follows.
(1.) 2D BMLJ: Binary Mixture (65:35) of Lennard-Jones particles is a well-
studied model glass forming system [65], whose mechanical properties are well
characterised. Two different system size were considered, N=10K and 125K, at a
reduced density $\rho=1.2$. The scales of A-A interactions within the model
provide the units for length ($a$), energy ($\epsilon$).
(2.) 3D Soft-rep: A system of 10% polydisperse soft spheres interacting via
the WCA potential [66], is a typical colloidal model in 3D. At a volume
fraction of $\phi=0.7$, system size of N=1K, 5K, 10K and 100K considered. the
particle has an average diameter of $a$ (unit length), mass $m$ (unit mass),
and energy $\epsilon$ (unit mass).
(3.) 3D SiO2: Most common molecular glass is modeled using the VSP modified
BKS potential of 1:2 binary mixture of silicon and oxygen atoms [67]. System
sizes of $N=13824,46656$ particles at a density of $\rho=2.8\,gm/cm^{3}$ are
considered. In this model unit length is $a_{Si-O}$, unit energy
$\epsilon_{Si-O}$ and unit mass $m_{Si}$.
(4.) 3D Metallic Glass: A metallic glass system $Cu_{64.5}Zr_{35.5}$ alloy
using the embedded atom model (EAM) [68]. System sizes $N=50000,100000$
particles are considered at zero pressure. The units of length, mass, and
energy are expressed in angstroms, grams/mole, and electron volts (eV),
respectively. The lengths are scaled in terms of the diameter of Cu (2.56
angstroms).

###  IX.1 Sample preparation protocol

For 2D BMLJ and 3D Soft-rep systems, equilibrium NVT MD simulations were
performed at $T=5.0\,\epsilon/k_{B}$ with $\rho=1.2$ and $\phi=0.7$
respectively and equilibrated configurations were quenched to $T=0.01$ at
cooling rates $\Gamma$ varying from $10^{-2}$ to $10^{-6}$. These
configurations were then subjected to potential energy minimization protocol
(FIRE or Conjugate gradient, CG). For 3D Silica, configurations were
equilibrated by NVT MD simulations at $T\approx 3100\,K$ and
$\rho=2.8gm/cm^{3}$ followed by instantaneous quench via potential energy
minimization (CG). For 3D metallic glass, NPT MD simulations were performed at
zero pressure and at temperatures $T=1100,1200,1300\,K$ to equilibrate the
system. Then configurations were quenched to $T=300K$ at cooling rates
$\Gamma=10^{10},10^{11},10^{12}\,K/s$ followed by potential energy
minimization by CG. In each model and cooling rate, a minimum of 10 to a
maximum of 250 samples were prepared. The shear stress $\sigma$ is computed
from the off-diagonal term of the virial stress tensor, using the Irving-
Kirkwood form [66].

###  IX.2 Shearing protocol

We probe the mechanical response using the Athermal Quasi-static Shear (AQS)
protocol [69]. In this procedure, an elementary shear strain
$\delta\gamma\sim\mathcal{O}(10^{-5}-10^{-6})$ is applied to a sample by the
uniform affine deformation rule:
$x_{i}^{{}^{\prime}}=x_{i}^{{}^{\prime}}+\delta\gamma\,y_{i},\,y_{i}^{{}^{\prime}}=y_{i},\,z_{i}^{{}^{\prime}}=z_{i}$
to all particles $i$. Following this step, the system is relaxed using an
energy minimization technique (FIRE or CG). This sequence of deformation
followed by minimization is thereafter continued. Here we chose the
$\delta\gamma$ such that in most of the samples we find only a single plastic
event within the incremental strain window. The maximum applied strain
$\gamma_{max}$ was chosen such that at least one plastic event occurs within
$\gamma_{max}$. In all our simulations, X is the shearing direction, Y is the
gradient direction and Z is the vorticity direction.

###  IX.3 Plastic event or Yielding

To detect the plastic event or the yielding event, we measure the largest non-
affine displacement in the system. The displacement is computed from two
minimized configurations (separated by the elementary strain $\delta\gamma$)
and defined as $\Delta
z^{2}_{max}(\gamma)=\max[z_{i}(\gamma+\delta\gamma)-z_{i}(\gamma)]^{2}$ where
$z_{i}(\gamma)$ is the coordinate of a particle $i$ in the direction
perpendicular to shear (gradient direction in 2D and vorticity direction in
3D). By studying the corresponding drops in stress and the potential energy,
we put a threshold on $\Delta z^{2}_{max}$, to identify the occurrence of the
event; $\Delta z^{2}_{max}$ is model dependent. We have also verified this
identification via the back-shear protocol [46] and we find that our method
provides an excellent balance between sensitivity of detection and
computational efficiency. More details are provided in SI.

## Acknowledgements

PC, SK, SS and VVV acknowledge financial support from NSM grant
DST/NSM/R&D_HPC_Applications/2021/29 as well as computational time on PARAM
Yukti computing facility. SK acknowledges the Swarna Jayanti Fellowship,
grants DST/SJF/PSA01/2018-19, and SB/SFJ/2019-20/05. SS acknowledge the
Faculty Initiation Grant from IIT Roorkee and PARAM Ganga computational
facility. VVV acknowledge the computational time on the Chandra cluster at IIT
Palakkad. Further we would like to thank Srikanth Sastry and Peter Sollich for
fruitful discussions.

## References

  * Bonn  _et al._ [2017] D. Bonn, M. M. Denn, L. Berthier, T. Divoux, and S. Manneville, Yield stress materials in soft condensed matter, Reviews of Modern Physics 89, 035005 (2017).
  * Varnik  _et al._ [2004] F. Varnik, L. Bocquet, and J.-L. Barrat, A study of the static yield stress in a binary lennard-jones glass, The Journal of chemical physics 120, 2788 (2004).
  * Schuh  _et al._ [2007] C. A. Schuh, T. C. Hufnagel, and U. Ramamurty, Mechanical behavior of amorphous alloys, Acta Materialia 55, 4067 (2007).
  * Rodney  _et al._ [2011] D. Rodney, A. Tanguy, and D. Vandembroucq, Modeling the mechanics of amorphous solids at different length scale and time scale, Modelling and Simulation in Materials Science and Engineering 19, 083001 (2011).
  * Barrat and Lemaitre [2011] J.-L. Barrat and A. Lemaitre, Heterogeneities in amorphous systems under shear, Dynamical heterogeneities in glasses, colloids, and granular media 150, 264 (2011).
  * Falk and Langer [2011] M. L. Falk and J. S. Langer, Deformation and failure of amorphous, solidlike materials, Annu. Rev. Condens. Matter Phys. 2, 353 (2011).
  * Fielding [2014] S. M. Fielding, Shear banding in soft glassy materials, Reports on Progress in Physics 77, 102601 (2014).
  * Nicolas  _et al._ [2018] A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat, Deformation and flow of amorphous solids: Insights from elastoplastic models, Rev. Mod. Phys. 90, 045006 (2018).
  * Argon [1979] A. Argon, Plastic deformation in metallic glasses, Acta metallurgica 27, 47 (1979).
  * Falk and Langer [1998] M. L. Falk and J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Physical Review E 57, 7192 (1998).
  * Maloney and Lemaitre [2006] C. E. Maloney and A. Lemaitre, Amorphous systems in athermal, quasistatic shear, Physical Review E 74, 016118 (2006).
  * Eshelby [1957] J. D. Eshelby, The determination of the elastic field of an ellipsoidal inclusion, and related problems, Proceedings of the royal society of London. Series A. Mathematical and physical sciences 241, 376 (1957).
  * Tanguy  _et al._ [2006] A. Tanguy, F. Leonforte, and J. L. Barrat, Plastic response of a 2d lennard-jones amorphous solid: Detailed analysis of the local rearrangements at very slow strain rate, The European Physical Journal E 20, 355 (2006).
  * Jensen  _et al._ [2014] K. Jensen, D. A. Weitz, and F. Spaepen, Local shear transformations in deformed and quiescent hard-sphere colloidal glasses, Physical Review E 90, 042305 (2014).
  * Karmakar  _et al._ [2015] S. Karmakar, C. Dasgupta, and S. Sastry, Length scales in glass-forming liquids and related systems: a review, Reports on Progress in Physics 79, 016601 (2015).
  * Jaiswal  _et al._ [2016] P. K. Jaiswal, I. Procaccia, C. Rainone, and M. Singh, Mechanical yield in amorphous solids: A first-order phase transition, Physical Review Letters 116, 085501 (2016).
  * Shrivastav  _et al._ [2016a] G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Yielding of glass under shear: A directed percolation transition precedes shear-band formation, Physical Review E 94, 042605 (2016a).
  * Divoux  _et al._ [2016] T. Divoux, M. A. Fardin, S. Manneville, and S. Lerouge, Shear banding of complex fluids, Annual Review of Fluid Mechanics 48, 81 (2016).
  * Popović  _et al._ [2018] M. Popović, T. W. de Geus, and M. Wyart, Elastoplastic description of sudden failure in athermal amorphous materials during quasistatic loading, Physical Review E 98, 040901 (2018).
  * Vasisht  _et al._ [2020] V. V. Vasisht, G. Roberts, and E. Del Gado, Emergence and persistence of flow inhomogeneities in the yielding and fluidization of dense soft solids, Physical Review E 102, 010604 (2020).
  * Ozawa  _et al._ [2018] M. Ozawa, L. Berthier, G. Biroli, A. Rosso, and G. Tarjus, Random critical point separates brittle and ductile yielding transitions in amorphous materials, Proceedings of the National Academy of Sciences 115, 6656 (2018).
  * Barlow  _et al._ [2020] H. J. Barlow, J. O. Cochran, and S. M. Fielding, Ductile and brittle yielding in thermal and athermal amorphous materials, Physical Review Letters 125, 168003 (2020).
  * Leishangthem  _et al._ [2017] P. Leishangthem, A. D. Parmar, and S. Sastry, The yielding transition in amorphous solids under oscillatory shear deformation, Nature communications 8, 14653 (2017).
  * Sastry [2021] S. Sastry, Models for the yielding behavior of amorphous solids, Physical Review Letters 126, 255501 (2021).
  * Dhont [1999] J. K. Dhont, A constitutive relation describing the shear-banding transition, Physical Review E 60, 4534 (1999).
  * Schall and Van Hecke [2010] P. Schall and M. Van Hecke, Shear bands in matter with granularity, Annual Review of Fluid Mechanics 42, 67 (2010).
  * Vasisht  _et al._ [2018a] V. V. Vasisht, M. L. Goff, K. Martens, and J.-L. Barrat, Permanent shear localization in dense disordered materials due to microscopic inertia, arXiv preprint arXiv:1812.03948 42 (2018a).
  * Weeks  _et al._ [2007] E. R. Weeks, J. C. Crocker, and D. A. Weitz, Short-and long-range correlated motion observed in colloidal glasses and liquids, Journal of Physics: Condensed Matter 19, 205131 (2007).
  * Mosayebi  _et al._ [2012] M. Mosayebi, E. Del Gado, P. Ilg, and H. C. Öttinger, Deformation of inherent structures to detect long-range correlations in supercooled liquids, The Journal of Chemical Physics 137, 024504 (2012), https://doi.org/10.1063/1.4732859 .
  * Mosayebi  _et al._ [2014] M. Mosayebi, P. Ilg, A. Widmer-Cooper, and E. Del Gado, Soft modes and nonaffine rearrangements in the inherent structures of supercooled liquids, Phys. Rev. Lett. 112, 105503 (2014).
  * Vasisht  _et al._ [2018b] V. V. Vasisht, S. K. Dutta, E. Del Gado, and D. L. Blair, Rate dependence of elementary rearrangements and spatiotemporal correlations in the 3d flow of soft solids, Physical review letters 120, 018001 (2018b).
  * Goyon  _et al._ [2008] J. Goyon, A. Colin, G. Ovarlez, A. Ajdari, and L. Bocquet, Spatial cooperativity in soft glassy flows, Nature 454, 84 (2008).
  * Bocquet  _et al._ [2009] L. Bocquet, A. Colin, and A. Ajdari, Kinetic theory of plastic flow in soft glassy materials, Physical review letters 103, 036001 (2009).
  * Jop  _et al._ [2012] P. Jop, V. Mansard, P. Chaudhuri, L. Bocquet, and A. Colin, Microscale rheology of a soft glassy material close to yielding, Physical review letters 108, 148301 (2012).
  * Biroli  _et al._ [2008] G. Biroli, J.-P. Bouchaud, A. Cavagna, T. S. Grigera, and P. Verrocchio, Thermodynamic signature of growing amorphous order in glass-forming liquids, Nature Physics 4, 771 (2008).
  * Karmakar  _et al._ [2009] S. Karmakar, C. Dasgupta, and S. Sastry, Growing length and time scales in glass-forming liquids, PNAS 106, 3675 (2009).
  * Karmakar  _et al._ [2012] S. Karmakar, E. Lerner, and I. Procaccia, Direct estimate of the static length-scale accompanying the glass transition, Physica A: Statistical Mechanics and its Applications 391, 1001 (2012).
  * Spaepen [1977] F. Spaepen, A microscopic mechanism for steady state inhomogeneous flow in metallic glasses, Acta metallurgica 25, 407 (1977).
  * Jack  _et al._ [2014] R. L. Jack, A. J. Dunleavy, and C. P. Royall, Information-theoretic measurements of coupling between structure and dynamics in glass formers, Physical review letters 113, 095703 (2014).
  * Tsamados  _et al._ [2009] M. Tsamados, A. Tanguy, C. Goldenberg, and J.-L. Barrat, Local elasticity map and plasticity in a model lennard-jones glass, Phys. Rev. E 80, 026112 (2009).
  * Rottler  _et al._ [2014] J. Rottler, S. S. Schoenholz, and A. J. Liu, Predicting plasticity with soft vibrational modes: From dislocations to glasses, Physical Review E 89, 042304 (2014).
  * Sollich [2010] P. Sollich, Lecture note, school on glass formers and glasses, bengaluru 2010 (2010).
  * Mizuno  _et al._ [2013] H. Mizuno, S. Mossa, and J.-L. Barrat, Measuring spatial distribution of the local elastic modulus in glasses, Phys. Rev. E 87, 042306 (2013).
  * Puosi  _et al._ [2015] F. Puosi, J. Olivier, and K. Martens, Probing relevant ingredients in mean-field approaches for the athermal rheology of yield stress materials, Soft Matter 11, 7639 (2015).
  * Patinet  _et al._ [2016] S. Patinet, D. Vandembroucq, and M. L. Falk, Connecting local yield stresses with plastic activity in amorphous solids, 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, Local yield stress statistics in model amorphous solids, Physical Review E 97, 033001 (2018).
  * 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, Predicting plasticity in disordered solids from structural indicators, Phys. Rev. Mater. 4, 113609 (2020).
  * Vasisht  _et al._ [2022] V. V. Vasisht, P. Chaudhuri, and K. Martens, Residual stress in athermal soft disordered solids: insights from microscopic and mesoscale models, Soft Matter 18, 6426 (2022).
  * Liu  _et al._ [2021] C. Liu, S. Dutta, P. Chaudhuri, and K. Martens, Elastoplastic approach based on microscopic insights for the steady state and transient dynamics of sheared disordered solids, Physical Review Letters 126, 138005 (2021).
  * Karmakar  _et al._ [2010] S. Karmakar, E. Lerner, and I. Procaccia, Statistical physics of the yielding transition in amorphous solids, Physical Review E 82, 055103 (2010).
  * Lerner  _et al._ [2018] E. Lerner, I. Procaccia, C. Rainone, and M. Singh, Protocol dependence of plasticity in ultrastable amorphous solids, Physical Review E 98, 063001 (2018).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales, Comp. Phys. Comm. 271, 108171 (2022).
  * Moorcroft  _et al._ [2011] R. L. Moorcroft, M. E. Cates, and S. M. Fielding, Age-dependent transient shear banding in soft glasses, Physical review letters 106, 055502 (2011).
  * Shrivastav  _et al._ [2016b] G. P. Shrivastav, P. Chaudhuri, and J. Horbach, Heterogeneous dynamics during yielding of glasses: Effect of aging, Journal of Rheology 60, 835 (2016b).
  * Bhaumik  _et al._ [2021] H. Bhaumik, G. Foffi, and S. Sastry, The role of annealing in determining the yielding behavior of glasses under cyclic shear deformation, Proceedings of the National Academy of Sciences 118 (2021).
  * Richard  _et al._ [2021a] D. Richard, E. Lerner, and E. Bouchbinder, Brittle-to-ductile transitions in glasses: Roles of soft defects and loading geometry, MRS Bulletin 46, 902 (2021a).
  * Note [1] Around $250$ samples were prepared at infinite to finite but higher $\Gamma$, $100$ samples at moderate $\Gamma$, and $50$ samples at the lowest $\Gamma$.
  * Note [2] We find that in 2D systems, the fluctuations dominate the local load curve and hence noisy for $N=10K$. Such issues are minimized when we study 3D systems.
  * Jiang and Murthy [2011] R. Jiang and D. Murthy, A study of weibull shape parameter: Properties and significance, Reliability Engineering & System Safety 96, 1619 (2011).
  * Ruan  _et al._ [2022] D. Ruan, S. Patinet, and M. L. Falk, Predicting plastic events and quantifying the local yield surface in 3d model glasses, Journal of the Mechanics and Physics of Solids 158, 104671 (2022).
  * Lemaître  _et al._ [2021] A. Lemaître, C. Mondal, M. Moshe, I. Procaccia, S. Roy, and K. Screiber-Re’em, Anomalous elasticity and plastic screening in amorphous solids, Physical Review E 104, 024904 (2021).
  * Biroli  _et al._ [2013] G. Biroli, S. Karmakar, and I. Procaccia, Comparison of static length scales characterizing the glass transition, Physical review letters 111, 165701 (2013).
  * Richard  _et al._ [2021b] D. Richard, G. Kapteijns, J. A. Giannini, M. L. Manning, and E. Lerner, Simple and broadly applicable definition of shear transformation zones, Physical Review Letters 126, 015501 (2021b).
  * Richard  _et al._ [2023] D. Richard, G. Kapteijns, and E. Lerner, Detecting low-energy quasilocalized excitations in computer glasses, arXiv preprint arXiv:2303.12887 (2023).
  * Brüning  _et al._ [2008] R. Brüning, D. A. St-Onge, S. Patterson, and W. Kob, Glass transitions in one-, two-, three-, and four-dimensional binary lennard-jones systems, Journal of Physics: Condensed Matter 21, 035117 (2008).
  * Vasisht and Del Gado [2020] V. V. Vasisht and E. Del Gado, Computational study of transient shear banding in soft jammed solids, Physical Review E 102, 012603 (2020).
  * Saika-Voivod  _et al._ [2004] I. Saika-Voivod, F. Sciortino, and P. H. Poole, Free energy and configurational entropy of liquid silica: Fragile-to-strong crossover and polyamorphism, Phys. Rev. E 69, 041503 (2004).
  * Mendelev  _et al._ [2019] M. Mendelev, Y. Sun, F. Zhang, C.-Z. Wang, and K.-M. Ho, Development of a semi-empirical potential suitable for molecular dynamics simulation of vitrification in cu-zr alloys, The Journal of Chemical Physics 151, 214502 (2019).
  * Maloney and Lemaître [2006] C. E. Maloney and A. Lemaître, Amorphous systems in athermal, quasistatic shear, Phys. Rev. E 74, 016118 (2006).

## Supplementary Information

Figure 6: Evolution with strain Stress $\sigma$ and potential energy U evolve
with applied strain. Jumps in these quantities observed as the system show
local yielding due to plastic events, hence non-affine displacement captured
by maximum displacement measure. This data is for 2D BMLJ model and poorely
aged sample

### Detecting plastic events

In Fig. 6 we show the evolution of the system’s shear stress ($\sigma$),
potential energy (U) as well as the Maximum displacement with applied strain
$\gamma$, for 2D BMLJ. Here we have considered the displacement in Y since the
shearing direction is along X. In 3D systems, we consider Z, the vorticity
direction. Clearly the jumps in the $\sigma$ and U is signed in term of kinks
in the displacement. In order to differentiate the signal from the background
noise, we compute the distribution of these displacements and is shown in Fig.
7. It is quite evident that we have a bimodal distribution and we take
conservative estimate for the threshold to recognize the plastic event as
$0.01$. This threshold is robust for varying ages of the sample. We carry out
this analysis for other models to obtain the appropriate thresholds.

Figure 7: Maximum displacement distribution Top plot show the raw data of
maximum displacement and bottom its distribution. This data is for 2D BMLJ
model and poorely aged sample

Plastic events in subsystem with elastic background :  We perform the above
analysis in the soft matrix as well wherein only a subsystem is allowed to
have plastic events. Even though the bimodality of the displacement
distribution is coarsened (see Fig. 8), the threshold for the plastic event is
quite evident and recognize the events accurately (see Fig. 9).

### Choice of k value in soft matrix

In the soft matrix method we allow the sub-system to relax completely without
any constraint and the background is allowed relax affinely. This is achieved
by addition of a spring potential with a sprint constant k. The choice of k is
made such that (1.) no large displacement occur in the background but also
(2.) do not alter the system’s mechanical properties. In the Fig. 10 we show
(for 2D BMLJ) the maximum displacement measured in the subsystem (left) and
the full system (right), including the subsystem. For k=2 and k=5, we find
that relative large displacement are occurring in the surround and not in the
subsystem. But with the increase the k, these events in the surrounding are
curbed. In the Fig. 11 we show curvature of the PE profile (vs. $\gamma$)
before the plastic event, the storage module (from the slope of $\sigma$ vs
$\gamma$ before the plastic event) as well as the pressure difference $\Delta
P=P_{k}-P_{k=0}$ for varying k values. These quantities are computed for the
full system with soft matrix background. We find that for a range of k values
(k¡=20), these measures do change and are simliar to bulk value (dashed line
in Fig. 11). Hence $k=10$ is a good estimate which satisfy both the criterion
defined above. Similar protocols are followed for other models.

Figure 8: Maximum displacement distribution in soft matrix Top plot show the
raw data of maximum displacement and bottom its distribution, within the
subsystem. This data is for 2D BMLJ model and poorely aged sample with a
subsystem size of $L_{s}=16$. Figure 9: Evolution with strain in the subsystem
Stress $\sigma$ and potential energy U evolve with applied strain. Jumps in
these quantities observed as the system show local yielding due to plastic
events in the sub-system. This data is for 2D BMLJ model (poorely aged) in a
soft matrix of k=10 and $L_{s}=16$. Figure 10: k dependence on events inside
and outside the sub-system Maximum displacement measure discussed in the text
is computed for various k values, in the subsystem (left) and the full system
(right) as a function of strain. The plastic event is signalled by the jump in
the measure. Increasing the value of k, events occur only within the
subsystem. This is for the 2D BMLJ system.

Figure 11: Various measures of the full system with a soft matrix background
(2D BMLJ) For $k\leq 20$, these measures do not vary substantially. The dotted
line is the measure for unconstrained (without soft matrix) full system. Here
$\Delta P=P_{k}-P_{k=0}$ and is measure at small strain of $10\delta\gamma$

.