On the critical exponents of the yielding transition of amorphous solids
I. Fernández Aguirre E. A. Jagla



 
Comisión Nacional de Energía Atómica, Instituto Balseiro (UNCu), and CONICET

Centro Atómico Bariloche, (8400) Bariloche, Argentina





 

We investigate numerically the yielding transition of a two dimensional model amorphous solid under external shear. 
We use a scalar model in terms of values of the total local strain, that we derive from the full (tensorial) description of the elastic interactions in the system, in which
plastic deformations are accounted for by introducing a stochastic “plastic disorder" potential.
This scalar model is seen to be equivalent to a collection of Prandtl-Tomlinson particles,  which are coupled through an Eshelby quadrupolar kernel. Numerical simulations of this scalar model reveal that the strain rate vs stress curve, close to the critical stress, is of the form γ̇∼ (σ-σ_c)^β. Remarkably, we find that the value of β depends on details of the microscopic plastic potential used, confirming and giving additional support to results previously obtained with the full tensorial model. 
To rationalize this result,  we argue that the Eshelby interaction in the scalar model can be treated to a good approximation in a sort of “dynamical" mean field, which corresponds to a Prandtl-Tomlinson particle that is driven by the applied strain rate in the presence of a stochastic noise generated by all other particles. The dynamics of this Prandtl-Tomlinson particle displays different values of the β exponent depending on the analytical properties of the microscopic potential, thus giving support to the results of the numerical simulations.
Moreover, we find that other critical exponents that depend on details of the dynamics show also a dependence with the form of the disorder,  while static exponents are independent of the details of the disorder.
Finally, we show how our scalar model relates to other elastoplastic models and to the widely used mean field version known as the Hébraud-Lequeux model.




On the critical exponents of the yielding transition of amorphous solids
    I. Fernández Aguirre and E. A. Jagla
    
========================================================================






§ INTRODUCTION



Amorphous solid materials are ubiquitous in every day life, and of great practical importance in many industrial processes <cit.>. They consist of a collection of elementary units that accommodate in space without a well defined ordering (contrary to what happens with crystals). The nature of the elementary units that form the material may span a wide range, roughly from ∼ 0.1 nm to ∼ 1  m<cit.>. The fact that in many cases these units are not microscopic leads to the fact that thermal fluctuations may be negligible in explaining the mechanical properties of these materials, which are then termed “athermal".

In recent years, there has been an increasing effort aimed at elucidating the mechanical properties of amorphous solids. 
One main piece of the phenomenology of amorphous materials is the existence of a yielding transition: In the absence of appreciable thermal activation effect, the material remains rigid if the applied stress is below some threshold, and it flows continuously if this threshold is exceeded. The properties of the material around this critical stress σ_c, or yield point, has attracted much attention. It is experimentally found  <cit.>  that the strain rate in the system γ̇ as a function of stress excess σ-σ_c follows in many cases a power law behavior of the form γ̇∼ (σ-σ_c)^β.  The flow exponent β is an important parameter characterizing the problem.
Other important critical exponents emerge when one considers the nature of the dynamics close to the transition. This dynamics proceeds through abrupt rearrangements in the system <cit.>, which share many features with the avalanches observed in the related model of depinning of an elastic interface<cit.>. This allowed to use in the analysis of the yielding transition the previously developed tools used in the depinning problem.<cit.>
In particular one can define for yielding additional critical exponents associated to the statistics of avalanches close to the transition.

Since the universal aspects of the depinning transitions are well known, the issue of the universality of the yielding transition has attracted much interest.  In particular, are the values of the critical exponents independent of details of the model and only dependent on some very general characteristics as dimensionality, for instance? And if this is not the case, what are the system features that determine the differences? 


In this work we argue that there are differences in the values of some critical exponents in the yielding problem, related to the form of the plastic yielding potential that is used to model the plastic rearrangements in the system. This is an interesting finding since it does not occur in the (short range) depinning problem. 
We find evidence that this result is related to the long range nature of the elastic interactions in the yielding problem, which leads to a sort of effective “dynamical mean field” description. In fact, in mean field depinning the same dependence of exponent β on the form of the pinning potential is well known.<cit.>

The paper is organized as follows. In section II we present the model, which is a reduction to a scalar problem of a tensorial model of the yielding transition that has been presented previously <cit.>. In section III the main numerical results are presented, showing the dependences of some critical exponents on the form of the plastic disorder potential. Section IV contains the arguments leading to a “mean field like" description of the problem and then to the justification of the different values of the critical exponents found numerically. In section V we discuss to what extent the present model is comparable to the elastoplastic models
discussed in the literature. Finally, in Section VI we summarize and conclude.



§ MODEL


We motivate here the model in an heuristic way, emphasizing the physical ingredients it incorporates. 
In Ref. <cit.> there is a derivation of the model from a full tensorial description of the elasticity of the material. In addition, in Appendix I, we present an alternative view in which a very similar model is deduced assuming the deformation field of the material is strictly one-dimensional under a single shear imposed deformation. 

The mesoscopic model we present describes the evolution of the system under a given external shear deformation of uniform symmetry. The goal is to predict the evolution of the corresponding local deformation e(r) compatible with the externally applied load.

The dynamics to be used is an over-damped dynamics in which the rate of change of e(r) is equalled to an effective force acting at r. There are two main parts of this force. One is a local term encoding the internal dynamics of the element at r. This part is derived from a potential function V_r(e). The form of V_r(e) takes into account both the local elasticity of the material and also the possibility of different locally stable configurations: V_r(e) has minima at a sequence of e values, corresponding to equilibrium configurations. The transition between consecutive minima correspond to plastic events in the system. Around each minima V_r(e) behaves quadratically, reflecting the elasticity at the actual configuration. 
In addition, there is a term in the evolution equation that reflects the elastic interaction between elements at different spatial positions. This term is written in terms of a kernel G(r-r') which is usually referred to as the Eshelby kernel. The model reads (using a discrete spatial representation):


    ηė_i=-dV_i(e_i)/de_i+∑_j G_ije_j+σ

where σ is the applied stress. From now on we will set the viscous damping coefficient η to η=1.
G_ij depends only on the distance between i and j, and is more compactly described by its Fourier transform G_ q:

    G_ q=-2μ B(q_x^2-q_y^2)^2/μ q^4+2B q_x^2q_y^2

and G_ q=0=0. B and μ are the bulk and shear modulus of the material.
These two equations define the model completely.
Since ∑_i G_ij∼ G_ q=0=0, spatially averaging Eq. (<ref>) we obtain

    ė=-dV_i(e)/de+σ,

that determines the instantaneous value of the deformation rate γ̇≡ė.
Alternatively, in an implementation that fixes the value of the deformation rate γ̇, Eq. (<ref>) defines the value of the instantaneous stress as

    σ=γ̇+dV_i(e)/de


Note that in the present formalism there is a single quantity e_i for each site describing the state of the system, and the separation 
between elastic strain and plastic strain usually done in elasto-plastic models is not made. We will come back to the relation with other elasto-plastic models later on.


The V(e) are stochastic potentials chosen to be uncorrelated among different spatial positions. The values of e at which V(e)
has local minima correspond to locally stable configurations of the system. The form of V(e) is quadratic around these minima to model an elastic material. To fully define the form of V(e) we must specify how the wells corresponding to different minima are connected.
We consider two qualitatively different forms of the V(e) potentials (see Fig. <ref>). 
In the first case the wells are connected sharply, at points in which dV(e)/de has jumps. In the second case the connection is made smoothly. 

In concrete, to define a potential V(e), the e axis is divided in intervals [a_n,a_n+1]  (n integer), in such a way that 
Δ_n≡a_n+1-a_n is stochastically chosen from a flat distribution between Δ_min=2 and Δ_max=4  (we have checked that the use of an exponential distribution does not affect the results). 

The intervals are centered at a_n≡(a_n+1+a_n)/2. At each interval n, V(n) is defined as 

    V_n(e)=1/2[(e-a_n)^2-Δ_n^2]

for the case in Fig. <ref>(a), and as

    V_n(e)=-5(Δ_n/2π)^2[1+cos(2π(e-a_n)/Δ_n)]

for the case in Fig. <ref>(b). The first case will be referred to as the “parabolic" potential, and the second case as the “smooth" potential.
The qualitative main difference between the two cases concerns the behavior at the transition points between different wells.
These are the point of maximum force, which are indicated as green dots in Fig. <ref>. 
In the parabolic case these points coincide with the potential maxima, where there is a discontinuity in the force.
In the smooth potential case, the maximum force occurs at points where the curvature of the potential changes sign continuously.







§ RESULTS OF NUMERICAL SIMULATIONS


In this Section we present results of simulations of the model, to elucidate the effect of the form of the potential on the critical exponents of the transition and the avalanche statistics. We set units such that B=1, and work in the case μ=B.

We focus first on the value of the flow exponent β. The value of β can be measured straightforwardly by driving the system at a constant strain rate, and measuring the stress as it is defined in Eq. (<ref>). As a result, we obtain the flow curves shown in the Figure <ref>. This graph displays clearly the existence of a critical stress and a monotonic growth for larger stress. The logarithmic plot in <ref>(b) clearly indicate that the values of β are dependent on the form of the potential. We obtain β_p≃1.51 for parabolic potentials and  β_s≃2.00 for smooth potentials.


  



Motivated by this difference between the two kinds of potentials, we moved to study the exponents characterizing the avalanche dynamics. In order to calculate these quantities and to see in particular if they depend on the kind of potential used, we ran quasistatic simulations in the following way (see Figure <ref>). In a simulation with a small γ̇, the maximum value of de/dt across the system is calculated: V_max≡max_i (de_i/dt) . This quantity stays lower than a well-chosen threshold as long the system is stable. However, when an avalanche is developing this quantity becomes order 1. When the avalanche finishes V_max becomes very small again. In this way we can identify individual avalanches in the system.
It is important to point out that to obtain more precise results we stop the driving while an avalanche is taking place. This avoids spurious additional avalanche triggering by the driving.

As indicated in Fig. <ref>, avalanche size S is proportional to the global stress drop that the avalanche causes, and its duration T is measured as the time between the first jump of any site from one minimum to another and the last one. 

Additionally, we also monitor the strain increases Δγ that have to be applied after one avalanche to trigger a second one.




Results for the avalanche size distribution (Figure <ref>) display a power law P(S)∼ S^-τ, with τ≃ 1.40 for parabolic potentials and τ≃ 1.38 for smooth potentials. The power laws are cut off at large avalanche size by the system size. This cut off defines the fractal dimension d_f of the avalanches which describes how the maximum size of avalanches grows with the linear size of the system as S_max∼ L^d_f. In order to determine S_max most reliable from the simulation we use a relation of S_max with the average size of S and S^2 that reads<cit.> S_max∼S^2/(2 S).


We obtain from the simulations that d_f≃ 1.01 for parabolic potentials and d_f≃ 1.01 for smooth potentials. Taking into account the numerical uncertainties, we conclude that both τ and d_f are independent on the potentials being of the smooth or parabolic type.


  

In order to calculate the dynamical exponent z, we first plot the relation between duration and size of avalanches.
This is done in Figure <ref>. The data show a wide dispersion, but averaging over avalanche size windows of logarithmic width, a well defined power law T ∼ S^p is obtained. The values of p that are obtained from fitting are 
p_p≃0.53 for smooth potentials and p_s ≃ 0.41 for the parabolic potential. From this values and the definition of the dynamical exponent z as  z = p d_f <cit.>, it is obtained that
 z_p ≃ 0.53 for smooth potentials and z_s ≃ 0.42 for parabolic potentials. We consider this difference to be significant, and our conclusion is that the dynamical exponent z is different for parabolic and smooth potentials.

  



The last exponent that was calculated is the θ exponent, measuring the distribution of distance-to-instability at different position of the sample. If x is the additional stress that has to be added to a given site to become unstable and jump to the next potential well, then θ is defined through the probability distribution of x for x close to zero, as P(x) ∼ x^θ. It is not straightforward (particularly in the smooth potential case) to calculate θ from a given equilibrium configuration in the system. However, the following trick can be used <cit.>: θ can be calculated by following the average strain increase Δγ that has to be applied in order to activate consecutive avalanches. The results is that Δγ∼ L^-d/1+θ. By calculating Δγ for different values of L, θ can be determined. The corresponding graph is presented in Figure <ref> and we obtain that θ≃ 0.44 for parabolic potentials and θ≃ 0.47 for smooth potentials. The two values coincide within the numerical precision.



The conclusion from the numerical simulations is that the “dynamical" exponents β and z (those  that crucially depend on the time that particles take to jump between consecutive potential wells) depend on the kind of potential used, whereas static exponents such as τ, d_f, and θ do not.<cit.> The analysis of the next section rationalizes this behavior.
Yet, an additional unexpected difference between smooth and parabolic potentials was observed. Figure <ref> shows curves of average avalanche duration vs. size for different system sizes. In the parabolic case as the system size increases we simply observe that the data extend to larger values of S and T. In the smooth potential case, we observe that the curves for different system sizes do not overlap
even for small avalanches. This indicates that there is a non-trivial dependence of the avalanche duration with L. If we suppose that the L dependence can be factorized as a power of L, then we can define a normalized time as

    T_n(S)=T(S,L)/L^ψ.
 
Using ψ=0.30 we obtain the curves for the normalized times seen in Figure <ref>. The collapse of all these curves is an indication that the Eq. <ref> is well satisfied.



  





The additional dependence of duration T on system size L that was observed in the case of smooth potentials allows an alternative definition of the dynamical exponent z. Instead of comparing the duration of avalanches with different sizes for a fixed L, we can compare the duration of the largest avalanche that occur for different system size. This allows to define an alternative exponent 
z_s^*=z_s+ψ. Consequently, we have z_s^* ≃ 0.72.
















































§ MEAN FIELD DESCRIPTION


We want to explore here the reasons why there are two different values of the dynamical exponents β and z for smooth and parabolic potentials, whereas those describing static properties, such as τ, d_f, θ are the same.
Let us first analyze in more detail the form of the equations of the model (Eqs. <ref>-<ref>) in real space. Since G_ q≤ 0, we see that G_ii∼∑_ q G_ q< 0. Noting k≡ -G_ii we can write  Eq. (<ref>) 
as


    ė_i=-dV_i/de_i-ke_i+∑_j iG_ije_j+σ

Now we separate G_ij in the last term as it average value, and its fluctuating part:

    G_ij= k/(N-1) +G_ij

The equations of the model are then written as

    ė_i=-dV_i/de_i+k(e_j-e_i)+∑_j iG_ije_j+σ


The kernel G_ij still has the (r_i-r_j)^-2 decay with distance, and the quadrupolar angular symmetry. But we emphasize that its spatial average vanishes: ∑_j iG_ij=0. 
Eq. <ref> is appropriate to consider approximate treatments of the model.







 §.§ Naive mean field: the Prandtl-Tomlinson problem



Neglecting the fluctuating term proportional to G_ij in Eq. (<ref>), it reduces to 

    ė=-dV/de+k(e-e)+σ

In order to obtain the flow curve (σ vs γ̇) in this limit, we note that e≡γ̇t. Defining also 

    w(t) ≡γ̇t +σ/k

it is obtained:

    ė=-dV/de+k(w(t)-e)

Written in this form, we see that e is driven by the applied w(t) on top of the potential V(e) through a spring of constant k.
According to Eq. (<ref>), the stress can be calculated as the average force on the driving spring: σ=k(w(t)-e).
This is just the Prandtl-Tomlinson (PT) model used to qualitatively describe the origin of a friction force between sliding solid bodies <cit.>. In the absence of thermal fluctuations –as it is the case here–, the PT model has a critical stress σ_c for γ̇→ 0 (as long as there are points at which -d^2 V(e_2)/d e_2^2>k ), and a power law increase of σ for finite γ̇, i.e, γ̇∼ (σ-σ_c)  ^β. The value of β turns out to be dependent of the kind of potential that is used. For smooth potentials β=3/2, whereas for parabolic potentials (with points at which the first derivative has jumps) the value β=1 is obtained <cit.>. Namely, Eq. <ref> provides a simple case in which the value of β depends on the form of the potential. 

In order to qualitatively consider the avalanche statistics and its possible dependence on the kind of potential in this mean field approach, we will go back to Eq. (<ref>), and replace the uniform force σ by driving at a constant speed γ̇ through a spring of a small stiffness k_0. 


    ė=f(e)+k(e -e)+k_0(w(t)-e)


The time fluctuations of e are of order 1/N (N being the total number of particles in the system). Although this fluctuation goes to zero in the thermodynamic limit (and then it reproduces the same average evolution than Eq. (<ref>) when k_0→ 0), it is enough to produce non trivial avalanches in the system. In fact, the statistics of the avalanches produced by a model like Eq. (<ref>) are well known.
Avalanches distribute with a cut-off power law P(S) ∼ S^-τ g(S/S_max), where g is a cut off function, τ=3/2, and the cut-off value S_max depends on k_0 as S_max∼ k_0^-2. As the value of k_0 is progressively reduced, avalanches with a critical size distribution P(S)∼ S^-3/2 are obtained. The 3/2 value of the τ exponent hold both for smooth and parabolic potentials. 








However, differences appear between smooth and parabolic potentials when considering the duration of the avalanches.
For the calculation of this time, it becomes crucial to take into account the time that an unstable site actually takes  to move to the new equilibrium position in the next potential well. In the case of parabolic potentials, the pushing force is finite as soon as the instability point is overpassed, and this implies that this time is independent of the stress excess Δσ. 
The situation is different for smooth potentials. A site that becomes unstable feels a pushing force that is 
strongly dependent on the stress excess Δσ, over the threshold for instability. It turns out that the time an unstable site takes to reach the new equilibrium position at the next potential well scales as ∼Δσ ^-1/2 <cit.>.

An analysis based on this difference between smooth and parabolic potentials (to be presented elsewhere <cit.>) leads to the conclusion that
for parabolic potentials, avalanche duration T scales with the avalanche size S as
T∼ S^1/2, whereas for smooth potentials T∼ S^1/4. 

Thus in addition to β, the dynamical exponent z is different for parabolic and smooth potentials in mean field.
This is a remarkable result. It shows that even in mean field, and in addition to the already discussed difference in the β exponent, there are differences in the dynamical exponent when comparing parabolic and smooth potentials. We remark that the static exponents 
are the same for both potentials in mean field. The exponents that are different are those
related to the dynamical characteristics of the avalanches, and the difference originates in the qualitatively different way in which a particle jumps from one potential well to the next, for smooth or parabolic potentials.






 §.§ A “dynamical" mean field approach: stochastically driven Prandtl-Tomlinson particles


By introducing the definition of w(t) into Eq. (<ref>) we obtain

    ė_i=-dV_i/de_i+k(w(t)-e_i)+∑_j iG_ije_j

This defines a set of coupled PT models, in which the variable e_i evolves under the external uniform driving w(t) on the potential V_i, and it is affected by all other e_j through the coupling term G_ij. 
We will now make a description in which this term is decoupled and treated as an external perturbation. 

To begin with, we start with a brief digression.
The accuracy of a mean field approximation depends essentially on the range of the interaction. Let us consider for the moment a standard, ferromagnetic Ising model in two spatial dimensions, with interactions decaying as 1/r^α. The values of the critical exponents depend continuously on the value of α and move towards mean field values as α is reduced. When α=2 (in general, when α is equal to space dimensionality), the model becomes mean field and the critical exponents are exactly given by their mean field values. 
A simple way to understand this result is the following. A spin in a given position interacts with a weighted sum of all other spins in the system. For α> 2 the influence of any individual spin on this sum has a finite (non-zero) weight. However for α≤ 2, the influence of any individual spin on the effective field seen by any other spin is infinitesimal (for an infinite size system). This means that fluctuation effects are unimportant, and mean field results are exact. Note that this implies not only that the exponents are mean field, but that the full solution to the problem is exactly given by the mean field approximation.

In our case, the interaction term in Eq. <ref> has precisely the ∼ 1/r^2 decay. However the sign is alternating with zero average, and this prevents the application of the arguments of the previous paragraph in a direct form. On average, the mean value of the last term is zero, it is its fluctuation in time what is relevant. In this sense, we note that although the sum of the last term is zero on average, the contribution of any individual e_i to the fluctuation is still infinitesimal when α≤ 2. This suggests that a “mean field" description should be rather accurate, if not exact, in the present case too. In this context the meaning of “mean field" is that the last term can be treated as an externally given fluctuating term, and in this way the evolution of each local variable becomes a one particle problem.


In other words, we will write formally Eq. <ref> as

    ė_i=f_i(e_i)+k(w(t) -e_i)+ξ_i(t)

where 

    ξ_i(t)=∑_j iG_ije_j

Now, ξ_i(t) will be taken to be an external noise. 
In the end, we should require this noise to be compatible with the evolution of the local variables, i.e., Eq. (<ref>) be satisfied. However, as a first step we will consider Eq. (<ref>) on its own, assuming some statistical properties of the stochastic noise ξ_i(t).

Taking into account that according to Eq. (<ref>) the time evolution of ξ_i depends on the variation rate of e_j, the statistical properties of the noise term must scale with the velocity at which the system is driven. 
We are interested mainly in the case in which driving is very slow.
In this limit, ξ_i can be considered to depend directly on the control variable in the system, that is, on the applied external strain γ̇t.
This means that the dependence of ξ_i on the strain rate can be explicitly incorporated by writing:

    ė_i=f_i(e_i)+k(w t -e_i)+ξ_i(γ̇t)

This equation defines what we call the stochastically driven Prandtl-Tomlinson model. The evolution of e_i will depend on the amplitude and correlations of the noise term ξ, as well as on the form of the force f_i(e_i).

We will consider the case (that will be shown is relevant in the yielding context) of a ξ(x) noise that has self-similar correlation properties characterized by the so called Hurst exponent H. This means that, statistically

    ξ_i(λ x) ∼λ^Hξ_i(x)


Note that a standard random walk has H=1/2.








The flow exponent of the model defined by Eqs. (<ref>) and (<ref>) was worked out in
Ref. <cit.>. There it was shown that 

    γ̇∼ (σ-σ_c)^β,

with a flow exponent

    β=1/H-1/α+1


where α is related to the analytic form of the potential at the transition point between consecutive potential wells: α=1 for parabolic potentials and α=2 for smooth potentials. Although we do not know for the moment what the appropriate value of H is, we note that from Eq. (<ref>), the difference between β values for smooth and parabolic potentials is 1/2, independently of the value of H. This is well satisfied by the results of the full simulations presented in Section III.

The present independent particle analysis gives a prediction also on the value of the θ exponent in the system.
We remind that this exponent characterizes the equilibrium distribution of distances x to the instability point. This distribution P(x) is expected to behave as P(x)∼ x^θ for small x. If there is no stochastic term in the driving (ξ=0 in (<ref>)) the value of x reduces linearly in time until destabilization, the distribution P(x) is flat and we obtain θ=0. If there is a stochastic term in the driving the value of θ is determined from the distribution P(x) of a Fractional Brownian Motion with an absorbing wall at x=0, which is <cit.> P(x)∼ x^1/H-1, i.e., 

    θ=1/H-1.

Note that this result is independent of the potential being of the parabolic or smooth type.
It is interesting to eliminate the (still undetermined) value of  H from the expressions of β and θ to obtain that in this mean field situation they are related by 

    β=θ+2-1/α

We notice again that the numerical values obtained in Section III for θ and β quite closely satisfy Eq. (<ref>).



 §.§ Finding the value of H



To check the consistency of our approach, and in particular the value of β predicted by Eq. (<ref>) we must calculate the H value of the signal ξ. In order to do this, we generate time series of ξ according to its definition as given by Eq. (<ref>) running a full simulation of the model as described in Section III.
The simulation is done in a quasistatic case, with γ̇→ 0. Examples of the ξ(t) signals that are obtained both in the parabolic and smooth cases are shown in Fig. <ref>. 






Overall, we generate five signals that correspond to points of different regions of the system, for each potential. The Hurst exponent of the signals is then obtained using the Detrended Fluctuation Analysis technique <cit.>. Essentially, the method studies how the signal growths within windows of different widths allowing to determine the value of H. The average result obtained is H ≃ 0.67 for parabolic potentials and H ≃ 0.64 for smooth potentials. Within the numerical errors, the two values coincide. This should come with no surprise at this point since we do not expect differences between the two kinds of potentials in the quasistatic limit. Moreover, when plugged into Eq. (<ref>) this value of H provides β values for parabolic and smooth potential that perfectly fit those obtained in the numerical simulations of Section III.




	
	
	
	



It remains to be understood why a value H≃ 0.65 shows up in the simulations.
In oder to address this point, we notice that according to its definition in Eq. (<ref>), the fluctuating term ξ(t) gets a cumulative contribution every time an avalanche occurs in the system. So we can try to make an estimation of the form of ξ(t) by assuming a random and uncorrelated distribution of avalanches in the system, with a size distribution characterized by some exponent τ.
Each avalanche will generate a contribution to ξ that we note δξ.
Under the assumption of uncorrelated avalanches, we can determine the distribution on increments P(δξ). If P(δξ) happens to have long tail, namely 

    P(δξ) ∼1/|δξ| ^ν+1

(with ν<2) for large |δξ|, then its random accumulation will produce a generalized random walk ξ(t) characterized by a non trivial Hurst exponent 
where H=1/ν <cit.>. 




We consider a square system of linear size L with periodic boundary conditions (Figure <ref>), and we focus on the effect of random avalanches in the system on the strain at the central point. Avalanches are assumed to occur with a size distribution S^-τ and in the two orthogonal easy directions in the system. The goal is to calculate the strain increment on the central site produced by each avalanche and mediated by the Eshelby propagator. 
In principle this problem reduces formally to the calculation of a (three dimensional) integral, but we have not been able to find a closed form of the result, so we first show the result obtained using a Monte Carlo method.
The numerical implementation of this process generates the form of P(δξ) observed in Fig. <ref>.
P(δξ) displays a power law for large δξ that becomes more robust as the system size is increased. The value of the decay exponent depends slightly on the value of τ from which the avalanches were chosen, but using the actual value of τ (τ∼ 1.4), we find P(δξ)∼δξ^-(ν + 1) with ν≃ 1.5 and thus H≃ 0.65,
which coincides with the value directly determined from the numerical simulations through (Eq. <ref>).

In Ref. <cit.>, Lin and Wyart also considered the strain fluctuations at a given site caused by the rest of the system (for the case of elasto-plastic models, and thus akin to our case of parabolic potentials, see below). Then, in a mean field approach, they 
were able to link the exponent ν in the distribution P(δξ) with the flow exponent β, finding a relation that is compatible with our results for parabolic potentials, namely β=ν for 1<ν<2 (if ν=1 they find the flow curve has logarithmic corrections).
However, they calculate the statistics of δξ assuming it is formed by random kicks from individual sites, with an intensity given by the Eshelby kernel, ∼± 1/r^2, finding P(δξ)∼δξ^-2, i.e, ν=1, which is not the result we obtain. 
The reason of the difference is  that the contributions δξ that generate the stochastic noise cannot be considered as generated in isolated points, since they are typically produced by avalanches, that are extended objects. 

The effect of avalanche size distribution on the value of ν can also be estimated using the following argument.
An avalanche of size S located at a distance D from the origin produces an increase δξ of strain at the origin with the following characteristics. If D ≫ S then δξ∼± S/D^2 (in this case the avalanche behaves as a point-like object, the ± sign is a short hand for the angular dependence of the Eshelby interaction). This behavior prevails until D∼ S. However for D≲ S the effect of the avalanche 
becomes proportional to  1/S, and independent of D, as a direct integration shows. 
Then the contribution to P(δξ) from avalanches with a fixed value S_0 takes the form:

    P(δξ|S_0)∼S_0/|δξ|^2     δξ≤1/S_0
    
    P(δξ|S_0)∼ S_0^2δ_𝖣(δξ-1/S_0)     δξ≃1/S_0

where δ_𝖣 notes a Dirac delta function.
Now the total P(δξ) is obtained integrating this result over S_0, considering the probability distribution of S_0:

    P(δξ)=∫  P(δξ|S_0)S_0^-τ dS_0

The result is 

    P(δξ)∼1/|δξ|^4-τ

which modifies the result obtained in <cit.> in the right direction: The value of H is H=1/(3-τ) which provides (using τ∼ 1.4 from the simulations in Section III) H≃ 0.62, quite close to the value H≃ 0.65 directly measured before.


All these verifications of self-consistency indicate that the treatment of the interaction term in Eq. (<ref>) as a mean field fluctuating noise is a consistent and quantitatively accurate approach.









§ RELATION WITH ELASTO-PLASTIC MODELS


It was already mentioned that one of the characteristics of the present model is that there is a single strain variable e, and no clear cut separation is made between elastic and plastic strains, contrary to what is usually done in EP models. However, for the case of parabolic potentials
this separation can in fact be proposed, and it is possible to discuss in detail the relation between EP models and the present one.

We consider our model with parabolic potentials. In this case, the central position of parabola at site i (to be noted γ^pl_i) can be identified with the plastic deformation at site i, and Eq. (<ref>) can be written as


    ė_i=μ(γ^pl_i-e_i)+∑_j G_ije_j+σ

where μ is the curvature of the potential. Note that this curvature is assumed to be equal at every potential well.
We will suppose that γ̇ is so small that it can always be assumed that e_i is in an equilibrium position. i.e., ė_i=0.
In this case we can write 

    σ_i≡-μ(γ^pl_i-e_i)=∑_j G_ije_j+σ

where the local stress σ_i has been introduced. Since all parabola have the same curvature, if the average strain increases at a rate γ̇, the value of σ_i increases uniformly in the system with the same rate, as long as no particle goes out of its local parabola. Namely

    δσ_i=μγ̇δ t

If σ_i becomes larger than the maximum stress that site i can sustain, the corresponding γ^pl_i changes to  a new value γ^pl_i+δγ^pl_i and the strains e_i will accommodate to new values satisfying Eq. (<ref>). Upon changes in γ^pl_i, the corresponding changes δσ_i in the stresses can be obtained from that equation. Working in Fourier space the result is

    δσ_ q=μ G_ q/μ-G_ qδ e_ q≡ H_ qe_ q

where G_ q is given in Eq. (<ref>). Since the denominator is strictly positive for all q, H_ q still has the same zero modes that the original G_ q and its cos(4θ) symmetry, and being independent of the norm of q (as G_ q itself), it has a decay in  real space as 1/r^2.

In this way, the previous equation gives the effect of an increase in plastic deformation on the stress in the sample. The kernel for this influence has  the  Eshelby structure ∼cos(4θ )/r^2. Such an influence of the plastic strain on the stress (Eq. <ref>), plus the linear increase of stress with applied strain (Eq. <ref>), are exactly the ingredients used for instance in the implementation of EP models given by <cit.>. 
Yet an additional consideration is necessary. In EP models it is typically (sometimes implicitly) assumed that there is a fixed time scale for a site that has overpassed its maximum stress, to move to a state with σ≃ 0 (this typically occurs at a constant rate, or in a single time step).
In our case, plastic strains change instantaneously when strain reaches the crossing between successive parabola, however, the starting Eq. (<ref>) has in fact a typical time scale τ∼ 1/μ for an unstable site to reach its new equilibrium position. It is thus clear that qualitatively, our model with parabolic potentials can be interpreted as an elastoplastic model, and then it is not surprising that we get the same critical exponents as found for instance in <cit.>.
The present comparison also suggests that the phenomenology of our model with smooth potentials might not be captured by usual elastoplastic models. But at the same time it suggests the appropriate modification in EP models to match this case too. In fact, the main difference between the dynamics of smooth and parabolic potentials seems to be the different time that it takes for a particle at a given potential well to reach the next one when it jumps over the barrier. For parabolic potentials, as we argued before, this time τ is roughly constant, independently of the stress excess over the critical value. This is what leads to consider a constant transition rate, and what makes possible the comparison with standard EP models. For smooth potentials however, the time τ that it takes to reach the new equilibrium position strongly depends on the stress excess σ_i-σ_i^c over the critical value σ_i^c, actually τ∼ (σ_i-σ_i^c)^-1/2<cit.>. In an implementation in terms of transition rates, smooth potential would require to consider stress dependent transition rates. We think that with this additional ingredient EP models can be used to reproduce also the results we obtain here with smooth potentials.

Finally, we note that there are other kinds of EP models (such as the model of Picard <cit.>, or that used by Barrat et al. <cit.>) that are directly defined in terms of its dynamics and cannot be derived from the minimization of a Hamiltonian function. It is possible, but not proved at present, that a modification of the transition rate in these models would also produce a change in the dynamical exponents β and z as we observed in our case changing from parabolic to smooth potentials.




 §.§ Relation to the Hébraud-Lequeux mean field


The Hébraud-Lequeux model <cit.>  is a further simplification on an elasto-plastic model, in which any plastic rearrangement is assumed to 
produce a random variation of stress on any other site. Note that the value of the random variation is renewed if the same site yields plastically a second time.

In our model as described by Eq. (<ref>) this random effect can be mimicked by replacing
in the last term the kernel G by a random coupling that is renewed every time e_i jumps to a new potential well. It is clear that this produces a noise term ξ(γ̇t) as in Eq. (<ref>) that is the accumulation of random contributions from all the strain jumps that occurred all across the system, i.e., a standard random walk, with a Hurst exponent H=1/2. According to our previous analysis, we know that this case provides (for parabolic potentials) β=2, θ=1 in fact, similar to the values that are obtained in the Hébraud-Lequeux model.
Note that this approximation for the case of smooth potentials produces a value β=5/2 instead.




§ CONCLUSIONS


In this paper we have investigated the critical properties of the athermal yielding transition in a two-dimensional model that includes structural disorder and long range elastic interactions as two main ingredients. Our results strongly suggest that some critical exponents depend on the form of the plastic disorder potential, finding differences between the cases of a “smooth” potential case (in which minima are smoothly connected) and a “parabolic” case (in which the potential ia a concatenation of parabolic pieces) with discontinuous forces at the transition points. The exponents that differ between the two cases are the flow exponent β and the dynamical exponent z. Other exponents are the same in the two cases. We interpret the differences as a consequence of the qualitatively different dynamics of the system around the transition points in the smooth and parabolic cases. We also claim that contrary to what happens in the depinning problem (where the two kind of potentials are known to produce no difference in the critical properties)
here the difference remains because of the long range nature of the elastic interaction. In fact, this long range nature of the interaction transforms the problem into an effective mean field one. We constructed explicitly the mean field theory describing the problem and showed it corresponds to a particle driven on top of the disordered plastic potential. The driving incorporates the mechanical noise of all other sites in the system as a stochastic contribution. We gave the values of most of the critical exponents in terms of the statistical properties of this noise, particularly its Hurst exponent H. As a consistency check we measured directly the value of H in the full simulation and also estimated it from a simplified analysis, finding H≃ 2/3. Overall, the values of the critical exponents found both for parabolic and smooth potentials, and the value of H are totally consistent, giving support to our mean field interpretation of the transition.



§ ACKNOWLEDGMENTS




We thank Ezequiel Ferrero, Jean-Louis Barrat, Vivien Lecomte and Alejandro Kolton for helpful discussions.






§ AN ALTERNATIVE DERIVATION OF A SCALAR MODEL UNDER THE ASSUMPTION OF SINGLE SHEAR


In addition to the heuristic presentation given in the main text, and the derivation from a full tensorial model given in <cit.>, we present here an alternative  derivation of the model in a case in which the deformations in the material are assumed from the beginning to be scalar.
This derivation may be applied to a case in which the material is submitted to an external single shear stress (instead of the deviatoric stress e assumed in the main text, which is composed of two orthogonal shears). Under these conditions (see Fig. <ref>) we will make the assumption that the local displacements u≡ (u_x,u_y) describing the sample deformation occur only along the x direction, namely u_y≡ 0. We will refer to u_x simply as u, and we will use two subindexes to indicate spatial positions in the sample along x and y directions. 



We assume that the interaction between different u_ij is local, i.e., it depends on the difference of u_ij on neighbor sites.
Along y, the interaction between u_i,j and u_i,j+1 must allow for “slips" between consecutive planes. 
We introduce corrugated potential energy functions V(u_i,j+1-u_i,j) with these properties. 
Along x the interaction between u_i,j and u_i+1,j is taken to be perfectly harmonic, i.e., described by an energy term μ(u_i+1,j-u_i,j)^2/2. However, and in order to get a final version as similar as possible to the one presented in the main text, we also add the corrugated potential along the x direction with a term V(u_i+1,j-u_i,j). This will not affect qualitatively the phenomenology, since different wells of the corrugated potential along x will not be explored because of the existence of the harmonic term.

The actual form of V is stochastic and depends on the actual values of i and j, but we do not indicate this in the notation. The elastic energy of the system is 


    E=∑_ijμ/2 (u_i+1,j-u_i,j)^2 + V(u_i,j+1-u_i,j)+ V(u_i+1,j-u_i,j)

We will write a first order evolution dynamics for the model. One may think to write something like



    u̇_ij= -δ E/δ u_ij

However, this equation relates the velocity of variation of u_ij to the total force acting on u_ij, and breaks Galilean invariance. It is more natural to postulate a dissipation mechanism in which viscous forces appear when there are relative motions between neighbor particles. In a mechanical analogy, instead of adding a dashpot between position u_ij and a reference position as Eq. (<ref>) implies, we add dashpots between neighbor sites on the sample. This leads to write the force balance at position i,j as

    4u̇_i,j-u̇_i,j+1-u̇_i,j-1-u̇_i+1,j-u̇_i-1,j=
    
    =μ (u_i+1,j+u_i-1,j-2u_i,j)+
    
    +F(u_i,j+1-u_i,j)+F(u_i,j-1-u_i,j) +
    
    +F(u_i+1,j-u_i,j)+F(u_i-1,j-u_i,j)

where F(x)≡ -∂ V(x)/∂ x.
This is already the model we are seeking for. However, we need to rearrange its terms in order to display its similarity
with the model presented in the text.
We define e_i,j≡ u_i,j+1-u_i,j, in such a way that combining the previous equations at sites i,j, and i,j+1 we obtain

    4ė_i,j-ė_i,j+1-ė_i,j-1-ė_i+1,j-ė_i-1,j=
    
    =μ (e_i+1,j+e_i-1,j-2e_i,j)+
    
    +F(e_i,j+1)+F(e_i,j-1)+F(e_i+1,j)+
    
    +F(e_i-1,j)-4F(e_i,j)

By introducing the notation 
∂^2_x U≡ U_i+1,j+U_i-1,j-2U_i,j, and ∂^2_y U≡ U_i,j+1+U_i,j-1-2U_i,j,  
we can write the previous equation in the compact form

    (∂^2_x+∂^2_y)ė= -μ∂^2_x e+(∂^2_x+∂^2_y) F(e)

Note that this equation does not fix the evolution of the mean value e, which must be determined according to the driving mechanism that is assumed to hold. 

Going to Fourier space and dividing by q_x^2+q_y^2, Eq. (<ref>) can be written (for q 0) as

    ė_ q= -μq_x^2/q_x^2+q_y^2  e_ q + .F(e) |_ q

In real space this equation reads (introducing the applied stress σ)

    ė_r=-dV/de_r +∑ _r'G_rr'e_r'+σ

where G_r,r is the real space form of 

    G_ q≡-μ q_x^2/q_x^2+q_y^2

In this form, the structure of the model is seen to be identical to that of Eqs. (<ref>) and (<ref>). Note that the kernel we find here has the property G_ q≤ 0, as it was the case for Eq. (<ref>). 
The only difference is in the symmetry of G, which is now dipolar instead of quadrupolar. This is naturally originated in the single shear geometry assumed in this restricted version.




Ferrero A. Nicolas, E. E. Ferrero, K. Martens, J-L Barrat, Deformation and flow of amorphous solids: a review of mesoscale elastoplastic models, arXiv:1708.09194.


FC_1 Höhler, R., Cohen-Addadv, S., Rheology of liquid foam, J. Phys.: Condens. Matter 17, R1041 (2005).


FC_2 D. Bonn,  M. M. Denn, L. Berthier, T. Divoux,  and S. Manneville, Yield stress materials
in soft condensed matter, Rev. Mod. Phys. 89, 035005 (2017).

SZ_Num_1 M. L. Falk,	 J. S. Langer, Dynamics of viscoplastic deformation in amorphous solids, Phys. Rev. E 57, 7192 (1998).

SZ_Num_2 C. E. Maloney, A. Lemaître, Subextensive Scaling in the Athermal, Quasistatic Limit of Amorphous Matter in Plastic Shear Flow, Phys. Rev. Lett. 93, 016001 (2004).

SZ_Num_3 C. E. Maloney, A. Lemaître, Amorphous systems in athermal, quasistatic shear, Phys. Rev. E 74, 016118 (2006).

SZ_Num_4 A. Tanguy, F. Leonforte, J. -L. Barrat, Plastic response of a 2D Lennard-Jones amorphous solid: Detailed analysis of the local rearrangements at very slow strain rate, The European Physical Journal E 20, pp 355–364 (2006).

SZ_Num_5 L. Gartner, E. Lerner, Nonlinear plastic modes in disordered solids, Phys. Rev. E 93, 011001(R) (2016).

SZ_Num_6 D. Rodney, C. Schuh, Distribution of Thermally Activated Plastic Events in a Flowing Glass, Phys. Rev. Lett. 102, 235503 (2009).

SZ_Num_7 D.Srolovitz, V.Vitek, T.Egami, An atomistic study of deformation of amorphous metals, Acta Metallurgica 31, pp 335-352 (1983).

SZ_Num_8 T. Albaret, A. Tanguy, F. Boioli, D. Rodney, Mapping between atomistic simulations and Eshelby inclusions in the shear deformation of an amorphous silicon model, Phys. Rev. E 93, 053002 (2016).

SZ_Num_9 C. Fusco, T. Albaret, A. Tanguy, Rheological properties vs. local dynamics in model disordered materials at low temperature, The European Physical Journal E 37, pp 37-43  (2014).

SZ_Num_10 G. J. Papakonstantopoulos, R. A. Riggleman, J-L. Barrat, J. J. de Pablo, Molecular plasticity of polymeric glasses in the elastic regime, Phys. Rev. E 77, 041502 (2008).

SZ_Num_11 A. Smessaert, J. Rottler, Distribution of local relaxation events in an aging three-dimensional glass: Spatiotemporal correlation and dynamical heterogeneity, Phys. Rev. E 88, 022314 (2013).


SZ_Exp_1 V. Chikkadi, P. Schall, Nonaffine measures of particle displacements in sheared colloidal glasses, Phys. Rev. E 85, 031402 (2012).

SZ_Exp_2 K. E. Jensen, D. A. Weitz, F. Spaepen, Local shear transformations in deformed and quiescent hard-sphere colloidal glasses, Phys. Rev. E 90, 042305 (2014).

SZ_Exp_3 P. Schall, D. A. Weitz, F. Spaepen, Structural Rearrangements That Govern Flow in Colloidal Glasses, Science 318, 5858 (2007).

SZ_Exp_4 K. W. Desmond, E. R. Weeks, Measurement of Stress Redistribution in Flowing Emulsions, Phys. Rev. Lett. 115, 098302 (2015).

SZ_Exp_5 A-L. Biance, S. Cohen-Addad, R. Hohler, Topological transition dynamics in a strained bubble cluster, Soft Matter 5, 4672-4679 (2009).

SZ_Exp_6 A-L. Biance, A. Delbos, O. Pitois, How Topological Rearrangements and Liquid Fraction Control Liquid Foam Stability, Phys. Rev. Lett. 106, 068301 (2011).

SZ_Exp_7 G. Debrégeas, H. Tabuteau, J.-M di Meglio, Deformation and Flow of a Two-Dimensional Foam under Continuous Shear, Phys. Rev. Lett. 87, 178305 (2001).

SZ_Exp_8 A. Amon, V. B. Nguyen, A. Bruand, J. Crassous, and E. Clément, Hot Spots in an Athermal System, Phys. Rev. Lett. 108, 135502 (2012).

SZ_Exp_9 A. Amon, R. Bertoni, J. Crassous, Experimental investigation of plastic deformations before a granular avalanche, Phys. Rev. E 87, 012204 (2013).

SZ_Exp_10 D. V. Denisov, K. A. Lörincz, J. T. Uhl, K. A. Dahmen, P. Schall, Universality of slip avalanches in flowing granular matter, Nature Communications  7, 10641 (2016).

SZ_Exp_11 A. Le Bouil, A. Amon, S. McNamara, J. Crassous, Emergence of Cooperativity in Plasticity of Soft Glassy Materials, Phys. Rev. Lett. 112, 246001 (2014).

fisherD. S. Fisher, 
Collective transport in random media: from superconductors to earthquakes,  
Phys. Rep. 301, 113 (1998).

kardar M. Kardar,
Nonequilibrium dynamics of interfaces and lines,
Phys. Rep. 301, 85 (1998).

Rosso_PNAS Lin, J., Lernera, E., Rosso, A., Wyart, M., Scaling description of the yielding transition in soft amorphous solids at zero temperature, Proceedings of the National academy of Sciences of the United States of America, 111, pp 14383-14387 (2014).

Narayan O. Narayan, D. S. Fisher, Critical behavior of sliding charge-density waves in 4-e dimensions, Phys. Rev. B, 46, 11520 (1992).

Jagla_Yiel E. A. Jagla, Different universality classes at the yielding transition of amorphous systems, Phys. Rev. E 96, 023006 (2017).


Rosso_2009 A. Rosso, P. Le Doussal, K. J. Wiese, Avalanche-size distribution at the depinning
transition: A numerical test of the theory, Phys. Rev. B 80, 144204 (2009).

footnote The correlation length exponent ν is also potential-independent. This follows as a consequence of the statistical tilt symmetry of the problem, that implies ν=1/(d-d_f), with d=2 in our case.

Prandtl L. Prandtl, Ein Gedankenmodell zur kinetischen Theorie der festen Körper, Z. Angew. Math. Mech. 8, 85 (1928).

Tomlinson G. A. Tomlinson, A molecular theory of friction, Philos. Mag. 7, 905 (1929).

Popov V. L. Popov and J. A. T. Gray, Prandtl-Tomlinson model: History and applications in friction, plasticity, and nanotechnologies, ZAMM. Z. Angew. Math. Mech. 92, 683 (2012).

Strogatz S. H. Strogatz, Nonlinear Dynamics And Chaos: With Applications To Physics, Biology, Chemistry, And Engineering (Studies in Nonlinearity) (Addison-Wesley Publishing Company, 2001), 1^st Edition, pp 100.

arj I. F. Aguirre, A. Rosso, and E. A. Jagla, unpublished.

Jagla_P-T E. A. Jagla, The Prandtl–Tomlinson model of friction with stochastic driving, Journal of Statistical Mechanics: Theory and Experiment 2018, 013401 (2018).



majumdar1 G. M. Molchan,
Maximum of a fractional Brownian motion: probabilities of small values, Comm. Math. Phys.
205, 97 (1999).

majumdar2F. Aurzada and C. Baumgarten Persistence of fractional Brownian motion with moving boundaries and applications J. Phys. A  46 125007 (2013).

DFA_1 R. Hardstone, S.-S. Poil, G. Schiavone, R. Jansen, V. Nikulin, H. Mansvelder, et al, Detrended 
fluctuation analysis: A scale-free view on neuronal oscillations, Frontiers in Physiology 3, 450 (2012).

DFA_2 C. Peng, S. Havlin, H. E. Stanley, A. L. Goldberger, Quantification of scaling exponents and crossover phenomena in nonstationary heartbeat time series, AIP Chaos 5, pp 82-87 (1994).
	
DFA_3 R. M. Bryce, K. B. Sprague, Revisiting detrended fluctuation analysis, Scientific Reports 2, 315 (2012).


Lin_Wyart J. Lin, M. Wyart, Some views on the Herschel-Bulkley exponent, arXiv:1708.00516 (2017).

picardG. Picard, A. Ajdari, F. Lequeux,  and L. Bocquet, Slow flows of yield stress fluids:  Complex spatiotemporal behavior within a simple elastoplastic model, Physical Review E 71, 010501 (2005).

barratC. Liu, E. E. Ferrero, F. Puosi, J-L. Barrat, and K. Martens,
Driving Rate Dependence of Avalanche Statistics and Shapes at the Yielding Transition, 
Phys. Rev. Lett. 116, 065501 (2016)

H-L P. Hébraud, F. Lequeux, Mode-coupling theory for the pasty rheology of soft glassy materials, Phys. Rev. Lett. 81, pp 2934-2937 (1998).


