Properties of the density of shear transformations in driven amorphous solids
PProperties of the density of shear transformations in driven amorphous solids
Ezequiel E. Ferrero and Eduardo A. Jagla Instituto de Nanociencia y Nanotecnolog´ıa, CNEA–CONICET,Centro At´omico Bariloche, (R8402AGP) San Carlos de Bariloche, R´ıo Negro, Argentina. Centro At´omico Bariloche, Instituto Balseiro, Comisi´on Nacional de Energ´ıa At´omica, CNEA,CONICET, UNCUYO, Av. E. Bustillo 9500 R8402AGP S. C. de Bariloche, R´ıo Negro, Argentina (Dated: November 25, 2020)The strain load ∆ γ that triggers consecutive avalanches is a key observable in the slow defor-mation of amorphous solids. Its temporally averaged value (cid:104) ∆ γ (cid:105) displays a non-trivial system-sizedependence that constitutes one of the distinguishing features of the yielding transition. Detailsof this dependence are not yet fully understood. We address this problem by means of theoreticalanalysis and simulations of elastoplastic models for amorphous solids. An accurate determinationof the size dependence of (cid:104) ∆ γ (cid:105) leads to a precise evaluation of the steady-state distribution of lo-cal distances to instability x . We find that the usually assumed form P ( x ) ∼ x θ (with θ beingthe so-called pseudo-gap exponent) is not accurate at low x and that in general P ( x ) tends to asystem-size-dependent finite limit as x →
0. We work out the consequences of this finite-size de-pendence standing on exact results for random-walks and disclosing an alternative interpretation ofthe mechanical noise felt by a reference site. We test our predictions in two- and three-dimensionalelastoplastic models, showing the crucial influence of the saturation of P ( x ) at small x on the sizedependence of (cid:104) ∆ γ (cid:105) and related scalings. Punctuated dynamics is inherent to many out of equi-librium driven systems. When energy is loaded at a smalland fixed rate, the nature of the system is such that thisenergy is dissipated in sudden bursts of activity typicallycalled slip events or avalanches. This kind of systems arereferred to as displaying a stick-slip dynamics. Examplesinclude the relative motion of tectonic plates giving riseto earthquakes [1], the sliding of charge density waves [2],the driven movement of a magnetic interface in thin mag-netic films [3], the intermittent motion of rain dropletson a windshield [4] and the plastic rearrangements oc-curring in amorphous solids under a slow and sustainedstrain increase [5]. In all these cases, a stationary situ-ation is established in which, on average, the stress (orenergy) increase during quiescence periods is equal to thestress (or energy) drop released during avalanches.Suppose that we drive a system with stick-slip dynam-ics on its steady state, and we are interested in the statis-tics of strain increases needed to produce a new slip event,for systems of different sizes. If the system consists on N ‘blocks’ that can be locally destabilized, one expects thatthe load needed to trigger the weakest block scales with1 /N . This is, if we double the system size, the closest in-stability will be halfway apart in terms of strain increaseneeded. Equivalently, if we drive the system at a smallfinite rate, the pace at which we observe slip events dou-bles when we double the system size. More rigorously,if avalanches have a maximum extent that does not di-verge as the system size goes to infinity, then the systemis extensive. The previously mentioned balance betweenaccumulation and release of energy then implies that ifthe system size is doubled, the average load increase thathas to be applied to generate a new avalanche is halved.While this is the case for most stick-slip phenomena (e.g.,friction, depinning, wetting, etc.), the behavior of amor-phous solids under deformation disobeys this logic. In the deformation of amorphous materials, if we doublethe system size, the rate at which we observe slip eventsdoes not double. It increases, but less; it is sub-extensive in the system size. In other words, to trigger the next slipone needs to load more than expected . As a consequence,when the system finally yields, the slip of a single blockis not enough to compensate the load excess and systemspanning avalanches of plastic events emerge. Therefore,the plastic activity is rarely confined to localized plasticevents and, instead, it is mostly originated in extendedstructures [6] This points clearly to the non-extensivenessof the problem. In fact, if we consider conversely that thedynamics of the problem produces system size spanningavalanches, then a doubling in the system size would notduplicate the number of avalanches.It is now well established that the statistics of the meanstrain load (cid:104) ∆ γ (cid:105) needed to trigger consecutive avalanchesin the steady state of quasistatically driven amorphoussolids has profound consequences on the criticality of theyielding transition [7–9]. In particular, its finite-size scal-ing is expected to be a manifestation of the distributionof putative shear transformation zones P ( x ) ( x ≡ Σ th − Σstands for the local ‘stress distance’ to the local yieldingthreshold Σ th ) and bounds through scaling relations thepossible exponents governing the avalanche size distri-bution [5, 10–14]. While consensus on this scaling beingsub-extensive prevails, i.e., (cid:104) ∆ γ (cid:105) ∼ N − α , with 0 < α < α andits justification. On one hand molecular dynamics (MD)simulations [7–9] of model glasses under quasistatic de-formation support an universal value of α (cid:39) /
3, validboth in d = 2 and d = 3 dimensions. On the other hand,elasto-plastic (EP) models for amorphous solids [10, 11]display dimension-dependent values of α . a r X i v : . [ c ond - m a t . d i s - nn ] N ov Despite large advances in the field, theoretical argu-ments have not taken yet in full account the statisti-cal significance of the (cid:104) ∆ γ (cid:105) sub-extensivity in the steadystate. Such a scaling has been well accounted as a justifi-cation for system-spanning avalanches of plastic activityin the system; but the fact that it also implies an inher-ent discrete evolution for the ‘local distances to threshold’ x was disregarded. In this work, we address this issue,presenting an alternative and consistent picture for thefinite-size scaling of x min ≡ min i x i . Standing on theground provided by previous works [15–18], we interpretthe evolution of the stress (and thus also x i ) in a genericregion of the system as an effective stochastically-drivenrandom walk. Working out the finite-size scaling of therelevant discrete jump in this walk we derive a genericscaling law (cid:104) x min (cid:105) ∼ N − α with α = 2 /
3. In doing so,we revisit the significance and shape of the distribution P ( x ), which shows a finite limit P = lim x → P ( x ) thatscales as P ∼ N − α in the thermodynamic limit. Thisscaling occurs independently of the eventual value of θ observed at intermediate values of x where P ( x ) ∼ x θ can be fitted.In Sec. I we refresh the subject under discussion in amini-review. In Sec. II we motivate and perform an anal-ysis in terms of simple random-walkers problems withexact solutions to understand the effect of discrete steps.In Sec. III we rationalize the collective effect of plasticevents during avalanches as an effective mechanical noisewith a discrete step effect on the ‘walks’ of local stresses.In Sec. IV we test our hypothesis in extensive simulationsof 2D and 3D elastoplastic models, presenting rigorous fi-nite size analysis for (cid:104) x min (cid:105) and P ( x ) in different cases.Finally, in Sec. V we summarize our results, that we be-lieve allow to construct a consistent scenario for whatseemed a priori contrasting results in literature. Figure 1. Steady-state stress-strain scheme in the quasistaticshear deformation of an amorphous solid.
I. OVERVIEW OF THE SUBJECT
Let us start by briefly reviewing the main conceptsand literature results on this topic. In essence, onecan think of a yielding material as a fully-connected setof elastoplastic blocks characterized by a local stress Σ i or, equivalently, a local distance to the stress threshold x i = Σ th i − Σ i . These blocks evolve according to a globalload that drives the x i values towards zero. When a par-ticular x j reaches zero, the block yields, reaching a newequilibrium position (at some new, positive value of x j )while at the same time producing (via elastic interac-tions) a modification of the values of other x i all acrossthe system. We say that these are ‘mechanical kicks’given to the blocks each time one of the blocks yields.The yielding of block j may produce (due to the me-chanical kicks) the yielding of other blocks in cascade.This is the origin of avalanches in the system that char-acterize the dynamics. Because of this avalanche domi-nated dynamics, the stress-strain evolution in the systemhas a qualitative form as depicted in Fig.1. Once in itssteady state, a driven amorphous solid performs an in-terspersed sequence of load periods and slip events whenthe relevant stress component is monitored. A station-ary average value is expected for the stress on a steadystate. On top of this average value, stress fluctuationscontain information on the physics of the problem. Inparticular, the average strain increase of the loading pe-riods (referred to as (cid:104) ∆ γ (cid:105) ) and the average stress-dropduring the slip events ( (cid:104) ∆Σ (cid:105) ) must be proportional in astationary situation (see Fig.1), namely (cid:104) ∆Σ (cid:105) = B (cid:104) ∆ γ (cid:105) (1)with B an elastic constant.In a quasistatic athermal dynamics, the stress in-crement that needs to be applied to trigger a newavalanche is nothing but the minimum x i across the sys-tem, namely x min . Then, (cid:104) x min (cid:105) = (cid:104) ∆Σ (cid:105) . (2)Further, energy drops, quantifying the energy dissipatedduring plastic avalanches, can be easily related to thestress-drops as (cid:104) ∆ U (cid:105) = Σ Y (cid:104) ∆ γ (cid:105) V = Σ Y B (cid:104) ∆Σ (cid:105) V ; where V = L d = N is the system volume and Σ Y is the globalyield stress [7, 8]. In fact, the starting point of the cur-rent discussion can be traced back to a series of MD qua-sistatic simulations [7–9] where the following system-sizescaling laws were verified (cid:104) ∆ U (cid:105) ∼ N δ ; (cid:104) ∆Σ (cid:105) ∼ N − α (3)holding δ + α = 1, with δ (cid:39) / α (cid:39) /
3, both in twoand three spatial dimensions. Also [19] independentlyshowed compatible results. In Ref. [8], when thinkingon the distribution of possible plastic events, an ansatzwas introduced. Arguing that the distribution of energybarriers felt in a quasistatic loading protocol should growas a power-law, it was proposed for the distribution oflocal distances to instabilities ( x is the generic value ofall x i across the system) the expression P ( x ) ∼ x θ , (4)for small x and with θ >
0. Then the distribution of x min follows a Weibull distribution P ( x min ) ∼ x min θ exp( − N x min θ ) , (5)with its mean value scaling as (cid:104) x min (cid:105) ∼ N − / (1+ θ ) (6)which provides a justification for the scaling of Eq. 3,linking α and θ : α = 11 + θ . (7)Yet, notice that Eq.6 does not imply P ( x ) ∼ x θ . Thesmall argument power-law form of the Weibull distribu-tion for P ( x min ) was verified in the statistics of the ‘as-quenched’ state or isotropic solid state, both in d = 2and d = 3 dimensions [7–9]. And it is in fact for this casethat the ansatz (4) was proposed [8]. Nevertheless, thesepioneer MD simulations couldn’t easily access the wholedistribution P ( x ), and results where only presented for P ( x min ) (or P (∆ γ )).Luckily, soon after the problem was addressed byEP model simulations measuring the full P ( x ) distribu-tion [10]. In there, a plausible law P ( x ) ∼ x θ was foundnot only in the ‘as-quenched’ state but also at the crit-ical stress. The P ( x ) ∼ x θ ansatz was subsequently ex-tended to describe in EP simulations not only the steadystate [11], but also the transient regime [20] where astatistics of extended avalanches was equally observed.It was concluded that θ , and therefore α according tothe construction, should be dimension and system pa-rameter dependent, which was formalized in an analyticmean-field approach [21]. This theory has the virtue offormally catching a strain-dependence of θ , a feature thatis observed in the transient regime both in EP [20, 21]and MD [22–25] simulations. In such transient, the val-ues of θ observed are highly non-universal, dependingon system preparation, system parameters and dimen-sion [21, 23, 26].In the construction summarized in [21], α is expectedto follow the same trend as θ all the way from the ‘as-quenched’ state to the steady-state, keeping the relation α = 1 / (1 + θ ), and binding α to be also highly non-universal. Nevertheless, one naturally expects θ and α to stop depending on strain in the steady-state, and in-deed the literature has collected from the beginning ev-idence for such expectation [8, 10, 11]. Moreover, wehave recently showed that in that limit those exponentsare model-independent [15] for a large set of EP modelrules; they do depend on dimension though. So, at somepoint the variation of α and θ with strain should vanish. How that happens, may be a matter of theoretical dis-cussion itself. For the time being, we will focus on thelimit of large strains where a self-consistent and station-ary stick-slip phenomenon is expected to occur.Interestingly, in contrast with the case of ‘as-quenched’systems, the relation α = 1 / (1 + θ ) does not seem to holdso well in the numerical results of the steady-state in EPmodels. For example, in [11] θ is reported to be ∼ . ∼ .
35 respectively in d = 2 and d = 3, while α results form the x min scaling in ∼ .
67 and ∼ .
79 forthose cases [27]. More recent EP simulations [28] show α (cid:39) .
675 combined with θ (cid:39) / d = 2. And in [15]we have observed α (cid:39) / θ (cid:39) .
75 for 6 different d = 2 EP models, pushing the relation α = 1 / (1 + θ )even further away from validity. The apparent violationof such relation in the steady state is accompanied bytwo related observations. First, it is well known from thebeginning of this discussions that P ( x min ) ∼ ( x min ) θ doesnot show up in the steady-state [8, 10]; in fact that law,valid for the ‘as-quenched’ state, is rapidly suppressedas soon as the applied stress is finite [22, 23]. Secondly,recent numerical results in both MD simulations [25] andEP models [15, 28] have consistently made evident thatin the steady state P ( x ) displays a plateau at smallvalues of x (a non-zero base value, unmistakable in adouble logarithmic plot P ( x ) vs x ), and suggested thatthe finite-size scaling of (cid:104) x min (cid:105) can be dominated by thebehavior with system size of the finite asymptotic valueof P ( x ) at vanishing x rather than by the exponent θ .What seems to be clear, at least, is that the plain as-sumption P ( x ) ∼ x θ , which was somehow inherited fromthe ‘as-quenched’ phenomenology and carried by for allvalues of strain, is insufficient in the steady state. Nev-ertheless, for instance, a scaling relation based on Eq. 7,and linking the exponents τ, d f that describe the distri-bution of avalanche sizes with the exponent θ (namely, τ = 2 − θ θ dd f ) has been largely adopted in the the EPmodels literature [5, 12, 14, 29], always accepted withoutfurther justification and relying sometimes on generouserror bars for the exponents. Something is missing inthe understanding of what controls α , which might causethat even the latter relation among exponents should berevised. In this work we address the issue, admittedlylimiting ourselves to the steady-state case, were we ex-pect universal values of α [8, 15, 22]. A mean-field approach to yielding
In [21], Lin&Wyart extending a work by Lemaitre andCaroli [30] presented a mean-field approach which isbased in the assumption that the mechanical kicks pro-duced by yielding sites on every other site can be takenfrom a given distribution defined once and for all, inde-pendently of the state of the system, and, more impor-tantly, that this distribution is heavy-tailed. The mean-field dynamics can be described as follows. If at time t the block j yields (reaches x j ≤ x j (e.g. x j = 1) and the rest of the absorbing boundary external drive
12 3
Figure 2. Schematic wandering of the x i values. The externaldrive pushes everyone towards zero 1 ○ . When x j yields, it isre-injected 2 ○ at a positive value of x and all the remaining x i receive kicks ξ i ○ . These kicks can be either positive ornegative and of different intensities (as illustrated by the ar-rows’ different colors and lengths), according to the prescribeddistribution w ( ξ ). The signed nature of the mechanical noiseallows the system to ‘forecast’ the boundary [11] and createsa density depletion of P ( x ) close to x = 0. blocks suffer mechanical kicks ξ i taken form a distribu-tion w ( ξ i ) (see Eq. 9) with zero mean (cid:104) ξ i (cid:105) = 0, x j ( t + 1) = 1 ,x i ( t + 1) = x i ( t ) + ξ i − − x j ( t ) N − , (8)where the last term grants stress conservation and in thiscase the re-injection point has been chosen to be x = 1.A highly distinctive feature of yielding phenomena laysin the fact that the mechanical noise distribution w ( ξ i )comprises both positive and negative values. This fea-ture has its roots in the Eshelby response observed uponplastic events in amorphous materials and described inApp. A. Once this is guaranteed, these mean-field mod-els behave qualitatively like elastoplastic models of amor-phous solids [5]: There is a global yield stress Σ Y suchthat for Σ < Σ Y the dynamics eventually stops, cor-responding to the solid phase. For, Σ > Σ Y the dy-namics does not stop and is characterized by a globalstrain rate. The dynamics therefore can be rational-ized as ‘random walks’ of the elastoplastic blocks in the x -coordinate, with an absorbing boundary at x = 0 [21],as qualitatively depicted in Fig. 2. This representation ofthe yielding phenomenon allows us to start by analyzingsimple random-walk processes and extrapolate conclu-sions from there, see Sec. II.The signed nature of the mechanical noise gives rise toa density depletion of P ( x ) close to the absorbing bound-ary [11, 15]. In contrast to the case of a purely posi-tive interaction among sites where each destabilized blocktend to destabilize the others, the signed kicks allow someblocks to escape the boundary and survive longer withoutyielding. This is the hand-waving argument for the exis-tence of a pseudo-gap P ( x ) ∼ x θ with θ >
0. Formallysolving the stochastic problem of Eq. 8, Lin&Wyart con-cluded that θ depends continuously on the applied shearstress, non monotonically and without signs of universal-ity at the yield stress. The latter observation leaves littleroom for the expectation of universal exponents among different EP models, not to talk about MD simulations.We will contrast this view. Alternative views
The mean field construction in [21] is based on theassumption of a mechanical noise w ( ξ ) of the form w ( ξ ) ∼ | ξ | µ +1 (9)in the particular case of µ = 1, that the authors claim isthe only value with ‘physical meaning’ to be expected tooccur.The necessity of the particular value µ = 1 has beenrecently questioned [15–17]. In particular, it was arguedthe assumption µ = 1 in Eq.9 is not in agreement withthe observation of sub-extensive avalanches dominatingthe plastic activity in the quasistatic limit [15, 16]. Othervalues of µ with 1 < µ < µ (cid:39) / < µ < θ = µ/ < µ <
2, in-dependent on other parameters. Yet, a uniquely-definedvalue of θ = µ/ α observed in both MD and EP simulations.Remarkably, it has been recently observed quite clearlyin both EP [15, 28] and MD [25] simulations that thetrue shape of P ( x ) at small values of x deviates from apure power-law ∼ x θ , and has a finite limit P ≡ lim x → P ( x ) (cid:54) = 0 . (10)Namely, for any finite system size, P ( x ) has a plateau (when P ( x ) is presented in a logarithmic plot) at smallenough x . As we will argue in the following, the plateauin P ( x ) is originated in the discrete nature of the me-chanical noise that produces the “kicks” felt by x i . Thesekicks (that we consider to be generated by extended plas-tic avalanches elsewhere in the system) push each x i toperform a (non-Gaussian) random walk. In this scenario,it is the system size scaling of P what dominates the scal-ing of (cid:104) x min (cid:105) and controls the values of α and δ in Eq. 3,which now turn to be compatible with the independentexponent θ = µ/ P and (cid:104) x min (cid:105) and conciliates them with the existence ofa well defined value of θ that indeed describes an inter-mediate region of x values where P ( x ) ∼ x θ . x P ( x ) -2 -3 -4 -6 -4 -2 x -4 -2 P ( x ) -4 -2 x/ σ P ( x ) / σ ~x~x σ (a)(b)(c) Figure 3. Numerically determined probability distributionfor a variable x performing a standard (Gaussian) RW in theinterval (0 ,
1) with absorbing boundary conditions and ran-dom reinjection. Different curves correspond to different val-ues of the width of the single step distribution, as indicated.(a) Linear scale. (b) Log scale to emphasize the behavior atlow x . The straight line shows the expected asymptotic limitfor σ →
0. In (c) the axis are rescaled with σ to show thatthe value P (0) scales as σ . II. SIMPLE RANDOM WALKS AND THE P ( x ) PLATEAU
We analyze first a simple case. Consider a variable x i performing a random walk in the interval [0 , x i moves out ofthe interval, it is “absorbed” and re-injected in some ran-dom way [32]. In the case of a continuous time randomwalk (a Wiener process), and when the reinjection is doneproportionally to the local value of the probability, theform of the distribution of x i values observed along timein the steady state can be analytically computed to be P ( x ) = π sin( πx ). For small x it behaves as P ( x ) ∼ x ,i.e. θ = 1. If we consider N variables ( N (cid:29)
1) perform-ing the same random walk, the minimum among themwill be in the region in which P ( x ) is linear, and we will x P ( x ) -2 -3 -4 -5 -8 -6 -4 -2 x -4 -2 P ( x ) -2 x/ σ P ( x ) / σ / ~(x / σ ) ~x σ (a)(b)(c) Figure 4. Probability distribution for a variable x performinga RW with Hurst exponent H = 2 / , x . The straight line showsthe expected asymptotic limit for σ →
0. In (c) the axis arerescaled with σ to show that the value P scales as σ − / ( σ − / H for a generic H ). have P ( x min ) ∼ x min exp( − N x min ), and (cid:104) x min (cid:105) ∼ N − / .The random walks that we introduced in the previ-ous section to describe the phenomenological dynamicsof yielding are inherently discrete, and one needs to in-vestigate the consequences of this fact on P ( x ). In fact,for a finite step random walk, although the overall formof P ( x ) is the same as before, there is a small correctionat small x that depends on the step size, and has a strongeffect on the value of x min . Let’s think for a moment of aparticle performing a discrete random-walk characterizedby a step that is Gaussian-distributed, with a dispersion σ . Assuming the particle is at some position x at agiven step, the next jump makes x to be distributed as P ( x ) ∼ Θ( x )Θ(1 − x ) exp[ − ( x − x ) / σ ], where theHeaviside functions Θ appear because of the absorbingboundary conditions. We note that the value of P (0 + )is finite. It turns out that this effect remains in the fullsolution for the stationary form of P ( x ). So, the discretenature of the steps taken by x i suffices to explain thefinite limit of P ( x ) as x →
0. In Fig.3 we see the distri-bution of P ( x ) for Gaussian random walks with differentmagnitudes of the average elementary step, namely, dif-ferent width σ of the Gaussian ‘kicks’. Fig.3(a) shows thestationary distributions in lin-lin scale, Fig.3(b) showsthem in log-log scale, and the scaling proposed in Fig.3(c)shows that the value of P is proportional to σ .The situation is conceptually identical in the case inwhich we consider generalized RWs with a non-trivialHurst exponent H ; this is, random walks generated byjumps ξ drawn from a heavy-tails distribution of the form w ( ξ ) ∼ | ξ | H +1 , (11)for large | ξ | with 1 / < H <
1. Note first of all thatin this case, the ‘typical jump’ or distribution width σ cannot be defined as being variance of the distributionbecause of its heavy tails, but it can be alternatively de-fined as σ ≡ (cid:104)| ξ |(cid:105) . As it was the case for a Gaussianvariable, in the limit of vanishingly small jumps (i.e., σ →
0) RWs, the form of P ( x ) for x close to zero is stillexpected to be P ( x ) ∼ x θ , where now θ = 1 / (2 H ) [21].Yet, for finite σ , a finite value for P appears, as shownin Fig.4 for H = 2 /
3. For concreteness, in this numer-ical example we have taken the distribution w ( ξ ) to begiven by Eq.11 if | ξ | > ξ , and w ( ξ ) = 0 if | ξ | < ξ .This distribution has a width σ = (cid:104)| ξ |(cid:105) = ξ / (1 − H ).We see that the limiting value of P as a function of σ scales as P ∼ σ / (Fig.4c) . In the generic case with1 / < H < P scales as P ∼ σ θ . (12)This can be justified by noticing that close to x = 0, σ isthe only possible scaling quantity with the same dimen-sion as x . Then we can write P ( x ) = P f ( x/σ ) (13)with f ( u ) (cid:39) u (cid:28)
1. On the other hand, for x (cid:29) σ (but still ‘small’) we must have P ( x ) (cid:39) Cx σ with C inde-pendent of σ , therefore implying Eq. 12. In other words, σ marks a scale crossover below which the distributionof x values tends to a constant [33].Finally, notice that everything we have said for thesteady state distribution P ( x ) populated along time isalso true if we populate the distribution with the x i valuesof many independent walks in their steady state. A. N random walks without or with drift Let’s consider then N independent random walkerssubject to the following protocol. Now, starting forma condition where every x i is in the interval (0 ,
2) we -6 -4 -2 x -3 -2 -1 P ( x ) -5 -3 -1 x/ σ -1 P ( x ) / σ / ~x ~(x/ σ ) N (a)(c) -6 -4 -2 x -2 -1 -4 -2 x/ σ ~x N (b)(d) Figure 5. Probability distribution populated by N variables x i subject to a dynamics of random kicks taken from a heavytailed distribution like the one in Eq. 14, with Hurst exponent H = 2 / A = 0 .
1. Sites absorbed at the boundaries arere-injected in x = 1. Different curves correspond to different σ ∼ M − /µN ∼ N − / Left:
Without drift: (a) raw-data, (c)rescaled with σ to show P ∼ σ − / H . Right:
With drift. (b)raw-data, (d) rescaled with σ . look for the minimum x i , that we indicate as x min . Everysite is shifted by an amount − x min . The site resultingwith x i = 0 is re-injected in the box at x i = 1 and every-one updated by a (randomly) signed random quantity ξ taken from a distribution similar to Eq. 11 ( µ = 1 /H ) w ( ξ ) = AM N | ξ | µ +1 , (14)but with upper and lower cutoffs set to ξ up = (2 A/µ ) µ and ξ lo = (2 A/µ ) µ M − µ N for it to be normalized [21].Importantly, here M N is an N -dependent parameter, fre-quently chosen as N itself (see [21]). In the simulations ofthis toy model we will use M N = 1 /x min for reasons thatwill be clearer later on [34]. Every site resulting in x i ≤ x i ≥
2) after the random kicks is alsore-injected at x i = 1 (but not producing further kicks).The N walkers feel these kicks independently, yet theyare drifted globally by − x min after each kick update. Inorder to clearly identify the effect of such global drift, wewill also analyze the case where we avoid the global driftstep and simply: re-inject the site with the minimum x i ,give random kicks to everyone and further re-inject thosethat go out of the box.In both protocols, with and without drift, a steadystate is established after a transient and the resulting P ( x ) distributions are shown in Fig.5. We can see thatthe drift couples the dynamics of the walkers and pro-duces the effect of a ‘belly’ on the curves that delays thedecrease of P ( x ) as we sense x decreasing. The choice ofthe parameter A now becomes relevant. If A is small, thedrift effect overtakes good part of the P ( x ) distributionand it masks the power-law regime which gets difficultto determine, forcing us to simulate very large systems(or very small σ ). If instead A is big enough (closer to1) the drift effect is much diminished (data not shown).In any case, when a reasonably large power-law region isgranted, the θ exponent is preserved for any A , θ = µ/ < µ <
2. Notice that, despite this ‘belly’effect, the existence of a plateau at small x is unchanged,and the predictions P ∼ σ θ still holds, as can be seen inthe data collapse of Fig. 5(b,d).The dynamics that we have just described can bethough as a mean-field model for a system of elasto-plastic blocks with local thresholds where each of themfeels an external drive and a noise represented in w ( ξ ).We will now analyze a spatially extended system of driveninteracting blocks in this context. III. EFFECTIVE MECHANICAL NOISE OF ANINTERACTING SYSTEM AND THE P ( x ) DISTRIBUTION
Let us imagine a coarse-grained representation of anamorphous material under deformation, represented by ascalar stress Σ i on each block and local yielding thresh-olds Σ th i . The variables of interest will be the local dis-tances to threshold x i = Σ th i − Σ i . Our argumentationline is based on the analysis of the mechanical noise feltby a given site of such a system, caused by the plasticactivity elsewhere and governing the ‘wandering’ of x i .For the results of the previous section to be applicableto the present case, this noise must consist ideally of in-dependent, uncorrelated kicks. As previously mentioned,Refs. [21, 26] present a mean-field model considering kicksof a mechanical noise generated by single Eshelby events.We will refer to these kicks generated by single sites as‘elementary’ kicks. The approximation of Ref. [21, 26]describes qualitatively well the overall phenomenologyobserved in numerical simulations, but fails in predict-ing the exponents observed, at least below d = 4. Thisdiscrepancy was indeed ascribed to the presence of “di-mensional effects” or correlations between the elemen-tary kicks produced in different positions of the system.We believe that the quantitative predictive power of thiskind of analysis can be improved, still keeping the “mean-field” character of the approach, by noticing and takinginto account that elementary kicks are not independent.Elementary kicks produced by sites that participate ofthe same avalanche are highly correlated among them,but those from different avalanches are not. This factallows us to build a mean-field approach based on inde-pendent non-elementary kicks. One possible choice is todefine them as the integrated kicks given by avalanches,that in the quasistatic limit are by definition uncorrelatedevents.The fact that the uncorrelated mechanical noise underconsideration is produced by avalanches as a whole isthe reason why now µ in Eq. 14 can be different from thevalue µ = 1 that was obtained considering the effect ofhypothetical uncorrelated elementary kicks instead [21]. Actually, this alternative approach of avalanche-levelnoise was already followed in [15, 17]. Simulations ofdifferent EP models in two dimensions produce in a testsite a noise characterized by a Hurst exponent H (cid:39) / µ (cid:39) /
2. With thatbeing proved to be effectively the case for a fully interact-ing system [15, 17], we cannot expect anything differentfor its full distribution P ( x ) than the features discussedin previous sections. Finite size scaling of the P ( x ) plateau The mechanical noise represented by Eq. 14 contains asa fundamental parameter the value of µ (or H ≡ /µ ). Asecond property of the distribution that has an importantphysical impact is its “width” σ . In particular, we areinterested in how it scales with system size N . The lowercutoff of the distribution ξ lo = (2 A/µ ) µ M − µ N is related tothe system size and fixed by normalization. If 1 < µ < σ can be shown to be proportional to ξ lo , andso σ ∼ M − /µN (15)It will be fact the finite-size behavior of the lower cutoff in w ( ξ ), the noise produced by the far away plastic activity,what will dominate the scaling of interest. There is alsoan upper cutoff for the kick distribution, ξ up , but thatis related with the strongest, nearest plastic events, andindependent on the system size [21, 26].We have shown in the previous section that any finitestep random walk process of a variable x with absorbingboundaries, subject to such a random noise with 1 < µ < P ( x ) ∼ x θ for x (cid:38) σ (16) P ( x ) ∼ σ θ for x → θ = µ/
2, and σ is the “width” of the distribu-tion w ( ξ ), as previously defined. For instance, possiblefunctional forms for P ( x ) at small x are P ( x ) (cid:39) σ θ + x θ or P ( x ) (cid:39) ( σ + x ) θ . Furthermore, we have shown that N random walkers, coupled by a common global driftgenerate the same limiting form of P ( x ) as x → M N in Eq. 14 on the system size N , and use it to calculate thescaling of σ (Eq. 15) and thus the N -dependence of (cid:104) x min (cid:105) .Note that the approach of [21, 26] uses M N = N whichimplicitly considers that each of the N sites produces in-dependent kicks on the generic block i , perturbing x i .We would like to stress here that this is clearly not real-istic. Furthermore, in careful consideration, it goes itselfagainst the basic feature of yielding phenomena display-ing size-spanning avalanches and sub-extensive scaling ξξ w Z (N ) −(µ+1) ∼ξ ∼ξ −(µ+1) /M N ξξ (N ) (a) (b) Figure 6. (a) A schematic plot of the number density of kicks Z of a given intensity ξ observed in systems of two differentsizes N , and N > N . The two curves differ below the smallsize threshold ξ but are coincident in the heavy tail part,for large ξ . (b) The two curves in (a) normalized to becomethe probability distribution w ( ξ ). The normalizing factor isthe number M N of avalanches that occur in the two systemsunder the same increase of external strain. for the rate of plastic events. Using M N = N and µ = 1in Eq.14 implies somehow extensivity if kicks are sup-posed to be independent. Instead, we think on the totalnoise produced by one avalanche. Among the marginalkicks that a site receives (the ones that it almost fail tocatch because of working in a finite system N ), the dom-inant one is not the kick coming from a single site atthe maximum possible distance, but the largest possiblekick coming from such a distance. That is, a kick com-ing from the largest avalanche at the largest distance.If Eq. 14 represents the distribution of kicks generatedby individual avalanches in the system, the value of M N must be chosen in accordance with this interpretation.The dependence of M N on system size N can beworked out as follows. Consider two systems with dif-ferent sizes N and N > N , and suppose that we wantto compare the number of kicks of intensity ξ producedonto some reference site when a fixed (long) deformationstrain is applied to the system. The Eshelby interactingkernel decays in space as ∼ /r d , and this implies thatincreasing the system size from N to N > N does notproduce new large kicks [35], but instead increases thenumber of small ones, those generated at large distancesin the system with N sites. This means that if we plotthe density number of kicks observed at a given site as afunction of the kick magnitude, we would obtain a plot asthe one qualitatively depicted in Fig. 6(a). The portionof these curves following the 1 / | ξ | µ +1 law will be mostlyindistinguishable for the two system sizes. Now, in or-der to plot the probability distribution w ( ξ ), as shown inFig. 6(b), it is clear that we have to divide by the totalnumber of avalanches (kicks) that occurred in each case.This is why M N in Eq. 14 must be considered to be pro-portional to such a number. In other words, M N and theaverage size of avalanches in the system, noted S , mustbe related through M N ∼ N S − ∼ (cid:104) ∆Σ (cid:105) − (18)(which, together with Eq. 2 justifies our choice for M N (cid:39) /x min in the toy model of the previous section II A).Now, collecting the results of Eqs. 2, 15, and 18 wearrive at the important result[36] σ ∼ (cid:104) x min (cid:105) /µ . (19)Introducing this into Eq. 17 we get P ∼ σ θ ∼ (cid:104) x min (cid:105) / , (20)since, for 1 < µ < θ/µ = 1 / µ in such range.We are now only one step away from our general scal-ing results. As mentioned before, recent results in simula-tions of different EP models [15, 28] and also in MD sim-ulations [25] have shown that (i) a plateau exists for P ( x )at vanishing x , but also that (ii) (cid:104) x min (cid:105) shifts towards theplateau region of P ( x ) as the system size N is increased.This can now be analytically justified: From 16 and 17the crossover between the plateau and the power-law re-gion is expected at x cross (cid:39) σ . Combined with Eq. 19,this provides x cross ∼ (cid:104) x min (cid:105) /µ . For any µ >
1, this tellsthat (cid:104) x min (cid:105) becomes lower than x cross for large N . Inpractice, crossovers can be very broad, yet, in the limit N → ∞ the following relation holds (cid:104) x min (cid:105) P (cid:39) /N (21)Using Eqs. 20 and 21 we finally obtain the two importantpredictions: (cid:104) x min (cid:105) ∼ N − / (22)and P ∼ N − / . (23)Notice further that if we assume P ( x ) (cid:39) P + x θ , usingEq. 20: P ( (cid:104) x min (cid:105) ) (cid:39) (cid:104) x min (cid:105) / + (cid:104) x min (cid:105) θ . And, provided θ = µ/ > /
2, the second term becomes negligible overthe first when (cid:104) x min (cid:105) is small enough. We then could alsoexpect a good ansatz to be: P ( (cid:104) x min (cid:105) ) (cid:39) (cid:104) x min (cid:105) / . (24)Followed up from Eq. 20, the latter would interchange P by P ( (cid:104) x min (cid:105) ) in every subsequent expression. In thelimit N → ∞ both formulations are equivalent, since weexpect P ( (cid:104) x min (cid:105) ) to be part of the plateau and identical to P . Notice nevertheless that Eq.24 (and the ones derivedfrom it) may work well even before reaching that limit.The scaling provided by Eqs. 22 and 23 (or alterna-tively P ( (cid:104) x min (cid:105) ) ∼ N − / ) is quite generic, as it does notdepend on the actual value of µ neither on the dimen-sion of the problem. Even more, it is highly stimulating,since it agrees with the original observations of the (cid:104) x min (cid:105) scaling in MD simulations [7, 8] both in d = 2 and d = 3.Yet, there are assumptions implicitly made in their de-duction that can limit their validity. For instance, ourconstruction does not account for anisotropy effects on -6 -5 -4 -3 -2 -1 x -2 -1 P ( x ) -4 -2 xN (1/3)/0.75 -4 -2 xN -1 P ( x ) N / x L (=N ) (a)(b) (c) Figure 7. Distribution of local distances to threshold P ( x ) inthe quasistatic driven steady state of Picard’s 2D model. (a) The P ( x ) distributions. Different linear system sizes L = √ N are represented with different colors/symbols as declared inthe label. Pink crosses indicate the location of (cid:104) x min (cid:105) for thedifferent system sizes. (b) P ( x ) N / vs xN / testing thescalings of Eqs. 22 and 23. (c) P ( x ) N / vs xN (1 / / . topreserve the power-law regime ∼ x θ with θ = 0 .
75 observedin the main plot at intermediate values of x . the dimensions composing the system, which could af-fect the scaling of any observable with the global systemsize N . Such an effect appears clearly when consideringthree dimensional systems, as we discuss bellow. In addi-tion, Eqs. 22 and 23 do not apply in the case of a modelwith a (quenched) random kernel, that we describe inAppendix B, mainly due to the failure of the argumentabout the scaling of σ with N . In the next section we testthe predictions of Eqs.22 and 23 in elasto-plastic modelsin dimensions d = 2 and d = 3. IV. ELASTO-PLASTIC MODELS IN 2 AND 3DIMENSIONS
We now present results of quasistatic simulations ofspatially extended elasto-plastic models. We will limitourselves in particular to the Picard’s model [37]. Detailsabout model definition and simulation protocols can befound in the Appendix A, and data was produced withessentially the same codes used in [15].
Two-dimensional systems (2D)
We start with the d = 2 case. Fig. 7 shows the dis-tribution P ( x ) for different system sizes N = L × L . Wehave collected the values of x ≡ Σ th − Σ (see App. A for N -6 -4 -2 〈 x min 〉 [d=2]P [d=2]P( 〈 x min 〉 ) [d=2] ~N -1/3 ~N -2/3 Figure 8. Dependence of (cid:104) x min (cid:105) , P and P ( (cid:104) x min (cid:105) ) with systemsize N for 2D Picard’s model . parameters definitions) from every block in the systemfor several configurations in the steady state right afteran avalanche has finished and before loading the systemto the next avalanche. As discussed in previous sections, P ( x ) displays –also in this fully spatial model– an excessof probability at x = 0, evidencing the occurrence of anaturally emerging discrete step for the wandering of the x values. Already from the upper panel (Fig.7(a)) it isevident the settling of a system-size dependent plateauat x →
0. This plateau occurs systematically at smallervalues of x as L increases. The form of P ( x ) has morestructure than in the random-walk experiments of Sec. II.Now the crossover region between the power-law regimeand the plateau is broader, the power law range is shrunkdue to the natural existence of a global drift, and forsmall systems P ( x ) even displays an “S” shape beforecutting-off when x becomes order 1. Yet, we can identifyfor the largest system size a power-law regime spanningtwo orders of magnitude in x ([ ∼ . − , ∼ . − ]) inexcellent agreement with x θ with θ = 0 .
75 (the value ex-pected when µ = 3 / P ( x ) N / vs. xN / . The magenta crosses indicatethe position ( x min , P ( x min )) on each P ( x ) curve [38]. Thecoincidence of the horizontal coordinate of these pointsis the indication that Eq. 22 is very well satisfied. Ac-cording to Eq. 23 we also expect that the plateaus of allcurves in Fig. 7(b) level up. We see that they do butnot perfectly. Instead, note that the values of P ( x ) at x = (cid:104) x min (cid:105) (i.e., the vertical coordinate of the crosses)do become coincident in Fig. 7(b), fulfilling better thecombination of Eqs. 22 and 24 (cid:104) x min (cid:105) P ( (cid:104) x min (cid:105) ) ∼ /N. (25)0 -6 -5 -4 -3 -2 -1 x -1 P ( x ) -4 -2 xN (2/9)/0.37 -4 -2 xN (7/9) -1 P ( x ) N / (a)(b) (c) L (=N ) Figure 9. Distribution of local distances to threshold P ( x ) inthe quasistatic driven steady state of Picard’s 3D model. (a) The P ( x ) distributions. Different linear system sizes L = √ N are represented with different colors/symbols as declared inthe label. Pink crosses indicate the location of (cid:104) x min (cid:105) for thedifferent system sizes. (b) P ( x ) N / vs xN / testing thescalings of Eqs. 22 and 23. (c) P ( x ) N / vs xN (2 / / . topreserve the power-law regime ∼ x θ with θ = 0 .
37 observedin the main plot at intermediate values of x . Fig. 7(b) is built to display the combined scaling of (cid:104) x min (cid:105) and P (or P ( (cid:104) x min (cid:105) )). If instead we want to get a col-lapse of the power-law range of the P ( x ) distribution fordifferent system sizes, we must preserve the power-lawexponent in the transformation. This is done in Fig. 7(c)where we plot P ( x ) N / vs. xN (1 / / . , according tothe observed θ (cid:39) .
75. Following our generalized mean-field picture the value θ (cid:39) .
75 observed in the 2D elasto-plastic model corresponds to a mechanical noise with aHurst exponent H = µ − (cid:39) / µ = 2 θ (cid:39) / H (cid:39) /
3. Further-more, very recently compatibility with µ (cid:39) / (cid:104) x min (cid:105) , P ( (cid:104) x min (cid:105) ) and P (estimated from the curves in Fig. 7) as a function of N = L . Dashed straight lines are displays of the exactpower-laws N − / and N − / , not fits. We can see thatthe prediction of Eq. 22 work remarkably well and Eq. 25accompanies it perfectly. The original prediction for thescaling of P (Eq. 23) is also good (as could be seen inthe collapses of Figs. 7(b)-(c)), but we can also noticethat P is slowly merging with P ( (cid:104) x min (cid:105) ) as system sizeincreases, and it is indeed when N → ∞ when we expectthem to be equal and Eq. 23 to hold. N -6 -4 -2 〈 x min 〉 [d=3]P [d=3]P( 〈 x min 〉 ) [d=3] ~N -2/9 ~N -7/9 Figure 10. Dependence of (cid:104) x min (cid:105) , P and P ( (cid:104) x min (cid:105) ) withsystem size N for (cubic box) 3D Picard’s model. . Three-dimensional systems (3D)
Now, let us discuss the three-dimensional case. Con-trary to the 2D case, where the few interaction kernelsthat one can choose (corresponding to the different kindof volume-preserving applied deformations) are symmet-ric under the exchange of q x and q y , in 3D the manydifferent possibilities for choosing the elastic kernel allare non-symmetric respect to the permutations of q x , q y and q z . The precise symmetry of the six independentdeviatoric modes in 3D can be seen for example in [39].The results we present here correspond exclusively to thekernel shown in Eq. A7, where the way in which the z dimension enters differs from that of x and y . In Fig. 9we show data similar to that in Fig. 7 but for the d = 3case. We can fist observe in the raw data of Fig. 7(a)that the determination of the θ exponent is more am-biguous than in d = 2. At intermediate values of x , say ∼ (0 . − . θ (cid:39) . − .
37, as reported inprevious works [11, 14]. Yet, such a value for θ wouldimply µ = 2 θ (cid:39) . − . < H > x close to x = 0 to hold. Noticenevertheless that, for the largest system sizes, anotherpower-law regime at smaller x ∼ (10 − − − ) is insinu-ated. We will come back on this when discussing systemswith different aspect ratios, but let us advance that suchpower-law with a steeper slope would represent a moreconsistent value for θ in d = 3.In any case, let us now discuss scalings for the datain Fig. 9. In Fig. 9(b) we see that the N dependenceof (cid:104) x min (cid:105) an P follows a power-law behavior like the onepredicted by Eqs. 22 and 23 but with clearly different ex-ponents. Actually, the observed scaling is (cid:104) x min (cid:105) ∼ N − / Figure 11. Dependence of x min with system size N for the 3DPicard’s model with different aspect ratios. The sketch of theinset represents the 3D simulation box with its dimensions L x , L y and L z , and arrows indicating the strain deformationthat gives rise to the propagator that we use for d = 3 in thiswork (Eq.A7). and P ∼ N − / . Using these values we rescale the P ( x )data to obtain Fig. 9(b). Again, notice that as in thecase of d = 2 the collapse of the points ( (cid:104) x min (cid:105) , P ( (cid:104) x min (cid:105) ))(Eqs. 22 and 24) is better than the scaling of the plateaus,which are even hard to define. If we further consider thepower-law regime with an exponent θ (cid:39) .
37 we can do asin the d = 2 case and produce Fig. 9(c), for completeness.In Fig. 10 we show the values of (cid:104) x min (cid:105) , P ( (cid:104) x min (cid:105) )and P (estimated from Fig. 9) as a function of N = L for d = 3. Dashed straight lines simply display thepower laws ∼ N − / and ∼ N − / , they are not fits. Themeasured values shown in Fig. 10 follow these trends verywell. These values do not coincide with the predictions ofEqs. 22 and 23. We believe the main reason is that ourargumentation in the previous section implicitly assumedthat all spatial dimensions of the system participate onthe same footing. As we already stressed it, while the d = 2 Eshelby propagator (Eq. A2) is in fact symmetricagainst exchange of axis, this is not the case for the d = 3propagator (Eq.A7).We can provide a partial explanation for the valuesfound for the N dependence of (cid:104) x min (cid:105) and P ( (cid:104) x min (cid:105) ) (or P ) in 3D in the following way First of all, notice thatfor q z = 0 the three dimensional kernel (Eq. A7) reducesto the two dimensional one (Eq. A2). We will make theassumption that the non-trivial scaling of (cid:104) x min (cid:105) is stillgoverned by the finite-kick walk analysis that we did inSec. III, but in which the z coordinate has to be treatedas a ‘dumb’ independent dimension. This is, let’s thinkon the d = 3 case as a collection of several d = 2 systemsstacked in the z direction, and evolving in parallel. Ifwe take, L z systems of size L × L and choose after eachavalanche the minimum x among all of them, we would -6 -5 -4 -3 -2 -1 x -1 P ( x )
64 64 864 64 3264 64 19264 64 38416 16 6432 32 6464 64 64128 128 64256 256 64512 512 64 ~x ~x L x L y L z Figure 12. P ( x ) for Picard’s 3D model with different aspectratios. For the largest system size of each size grow scheme(varying L x = L y or varying L z ), the power-law regime visibleat the smallest values of x is marked by a straight line to guidethe eye. have a (cid:104) x min (cid:105) scaling as (cid:104) x min (cid:105) ∼ L − / L − z (26)and P ∼ L − / (27)(note that P turns out to be independent of L z ). When L z = L this leads to the scaling (cid:104) x min (cid:105) ∼ N − / and P ∼ N − / (with N = L ) that we observe in Fig. 8. Infact, simulations in systems with different L x = L y = L and L z show that Eqs. 26 and 27 are very well satisfied,as we will see in the following.Fig. 11 shows the scaling of (cid:104) x min (cid:105) for different cases.First, the N = L case is reproduced from Fig. 10 forcomparison. Then, we increase the system size while fix-ing L x = L y and varying only L z (the dimension perpen-dicular to the shear plane, that enters in a ‘different’ waythan the other two in the propagator A7). This yields ascaling (cid:104) x min (cid:105) ∼ N − controlled by (cid:104) x min (cid:105) ∼ L − z since thesystem size in the other two dimensions is fixed. Finally,we do inversely and we increase the system size by grow-ing L x = L y and keeping L z fixed. This yields a scaling (cid:104) x min (cid:105) ∼ N − / controlled by a scaling of (cid:104) x min (cid:105) ∼ L − / (eq. 26) for both L x and L y . Notice that when the systemsize is increased in this way (at a fix perpendicular direc-tion to the shear plane) we recover the scaling observedin the MD simulations of [7–9], that shows no exponentdifference between d = 2 and d = 3.In Fig. 12 we take a look to the P ( x ) distributions inthese asymmetric boxes for different aspect ratios. Onone hand, we have fixed L x = L y = 64 and vary L z between 8 and 384. On the other hand, we have fixed L z = 64 and vary L x = L y between 16 and 512. No-tice first that, when L z is the only changing dimension,2the plateau level actually increases, with a small positivepower, and it seems to saturate for large sizes around P ∼ .
12. So, the strong (cid:104) x min (cid:105) scaling decreasing as1 /N , is accompanied by a barely changing P with N ,as we could have expected from Eq. 20. These curvesfor P ( x ) have the particularity that they only show apower-law regime at ‘large’ values of x , and they cor-respond to an ‘abnormally small’ value of θ , coincidentwith the many times reported [11, 14] but never trulyjustified θ (cid:39) . − .
37 in 3D. This θ value would pointto µ <
1, beyond the assumptions used for the derivationof our scaling arguments.Now, if we analyze the curves when varying L x = L y at fix L z things change dramatically. First, the plateaulevel is ‘well behaved” decreasing as N increases. In fact,a reasonable P ∼ N − . accompanies the scaling of x min shown in Fig.11 for this case. Secondly, the larger systemsizes clearly display a different power-law at intermediate x values, with P ( x ) ∼ x . in such range. As it standscloser to the boundary x = 0, that will be the power-lawdominating the system’s dynamics close to the transition(e.g., the value of the flowcurve exponent β [15, 16]). Infact, in d = 3, H = 1 / (2 θ ) (cid:39) . − .
77 is expected [16],consistent with θ (cid:39) . − .
66. Moreover, θ (cid:39) . µ = 2 θ (cid:39) .
33 in d = 3, which bringsthe problem back into the range of validity of our generalassumptions for the derivation of the scalings (22 and 23).It needs to be stressed that the occurrence of these twoclearly different finite size scalings in 3D – (i) growing thesystem in the direction perpendicular to the shear planeor (ii) growing the system in the directions of the shearplane – remains as an open issue (see discussion below). V. SUMMARY AND DISCUSSION
In this paper we have considered the problem of thestrain load ∆ γ needed to trigger consecutive avalanchesin the steady state of quasistatically deformed amorphoussolids. In particular, we studied the finite-size scaling ofits mean value (cid:104) ∆ γ (cid:105) . The values of ∆ γ are intimatelyrelated to the distribution P ( x ) of local distances toinstability x ; (cid:104) ∆ γ (cid:105) is simply proportional to the aver-age value of the minimum x across the system, namely (cid:104) ∆ γ (cid:105) ∼ (cid:104) x min (cid:105) . We have built a theoretical argumentstarting by simple random walks of x with an absorbingboundary to show how the effect of a discrete step in-duces a finite value of P ( x ) at the boundary. Then westood on an alternative mean-field modeling approach forthe yielding phenomena [15, 17], considering as the physi-cally relevant case the one in which the mechanical noiseis generated by extended and collective plastic events,leading to a fat tail noise distribution with 1 < µ < P ( x ) is expected to acquire a finite value as x → P ( x →
0) = P (cid:54) = 0. More importantly, the dis-creetness in the mechanical kicks is not the trivial ∼ /N finite-system effect, but one that has to do also with the mean avalanche size in the system (e.g., see Eq. 18). Thisholds for any 1 < µ < P vs N scaling in that case. The scenario is confirmed by ex-tensive numerical simulations of a classical elastoplasticmodel in 2 and 3 dimensions.Even though the value of P decreases to zero as N → ∞ , and therefore it could be naively considereda finite-size effect, its behavior with system size hap-pens to be precisely what governs the scaling of (cid:104) x min (cid:105) ,and thus of (cid:104) ∆ γ (cid:105) , our quantity of interest. Our theo-retical analysis is able to justify a universal dependence (cid:104) ∆ γ (cid:105) ∼ N − α , with α = 2 /
3, independent of spatial di-mension and system parameters, as is actually found inMD simulations [7, 8]. Moreover, we have no need to as-sume a particular shape for the energy barriers [8] in do-ing so. It is worth mentioning nevertheless, that, as mostof the numerical literature on the field, our constructionassumes so far an athermal system. In this case the dy-namics is dominated by the minimal value of distance toinstability, x min , at every loading step. A finite tempera-ture in a thermodynamic system ( N → ∞ ) may blur this(otherwise strictly) extremal dynamics. It might be aninteresting problem for future works to analyze how ourpredictions are impacted by a finite temperature.In the numerical results presented here for EP modelsin d = 2 the value α = 2 / d = 3systems display a different value α (cid:39) /
9. We have iden-tified a possible reason for this discrepancy in d = 3 inan unforeseen (cid:104) x min (cid:105) scaling dependence with the linearsize of the sample along different directions relative tothe externally applied shear. In contrast with the d = 2case, the interaction kernel in d = 3 is not symmetricin all coordinates. Growing the system in the directionperpendicular to the shear plane has an effect markedlydifferent on (cid:104) x min (cid:105) , than growing it in the other direc-tions. By studying d = 3 systems of different aspect ratiowe addressed these multiple scalings, showing that in thecase in which the dimension perpendicular to the shearplane is kept fixed, α (cid:39) / d = 3elasto-plastic systems. The different scalings of (cid:104) x min (cid:105) and P with the different linear dimensions of a 3D sys-tem can be rationalized by a gedanken problem in whichthe 3D system behaves as a collection of independent 2Dsystems. However, there is no basis to expect that this isactually the way in which a 3D system behaves and weknow that the dynamics of interactions is more complexthan that. The kernel asymmetry might be a weaknessof the EP simplification and non-physically dominant atthe end of the day. This is what one could interpretfrom the fact that classical MD results [7–9] maintainthe (cid:104) x min (cid:105) ∼ N − / scaling. Moreover, a recent proposalof ‘augmented’ elastoplastic models tries to incorporate(among other things) the fact that shear strain in any direction due to a rearrangement can trigger the nextrearrangement equally well. Successive rearrangementsobserved in MD are “isotropically distributed” and notconcentrated in the strain direction prescribed by the im-3posed deformation [40]. If the interaction kernel is sym-metrized somehow our predictions Eqs. 22 and 24, turnto be valid in the elastoplastic 3D case as well. We havechecked this so far for synthetic, non-physical, kernelsonly (data not shown).In any case, a definite value of α implies additionalpredictions on other critical exponents of the yieldingtransition. For instance, the avalanche distribution ex-ponent τ and the fractal dimension of avalanches d f arelinked to α through [11, 15] dα = d − d f (2 − τ ) . (28)(note that this relation is usually written using θ insteadof α , by applying the extra assumption α = 1 / (1 + θ ),that we consider not justified in the steady state). Aunique value α = 2 / d = 2 implies d f (2 − τ ) = 2 / d f (cid:39) τ (cid:39) .
33. For d = 3, we muststill understand which is the value of α that we shouldexpect, but d f and τ could also suffer from an asymmetryeffect if Eq. 28 is expected to hold.Finally, all this picture should be compatible withknown results for the as-quenched state; with rigorouspower-laws for P ( x min ) and P ( x ) at small arguments. Webelieve that the effective mechanical noise governing thedistribution P ( x ) and its properties, like the one thatdefines the finite-size scaling of (cid:104) x min (cid:105) , must display sys-tematic biases in the non-universal transient. While wewill not venture to link transient values of µ (or H ) with θ in such a regime (which, furthermore, is only measurableon a given system size for certain ranges of initial anneal-ing), our guess is that avalanches progressively build upand their geometry -encoded in d f [11, 20]- varies withstrain, therefore modifying the effective noise, until itreaches a steady distribution governed by 1 < µ < Conclusion
In conclusion, we have provided a novel interpretationof the finite size scaling of (cid:104) ∆ γ (cid:105) in the steady state ofamorphous systems under deformation. This interpreta-tion seems to conciliate MD simulation results and EPconstructions, otherwise in contradiction in this limit.While the hypothesis of a marginal stability behavior,rooted in the celebrated P ( x ) ∼ x θ pseudo-gap, has beenproved to hold in the as-quenched isotropic state of modelglasses [7, 8] and still renders important outcomes in thetransient [20, 24], it does not seem to apply ‘as-is’ tothe steady state case. There, at least, the system dy-namics is correlated at the level of avalanches and thisnaturally produces a finite value of P ( x ) as x →
0, whenobserving the P ( x ) distribution in the quasistatic limit,justified on a discrete step for the effective dynamics ofthe x values. This behavior of P ( x ) does not invalidatethe essence of the yielding transition, anchored in thesub-extensive scaling of (cid:104) ∆ γ (cid:105) − ; since the level of such asymptotic plateau at small x is itself dependent on N and is shown to govern the behavior of (cid:104) ∆ γ (cid:105) ∼ N − α ,independently on θ .Some questions remain open, and we hope they willmotivate further endeavors on the subject. But we be-lieve that this is a first step in shedding light on a prob-able misconstruction in the field, based in a wrong ex-trapolation of arguments valid in the early deformationregime to the steady state case. ACKNOWLEDGMENTS
We are indebted to D. Vandembroucq and C. Mal-oney for illuminating discussions on an early version ofthis manuscript. We sincerely thank the critical feedbackprovided by M. Wyart on a first draft of this work. Wealso acknowledge exchanges with J.-L. Barrat, E. Lerner,J. Rottler and B. Tyukodi. EEF acknowledges supportfrom PICT-2017-1202.
Appendix A: Elastoplastic model and simulationprotocol
EP models are intended to describe amorphous ma-terials at a coarse-grained-level, laying in between theparticle-based simulations and the continuum-level de-scription [5]. In short, the amorphous solid is representedby a coarse-grained scalar stress field Σ( r , t ), at spatialposition r and time t under an externally applied shearstrain. Space is discretized in blocks (e.g., square lattice).At a given time, each block can be “inactive” or “active”(i.e., yielding). This state is defined by the value of anadditional variable: n ( r , t ) = 0 (inactive), or n ( r , t ) = 1(active). An over-damped dynamics is imposed for thestress on each block, following some basic rules: (i) Thestress loads locally in an elastic manner while the blockis inactive. (ii) When the local stress overcomes a localyield stress, a plastic event occurs with a given probabil-ity, and the block becomes “active” ( n ( r ) is set to one).Upon activation, dissipation occurs locally, and this is ex-pressed as a progressive drop of the local stress, togetherwith a redistribution of the stresses in the rest of the sys-tem in the form of a long-range elastic perturbation. Ablock ceases to be active when a prescribed criterion ismet. The auxiliary binary field n ( r , t ) shows up in theequation of motion for the local stress Σ( r , t ), defininga dynamics that is typically non-Markovian. While thestructure of the equation of motion for the local stressesis almost unique in the literature, both its parametersand the rules governing the transitions of n ( r ) (0 (cid:10) d -dimensional scalar fieldΣ( r , t ), with tipically d = 2 or 3, and r discretized ona square/cubic lattice and each block Σ i subject to the4following evolution in real space ∂ Σ i ( t ) ∂t = µ ˙ γ ext + (cid:88) j G ij n j ( t ) Σ j ( t ) τ ; (A1)where ˙ γ ext is the externally applied strain rate, and thekernel G ij is the Eshelby stress propagator [41].It is sometimes convenient to explicitly separate the i = j term in the previous sum, as ∂ Σ i ( t ) ∂t = µ ˙ γ ext − g n i ( t ) Σ i ( t ) τ + (cid:88) j (cid:54) = i G ij n j ( t ) Σ j ( t ) τ ;(A2)where g ≡ − G ii > G is G ( r , r (cid:48) ) ≡ G ( r, ϕ ) ∼ πr cos(4 ϕ ) in polar coordinates,where ϕ ≡ arccos(( r − r (cid:48) ) · r ˙ γ ( ext ) ) and r ≡ | r − r (cid:48) | . Forour simulations we obtain G ij from the values of the prop-agator in Fourier space G q , defined as G q = − q x q y ( q x + q y ) . (A3)for q (cid:54) = and G q = = − κ (A4)with κ a numerical constant (see below). Note that inour square numerical mesh of size L × L , q x , q y must beunderstood as q x,y ≡ − (cid:16) πm x,y L (cid:17) (A5)with m x,y = 0 , ..., L − µ = 1 defines thestress unit, and the mechanical relaxation time τ = 1,the time unit of the problem. The last term of (A2)constitutes a mechanical noise acting on Σ i due to theinstantaneous integrated plastic activity over all otherblocks ( j (cid:54) = i ) in the system.The picture is completed by a dynamical law for thelocal state variable n i = { , } . We define hereafter therule corresponding to the Picard’s model [37] that we use: n i : (cid:40) → τ − on if Σ i > Σ th i ← τ − off (A6)where τ on and τ off are parameters and P (Σ th i ) = δ (Σ th i − d = 3, the Eshelby kernel for one scalar componentof the deviatoric strain in Fourier space can be writtenas G q = − q x q y + q z ( q x + q y + q z )( q x + q y + q z ) (A7)and the dimensional extension of the dynamics isstraightforward.
1. Quasistatic protocol
For the analysis of avalanche statistics, it is convenientto have a protocol that allows for the triggering and un-perturbed evolution (no driving) of avalanches until theystop, guaranteed by a degree of stress non-conservation κ > κ = 1, as in previous strain-controlledEP models implementations [14, 42, 43], unless otherwisespecified). This is the quasi-static protocol describedhere.Starting from any stable configuration, i.e., no siteis active and no site stress is above its local threshold( n i = 0 and Σ i < Σ th i for all sites), the next avalancheof plastic activity is triggered by globally increasing thestress by the minimum amount necessary for a site toreach its local threshold. That site (the weakest) is ac-tivated at threshold with no stochastic delays; it per-turbs the stress values of other sites and the rest of theavalanche evolves without any external drive followingthe dynamics prescribed by Eq. (A2) (and the corre-sponding activation rule) with ˙ γ = 0. The avalanchestops once there are no more active sites and all stressesare below their corresponding thresholds again. At thispoint the loading process is repeated. For each simula-tion run, data is collected only in the steady-state. Appendix B: Model with a quenched random kernel
In this Section we analyze the properties of a modelwith a different form of the interaction kernel. Instead ofusing the appropriate interaction to describe the proper-ties of yielding, namely the Eshelby kernel presented inEq. A3, we consider a model in which the G q kernel takesrandom values. In concrete, we use G q = − RND( q ) . (B1)where RND( q ) stands for an independent random num-ber chosen from a flat distribution between 0 and 1 foreach value of q . Note that this is a “quenched” randomkernel, since the form of G q is chosen once and for all atthe beginning of the simulation.Although this is probably not a realistic model to de-scribe any physical situation, there are a few reasonsthat make the study of this model interesting. The firstone concerns its relation with another version of a “ran-dom” yielding model, namely the H´ebraud-Lequeux (HL)model [44, 45]. In its essence, the HL model for a systemwith N sites considers that every time a single site per-forms a plastic re-accommodation, it produces a randomkick of finite variance σ (with σ ∼ N − / ) on every othersite. Note however that in this case the values of the ran-dom kicks are refreshed at every plastic event[46]. Fromits very definition the mechanical noise in the HL modelis a standard random walk, corresponding to a value of µ = 2. In the quenched random case we are examin-ing, we must first understand what are the properties5 -7 -6 -5 -4 -3 -2 -1 x -3 -2 -1 P(x) N -4 -3 -2 -1 P < x min > N Figure 13. (a) The form of P ( x ) for system with differentnumber of sites N in the quenched random kernel case. Thedashed line displays the expected behavior P ( x ) = C + C x on the N = 256 data. (b) The scaling of P and (cid:104) x min (cid:105) with N . Symbols are the result of simulations. Straight linesindicate the expected ∼ N − / dependence. of the uncorrelated mechanical noise felt by a particulartarget site. The quenched random kernel G q generatesvalues G r that are mostly uncorrelated spatially, and dis- tributed with a finite variance σ . This is enough to guar-antee that we will find a value µ = 2 (and therefore [21] θ = 1) as in the HL model. In addition, the dependenceof σ on the number of sites N in the system is σ ∼ / √ N ,as in the HL model. Then we can write down the scalingof P with N from (the limit of validity of) Eq. 12, whichis independent of the details of the kernel, as P ∼ N − / (B2)and also (cid:104) x min (cid:105) ∼ N − / (B3)thus finding in the present case a different scaling thatthe one given by Eqs. 22 and 23. The arguments that ledto Eqs. 22 and 23 fail here because the scaling of σ with N obtained in Eq. 19 was based in the conservation ofthe number of large kicks when system size is increased(see Fig. 6), something that does not occur here becauseof the assumed non-decaying nature of the interactions.We performed simulations with a quenched randomkernel and evaluated the distribution P ( x ), and the valueof (cid:104) x min (cid:105) . The simulations shown here were done in twospatial dimensions, but we verified that exactly the sameresults are obtained in three dimensions if the numberof sites in the system is maintained. This is of courserelated to the fact that in a randomly interacting modeldimensionality plays no relevant role.Figure 13 shows the results of simulations with therandom kernel. The upper panel shows the form of P ( x )for different values of N and we can see that the value θ = 1 for P ( x ) ∼ x θ in an intermediate range of x iswell established as N increases. The lower panel showsthe scaling of P and (cid:104) x min (cid:105) with N , where we observeclearly the expected ∼ N − / dependence. [1] J. M. Carlson, J. S. Langer, and B. E. Shaw, Reviews ofModern Physics , 657 (1994).[2] D. S. Fisher, Phys. Rep. , 113 (1998).[3] J. Ferr´e, P. J. Metaxas, A. Mougin, J.-P. Jamet, J. Gor-chon, and V. Jeudy, Comptes Rendus Physique , 651(2013), disordered systems / Syst`emes d´esordonn´es.[4] P. G. de Gennes, Rev. Mod. Phys. , 827 (1985).[5] A. Nicolas, E. E. Ferrero, K. Martens, and J.-L. Barrat,Rev. Mod. Phys. , 045006 (2018).[6] C. Maloney and A. Lemaˆıtre, Physical Review Letters , 016001 (2004).[7] E. Lerner and I. Procaccia, Phys. Rev. E , 066109(2009).[8] S. Karmakar, E. Lerner, and I. Procaccia, Phys. Rev. E , 055103 (2010).[9] S. Karmakar, E. Lerner, I. Procaccia, and J. Zylberg,Phys. Rev. E , 031301 (2010).[10] J. Lin, A. Saade, E. Lerner, A. Rosso, and M. Wyart,Europhysics Letters (EPL) , 26003 (2014).[11] J. Lin, E. Lerner, A. Rosso, and M. Wyart, Proceedingsof the National Academy of Sciences , 14382 (2014). [12] Z. Budrikis, D. F. Castellanos, S. Sandfeld, M. Zaiser,and S. Zapperi, Nat. Comm. , 15928 (2017).[13] B. Tyukodi, S. Patinet, S. Roux, and D. Vandembroucq,Phys. Rev. E , 063005 (2016).[14] C. Liu, E. E. Ferrero, F. Puosi, J.-L. Barrat, andK. Martens, Phys. Rev. Lett. , 065501 (2016).[15] E. E. Ferrero and E. A. Jagla, Soft Matter , 9041(2019).[16] E. E. Ferrero and E. A. Jagla, Phys. Rev. Lett. ,218002 (2019).[17] I. Fern´andez Aguirre and E. A. Jagla, Phys. Rev. E ,013002 (2018).[18] E. A. Jagla, Journal of Statistical Mechanics: Theoryand Experiment , 013401 (2018).[19] K. M. Salerno and M. O. Robbins, Phys. Rev. E ,062206 (2013).[20] J. Lin, T. Gueudr´e, A. Rosso, and M. Wyart, Physicalreview letters , 168001 (2015).[21] J. Lin and M. Wyart, Physical Review X , 011005(2016).[22] H. G. E. Hentschel, P. K. Jaiswal, I. Procaccia, and S. Sastry, Phys. Rev. E , 062302 (2015).[23] W. Ji, M. Popovi´c, T. W. J. de Geus, E. Lerner, andM. Wyart, Phys. Rev. E , 023003 (2019).[24] B. Shang, P. Guan, and J.-L. Barrat, Proceedings of theNational Academy of Sciences , 86 (2020).[25] C. Ruscher and J. Rottler, Soft Matter , 8940 (2020).[26] J. Lin and M. Wyart, Phys. Rev. E , 012603 (2018).[27] As a matter of fact, different values for the exponent θ are presented in Ref. [11] when either fitted from the P ( x )distribution or computed from the ‘extremal dynamics’(the scaling of (cid:104) x min (cid:105) ) through Eq.7; “a difference pre-sumably resulting from corrections to scaling” accordingto the authors.[28] B. Tyukodi, D. Vandembroucq, and C. E. Maloney,Phys. Rev. E , 043003 (2019).[29] K. Karimi, E. E. Ferrero, and J.-L. Barrat, Phys. Rev.E , 013003 (2017).[30] A. Lemaˆıtre and C. Caroli, arXiv preprintarXiv:0705.3122 (2007).[31] J. Parley, S. Fielding, and P. Sollich, arXiv preprintarXiv:2010.02593 (2020).[32] In the simulations presented here the reinjection is maderandomly and uniformly in the full interval (0 , x min todefine the distribution w ( ξ ), one could use the self-tunedmean value (cid:104) x min (cid:105) and the conclusions are identical.[35] In fact the largest kicks are produced by neighboravalanches, the coarse-grained lattice description imposesthe upper cutoff of the kick distribution, the minimal dis- tance.[36] See Appendix B for the discussion of a case in which theassumptions made to derive this result do not apply, andthen Eq. 19 does not hold.[37] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, Phys-ical Review E , 010501 (2005).[38] x min is independently computed for each system size asthe arithmetic average of the minimum x values (in the L × L system) for each after-avalanche configuration inthe steady state.[39] E. A. Jagla, Phys. Rev. E , 043004 (2020).[40] G. Zhang, S. Ridout, and A. J. Liu, arXiv preprintarXiv:2009.11414 (2020).[41] G. Picard, A. Ajdari, F. Lequeux, and L. Bocquet, TheEuropean physical journal. E, Soft matter , 371 (2004).[42] K. Martens, L. Bocquet, and J.-L. Barrat, Soft Matter , 4197 (2012).[43] A. Nicolas, K. Martens, and J.-L. Barrat, EPL (Euro-physics Letters) , 44003 (2014).[44] P. H´ebraud and F. Lequeux, Physical Review Letters ,2934 (1998).[45] E. Agoritsas, E. Bertin, K. Martens, and J.-L. Barrat,Eur. Phys. J. E , 71 (2015).[46] It has to be emphasized that when using a quenched ker-nel as in this case, there is a stability condition expressedin the fact that G q has to be non-positive, otherwise wewould obtain exponentially growing modes. This is whywe define the random kernel in q space. If we define arandom kernel in real space instead, the negativity of qq