## Residual stress distributions in athermally deformed amorphous solids from atomistic simulations

C´ eline Ruscher 1,2 and J¨rg Rottler 1 o

1

Department of Physics and Astronomy and Quantum Matter Institute, University of British Columbia, Vancouver BC V6T 1Z1, Canada 2 Institut Charles Sadron, 23 rue du Loess, F-67034 Strasbourg, France

The distribution of local residual stresses (threshold to instability) that controls the statistical properties of plastic flow in athermal amorphous solids is examined with an atomistic simulation technique. For quiescent configurations, the distribution has a pseudogap (power-law) form with an exponent that agrees well with global yielding statistics. As soon as deformation sets in, the pseudogap region gives way to a system size dependent plateau at small residual stresses that can be understood from the statistics of local residual stress differences between plastic events. Results further suggest that the local yield stress in amorphous solids changes even if the given region does not participate in plastic activity.

## I. INTRODUCTION

Linking atomic scale observables quantitatively to parameters in effective, coarse-grained descriptions is an important challenge in materials physics. Amorphous materials might appear homogeneous beyond local packing effects, but exhibit structural heterogeneity at the nanoscale [1]. The shear modulus, for instance, is not uniform everywhere but exhibits Gaussian fluctuations on the scale of ∼ 10 particle diameters that have their origin in the strong nonaffine atomic displacement field of disordered packings already in the elastic regime [2].

gions that exhibit a particularly low local yield stress in the given shear direction. Despite the artefacts potentially induced by the frozen constraint, key information is thus revealed about the distribution of energy barriers that control shear rearrangements. This situates the FM method well within the broader effort to predict the location of irreversible (dynamical) rearrangements from structural (static) features of amorphous packings [1114].

While local linear properties such as moduli can be probed without an external perturbation through the use of fluctuation formulas [3], no such options exists for nonlinear, plastic properties. Ideally one would like to induce plastic yielding in a small local region of interest while eliminating plasticity everywhere else. The only computational approach that has been proposed to date to accomplish this task consists in the (rather drastic) step of enforcing fully affine deformation everywhere but in the region of interest, where particles can move freely [4]. As a result, the nonaffine displacement field (with correlations extending over many particle diameters) is truncated at the boundary between the probe region and the 'frozen matrix' (FM). At the linear response level, it is well documented that this truncation shifts the average shear modulus to higher values and reduces the variance of the Gaussian fluctuations in comparison to a fully unconstrained approach [5]. However, the method still faithfully detects about the distribution of soft vs. hard regions in the amorphous solid.

For this reason, the FM technique has recently been applied to study the distribution of local yield stresses σ Y in amorphous solids [6, 7]. Patinet et al. have constructed local yield stress maps of two-dimensional amorphous mixtures after a quench from the liquid phase [8-10]. They found that when these solids are sheared quasistatically , the first plastic events are indeed occurring in re-

Motivated by this success, we seek to explore in this contribution the potential of the FM method to reveal robust statistical information of local mechanical observables in quiescent and flowing amorphous solids. Our focus here is not the yield stress per se, but instead the residual stress x = σ Y -σ 0 , which is the difference between the local yield stress σ Y and the local stress σ 0 and is a measure of how far a local region is from instability . The distribution P x ( ) plays a central role in contemporary descriptions of the statistical properties of the yielding transition, because its behavior at small arguments controls the statistics of macroscopic slip events (avalanches) under athermal quasistatic deformation. The distribution is assumed to be scale-free (i.e has a power-law form) and to vanish at zero, i.e. P x ( ) ∼ x θ as x → 0 and the pseudogap exponent θ enters various scaling relations linking critical exponents of the yielding transition.

In the following, we first present results for P x ( ) from the FM method applied to the quenched state of a 2D model amorphous solid. We show that a pseudogap form is indeed obtained with an exponent that agrees with one inferred from macroscopic deformation. We then extend the analysis to the transient and steady state of quasistatic deformation, where the distribution P x ( ) develops a system size dependent plateau in the small x region. By studying the distribution of stress increments P (∆ ) x between slip events, where ∆ x = ∆ σ Y -∆ σ 0 , we show that this plateau can be attributed to the discrete increments of the underlying mechanical noise. We argue that the FM method reveals correct generic trends

in the residual stress distribution P x ( ) , but a quantitative prediction of the pseudogap exponent under deformation requires larger system sizes that are currently computationally not accessible due to combined effects of frozen boundary and finite size artefacts.

## II. SIMULATION METHODS

## A. System

We consider a 2D Lennard-Jones (LJ) glass-forming binary mixture which has been introduced by Lancon ¸ et al. [15] originally to investigate the properties of 2D quasicrystals. Following refs. [9, 16], the N L large and N S small particles interact through the potential:

$$\text{man requires unrecognized} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{out} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{In} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{ing} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{is} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{ini} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{ in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{ine} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{\text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in} \text{in}

\text{.\tilde { \Delta } \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Pho \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi
\text{.\tilde { \Delta } \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Psi \Phi \Phi \Phi \Psi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi\Phi \Phi \Phi \Phi\Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phis \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \ Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \Phi \$$

where { a, b } = { L, S } and r is the distance between two particles. The potential is shifted at the cutoff distance r cut = 2 5 . σ LS and smoothed for r in &lt; r ≤ r cut where r in = 2 0 . σ LS in order to ensure that U ab ( r ) is twice differentiable. The shift in energy U S and the coefficients C k are:

$$C _ { k } & \text{ are:} \\ U _ { S } & = C _ { 0 } - 4 \varepsilon _ { a b } \left [ \left ( \frac { \sigma _ { a b } } { r _ { i n } } \right ) ^ { 1 2 } - \left ( \frac { \sigma _ { a b } } { r _ { i n } } \right ) ^ { 6 } \right ] _ { \quad \text{(2)} } \\ C _ { 0 } & = - ( r _ { c u t } - r _ { i n } ) [ 3 C _ { 1 } + C _ { 2 } ( r _ { c u t } - r _ { i n } ) ] / 6 \quad & \text{cular} \\ C _ { 1 } & = 2 4 \varepsilon _ { a b } \sigma _ { a b } ^ { 6 } ( r _ { i n } ^ { 6 } - 2 \sigma ^ { 6 } ) / r _ { i n } ^ { 1 3 } \\ C _ { 2 } & = - 1 2 \varepsilon _ { a b } \sigma _ { a b } ^ { 6 } ( 7 r _ { i n } ^ { 6 } - 2 6 \sigma ^ { 6 } ) / r _ { i n } ^ { 1 4 } \quad & \text{config} \\ C _ { 3 } & = - [ 3 C _ { 1 } + 4 C _ { 2 } ( r _ { c u t } - r _ { i n } ) ] / [ 3 ( r _ { c u t } - r _ { i n } ) ^ { 2 } ] \quad & \text{there} \\ C _ { 4 } & = [ C _ { 1 } + C _ { 2 } ( r _ { c u t } - r _ { i n } ) ] / [ 3 ( r _ { c u t } - r _ { i n } ) ^ { 3 } ] \quad & \text{done} \\ \text{The different } L J \text{ parameters are } \sigma _ { L L } \ = \ 2 \sin ( \pi / 5 ),$$

The different LJ parameters are σ LL = 2sin( π/ 5) , σ LS = 1 , σ SS = 2 sin( π/ 10) , ε LL = ε SS = 0 5 . , ε LS = 1 , and all masses are set to m = 1 . The ratio between large and small particles is chosen such as N /N L S = (1 + √ 5) / 4 , and we work at constant density N/V = 1 0206 . . In what follows, the length, mass, energy and time units are expressed in term of σ LS , m ε , LS and σ LS √ m/ε LS respectively. For this system, the glass transition temperature is T g = 0 325 . ε LS /k B where k B is Boltzmann's constant.

All simulations have been carried out using the LAMMPS software [17]. We consider a 2D triclinic simulation box of size L under periodic boundary conditions where L ∈ [53 , 300] in order to probe system size scaling. To generate the different glass configurations, we first equilibrate systems in the liquid phase at T = 2 T g using the Langevin thermostat with a damping parameter

FIG. 1. Top: Illustration of the frozen matrix method. A circular region of size R is selected from a global configuration and a shell of frozen particles of size 2 r cut is kept in addition. Bottom: Example of a quenched (applied strain γ = 0) global configuration of size L = 100 . The circular regions onto which the frozen matrix is applied, are represented with darker colors.The regions do not interpenetrate to ensure their independence during FM analysis.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-30b279aa-087f-4176-8bb4-461fea037082.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-f8d0169b-25e8-4c7f-a718-044a28b87095.jpg)

T damp = 1 0 . . The timestep is chosen as δt = 0 005 . . After equilibration, the different configurations are cooled down to T = 0 at dT/dt = 2 · 10 -3 . Finally an energy minimization is perfomed to ensure that the system is in its local minimum.

## B. Athermal quasistatic shear

Subsequent to cooling, the systems are deformed following the athermal quasistatic shear (AQS) protocol []. Simple shear is applied in the following way: an affine deformation is first performed by tilting the simulation box in the x direction by an amount δγ xy L , where δγ xy = 5 · 10 -5 is the strain increment, and then remapping the position of the particle inside the deformed box.

In a second step, we allow the system to relax through an energy minimization using the conjugate gradient method.

## C. Frozen matrix method

To probe the local properties of our samples, we consider circular regions of size R embedded into a frozen shell of size 2 r cut as illustrated in Figure 1 (top). During an AQS step, only the inner region can relax implying that if a plastic event occurs, it is necessarily located inside the circular region. In the following, we refer to deformations of the whole periodic simulation box as global and to deformations of the circular regions as local .

In practice, starting from a global configuration obtained at a given applied strain γ , we select regions, as shown in Figure 1 (Bottom), that we deform under simple shear. By monitoring for each site region the evolution of stress as a function of the applied strain, we can access σ Y and γ Y , the local yield stress and yield strain, respectively. Knowing the initial local stress σ 0 of the circular region, we can determine the residual stresses x associated with each local region.

## D. Detection of plastic events

To detect plastic events, we revisit an energy-based criterion introduced in ref. [18] that suits perfectly the AQS protocol. The observable κ = ( U aff -U 0 ) / Nδγ ( 2 ) measures the mismatch between the energy associated with the affine displacement U aff and the inherent structure U 0 . As we show in the Appendix, it is possible to determine analytically an upper bound to κ ≤ G / B (2 ρ ) below which system behaves only elastically. G B is the Born shear modulus. As it is associated to affine deformation, G B has the same value for both global and local deformations and it is therefore a convenient criterion that can be applied to determine plastic activity in both cases. We numerically found that G B ≈ 26 , therefore in this work, we consider events that have κ ≥ 30 as plastic.

## III. DISTRIBUTION OF RESIDUAL STRESSES

In order to assess the ability of the FM method to accurately obtain the distribution P x ( ) , we focus first on the quenched state. We prepare N g = 5 · 10 4 quenched global configurations of size L = 100 , and from these N g configurations, we select N glyph[lscript] ≥ 2 · 10 5 independent local sites onto which we apply the FM method by deforming these local sites up to the first plastic event. In order to probe the influence of the size R of the local region, we consider R ∈ [5 0; 20 0] . . .

Additionally, we deform all N g configurations globally without any constraints and record the global residual stress x min of the first time failure is observed as well

FIG. 2. Comparison between the distribution of weakest sites P y ( min ) ( glyph[star] ) computed from the global configurations and the distribution of P y ( ) (open symbols) obtained by the local FM method for circular regions of different sizes R . The top panel shows results for the residual stresses y = x while the middle panel shows the yield strains y = γ Y . In both panels, the dashed lines correspond to Weibull fits to the global data. The bottom panel shows the scaling of 〈 x min 〉 and 〈 γ Y,min 〉 with system size L .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-e0c75602-2e1b-44d8-b48a-89f1b3551c8e.jpg)

as the associated global yield strain γ Y . In the AQS protocol, the first failure event represents by definition the weakest site in the system. According to extreme value statistics, if the underlying distribution of the observable y = { x, γ Y } is power-law distributed at small arguments, P y ( ) ∼ y θ , then the distribution of independent minimal values y min sampled from P y ( ) is expected to follow a Weibull distribution: [19]

$$\intertext { a \, \omega } _ { \, } \, & P ( y _ { m i n } ) = \frac { 1 + \theta } { \langle y _ { m i n } \rangle } \left ( \frac { y _ { m i n } } { \langle y _ { m i n } \rangle } \right ) ^ { \theta } \exp \left [ - \left ( \frac { y _ { m i n } } { \langle y _ { m i n } \rangle } \right ) ^ { 1 + \theta } \right ] \\ \intertext { l u a l } _ { \, \dots } \,.$$

FIG. 3. Maps of residual stress x , local yield stress σ Y , and yield strain γ Y upon deformation for an applied strain γ = 0 (Top), γ = 0 03 . (Middle) and γ = 0 18 . (Bottom).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-ca5329ba-84bd-4cf9-b0a9-0d65fc9e62f4.jpg)

In Figure 2 we show P x ( ) and P x ( min ) (Top) and P γ ( Y ) and P γ ( Y,min ) (Middle). We first observe that the Weibull distribution, when fitted to our data, describes well the distributions P x ( min ) and P γ ( Y,min ) coming from the global configurations. However, we notice that the Weibull exponent θ = 0 37 . is smaller in the case of weakest residual stresses than in the case of smallest yield strain, where we obtain θ = 0 52 . . We attribute the difference to the limited statistics obtained for both quantities in the very small x min and γ Y,min regions.

important for small local regions but, as R increases, relaxation becomes easier and the two power law regimes merge into one unique regime when R ≥ 10 .

Inspection of the distributions P y R ( ; ) reveals the effects of size R and of the constraint imposed by the frozen shell. For both residual stresses and yield strains, increasing R induces a shift of the distribution toward smaller values of y . However, the distributions show a pseudogap form at small x and rapidly decaying part at larger x . For R &lt; 10 , we notice the presence of two power-law regimes. One regime at very small values exhibits an exponent θ 1 = 0 58 . for P x ( ) which is larger than the exponent θ = 0 37 . extracted from P x ( min ) . For the yield strain, however, θ 1 = 0 52 . is in perfect agreement with θ extracted from the Weibull fit of P γ ( Y,min ) and also in good agreement with θ 1 from the residual stress.

For larger values of y and R = 5 00 . , there appears to be a transition to a steeper power-law regime where θ 2 = 1 3 . . This second regime is a consequence of the frozen boundary, which prevents nonaffine relaxations outside the circular region and therefore makes the rearrangements more difficult. This effect is particularly

For the weakest sites which are on the verge of yielding, the rearrangements do not require a substantial deformation and the constraints induced by FM do not affect their relaxation. Therefore, the distributions obtained in small x or small γ Y regions show what we believe to be the right pseudogap form. To confirm our hypothesis, we look at the scaling of 〈 x min 〉 and 〈 γ Y,min 〉 with the system size in the global unconstrained configurations. Results are shown in the bottom panel of Figure 2, where we see that 〈 y min 〉 ∼ L -α with α = 1 24 . in the case of 〈 x min 〉 and α = 1 21 . in the case of 〈 γ Y,min 〉 . Using extreme value statistics, θ = 2 /α -1 (see also eq. (5)), we thus find θ = 0 61 . (resp. θ = 0 65 . ) for residual stresses (resp. yield strains). These values are in good agreement with θ 1 extracted from the distributions obtained with FM. The results obtained in the quenched state are in agreement with what has been reported by other atomistic simulations that found θ ≈ 0 6 . for deformations on global systems [19, 20]. They are also coherent with results obtained from elastoplastic mesoscopic and meanfield models, for which P x ( ) can be determined easily, which find θ ≈ 0 6 . [21] and θ ≈ 0 5 . [22] respectively.

As the FM method is able to capture correctly the pseudogap exponent in the quenched state, we can use it to probe what occurs during deformation. As shown before, the shape of P x ( ) is strongly influenced by the size of the circular region. One possibility is to consider local re-

gions of R ≥ 10 to avoid dealing with the two power law regimes, but then each local region would include more than 300 particles. Recent work from Barbot et al. shows a good correlation between local yield stress and plastic activity for R = 5 0 ( . ∼ 80 particles) and thus argues for using a smaller size of system to correctly capture plastic activity . This value also coincides with the lower bound for which Hooke's law still applies [1]. For these reasons, in what follows, we set the size of the frozen region to R = 5 0 . .

## IV. EVOLUTION UPON DEFORMATION

## A. Local plastic observables

In order to probe the effect of deformation on the local plastic observables x σ , Y and γ Y , we consider N g different global configurations of size L that we deform up to 20% . Snapshots are saved every 1% of deformation in the transient regime and every 2% of deformation in the stationary regime. As for the quenched state, we consider N glyph[lscript] independent local regions of size R = 5 0 . selected from global configurations saved at given applied strain

γ

.

Examples of maps of these quantities for several values of applied strain can be seen in Figure 3. As also found by Patinet et al. , [8, 9], heterogeneities are well-visible in the quenched state for residual stress, yield stress and yield strain. We see the presence of weak spots embedded into an harder medium. The maps of the residual stress and yield strain are very similar due to the nature of the loading protocol. Indeed less deformation is required to make weak sites fails, so sites with small x are also sites with small γ Y .

Upon deformation, we notice that the number of sites with very weak values of σ Y decreases significantly pointing toward a homogenization of the system in the stationary regime. This observation is also visible in the top panel of Figure 4 where the probability of the local yield stress, P σ ( Y ) , is smaller at small σ Y in the transient and stationary regimes than in the quenched state. Moreover, there appears to be a tendency of hard spots to become softer in the transient regime. Maps in the transient regime show the larger number of weak spots as we also observe in the probability distribution function (pdf) of the yield strain P γ ( Y ) displayed in the bottom panel of Figure 4. Indeed the fraction of small γ Y is larger in the transient regime than in the quenched and stationary regime which have almost the same fraction of small γ Y .

## B. Plateaus and mechanical noise

The evolution of residual stresses upon deformation is shown in Figure 5. At the earliest stages of deformation, we find P x ( ) ∼ x θ in the small x region and the pseudogap exponent θ is decreasing. After only 3%

FIG. 4. Probability distributions of the yield stress, σ Y (Top) and of the yield strain, γ Y (Bottom) for three applied strains in the quench, transient and stationary regimes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-99fdb6ef-b006-469a-bde6-c0c32921ffc4.jpg)

FIG. 5. Evolution of the distribution of residual stresses P x ( ) with increasing applied strain γ = 0 0 01 0 02 0 04 0 06 0 08 0 12 0 18 , , , . , . , . , . , . , . (from top to bottom) in systems of size L = 300 . The triangles represent the location of 〈 x min 〉 obtained from global deformation of the configurations. The different curves have been shifted in the upper direction for sake of clarity.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-d21d36e4-40c1-43da-b4c7-f2c6a0e8e9ab.jpg)

strain, we notice the appearance of a plateau which remains visible for larger applied strains. Such deviations from the pure pseudogap form of P x ( ) in the steady state have recently been reported in two different works that employed coarse-grained idealized elastoplastic models (EPM) [23, 24]. These observations could be described by rewriting the pdf of residual stresses as P x ( ) ∼ p 0 + x θ , where p 0 ∼ L -p . A system size dependence of the plateau value of P x ( ) also appears in our atomistic simulations as reported in Figure 6, where P x ( ) are shown

FIG. 6. Effect of the system size L on P x ( ) in the transient (Top) and stationary regime (Bottom). The triangles indicate the location of 〈 x min 〉 while the circles stand for the location of the lower cutoff of the mechanical noise ∆ x , see Fig. 7. The insets shows the scaling with L of the plateau p 0 and the crossover x c from the power law regime to the plateau.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-514bad44-daf3-41fa-bdf8-a77ea1e84fd0.jpg)

for different sizes L and for two given applied strains, γ = 0 03 . (transient) and γ = 0 18 . (stationary). In order to investigate the size scaling of the plateau value p 0 (see inset), we average the value of P x ( ) for x ≤ 10 -3 and find the plateau p 0 ∼ L -0 15 . in the transient regime and p 0 ∼ L -0 27 . in the steady-state regime. The values of these exponents are significantly lower than those reported in the EPM studies [23, 24]. In qualitative agreement with these studies, however, we observe that when the system size L increases, the crossover from the power law region to the plateau occurs at smaller values of x . We estimate the value x c at which the crossover takes place by looking at the intersection between the plateau and the power law region. As shown in the inset of Figure 6, in the stationary regime we find that x c ∼ L -0 78 . . This value is in agreement with EPM results for which the exponent has been estimated between 0.73-0.95 [23, 24].

The origin of the plateau can be attributed to the discreteness of the underlying stochastic process as suggested by Zoia et al. [25] in the context of fractional Brownian motion. The dynamics of the variable x can be viewed as a random walk in the presence of an absorbing boundary at x = 0 . One expects deviations from scalefree behavior of the steady state distribution P x ( ) near that boundary when x becomes of the order of the typ-

FIG. 7. Distribution of the changes in local stress ∆ σ 0 , yield stress ∆ σ Y and residual stress ∆ x between consecutive avalanches in the transient (left) and stationary (right) regimes. The bottom panel shows the system size scaling of the lower cutoffs of these distributions.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-4e7eafec-818b-4e21-b6ea-4602c5006fee.jpg)

ical increment of the mechanical noise coming from the stochastic redistribution of the stress during plastic rearrangements. The plateau is therefore a finite size effect that vanishes in the thermodynamic limit. The distribution of stress kicks ∆ σ 0 , which monitors the variation of stress between two consecutive plastic rearrangements, is known to be broadly distributed and can be described

FIG. 8. Upper panel: System size scaling of 〈 x min 〉 and ∆ σ in the transient (Top) and stationary (Bottom) regimes. Lower panel: Examples of fits performed on P x ( ) and cd f ( x ) in the transient and stationary state for L = 300 . The filled symbols show the range of x values for which x &gt; ∆ x c .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-dd27f004-e302-44bf-9cb8-79f2e551df98.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-ce92743c-f3c7-40ff-85cf-ea217f4e3cbd.jpg)

by a truncated Pareto distribution [22],

$$\rho ( \Delta \sigma _ { 0 } ) = \frac { \mu \Delta \sigma _ { 0, c } ^ { \mu } } { | \Delta \sigma _ { 0 } | ^ { 1 + \mu } } \, \quad \ \ \ \ ( 4 ) \, \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

where ∆ σ 0 ,c ∼ L -2 /µ is the lower cutoff of the distribution which depends on the system size and sets a characteristic scale. The distribution also has an upper (system size independent) cutoff which comes from the variation of stress felt by adjacent sites. The exponent µ is related to the elastic interaction kernel G ( r ) ∼ r -d/µ . For µ = 1 , we recover the Eshelby kernel [22] that can be well fit to elementary shear transformations observed in atomistic glass models [26].

The AQS scheme allows to resolve individual avalanches but cannot discriminate between single and multiple plastic events inside a given avalanche. Therefore, noise kicks are defined between the end points of consecutive avalanches. In the transient and the stationary regimes, we record all stress drops in an applied strain interval of 2% . From our data we can determine not only the distribution of stress kicks P (∆ σ 0 ) between two consecutive avalanches, but also the distribution of the variation of yield stress P (∆ σ Y ) and the distribution of the variation of the residual stress P (∆ ) x for the sites that did not fail during theses avalanches. The results are shown in Figure 7, where we observe that the stress kicks obey P ( ∆ | σ 0 | ) ∼ | ∆ σ 0 | -2 5 . with a system size dependent lower cutoff (defined by the intersection of the flat region of the distributions with the power law region) that varies as L -1 2 . (Bottom panel). The scaling of the lower cutoff implies a slightly overestimated value of µ = 2 1 2 / . ≈ 1 67 . , that could be related to the fact that we estimate the scaling on small interval [27]. The small peaks in the flat region come from the fixed strain increment as a typical step generates a stress ∼ Gδγ ≈ 1 2 . · 10 -3 . The data is thus consistent with eq. (4) with µ = 1 5 . , which is a larger exponent than what is expected for an Eshelby kernel. Within the AQS protocol, we can only probe stress changes between entire avalanches and not single plastic rearrangements. Therefore, the stress kicks that we measure are not coming from spatially localized events but rather from extended ones, for which µ &gt; 1 is expected [28].

Interestingly, the middle panel of Fig. 7 shows that the local yield stress also changes even if the region did not undergo a shear rearrangement. The atomistic simulations are thus at variance with EPMs, which usually assume ∆ σ Y = 0 , and the distribution of kicks P (∆ σ 0 ) is equivalent to the distribution P (∆ ) x . In our system, the yield strain on stable sites varies between two consecutive rearrangements as P ( ∆ | σ Y | ) ∼ | ∆ σ Y | -1 8 . . Even if the local environment of particles in terms of nearest neighbors does not change, long ranged nonaffine displacements that occur (even during the elastic branches of loading) perturb the local environment. Already small nonaffine changes in the atomic positions are likely to be sufficient to slightly modify the local yield stress.

Coming now to the distribution of the variation of the residual stress itself, the bottom panel of Fig. 7 shows that P ( ∆ | x | ) ∼ | ∆ x | -2 , implying µ = 1 and therefore a compatibility with the Eshelby kernel. However, we believe this value to be coincidental. Moreover, the system size scaling of the lower cutoff ∆ x c ∼ L -1 05 . is not consistent with eq. (4). We attribute these differences to the contributions ∆ σ Y from the changing yield stress to ∆ x . We also note that the scaling of ∆ x c is reasoably close to the scaling that we found for x c .

The location of the lower cutoff ∆ x c on the distribution P x ( ) is indicated with open circles in Figure 6. We observe that for the largest system sizes L = 200 and L = 300 , ∆ x c marks the entrance of the plateau region, while for L = 100 , it appears to be still in the power law region. Given the limitations of accessible system sizes with the atomistic simulations, these observations are consistent with the notion expressed above that ∆ x c is the relevant discretization scale below which deviations from the pseudogap form of P x ( ) become visible.

## C. Scaling relations

The weakest residual stresses are controlling of the flow of amorphous solids, which in the athermal quasistatic limit consists of periods of elastic loading punc-

FIG. 9. Evolution of the exponent α with applied strain γ computed from 4 different methods: scaling of the stress drops, α S , scaling of the weakest sites, α 〈 x min 〉 , from the pseudogap exponent θ extracted from P x ( ) computed from the FM, α FM , and from the scaling of 〈 x min 〉 obtained from the inversion of the cdf, α I .

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-c8b03291-d048-44da-9307-c53ea7a7b539.jpg)

tuated by sudden stress drops. The magnitude of these stress drops decreases with increasing system size as 〈 ∆ σ 〉 ∼ L -α S . In the upper panel of Fig. 8 we compute the average size of the stress drops induced by global deformation vs. L and find a steady state value α S = 1 27 . in very good agreement with previous atomistic simulations [29, 30]. In the transient regime, a larger value α S = 1 55 . is found. Also shown is the behavior of 〈 x min 〉 in these simulations, which also follows a scaling form 〈 x min 〉 ∼ L -α 〈 x min 〉 . If one assumes that on average loads and releases compensate, one expects 〈 ∆ σ 〉 ∼ 〈 x min 〉 . Figure 9 reports that these two exponents agree well with each other for the entire range of deformation strain. As soon as deformation sets in, the exponents jump from a (preparation dependent) initial value near 1.2 to a larger value of 1.7 and then decrease during the transient loading phase to a steady state value 1.27. Similar trends have been observed in prior atomistic simulation studies [20, 31].

The mean value of the weakest site 〈 x min 〉 can be computed from the underlying distribution of residual stresses via a well-known expression from extreme value statistics [19],

$$\int _ { 0 } ^ { \langle x _ { m i n } \rangle } P ( x ) d x = \frac { c } { L ^ { d } } \ \Rightarrow \langle x _ { m i n } \rangle \sim L ^ { - \alpha _ { \langle x _ { m i n } } \rangle } \quad ( 5 ) \quad \text{ for tb} \\ \text{cover} \\ 9. \ w \epsilon$$

From eq. (5) we see that depending on the form of P x ( ) , α 〈 x min 〉 can have different expressions. For instance, if we assume that x ≤ 〈 x min 〉 is inside the plateau region then P x ( ) ∼ L -p and 〈 x min 〉 ∼ L -( d -p ) . However, if we assume that the pseudogap description P x ( ) ∼ x θ is valid then 〈 x min 〉 ∼ L -d/ (1+ ) θ .

In what follows, we seek to ascertain if the distribution P x ( ) that we obtained with the FM method is compatible with the exponents α S and α 〈 x min 〉 that were obtained from unconstrained global AQS deformations. In order to distinguish from these exponents, we denote by α FM

the value of the exponents obtained from FM. From the observations made in Figure 6, we see that 〈 x min 〉 is located in the plateau region, which suggests α FM = d -p . In the transient regime α FM = 1 85 . while in the stationary regime α FM = 1 73 . . These two values are obviously too large to be in agreement with what has been measured from global configurations.

Alternatively, we might assume that the global scaling on 〈 x min 〉 reflects the thermodynamic limit, where the plateau is irrelevant. A careful inspection of P x ( ) reveals that the crossover between the two power law regimes observed in Figure 2 remains upon deformation. We can therefore locate the transition from a region where the power law is dominated by the FM artefacts coming from the truncation of nonaffine displacements to a region where we assume the pseudogap exponent to be valid. As 〈 x min 〉 &lt; ∆ x c and is located either in the plateau or near the crossover region, we focus on the power law region for x ≥ ∆ x c . Moreover, as noted before, the crossover toward the plateau decreases with increasing L . Therefore, we compute the cumulative distribution function (cdf) of P x ( ) for L = 300 and consider only the range of x values between the crossover to the steeper power law and the entrance to the plateau ( x ≥ ∆ x c ) as shown in the bottom panel of Figure 8. Once the pseudogap exponent θ FM is found, we determine the scaling exponent as α FM = d/ (1 + θ FM ) . In Figure 9 we see that the agreement between α FM , α S and α 〈 x min 〉 is good in the transient regime but α FM is slightly overestimated in the stationary regime. As mentioned above, the discreteness of the mechanical noise and the constraints imposed by the FM reduce significantly the range where P x ( ) ∼ x θ FM can be observed.

Given the difficulties to extract a reliable value of θ FM directly from fits to P x ( ) , we explore another approach, which considers eq. (5) and extracts 〈 x min 〉 from P x ( ) by inverting its cdf, i. e. 〈 x min 〉 = cd f -1 ( cL d ) . The resulting exponent is denoted α I . The only free parameter in this procedure is the proportionality constant c . The choice of c allows to explore regions of different values of x . We selected two values: c = 30 , which explores regions where x glyph[lessorsimilar] 0 01 . and c = 1000 , for which x glyph[lessorsimilar] 0 1 . for the smallest system size L = 53 . The former range covers the values of x in the plateau regions. In Figure 9, we show the evolution of α I upon deformation. For c = 30 , α I is much larger than α S , α 〈 x min 〉 and even α FM . For γ ≥ 0 03 . , the scaling only reflects the presence of the plateau as α I for γ = 0 03 . and γ = 0 18 . are in agreement with the values α FM = d -p computed above. The values for γ &lt; 0 03 . are also overestimated despite the absence of plateaus, possibly due to poorer sampling of small x . For c = 1000 , however, we observe a good agreement between α I and α S and α 〈 x min 〉 for the whole range of deformation. This suggests that relevant information is contained in the initial power law region beyond the plateau.

## V. CONCLUSIONS

Using the FM method, we have examined the distribution P x ( ) of residual stresses (or local thresholds to mechanical instability) in athermal 2D amorphous solids. In the quiescent (freshly quenched state), the FM method reveals a power law form with two different exponents. In the limit of small x , P x ( ) ∼ x θ 1 where θ 1 &lt; 1 agrees well with an independent estimate of the pseudogap exponent θ from extreme value statitics. The FM method then shows a second power law regime with exponent θ 2 &gt; 1 . In this regime we believe the local yield stress to be overestimated due to incomplete relaxation in the probed region, and this effect could be more pronounced for larger yield stresses.

As soon as deformation sets in, P x ( ) becomes analytic and develops a plateau as x → 0 . The FM results thus show that similar observations made with mesoscale EPMs are generic and extend to a more detailed atomistic model. In order to elucidate the origin of this plateau, we have examined P (∆ ) x , the distribution of the residual stress differences between two consecutive avalanches. This distribution has a power-law form with a system size dependent lower cutoff ∆ x c ( L ) that endows the stochastic process with a characteristic scale. We believe that the plateau appears for x &lt; ∆ x c ( L ) , which is consistent with our data in the limited range of system sizes L accessible to us in atomistic simulations.

The computation of the pseudogap exponent θ from the obtained distributions P x ( ) presents challenges. The crossover into the plateau region severely limits the region in which θ can be determined by direct fits to P x ( ) . Only larger system sizes can help. Moreover, even though the values of 〈 x min 〉 obtained from unconstrained global simulations fall into the plateau region for our largest system sizes, the observed system size scaling is not consistent with that predicted from the plateau itself. One possibility is that in the FM calculation, the contribution from the changing local yield stress ∆ σ Y to the total change ∆ x = ∆ σ Y -∆ σ 0 is somehow overestimated. The plateau thus appears sooner than in the actually sampled residual stress distributions.

One of the most intriguing observations in the present study is a changing local yield stress in regions not experiencing a plastic event. This possibility is not conceptualized in current mesocale elastoplastic [32] or meanfield [22] models of amorphous plasticity. Further improvements to the FM method that reduce the boundary artefacts and allow a better relaxation are needed to further explore this additional physics at the atomistic level.

A lingering question remains: what is the significance of the existence of this plateau in P x ( ) for the statistical properties of the yielding transition? We suggest that the answer depends on the precise form of the system size scaling of the characteristic scale of the mechanical noise increments ∆ x c ( L ) . For ideal 'Eshelby sources' and no changes of the local yield stress except upon yielding (as assumed in EPM and mean-field treatments), we have ∆ x c ( L ) ∼ L -2 from eq. (4). In this case, the plateau crossover decreases faster with L than 〈 x min 〉 ∼ L -α 〈 x min 〉 ( α 〈 x min 〉 &lt; 2) and 〈 x min 〉 must remain in the pseudogap region. The resultant scaling relations that connect the exponent θ with the exponents τ and d f that describe the statistics of macroscopic avalanches remain unaltered [33]. The present results, however, suggest at least the possibility that ∆ x c ( L ) vanishes slower than 〈 x min 〉 and thus the scaling of 〈 x min 〉 becomes dominated by the plateau regime for large enough system sizes. In this case, the scaling relationship that links the pseudogap exponent θ to the avalanche exponents τ and d f will be altered. Our present computational capabilities are insufficient to settle this question definitively .

## VI. ACKNOWLEDGEMENTS

We thank Peter Sollich for discussion and useful comments on our manuscript. J.R. thanks the Alexander von Humboldt Foundation for financial support. High performance computing resources were provided by ComputeCanada and the Quantum Matter Institute at the University of British Columbia. This research was undertaken thanks, in part, to funding from the Canada First Research Excellence Fund, Quantum Materials and Future Technologies Program. C.R. is part of the ANR LatexDry project, grant ANR-18-CE06-0001 of the French Agence Nationale de la Recherche.

## APPENDIX: DETERMINATION OF PLASTIC EVENTS

To determine plastic events, we revisit a criterion proposed by Lerner and Procaccia [18], which relies on the difference in potential energy between the affine deformation U aff and the underlying inherent structure U 0 after deformation:

$$\kappa = \frac { U _ { \text{aff} } - U _ { 0 } } { N \delta \gamma ^ { 2 } }$$

where δγ is the strain increment. While in ref. [18] a reasonable but arbitrary value of κ was selected, we demonstrate here that there exists a range of κ associated with purely elastic deformation and therefore a lower bound above which κ captures plastic events only.

The demonstration relies on the hypothesis that we work in the elastic regime. We look at the energy variation after an AQS step. The two stages in the AQS protocol lead to variation in the strain energy density that can be described as follow:

- 1. During the affine deformation the strain energy density evolves with respect to the starting configuration as:

$$\rho U _ { \text{aff} } = \rho U _ { I S } + \sigma _ { 0 } \delta \gamma + \frac { G _ { B } } { 2 } ( \delta \gamma ) ^ { 2 } \text{ \quad \ \ } ( 7 )$$

FIG. 10. Probability distribution function of κ where κ has been computed for each strain increment δγ . The vertical solid line represents the lower bound of κ to detect plastic events.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-0ad0c6a1-f4f7-4130-91ea-cabcedb9d751.jpg)

where σ 0 and U IS are respectively the residual stress and the inherent structure energy associated to the configuration before a strain increment is applied. G B is the affine shear modulus.

- 2. After the relaxation the strain energy density with respect to starting configuration is given by:

$$\rho U _ { 0 } = \rho U _ { I S } + \sigma _ { 0 } \delta \gamma + \frac { G } { 2 } ( \delta \gamma ) ^ { 2 } \quad \quad ( 8 ) \quad \text{ing} \quad \text{deal}$$

where G is the shear modulus.

Consequently the difference in strain energy density gives:

$$\rho ( U _ { \text{aff} } - U _ { 0 } ) = \frac { G _ { B } - G } { 2 } ( \delta \gamma ) ^ { 2 }$$

Finally, we can estimate an upper bound for κ in the elastic regime:

$$\kappa = \frac { G _ { B } - G } { 2 \rho } \leq \frac { G _ { B } } { 2 \rho } \, \quad \quad \ \ ( 1 0 ) \, \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

We determine that G B ≈ 26 , meaning that the lower bound for plastic activity is κ ≈ 13 . To verify our analytical description, we computed κ for each strain increment

[1] Michel Tsamados, Anne Tanguy, Chay Goldenberg, and Jean-Louis Barrat. Local elasticity map and plasticity in a model lennard-jones glass. Phys. Rev. E , 80:026112, Aug 2009.

[2] J. P Wittmer, A Tanguy, J.-L Barrat, and L Lewis. Vibrations of amorphous, nanometric structures: When does continuum theory apply? Europhysics Letters (EPL) ,

δγ = 5 · 10 -5 upon deformation (up to 30% ). Its probability distribution function is shown in Figure 10, where we observe that the distribution is peaked around κ ≈ 8 and a long tail persists for larger value of κ . The peak is associated with elastic contributions. Indeed we measured that G ≈ 12 for this system so κ is more likely to be (26 -12) / 2 ≈ 7 .

We decided to chose κ ≈ 30 as criterion to ensure to have only plastic contributions. In Figure 11, we observe that κ is matching with stress release in the stationary regime but more interestingly, as suggested in the inset,

FIG. 11. Main panel: Typical stress-strain curve obtained during the deformation of a global system where L = 100 . The peaks in the κ observable indicate stress drops. Inset: Zoom in the initial stage of the deformation.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.01081v4-figure-0a47c92b-14f8-44be-a720-c7632653fd68.jpg)

this criterion allows one to determine plastic events for which no stress drop is measured. This is due to the higher sensitivity of the potential energy to structural rearrangement. The detection in terms of stress drops is constrained by the choice of the strain increment as a stress release would have been certainly noticed with a smaller δγ but this would also imply longer simulation times. Therefore, κ is an efficient criterion to monitor plasticity while working with a finite value of δγ .

57(3):423-429, feb 2002.

[3] Kenji Yoshimoto, Tushar S. Jain, Kevin Van Workum, Paul F. Nealey, and Juan J. de Pablo. Mechanical heterogeneities in model polymer glasses at small length scales. Phys. Rev. Lett. , 93:175501, Oct 2004.

[4] P. Sollich. CECAM meeting , 2011.

- [5] Hideyuki Mizuno, Stefano Mossa, and Jean-Louis Barrat.

- Measuring spatial distribution of the local elastic modulus in glasses. Phys. Rev. E , 87:042306, Apr 2013.
- [6] Francesco Puosi, Julien Olivier, and Kirsten Martens. Probing relevant ingredients in mean-field approaches for the athermal rheology of yield stress materials. Soft Matter , 11:7639-7647, 2015.
- [7] Baoshuang Shang, Pengfei Guan, and Jean-Louis Barrat. Role of thermal expansion heterogeneity in the cryogenic rejuvenation of metallic glasses. Journal of Physics: Materials , 1(1):015001, sep 2018.
- [8] Sylvain Patinet, Damien Vandembroucq, and Michael L. Falk. Connecting local yield stresses with plastic activity in amorphous solids. Phys. Rev. Lett. , 117:045501, Jul 2016.
- [9] Armand Barbot, Matthias Lerbinger, Anier HernandezGarcia, Reinaldo Garc´ ıa-Garc´ ıa, Michael L. Falk, Damien Vandembroucq, and Sylvain Patinet. Local yield stress statistics in model amorphous solids. Phys. Rev. E , 97:033001, Mar 2018.
- [10] Armand Barbot, Matthias Lerbinger, Ana¨ el Lemaˆ ıtre, Damien Vandembroucq, and Sylvain Patinet. Rejuvenation and shear banding in model amorphous solids. Phys. Rev. E , 101:033001, Mar 2020.
- [11] M. L. Manning and A. J. Liu. Vibrational modes identify soft spots in a sheared disordered packing. Phys. Rev. Lett. , 107:108302, Aug 2011.
- [12] Majid Mosayebi, Patrick Ilg, Asaph Widmer-Cooper, and Emanuela Del Gado. Soft modes and nonaffine rearrangements in the inherent structures of supercooled liquids. Phys. Rev. Lett. , 112:105503, Mar 2014.
- [13] Jun Ding, Sylvain Patinet, Michael L. Falk, Yongqiang Cheng, and Evan Ma. Soft spots and their structural signature in a metallic glass. Proceedings of the National Academy of Sciences , 111(39):14052-14056, 2014.
- [14] C. Patrick Royall and Stephen R. Williams. The role of local structure in dynamical arrest. Physics Reports , 560:1 - 75, 2015. The role of local structure in dynamical arrest.
- [15] F. Lancon and L. Billard. ¸ Two-dimensional system with a quasi-crystalline ground state. J. Phys. France , 49(2):249256, 1988.
- [16] M. L. Falk and J. S. Langer. Dynamics of viscoplastic deformation in amorphous solids. Phys. Rev. E , 57:7192-7205, Jun 1998.
- [17] S. J. Plimpton. Fast parallel algorithms for short-range molecular dynamics. J Comp Phys , 117, 1995.
- [18] Edan Lerner and Itamar Procaccia. Locality and nonlocality in elastoplastic responses of amorphous solids. Phys. Rev. E , 79:066109, Jun 2009.
- [19] Smarajit Karmakar, Edan Lerner, and Itamar Procaccia. Statistical physics of the yielding transition in amorphous solids. Phys. Rev. E , 82:055103, Nov 2010.
- [20] H. G. E. Hentschel, Prabhat K. Jaiswal, Itamar Procaccia,
- and Srikanth Sastry. Stochastic approach to plasticity and yield in amorphous solids. Phys. Rev. E , 92:062302, Dec 2015.
- [21] Jie Lin, Alaa Saade, Edan Lerner, Alberto Rosso, and Matthieu Wyart. On the density of shear transformations in amorphous solids. EPL (Europhysics Letters) , 105(2):26003, jan 2014.
- [22] Jie Lin and Matthieu Wyart. Mean-field description of plastic flow in amorphous solids. Phys. Rev. X , 6:011005, Jan 2016.
- [23] Botond Tyukodi, Damien Vandembroucq, and Craig E. Maloney. Avalanches, thresholds, and diffusion in mesoscale amorphous plasticity . Phys. Rev. E , 100:043003, Oct 2019.
- [24] E. E. Ferrero and E. A. Jagla. Criticality in elastoplastic models of amorphous solids with stress-dependent yielding rates. Soft Matter , 15:9041-9055, 2019.
- [25] Andrea Zoia, Alberto Rosso, and Satya N. Majumdar. Asymptotic behavior of self-affine processes in semiinfinite domains. Phys. Rev. Lett. , 102:120602, Mar 2009.
- [26] T. Albaret, A. Tanguy, F. Boioli, and D. Rodney. Mapping between atomistic simulations and eshelby inclusions in the shear deformation of an amorphous silicon model. Phys. Rev. E , 93:053002, May 2016.
- [27] The statistics for the noise obtained for the smallest system L = 53 is not good enough to extract any reliable information. Indeed as reducing the size implies larger γ Y between events, we did not record enough stress drops on a strain interval of 2% .
- [28] I. Fern´ndez Aguirre and E. A. Jagla. a Critical exponents of the yielding transition of amorphous solids. Phys. Rev. E , 98:013002, Jul 2018.
- [29] K. Michael Salerno and Mark O. Robbins. Effect of inertia on sheared disordered solids: Critical scaling of avalanches in two and three dimensions. Phys. Rev. E , 88:062206, Dec 2013.
- [30] Chen Liu, Ezequiel E. Ferrero, Francesco Puosi, JeanLouis Barrat, and Kirsten Martens. Driving rate dependence of avalanche statistics and shapes at the yielding transition. Phys. Rev. Lett. , 116:065501, Feb 2016.
- [31] Misaki Ozawa, Ludovic Berthier, Giulio Biroli, Alberto Rosso, and Gilles Tarjus. Random critical point separates brittle and ductile yielding transitions in amorphous materials. Proceedings of the National Academy of Sciences , 115(26):6656-6661, 2018.
- [32] Alexandre Nicolas, Ezequiel E Ferrero, Kirsten Martens, and Jean-Louis Barrat. Deformation and flow of amorphous solids: Insights from elastoplastic models. Reviews of Modern Physics , 90(4):045006, 2018.
- [33] Jie Lin, Edan Lerner, Alberto Rosso, and Matthieu Wyart. Scaling description of the yielding transition in soft amorphous solids at zero temperature. Proceedings of the National Academy of Sciences , 111(40):14382-14387, 2014.
Céline Ruscher Jörg Rottler
Residual stress distributions in athermally deformed amorphous solids from atomistic simulations
The distribution of local residual stresses (threshold to instability) that
controls the statistical properties of plastic flow in athermal amorphous
solids is examined with an atomistic simulation technique. For quiescent
configurations, the distribution has a pseudogap (power-law) form with an
exponent that agrees well with global yielding statistics. As soon as
deformation sets in, the pseudogap region gives way to a system size dependent
plateau at small residual stresses that can be understood from the statistics
of local residual stress {\em differences} between plastic events. Results
further suggest that the local yield stress in amorphous solids changes even if
the given region does not participate in plastic activity.