Insights from the quantitative calibration of an elasto-plastic model from a Lennard-Jones atomic glass

We compare the macroscopic and the local plastic behavior of a model amorphous solid based on two radically different numerical descriptions. On the one hand, we simulate glass samples by atomistic simulations. On the other, we implement a mesoscale elasto-plastic model based on a solid-mechanics description. The latter is extended to consider the anisotropy of the yield surface via statistically distributed local and discrete weak planes on which shear transformations can be activated. To make the comparison as quantitative as possible, we consider the simple case of a quasistatically driven two-dimensional system in the stationary flow state and compare mechanical observables measured on both models over the same length scales. We show that the macroscale response, including its fluctuations, can be quantitatively recovered for a range of elasto-plastic mesoscale parameters. Using a newly developed method that makes it possible to probe the local yield stresses in atomistic simulations, we calibrate the local mechanical response of the elasto-plastic model at different coarse-graining scales. In this case, the calibration shows a qualitative agreement only for an optimized subset of mesoscale parameters and for sufficiently coarse probing length scales. This calibration allows us to establish a length scale for the mesoscopic elements that corresponds to an upper bound of the shear transformation size, a key physical parameter in elasto-plastic models. We find that certain properties naturally emerge from the elasto-plastic model. In particular, we show that the elasto-plastic model reproduces the Bauschinger effect, namely the plasticity-induced anisotropy in the stress-strain response. We discuss the successes and failures of our approach, the impact of different model ingredients and propose future research directions for quantitative multi-scale models of amorphous plasticity.


Introduction
Multi-scale models of plastic deformation aim to gain a qualitative and quantitative understanding of the relation between microstructural properties and the dynamics of plastic activity.Ultimately, such models can help design new microstructures tailored to meet the demands of specific engineering applications and industries [1,2].An essential ingredient for a successful multiscale approach is establishing links between atomistic and macroscale continuum descriptions in a physically grounded manner.In this perspective, an intermediate mesoscale description is a desirable conceptual step, as illustrated, for instance, in crack propagation [3] and crystal plasticity [4].Mesoscale approaches rely on a statistical description of the relevant phenomena, enhancing our understanding of the collective processes at play.By doing so, it overcomes the limitations of atomistic and continuum approaches.Namely, atomistic models are limited in terms of system sizes and time scales.On the other hand, continuum approaches have difficulties in capturing complex heterogeneous flow patterns or spatio-temporal correlations in plastic activity [5].Specifically, mesoscale models circumvent these limitations by replacing the quasi-infinite microscopic degrees of freedom at the atomistic scale with more manageable, discrete, and coarser ones, based on a continuum description [6].
In amorphous materials, plastic deformation occurs at the lowest scales via atomic rearrangements leading to localized shear transformations (STs) [7][8][9] and, in some cases, to a permanent dilation or contraction associated with local changes of free volume [10].The properties of STs have been studied in detail and present large distributions of activation energy barriers [11], directions [12], sizes [13] and shapes [8] and lead to a redistribution of elastic stresses within the system [14].Thus, in these materials, plastic deformation is made up of "quanta" corresponding to local topological changes, which, unlike dislocations [15], are localized in space and time [16].
This view of plasticity has been successfully implemented in discrete elasto-plastic models based on a solid-mechanics description [17].The elasto-plastic approach considers a discretecontinuum medium endowed with stochastic local evolution rules for plastic activity and structural properties [12,[18][19][20][21][22].One advantage of elasto-plastic models is that they rely on physically meaningful quantities that are easy to interpret.At the same time, they can qualitatively reproduce the phenomenology associated with amorphous plasticity in many conditions such as, e.g., under athermal quasistatic shear [20][21][22][23][24], at a finite strain-rate [25,26], under creep loading [27] or near mechanical failure [28,29].This simplicity and versatility can help establish links between statistical physics and engineering formulations of plasticity.
Nonetheless, the lack of well-defined plastic deformation units in amorphous solids has hindered the development of multi-scale strategies.Elasto-plastic models have thus remained widely phenomenological or, at best, based on educated guesses.In the shear transformation zone theory framework [30], Hinkle and Falk have for instance employed the potential energies of the atomistic model coarse-grained at different length scales as a surrogate of the effective temperature of the continuum-level field description [31].However, elasto-plastic mesoscale parameters need to be material specific and realistic to make quantitatively accurate predictions.For this reason, there have been many attempts in recent years to measure mesoscale parameters from atomistic simulations [11][12][13][32][33][34][35].
A fundamental ingredient of elasto-plastic models corresponds to a local measure of shear susceptibility, able to predict in principle the plastic activity [36][37][38][39][40][41][42][43][44][45].Among the proposed alternatives, the local yield stress exhibits an excellent correlation with plastic event locations in comparison with other indicators [46].Direct measurements of local yield stress in atomistic samples C. R. Physique -2021, 22, n S3, 135-162 Figure 1.Multi-scale approach: the macroscopic (top) and local (bottom) quasi-static mechanical responses are investigated in both elasto-plastic and atomistic bi-dimensional models.To measure the local yield stress, patches of equal area are isolated from the rest of the system.have become recently possible using a method [41] based on straining subsets of atoms, as illustrated in Figure 1.The method can characterize the local resistance to irreversible deformation, which can be interpreted as a mesoscale generalization of the macroscale yield stress signaling the onset of plastic deformation.Therefore, the local yield stress is well suited for elasto-plastic models, rooted in the framework of solid mechanics.Moreover, the method makes it straightforward to gather statistics of diverse local properties that traditionally have been challenging to measure, such as rearrangement amplitudes and directions [47], structural evolution [48] or local yield surface anisotropy [49].However, until recently [50], these promising properties have not yet been explored to tentatively fill the gaps necessary for quantitatively accurate elasto-plastic models.
Filling such gap is the aim of the present work in which we compare, as quantitatively as possible, atomistic and elasto-plastic models both at the macro and local scales by leveraging the advantages of local yield stress measurements.We prepare atomistic glass samples driven in the athermal quasistatic limit until the stationary flow state.The local mechanical properties of the as-quenched glasses and the stationary states are measured with the local yield stress method described in [41,47] at different length scales and for different shear orientations.On the other hand, we set up a mesoscale elasto-plastic model with statistically distributed structural properties and anisotropic local yield functions and mimic the atomistic glass preparation and loading protocol.Contrary to [50], we focus on the athermal quasistatic regime and consider a generalization of the local yield stress method of [41,47] for different length scales.Motivated by local atomistic observations showing anisotropic yield surfaces and non-normal plastic flow [47], we make a step forward by modeling the possibility of plastic rearrangements according to a set of various local discrete weak planes for each element of the mesoscale elasto-plastic model.We investigate the consequences of this feature in terms of model calibration and emergent plastic anisotropy upon plastic deformation.The comparison of behaviors at local scales is carried out by measuring the response of different subsets of the elasto-plastic matrix by mimicking the atomistic procedure (see Figure 1).As a result, we do not establish a hierarchical dependence between models, in the sense that atomistic data does not feed the elasto-plastic model.Instead, we consider atomistic and elasto-plastic descriptions on an equal footing and calibrate the elastoplastic model by requiring the same mechanical properties on both models when measured on the same length scale.In this way, it becomes possible to not only calibrate the model but to assess its predictions at different scales simultaneously, which allows us to carefully analyze agreements and discrepancies with the reference atomistic glass.
The article is organized as follows: in the second and third sections, we introduce the atomistic and elasto-plastic models and methods employed to characterize both the global and local mechanical response of glasses.Section 4 deals with the calibration of the elasto-plastic model and the comparison of atomistic and mesoscale approaches at both macroscopic and local scales.The various emergent phenomenologies obtained from the elasto-plastic model are presented in Section 5.In Section 6, we discuss the successes and failures of our approach, different key ingredients, and suggest future research directions for quantitative multi-scale models of amorphous plasticity.

Preparation and loading protocols
We use the two-dimensional binary atomic system studied extensively in [41,47].This system has been widely used as a model glass for its good glass-formability and to study the plasticity of glasses from the atomic scale [16].It makes it possible to reproduce the main features of the phenomenology observed in the mechanical response of amorphous solids while remaining as simple as possible.The composition is chosen so that the ratio between the number of large (L) and small (S) particles is equal to N L /N = (1 + 5)/4.The atoms, all with a mass m = 1, interact via a Lennard-Jones (LJ) interatomic potential of parameters σ AB and AB , where AB corresponds to the interacting species.The inter-species parameter of the potential defines the units of length σ SL , energy SL , time t 0 = σ SL m/ SL and stress Σ = SL /σ 2 SL , used in the rest of the article.For distances between R in = 2 and the cutoff radius R out = 2.5, the LJ potential is replaced by a polynomial function to be twice differentiable.
The systems are simulated at a constant density ρ = 1.02 using periodic boundary conditions.Except for the simulations exploring finite-size effects, where larger systems are employed (see Figure 5b), the systems contain 10 4 atoms with a linar size equal to L = 98.8045.The glasses are obtained by instantaneously quenching at constant volume a supercooled liquid at thermodynamic equilibrium from a temperature T = 1.132T sim g , where T sim g = 0.31 SL /k B .The quench is followed by a static relaxation to cancel forces via a conjugate gradient method down to machine precision.This preparation protocol, named equilibrated supercooled liquid (ESL) in [47], has the advantage of producing glass samples which exhibit no stress overshoot when strain is applied and thus without sample-scale shear banding.This is specially convenient for studying the stationary flow regime.
We consider 100 independent glass samples deformed by the athermal quasi-static shear (AQS) deformation method [51] in simple shear along the xy direction until a deformation γ ext = 5.This iterative method consists of deforming the system by small increments of affine strain ∆γ ext = 10 −4 , followed by a static relaxation.The trajectories obtained consist of a series of  x y induced by discrete steps of ∆γ ext = 10 −4 .The results shown correspond to the optimal parameters with l = 6.6.
states at mechanical equilibrium corresponding to the limit of zero temperature and strain rate.In the large strain regime, the typical mechanical response reported in Figure 2a is characterized by the fluctuation of the measured stress Σ ext x y around a plateau value.As shown schematically in Figure 1, it consists of a succession of quasi-linear elastic branches, where the external stress increases by increments ∆Σ ext x y , intercepted by stress drops δΣ ext x y that corresponds to the occurrence of irreversible deformation events dissipating elastic energy.The steady-state, signaled by the convergence of stresses and potential energy, is reached for a strain γ ext > 0.5.

Atomistic implementation of the local shear test
To measure the local mechanical response, we use the recently developed local yield stress method [33,41,47], as illustrated in Figure 1.This method consists of deforming in pure shear along direction α a zone of radius R, known as a patch, by applying a purely affine deformation to the surrounding medium.The central area is deformed using the same AQS method as for deforming the overall system.Therefore, plastic events are necessarily triggered within the relaxed zone, which makes it possible to measure the local coarse-grained yield stress τc in several directions independently of other plastic rearrangements.Since glasses are frustrated systems, they are locally decorated with a spatially heterogeneous local stress τi .The residual plastic strength, or distance to threshold ∆ τc , corresponds to the amount of stress added at the coarse-graining scale R to trigger the instability locally, thus ∆ τc = τc − τi .This method also makes it possible to estimate the relaxation amplitude δ τc by measuring the stress drop associated with the plastic instability.
The sampling of τi , τc , ∆ τc , and δ τc is performed in the as-quenched state and the stationary state on a regular square grid with lattice parameter 2.5.These local quantities are computed for different patch sizes R and orientations α.We note that the local yield stresses and distances to threshold (stress drops) are expected to be overestimated (underestimated) for too small patch radii R as a consequence of the hard boundary conditions imposed by the frozen matrix of the local yield stress method.A more accurate measure is thus anticipated as R grows.This aspect is discussed further when we compare the results obtained with the elasto-plastic model.

Elasto-plastic model
Elasto-plastic models define the state of a material subvolume in terms of continuum mechanics fields.Here, we discretize the material domain into a 2D regular lattice of mesoscale elasto-plastic quadrilateral elements.Each element corresponds to a material sub-volume where localized atomistic rearrangements can occur.Such rearrangements are represented, at the mesoscale, as localized slip events.This coarse-grained description aims at capturing the essential effects of the rearrangements while neglecting unnecessary details.As illustrated in Figure 1, the elements have a length l that is above the typical length scale of the atomistic rearrangements.Below this scale, the elasto-plastic model cannot resolve microscopic details.Consequently, we consider strain fields which are element-wise constant, and associate each mesoscale element with a single elastic and plastic strain tensor, ε el and ε pl respectively.
The elastic strain field ε el is the consequence of externally imposed boundary conditions and the presence of the plastic strain field ε pl .We solve the stress equilibrium equation assuming linear elasticity applying the Finite Element Method (FEM) [21,26,27,29,52,53].To this end, we consider a 2D regular grid of quadrilateral finite elements (FEs) with linear shape functions.Each FE exactly matches the shape of a single mesoscale element.By applying the FEM, we obtain the displacement field from where we compute the strain field ε as the symmetric gradient.The strain associated with a mesoscale element is defined as the strain averaged within the corresponding FE.The elastic strain is obtained as ε el = ε − ε pl and the stress is computed from linear elasticity as Σ = C : ε el + Σ 0 , where C denotes the rank-4 stiffness tensor and Σ 0 is the pre-stress present in the system prior to the initiation of the driving protocol.To establish as much contact as possible with the atomistic model, we compute the solution to the stress equilibrium equation under biperiodic boundary conditions with an externally applied shear strain γ ext along the xy direction.The external shear stress Σ ext x y is computed as the average shear stress over all the elasto-plastic elements composing the system.
In the next subsections, we detail the structural properties of the mesoscale elements and their stochastic evolution rules.Plastic activity occurs as a sequence of slip events coupled by stress redistribution competing against structural disorder.The specific choices we make regarding these rules are motivated by measurements on the atomistic systems, against which the elastoplastic model results will be bench-marked later.

Local slip systems
Plastic activity in glasses proceeds by atomistic rearrangements, which lead, in general, to a permanent localized deformation with a shear component and a hydro-static one, the latter associated with the creation or annihilation of free volume.Nonetheless, following the strong correlation between local shear thresholds and plastic activity [46,47] as well as the weak correlation between local free volume and softness [48] in the atomistic system considered here, we simplify the description and consider only the shear component.
Even when the tensorial formulation of plasticity is accounted for, the vast majority of the elasto-plastic models assume an isotropic yield criterion to simulate the activation of plastic deformation [6,21,29].However, atomistic measurements of the local response to shear of glasses reveal anisotropic yield surfaces and the presence of weak planes [47] with spatially-fluctuating properties.To recreate such a complex response to locally applied shear, we consider that each mesoscale element contains several slip systems.We define each slip system by a plane of normal unit vector n and a direction s, contained within the plane.A slip of amplitude ∆γ pl is thus represented by the plastic strain increment with Since we consider a 2D scenario, each slip system contains two possible slip senses +s and −s.However, we associate each slip system with a single sense.The reason for this choice is that, in the presence of strong structural disorder, planes that qualify as weak under the action of local shear stress may not qualify as such when the shear direction is reversed [47].We can write the tensor M in terms of the angle θ between the slip plane and the horizontal axis as with θ ∈ (−π/2, π/2] due to symmetry.The resolved shear stress τ on a slip plane is given by Each slip system has a critical resolved shear stress τ c > 0, in the following referred to as slip threshold.The distance to threshold is defined as ∆τ c = τ c − τ.Whenever a slip system fulfills ∆τ c < 0, it is denoted as active, and a slip event is performed as described later.

Structural properties
To set the density of slip systems, we consider that a zone of l ≈ 2 would contain an average of 4 atoms and a strict maximum number of slip systems equal to N = 4. Consequently, a reasonable upper bound for the slip system density is 1.In our simulations, we set the density of slip systems to ρ s = 0.8, below the estimated upper bound (we assess the robustness of our results upon variations of this parameter in Appendix A.1).We relate the number N of slip systems to the discretization length scale l as N = ρ s l 2 .Both the slip angles and thresholds are statistically distributed to represent structural heterogeneity.The orientation of the local slip systems defined within an element, although statistically distributed, must fulfill a certain constraint.Namely, the element must have a finite critical resolved shear stress τ c (α), i.e. defined for any shear orientation α.The simplest way to achieve this is to consider at least four slip systems with θ, θ + π/2, θ + π/4 and θ + 3π/4 which ensures the fulfillment of the constraint.Due to the a priori lack of privileged orientation, we consider θ uniformly distributed in the interval (−π/2, π/2].We populate each mesoscale element with N multiple of 4 to apply the procedure defined above.
Slip thresholds are independently drawn, with the only requirement of being positive-definite.For simplicity, we use the Weibull distribution where the parameter λ defines the scale and the exponent k the shape of the distribution.

Performing a slip event
When a slip event occurs within an element, we update the plastic strain field ε pl by adding a tensorial increment ∆ε pl defined by (1), which is homogeneous through the element and zero everywhere else.We consider the amplitude ∆γ pl of the strain increment statistically distributed to represent the effects of the heterogeneous microstructure on an initiated slip event.We impose the constraint of avoiding negative dissipation [54].Within our model (Appendix A.7), this constraint is fulfilled if the slip amplitude verifies ∆γ pl < γ max with where S is the Eshelby tensor of the elasto-plastic elements [55].We consider a bounded distribution in order to explicitly verify the constraint.Specifically, we choose a bounded power-law distribution of the form with 0 < ∆γ pl < γ max and χ > 0. This distribution has the advantage of having a single free parameter, which allows to vary its shape while limiting the complexity of the model calibration.We note, however, that its specific functional form does not significantly affect the results, as will be shown in Section 6.After a slip event has occurred in a certain element, we renew the orientations and thresholds of the N slip systems within the element from their respective probability distributions.The value of N remains fixed.This process aims to represent stochastic changes in the local microstructural properties induced by plastic activity.

Loading protocol
We implement an athermal quasistatic protocol to represent as faithfully as possible the AQS driving protocol of the atomistic model (Section 2.1).In the absence of temperature effects, we consider that whenever one or more slip systems are active, slip events are simultaneously performed in all those systems (Section 3.3).Due to stress redistribution, new slip systems might become active.This process is repeated in a series of relaxation steps until no slip system is active.During this process, the external strain γ ext is kept fixed.In the case where several slip systems within the same mesoscale element simultaneously become active, only the slip system with the lowest distance to threshold ∆τ c is considered active.
When there is no active slip system, we perform discrete external shear strain increments of ∆γ ext = 10 −4 .Slip events taking place between two external strain increments are known as an avalanche.We apply external strain increments until reaching a maximum external strain target of γ ext = 2.

Elasto-plastic implementation of the local shear test
To reproduce as faithfully as possible the results measured at local scales in atomistic simulations, we mimic the method described in Section 2.2 with the elasto-plastic model.Elasto-plastic patches are defined respecting the geometrical restrictions imposed by the lattice discretization.Thus, we define quadrilateral patches formed by n × n contiguous elasto-plastic elements selected from the system (Figure 1).We associate to each elasto-plastic patch an effective radius R = (l / π)n by considering the radius of a circular patch of equal area.This allows us to compare the results with atomistic patches of a specific radius.
To measure the mechanical response of a patch, we consider it in isolation from the rest of the system (Figure 1, left).The stress Σ of a patch is defined as the average stress of the n × n elements forming the patch.Before the initiation of a shear test, a patch has an initial stress Σi .
We apply to its free boundary a pure shear strain with orientation α and amplitude γ.When one slip system becomes active, the patch has reached a critical stress Σc .At this point, an avalanche of slip events occurs at a constant external strain.Afterward, the patch has a final stress Σf .This process is illustrated in Figure 1 (right).Except for the non-periodic boundary conditions, elastoplastic patches are driven, and slip events occur, according to the same rules described for the full system (Section 3).We note that, although the patch is led to mechanical instability by raising the stress till Σc , this value corresponds to the total amount of stress at the moment of instability.Thus, it depends only on the patch structural properties.We can use it to establish definitions which closely mimic the atomistic ones.Thus, we define the coarse-grained yield stress for the orientation α as τc = M(α) : Σc that is, the resolved shear stress on the shearing plane α at the moment of the patch mechanical instability.In order to characterize local stability in a manner analogous to the atomistic method, we define the coarse-grained distance to threshold ∆ τc as ∆ τc = τc − τi where τi = M(α) : Σi .Moreover, we define the coarse-grained stress drop δ τc from the yield stress to the final stress state as δ τc = τc − τf .
We remark the difference between the coarse-grained patch-scale yield stress τc and the slip thresholds τ c defined in Section 3.1.The former characterizes the measured resistance to plastic deformation at the coarse-graining scale defined by the patch radius R, while the latter refers to the resistance to slip of individual weak planes within the patch, with values drawn from (5).

Calibration of the elasto-plastic model
In this section, we calibrate the mesoscale resolution l and the elasto-plastic parameters λ, k, and χ associated with the threshold scale, threshold disorder, and slip event amplitude, respectively.The elasto-plastic model considers isotropic and homogeneous elastic properties.To this end, we use the system-scale effective shear modulus G = 13.2 and bulk modulus B = 59 (in LJ units) measured in the atomistic model.
First, we focus our attention on the macroscale behavior in the stationary state, which does not depend on the initial conditions.Macroscale behavior provides us with an initial general understanding of the system dynamics.Moreover, it introduces constraints in the mesoscale parameter space, establishing a strategy to calibrate the local spatially-fluctuating scale-dependent mesoscale properties.We perform the calibration at different element length scales l and establish an optimum range of values.
In addition, we also take advantage of this calibration procedure to generate reasonable quench states by assuming that the thermally-activated structural changes occurring in glassforming liquids are statistically close to those induced by plastic activity [56].The optimal set of parameters are thus resorted to generating the quench states from a Kinetic Monte Carlo approach where only the temperature needs to be calibrated, as detailed in Appendix A.3.

Macroscale behavior
As the externally applied strain γ ext increases, the system exhibits transient behavior related to the initial system configuration.Eventually, the system reaches a stationary regime in which its statistical properties are independent of the external strain (Figure 2a).This regime can be identified by a plateau in the external stress Σ ext x y .For our analysis, we focus on the region γ ext > 0.5.To quantitatively calibrate the macroscale behavior, we compute the overlap between the elasto-plastic and the atomistic probability densities of external stress values Σ ext x y , increments ∆Σ ext x y induced by the loading mechanism and drops δΣ ext x y caused by avalanches of plastic activity.
Based on the overlap, we define an error function L[P] for each probability density as where P EP and P MD respectively refer to the elasto-plastic and atomistic versions of the probability densities.When the densities perfectly overlap, L[P] = 0 and L[P] → 1 as they differ.We define the overall macroscale error as the average L macro = 1/3 P L[P] of the errors computed for each distribution.
We explore the 3-dimensional space of values for k, λ and χ.To do so efficiently, we consider discrete values of k from k = 1 to k = 3 in intervals of 0.05.We ensure that this is a range of plausible values for reproducing the atomistic measurements.For each k, we explore the 2dimensional space of values for λ and χ > 0, repeating each combination 48 times with a different realization of the initial conditions.Among the studied combinations of parameters, for each value of k we consider only the pair λ and χ with the lowest L macro .We repeat this procedure for different mesoscale element lengths l .For each l , the size of the elasto-plastic lattice and the number N of slip systems per element are set as detailed in Appendix A.2.We find that an excellent simultaneous agreement in all the compared macroscale magnitudes is possible for the different element lengths l considered (see Figure 2 for the specific case of l = 6.6).
Moreover, a fit of similar quality is possible for a wide range of k (Figure 3a).Consequently, macroscale behavior can be regarded as a constraint that must be fulfilled by any acceptable set of elasto-plastic mesoscale parameters but does not provide us with enough information for the calibration of the model, for which knowledge of local mesoscale properties becomes mandatory.

The error L[P(δΣ ext
x y )] in the fit to the distribution of external stress drops increases as the model resolution l is reduced (Figure 3a).External stress drops are the consequence of avalanches of slip events.The bigger the value of l , the bigger the portion that each mesoscale element covers, and therefore in this case avalanches become composed of fewer slip events.Eventually, for big enough l avalanches become single events, with amplitudes that are fit as part of the model optimization.On the other hand, for smaller values of l , avalanches become collective events whose sizes are not explicitly part of the model optimization and emerge from the selforganization of plastic activity.

Mesoscale local properties
For each simulated sample, an ensemble of patches representative of the system is defined, and for each patch, shear tests are performed with discrete orientations.Special attention is paid to the forward direction α = 0°aligned with the external load, and the backward direction α = 90°, which will be used for the quantitative calibration of the mesoscale properties.
We define a set of mesoscale magnitudes that allow us to quantitatively compare different aspects of the local mechanical response of both models when measured at the same scale R. Specifically, we consider the average and standard deviations of the local resolved shear stress τi (R, α), yield stress τc (R, α), distance to threshold ∆ τc (R, α) and stress drop δ τc (R, α).By generically denoting the elasto-plastic measurements as F EP (R, α) and the atomistic ones as F MD (R, α), we define for these 8 magnitudes the errors L[F ] as Although we will explore the impact of R in the statistics of local properties, we compute the error functions using the largest available value of R = 30, since presumably it is less affected by the  rigid boundary bias introduced by the patch measuring method.We define the mesoscale error in the forward direction α = 0°as the average We define the error L − meso in the backward direction α = 90°in a similar fashion.The overall stationary state error, L stat , is defined as the average between the macroscale and the mesoscale errors, L stat = 1/2(L macro +1/2(L + meso + L − meso )).We optimize L stat by exploring the neighborhood of the values that optimized L macro (Section 4.1).This allows us to optimize the fit to the mesoscale properties while fulfilling the constraint imposed by the macroscale behavior.
As discussed previously, macroscale properties are not enough to establish the optimum set of mesoscale parameters.However, thanks to the addition of the mesocale properties, L stat exhibits a minimum Lstat at a specific value of k (Figure 3b).We fit a parabola a + bx + c x 2 to L stat to estimate the optimum value k.Since to each value of k we associate a pair of λ and χ, we can by interpolation obtain the optimum values λ and χ.As shown in Figure 3c, Lstat depends on the element length l .We find an optimum value in the range of l = 4.4-6.6.The error minimum is however rather shallow and the error starts to increase more appreciably only for larger element length.
In order to further set the value of the element length l we take into account several considerations: first, the model loses spatial resolution for the coarser element lengths (but remarkably can reproduce macroscale properties with excellent quality).Second, we also expect that the representation of STs requires a larger plastic strain amplitude for a lower element length l .Indeed, lowering l leads to best-fit parameters with a decreasing value of χ, which implies larger local stress drops.We found that, below l ≈ 5, χ < 1.That is, the limit at which changing the sign of the resolved shear stress in the active slip plane becomes the dominant behavior.If this does not necessarily violate the non-negative dissipation criterion (Section 3.3), it seems unlikely regarding the expected physics of STs.Third, the range l = 4.4-6.6 is comparable to the lower limit established in [36] for the validity of linear elasticity.Moreover, when employed as the patch size in the local yield stress method, it has been shown that this length scale optimizes the correlation between plastic activity and ∆ τc [47].Since we do not observe a significant difference in the error computed within that range, in the following, we consider the upper bound l = 6.6.This scale is presumably above the typical size of atomistic rearrangements and larger than spatial correlations in the structural renewal process, both of which are fundamental assumptions of the elasto-plastic model.
We can have more insight into the origin of the error basin of L stat by looking at its components L macro , L + meso and L − meso separately.As shown in Figure 3d with l = 6.6, the error L + meso of the forward-oriented mesoscale properties exhibits qualitatively similar behavior as the macroscale error, with a wide range of parameters leading to quantitatively similar values without a clear minimum.However, the error L − meso of the backward-oriented mesoscale properties exhibits a minimum at a specific combination of parameters.Consequently, among the parameters that reproduce the macroscale and the mesoscale properties aligned with the external load, only a specific combination k, λ and χ can simultaneously reproduce the behavior of the patches when unloaded from the stationary state.
The estimated optimal parameter combination for the stationary state at the mesoscale l = 6.6 is k = 2.18, λ = 2.05 and χ = 2.25.The overall relative fit error L stat is 16%.With the optimal parameters, we perform 1024 simulation runs.We obtain a qualitative agreement in the quench and stationary states for all the investigated properties that progressively improves both on average and standard deviation as R increases.At R = 30, we find a good general quantitative agreement.In the quench state, the results are independent of α due to the lack of privileged orientation.On the other hand, we observe in the stationary state a dependence on shear orientation related to the mere static equilibrium ( τi ) but also to the induced anisotropy ( τc and ∆ τc ) known as the Bauschinger effect [49].
The local shear stress 〈 τi 〉 and the yield stress exhibit a remarkably good match for every orientation α.However, the distances to threshold 〈∆ τc 〉 and the stress drops 〈δ τc 〉 exhibit qualitative differences, suggesting that the elasto-plastic model captures the emergent anisotropy with a different degree of accuracy depending on the orientation α.Specifically, we observe a systematically higher accuracy in the forward direction (α = 0°) than in the backward direction (α = 90°).Moreover, the distance to threshold in the backward direction exhibits an incompatible scaling with R. Another qualitative discrepancy between models are the small local minima (Figure 4l, Figure 8d and f ), presumably an artifact induced by the FEM quadrilateral structured mesh.
We note that Figure 4(e,f ) reports negative yield stress 〈 τc 〉 < 0. Although slip systems have positive-definite slip thresholds τ c > 0 the definition of coarse-grained yield stress (Section 3.5) allows the measurement of negative values.A detailed discussion can be found in Appendix A.  Comparison of the local properties of the atomistic (blue) and elasto-plastic (red) models in the quench and stationary regimes: von Mises stress Σvm and resolved shear stress τi on the shear test plane (first row), yield stress τc (second row), distance to threshold ∆ τc (third row) and stress drop δ τc (forth row).The columns correspond to the probability distribution functions for R = 30 (left), the averages as a function of R (center) and the averages as a function of the shear orientation α for R = 30 (right).The results shown correspond to the optimal parameters with l = 6.6.

Emergent properties
The optimal elasto-plastic parameters are estimated by minimizing the discrepancies between models at the macroscale, and at a local scale only for the forward and backward orientations (see Section 4).We find that, for the optimal parameters found, a great deal of non-trivial and quantitatively accurate phenomenology naturally emerges without the need to explicitly include it into the optimization process.
An example of such finding is the root-mean-square deviation r of the external stress and the system size dependence of stress fluctuation std(Σ ext x y ), reported in Figure 5a and b, respectively.We compute the root-mean-square deviation over external strain intervals as where the average is performed over intervals of width ∆γ ext , centered at different positions γ ext of the stress-strain curve in the stationary regime.r (∆γ ext ) increases and saturates to the std(Σ ext x y ), the latter decreasing with system size as ∼L −0.94 .We find that both stress correlation and standard deviation show excellent quantitative agreement with molecular statics simulations.
Furthermore, in the elasto-plastic model specification we have defined a statistical relation between the amplitude of a slip event and the resolved shear stress τ on the active slip system at the moment of activation (7).To keep the model simple, we avoided a detailed description of the slip events and chose (7) based solely on fulfilling the physical constraint of avoiding negative dissipation.Nonetheless, when coarse-grained stress drops δ τc are measured, we observe an emergent non-linear relation between 〈δ τc 〉 and 〈 τc 〉.This relation is shown in Figure 5c for the quenched state and the stationary forward and backward directions, measured with R = 30.The elasto-plastic model reproduces accurately the relation found in the atomistic system.
On the other hand, the elasto-plastic model consistently reproduces the anisotropy in the stationary state observed in the atomistic measurements for the different local properties studied, as discussed in Section 4.2.Note that such anisotropy is not a built-in ingredient of the model.Neither the laws of slip system renewal nor slip event performance consider a privileged orientation.Rather, the anisotropic response emerges naturally from the system's dynamics in the presence of the external load.This effect can be understood as the result of statistical hardening, i.e., the growth of the yield stress that occur due to a statistically biased slip activation.For simplicity, but without loss of generality, we consider the anisotropic response in terms of only forward and backward responses.In steady state, forward-oriented slip systems (i.e., θ ≈ α ext = 0°) experience a higher resolved shear stress.Consequently, forward-oriented slip systems suffer a statistical rise of slip threshold while backward-oriented (i.e., θ ≈ α ext + 90°) systems do not.
To further assess this plasticity-induced anisotropy, we perform a reloading test for different orientations.To this end, systems are driven until the stationary state and from there are unloaded until they bear no external stress, Σ ext x y = 0.While in the stationary state the behavior of the system is history-independent, once unloaded the future response depends strongly on the past deformation history, since the local yield stresses have been subjected to an anisotropic bias due to plastic deformation.When the unloaded systems are reloaded with α ext = 0°, we observe a nearly elastic-perfectly plastic response (Figure 5d).However, when the system is loaded in the reverse direction α ext = 90°, we observe a much softer behavior showing a slow strain-hardening.The origin of this Bauschinger effect has been shown to come from the inversion of the polarization in the low distance to threshold population during unloading [49].As reported in Figure 5d, we find a remarkable quantitative agreement between the atomistic measurements and the predictions of the elasto-plastic model regarding the plasticity-induced asymmetrical response to an external load.

Successes, failures and key ingredients
The elasto-plastic model reproduces all the phenomenology investigated in the stationary state with a relative aggregated error of L stat ≈ 16%.Specifically, the model reproduces the macroscale behavior with an excellent agreement for all the element lengths l investigated, despite the diminishing spatial resolution as l increases.We find a good quantitative agreement at the mesoscale with R = 30.However, when mesoscale properties are measured in the forward direction, i.e. aligned with the external load, such measures do not restrict the range of optimum mesoscale parameters with respect to the macroscale estimation.To further restrict the range, we measure the mesoscale properties oriented against the external load, which allows us to establish a specific combination of optimum mesoscale parameters.Moreover, the quenched state results generated by the Kinetic Monte Carlo method are of similar quality as the measures performed in the stationary state.These successes have been thoroughly discussed in the previous sections.
Despite the success of the elasto-plastic model, some disagreements are worth mentioning.We find that for the acceptable combinations of mesoscale parameters (i.e., that fulfill the constraint of matching the macroscale behavior), the elasto-plastic model systematically underestimates the distances to threshold in the forward direction in the stationary state (Figure 4h).Moreover, we find a wrong scaling in the left tail of the probability density, with a negative exponent in contrast with the marginal stability expected in elasto-plastic models and observed in our atomistic results.This discrepancy implies that, in the stationary state, elasto-plastic elements are closer to mechanical instability than atomistic regions are.We note that the non-coarse-grained single element response exhibits a scaling which is qualitatively right, with a positive exponent.
On the other hand, measurements in the backward direction α = 90°are, in general, less satisfactory than in the forward direction α = 0°.A backward behavior with a lower fit quality is especially appreciable in the scaling of the distance to threshold 〈∆ τc 〉 (Figure 4h).Since the distances to threshold and the forward-backward anisotropy result from self-organization, the mentioned discrepancies suggest that a new ingredient that modifies the dynamical evolution of the system might be required in the elasto-plastic description.
A general trend in our results is the deviation of the elasto-plastic results from the atomistic ones at the lower scales.In the local properties, a convergence between measurements at coarser scales is expected since model-specific details and possible measuring artifacts related to local shear tests become less relevant.Nonetheless, this trend is also reflected in the macroscale behavior when looking at fine details.Specifically, Figure 2d reveals a discrepancy in the low-stress increment tail.While the elasto-plastic distribution exhibits discrete stress increments induced by the driving protocol, the atomistic distribution is smoothed.Such discrepancy in the low-stress increments results from elastic heterogeneity at the smallest scales, which leads to external stress increments that fluctuate around discrete values imposed by the discrete external strain increments.Moreover, the root-mean-square deviation of external stress (Figure 5a) deviates below strain windows of ∆γ ext < 10 −3 .This value is the lower bound above which both stress-strain curves exhibit a similar correlation structure between stress increments and drops.
Even before questioning the validity and the ingredients of the elasto-plastic model, other sources of error can be discussed, first of which is the accuracy of the atomistic local yield stress method itself.Since the method uses rigid boundary conditions, it is likely to overestimate the yield stress and undereste the stress drops.This is in line with the underestimation of the local yield stress below R < 15 in the elasto-plastic model, and consistent with the recent work done by Liu et al. [50].Moreover, the measurement is purely local and does not take into account the elastic heterogeneities influencing the effective mechanical loading of each zone.All of these effects are expected to vanish asymptotically as the patch size increases, which is consistent with our results.These limits highlight the interest of developing flexible loading conditions for the atomistic local yield stress method, which will be the subject of future works.
In order to establish the impact of each ingredient employed in the elasto-plastic model, we finally quantify the different relative errors for various model flavors (see Appendix A.8 for a detailed description).As reported in Figure 6, these errors are classified into three categories corresponding to macroscale and mesoscale properties in the forward and backward directions.
Elastic heterogeneity has been shown to become progressively more important for the correct representation of elasticity in glasses as the coarse-graining scale decreases [36].Therefore, it is the first candidate ingredient to enhance the agreement between the magnitudes mentioned above and their respective atomistic measurements.As discussed before, elastic heterogeneity leads to small disagreements between the atomistic and the elasto-plastic model at the lowest scales.Thanks to our FEM computation of the elastic fields, we study the impact of inhomogeneous but non-evolving elastic properties and conclude that it has a small effect on the measurements performed on the elasto-plastic model and specifically does not improve systematic disagreements such as the distances to threshold.
In our mesoscale description, quenched elastic inhomogeneity leads only to subtle variations of the optimal parameters and a very slight improvement (Figure 6A), but the global picture is not altered with respect to the homogeneous case.An explanation might be that the model is sensitive to very few fundamental ingredients, namely disorder in slip activation and amplitude of elastic coupling, while the physical origin of such ingredients does not alter the model dynamics.Since the leading effect of elastic inhomogeneity is an increment in the local stress fluctuations, the right fluctuation might already be captured, in an effective sense, by the calibrated homogeneous model.Along similar lines, a more accurate description based on evolving elastic properties induced by structural evolution might, through the correspondence between elastic inhomogeneity and eigenstrain [57], have a similar impact as increasing the fluctuation of slip event amplitudes.Considering correlations between slip thresholds and elastic properties might also lead to further improvements in the model results.
We introduced pressure sensitivity by considering a Drucker-Prager-like criterion.As shown in Figure 6B  direction.This finding is surprising since pressure sensitivity is known to be present in the atomistic system [48], and suggests that pressure effects are not correctly captured by a Drucker-Prager criterion alone as discussed in [10,58].An improvement would involve considering a local criterion for permanent local dilation and contraction based on free volume dynamics [59].
Along similar lines, we studied the impact of diverse extra model features that do not bring new fundamental ingredients.Namely, we considered (i) a fluctuating number N of slip system per elasto-plastic element; (ii) a typical plastic strain amplitude controlling correlations between the previous and renewed thresholds; (iii) an explicit built-in anisotropy in the slip thresholds as a function of orientation, understood as a bias in the structural renewal process induced by external stress.These ingredients induce minimal changes in the overall picture presented in this paper, and we find that the mentioned deficiencies in the comparison between both models do not improve.
We finally follow the opposite approach by assessing the quality of simpler models.To this end, we implement a scalar model with only two slip systems of orientations 0°and 90°, fixed slip amplitudes, and spatially homogeneous non-evolving slip thresholds (Figure 6H).In this case, the only source of disorder is the quenched pre-stress field.Consequently, the stationary regime exhibits quasi-deterministic features, and the atomistic macroscale behavior cannot be reproduced accurately.If we restore the original stochastic renewal of the slip thresholds (Figure 6F), a calibration of quality only slightly below the original model C is possible (note that this scalar model cannot, by construction, reproduce the behavior for all the orientations shown in the left columns of Figures 4 and 8).However, if instead, we restored the stochastic slip orientations and amplitudes (Figure 6E), the emergent plasticity-induced anisotropy is lost since the model fails to reproduce the backward response.Therefore, stochastic renewal of slip thresholds is the most physically relevant source of disorder in the system for reproducing the emergent anisotropy observed in the atomistic glasses.
The original model C uses an Eshelby-type elastic interaction obtained by solving the stress balance equation.In this case, slip events induce anisotropic stress variations with alternating signs as a function of orientation, which lead to spatial correlations in the stress field and to strain localization [23,60].To assess the role of spatial correlations, we consider a random interaction kernel obtained by element-wise randomly shuffling the stress variations induced by slip events in the rest of the system.This kernel maintains the same statistical properties as the original one except for spatial correlations.As shown by Figure 6D, a calibration of quality only slightly below the original model C is possible.We take a step further and consider a homogeneous mean-field interaction by averaging the stress variations induced by slip events in the rest of the system (Figure 6I).In this case, the atomistic macroscale behavior cannot be reproduced accurately.We conclude that stress fluctuations induced by slip events play a fundamental role in the steady state flow stress response, while spatial correlations do not.We note, however, that spatial correlations are mandatory to reproduce the strain localization observed in the transient regime observed in more stable glasses [47].
For all the model versions, the emergent relations depicted in Figure 5 remain qualitatively valid.The mean-field and the simplest model (Figure 6I and H, respectively), however, fail to reproduce the root-mean-square deviation of the external stress (Figure 5a) and thus the correlations in the stress-strain curve.

Conclusions
We have presented a comparison between two very distinct approaches to simulating the plastic deformation and structural evolution of glasses.On the one hand, we consider an atomistic model (Section 2), based on the detailed knowledge of individual atomic trajectories.On the other hand, we present an elasto-plastic model (Section 3), which coarse-grains details into a continuum mechanics description with local slip systems and random structural properties.We considered strain-driven systems in athermal quasistatic shear (AQS) conditions, which reduces the complexity of the process by, e.g., ruling out temperature and inertial effects, and therefore the number of ingredients necessary for a successful elasto-plastic simulation.
The method for measuring local properties proposed by [41] was extended here to an elastoplastic implementation.We compared the results obtained with the same method but implemented with two different models describing the same system.Despite the differences in fundamental building-blocks and time and length scales, the method allowed us to compare the mechanical behavior independently of microscopic details.Thus, its applicability is wide and beyond atomistic simulations, and could be generalized to other models that can resolve microstructural details.Consequently, it is well-suited for the calibration of coarse-grained models based on data obtained at lower scales.
The elasto-plastic model reproduces the macroscale behavior in AQS conditions with a high degree of accuracy, although macroscale behavior alone is insufficient for determining a unique set of optimum mesoscale elasto-plastic parameters.To this end, we probed portions of the system and measured their mechanical response in isolation.By comparing observations at different scales with the atomistic system, we established the optimal values of the parameters.The model can reproduce the main phenomenology associated with plastic activity and structural evolution observed in the atomistic samples, such as the emergent anisotropy in the yield response induced by previous plastic deformation history.By calibrating the model, we have found C. R. Physique -2021, 22, n S3, 135-162 an optimum mesoscale length within the range l = 4.4-6.6.For coarser element lengths, the model loses spatial resolution but still reproduces macroscale behavior with an excellent quality.
Despite the generally good quantitative agreement between models, we have shown that some discrepancies at the local scales cannot be solved within the proposed model even for the coarser scales investigated.Moreover, simplified versions of the model are still able to reproduce the rich emergent phenomenology investigated.Therefore, it seems sensible that alternative novel ingredients are necessary to enhance the agreement between elasto-plastic descriptions and atomistic measurements at the mesoscale.Possible extra ingredients include using a more complex structural renewal process or the evolution of local elastic properties and their correlations with local yield thresholds and slip amplitudes.Moreover, the assumption of small deformations is only justified by model simplicity and computational performance, but is at odds with the large deformations undergone by the atomistic system when measurements are performed in the stationary state.Along similar lines, the consideration of advection [19] would capture material flows and rotations, a kind of structural evolution that is entirely missing in the presented elasto-plastic approach.Extra ingredients used in the literature that are missing in the present work are a distribution of element sizes and shapes in an unstructured mesh [26,53,61], rounded potentials [62] for the activation of slip events, free volume creation/annihilation dynamics [59], the effects of viscosity and finite plastic relaxation times [25,50] or the effects of shear wave propagation [34,53,63].
Along similar lines, future research is expected to be based on quantitative optimization of coarse-grained models that include diverse physically-motivated ingredients.Moreover, designing new optimization techniques, definitions of error between models, and improved methods for coarse-graining microscopic observations are natural steps forward to build upon the present work.The optimization of coarse-grained models through quantitative comparison with lowerscale reference data is a wide-open area, with the potential to impact the multi-scale modeling of mechanical properties of amorphous materials and enhance our understanding of the complexity emerging from the mechanisms at play.

A.2. System sizes, number of slip systems and patch sizes
In our analyses, we consider mesoscale elements of different length l .The physical domain represented by the elasto-plastic model has, in atomistic units, a side of length L = 98.8045.
To match the domain, we construct elasto-plastic lattices with a different number of elements L EP × L EP .We compute the number of slip systems as N = ρ s l 2 , where ρ s = 0.8.However, we round N to its closest value multiple of 4 to distribute the slip orientations as explained in Section 3.2.Elasto-plastic patches are formed by n × n elements (Section 3.5), and each patch has an associated equivalent radius.Especially important is the patch with an equivalent radius of R = 30, which we use for calibrating the elasto-plastic model.The values used for these parameters are reported in Table 1.

A.3. Initial configuration-quenched super-cooled liquid
Before applying the strain-driven protocol detailed in Section 3.4, we must define the initial system configuration.To mimic the atomistic protocol as faithfully as possible, we consider an equilibrated super-cooled liquid (ESL) at a finite temperature T .The liquid is instantaneously quenched, and the resulting state is used as the initial system configuration.To this end, we simulate the atomistic rearrangements taking place in the ESL state within our elasto-plastic framework as a series of slip events that can spontaneously occur via thermal activation [64].
Since we aim at a coarse-grained description of the process, we simulate thermal slip activation using the Kinetic Monte Carlo (KMC) method [26,27,29,56].
The KMC method requires only knowledge of all the possible transitions to a next state and the energy cost associated with each one.Intermediate details of the transition process are not resolved.In the elasto-plastic model, a transition corresponds to the activation of a slip event among all the existing slip systems.Thus, one is probabilistically chosen from all the possible slip events based on the energy cost ∆E associated with their activation.We consider a thermal slip activation rate where the index n refers to each slip system present in the system in a sequential manner, i.e., n = 1, 2, . . ., N L 2 where N is the number of slip systems per element and L 2 is the number of elasto-plastic elements.We write the activation energy in terms of stress as ∆E n = ∆τ c n V a , where V a is a microscopic activation volume estimated to be ≈1 from molecular dynamics.Assuming thermal activation as a Poisson process, the total activation rate is given by ν tot = N L 2 n=1 ν n .We aim to select a slip event with a probability proportional to its activation rate.To this end, first we compute the activation probability as P n = ν n /ν tot , from where we create a list of partial sums S n = n i =1 P i .Then, a random number r from a uniform distribution in (0, 1] is drawn.We look in the list S n for the first index k such that S k > r .Slip event k is then chosen by the KMC method.Since the probability for r of landing in the portion of the partial sum S n is proportional to ν n , the algorithm chooses the slip event n with probability P n [65].The chosen event is then performed (Section 3.3), and the elastic fields are updated.After each thermally-activated event, an athermal adiabatic avalanche might be triggered, as described in Section 3.4.For thermallyactivated events, we set the maximum amplitude as γ max (τ c ) (6), since the local stress τ must reach, at the moment of activation, a value of at least τ c as result of thermal fluctuations.
Before starting the KMC simulation, the system is initially given random slip angles θ and thresholds τ c drawn from their respective distribution (see Section 3.2).We consider no prestress Σ 0 field and an initially homogeneous zero elastic strain ε el = 0. Slip events are thermally activated under a fixed external strain γ ext = 0 until the slip thresholds, and the local stress field become statistically stationary, at which point the system has been equilibrated at temperature T .
We mimic the atomistic instantaneous quench and subsequent relaxation by setting T = 0 in the elasto-plastic system and performing an athermal mechanical relaxation, as described in Section 3.4.The resulting state is used to create the initial configuration for the athermal quasistatic strain-driven protocol.Specifically, the local stress field of the quenched liquid is used as the pre-stress field Σ → Σ 0 for the athermal solid, and the elastic strain is homogeneously reset to zero ε el = 0. Similarly, the slip angles and slip thresholds are used as the initial ones for the AQS simulation.
During the calibration process detailed in Section 4, we find that the properties of the stress field in the quenched state depend almost exclusively on the equilibration temperature T .On the other hand, the stationary state properties are independent of T , since they do not depend on the initial conditions.Consequently, we calibrate the value of T in a straightforward manner based only on the quenched state stress field properties.For the optimal parameters reported in the main text, we find for l = 6.6 a value T = 0.25 V a , a very reasonable value when compared with the molecular dynamics equilibration temperature T = 0.351.

A.4. Negative threshold
We note that Figure 4(e,f) report negative average yield stress, 〈 τc 〉 < 0. Although slip systems have positive-definite slip thresholds τ c > 0 the definition of coarse-grained yield stress τc = M(α) : Σc (Section 3.5) allows the measurement of negative values.To understand this, let us consider a macroscale system in its stationary state, loaded along α ext = 0°.If we unload the system until the first instability occurs, we need to perform, given its macroscopic size, a very small load increment with orientation α ext = 90°.Thus, when the instability occurs, the system's overall stress state remains oriented with α ext = 0°.Consequently, the resolved shear stress at the moment of instability on the plane α ext = 90°is negative, i.e., the yield stress is negative.
We can find a criterion for negative yield stress more rigorously by considering a shear test of orientation α performed on a patch belonging to a system at the stationary state.The initial resolved shear stress along the tested orientation is τi (α) = M(α) : Σi .As a function of the von Mises stress Σvm = (1/2 Σ : Σ ) 1/2 it can be written as τi (α) = M(α) : M( β) : 2 Σvm .Taking into account that ∆ τc = τc − τi and M(α) : M( β) = cos(2(α − β))/2, the condition for a negative yield stress τc < 0 is Therefore, it becomes increasingly more likely to measure a negative yield stress τc (α) < 0 the more the shear tests orientation α differs from the orientation β of the patch initial state.On the other hand, a negative threshold becomes more likely to occur when the distance to threshold ∆ τc (α) is low compared to the patch's initial shear stress amplitude Σvm .Moreover, since the distances to threshold decrease, on average, with R, negative thresholds are more likely measured at bigger scales.

A.5. Local fluctuations
We reach similar conclusions for the standard deviation of the explored local properties reported in Figure 8 as we did for the averages in the main text (see Figure 4).Namely, the elastoplastic and atomistic model converge for big coarse-graining scales R. Nonetheless, we observed unexpected local minima in the vicinity of α ≈ ±45°, which we attribute to artifacts induced by the FEM quadrilateral structured mesh.The quench state is successfully reproduced from our KMC approach.

A.6. Backward local distributions
The probability densities of the explored local properties measured in the backward direction reported in Figure 9 are in good agreement in the quenched state and a qualitative one in the stationary state.The local drops exhibit an excellent agreement in both cases.

A.7. Limit to the amplitude of the slip events
When a local slip event takes place an elasto-plastic element, a uniform eigenstrain is prescribed to the element.Therefore, a deforming element can be interpreted as an Eshelby inclusion with uniform eigenstrain embedded within an elastic matrix.As shown by [54], there exists a limit to the amount of eigenstrain an Eshelby inclusion can undergo to prevent negative dissipation.
In the case of uniform eigenstrain, such limit is given by C. R. Physique -2021, 22, n S3, 135-162 To increase the numerical accuracy of our approach, we compute numerically the Eshelby tensor S of our quadrilateral finite elements.To do so, we prescribe local eigenstrains, solve the stress balance equation with the FEM and compute the induced strain response.The same operation carried out for different eigenstrain tensors allows us to obtain the components of S from a system of linear equations.

A.8. Alternative models
We have quantified the fitting error of different models with l = 6.6.Each model differs in certain implementation details and, therefore, in the optimal mesoscale parameters.This section details the differences between models and gives the values of the estimated optimum parameters.The fitting error for each compared quantity is reported in Figure 10.

Figure 2 .
Figure 2. Comparison of the macroscale behavior of the atomistic (blue) and elasto-plastic (red) models in the stationary regime (γ ext > 0.5).(a) Single simulation stress-strain curve.Probability densities of: (b) external stress values Σ ext x y , (c) drops δΣ ext x y and (d) increments ∆Σ extx y induced by discrete steps of ∆γ ext = 10 −4 .The results shown correspond to the optimal parameters with l = 6.6.

Figure 3 .
Figure 3. Relative errors in the fit as a function of the exponent k of the Weibull distribution of slip thresholds for different mesoscale element lengths l for: (a) the macroscale and (b) the aggregated macroscale and mesoscale properties.(c) Minimum of the aggregated error L stat as a function of l , where error bars indicate one standard deviation.(d) Individual components of L stat for l = 6.6.

Figure 4
shows von Mises stress Σvm = (1/2 Σ : Σ ) 1/2 (Figure 4 first row), yield stress τc (Figure 4 second row), distance to the threshold ∆ τc (Figure 4 third row) and stress drop δ τc (Figure 4 fourth row) in the quench and stationary states.Their probability distributions and dependence on orientation measured with R = 30 are shown in the first and third column, respectively.In the second column, we plot the averages measured for different scales R. The same analysis carried out for the standard deviations is reported in Appendix A.5.A detailed overview of the fitting errors for the individual magnitudes is shown in Figure 10(model C).

Figure 4 .
Figure 4. Comparison of the local properties of the atomistic (blue) and elasto-plastic (red) models in the quench and stationary regimes: von Mises stress Σvm and resolved shear stress τi on the shear test plane (first row), yield stress τc (second row), distance to threshold ∆ τc (third row) and stress drop δ τc (forth row).The columns correspond to the probability distribution functions for R = 30 (left), the averages as a function of R (center) and the averages as a function of the shear orientation α for R = 30 (right).The results shown correspond to the optimal parameters with l = 6.6.

Figure 5 .
Figure 5.Comparison of the behavior of the atomistic (blue) and elasto-plastic (red) models.(a) Root-mean-square deviation r of the external stress.(b) Finite-size scaling of the external stress standard deviation std(Σ ext x y ).(c) Average local stress drop versus average yield stress with R = 30.(d) Bauschinger tests: reloading (α ext = 0°) and reverse loading (α ext = 90°) for systems previously unloaded from the stationary state.The results shown correspond to the optimal parameters with l = 6.6.
, the error increases with respect to the original model, especially in the backward C. R. Physique -2021, 22, n S3, 135-162

Figure 6 .
Figure 6.Relative errors between the atomistic simulations and different versions of the elasto-plastic model (A-I) computed with their respective optimum parameters with R = 30 and l = 6.6 in the stationary state.Shown are the errors in the macroscale behavior ( Lmacro ) and mesoscale properties aligned with the external load ( L+ meso ) and against it ( L− meso ).The different model ingredients are: (i) randomly oriented slip systems with their own threshold versus scalar model with a single threshold value per element.(ii) Stochastic threshold renewal versus fixed value.(iii) Stochastic slip amplitudes versus fixed value.(iv) Eshelby, randomly shuffled Eshelby or mean-field interaction kernels.(v) Non-evolving and inhomogeneous versus homogeneous shear modulus.(vi) Drucker-Prager versus pressureindependent yield criterion.

Figure 7 .
Figure 7. Optimal slip threshold scale parameter λ as a function of the number N of slip systems in mesocale elements of length l = 6.6.

Figure 8 .
Figure 8.Comparison of the local property standard deviations of the atomistic (blue) and elasto-plastic (red) models in the quench and stationary regimes: von Mises stress Σvm and resolved shear stress τi (first row), yield stress τc (second row), distance to threshold ∆ τc (third row) and stress drop δ τc (forth row).The columns correspond to the averages as a function of R (left) and the averages as a function of shear orientation α for R = 30 (right).The results shown correspond to the optimal parameters with l = 6.6.

Figure 9 .
Figure 9.Comparison of the probability density of the local stress drop for the backward direction α = 90°of the atomistic (blue) and elasto-plastic (red) models in the quench and stationary regimes for R = 30: (a) local yield stress, (b) distance to threshold and (c) and stress drops.The results shown correspond to the optimal parameters with l = 6.6.

Figure 10 .
Figure 10.Relative errors in the stationary state between the atomistic simulations and different versions of the elasto-plastic model computed with R = 30 and l = 6.6.

Table 1 .
Element lenght l , elasto-plastic lattice size L EP , number N of slip systems per element and elasto-plastic patch size n with equivalent radius of R = 30