Avalanche size distributions in mean field plastic yielding models

I discuss the size distribution ${\cal N}(S)$ of avalanches occurring at the yielding transition of mean field (i.e., Hebraud-Lequeux) models of amorphous solids. The size distribution follows a power law dependence of the form: ${\cal N}(S)\sim S^{-\tau}$. However (contrary to what is found in its depinning counterpart) the value of $\tau$ depends on details of the dynamic protocol used. For random triggering of avalanches I recover the $\tau=3/2$ exponent typical of mean field models, which in particular is valid for the depinning case. However, for the physically relevant case of external loading through a quasistatic increase of applied strain, a smaller exponent (close to 1) is obtained. This result is rationalized by mapping the problem to an effective random walk in the presence of a moving absorbing boundary.

I discuss the size distribution N (S) of avalanches occurring at the yielding transition of mean field (i.e., Hebraud-Lequeux) models of amorphous solids. The size distribution follows a power law dependence of the form: N (S) ∼ S −τ . However (contrary to what is found in its depinning counterpart) the value of τ depends on details of the dynamic protocol used. For random triggering of avalanches I recover the τ = 3/2 exponent typical of mean field models, which in particular is valid for the depinning case. However, for the physically relevant case of external loading through a quasistatic increase of applied strain, a smaller exponent (close to 1) is obtained. This result is rationalized by mapping the problem to an effective random walk in the presence of a moving absorbing boundary.

I. INTRODUCTION
The yielding transition (at zero temperature) of an amorphous solid material occurs when the externally applied shear stress σ overpasses a critical value σ c . For σ < σ c the system is blocked and remains in this state indefinitely. As σ > σ c the system enters a plastic flow regime in which the strain in the system increases linearly with time. This is referred to as a plastic flow regime. The transition at σ = σ c is the yielding transition. At σ c the system is in a critical state, and the dynamics proceeds via avalanches of all sizes, characterized by its size distribution N (S).
The previous scenario has been confirmed by a variety of experimental measurements, 1-5 as well as analytical and numerical techniques. [6][7][8][9][10][11][12][13][14][15][16][17][18] From all these sources, the idea has emerged that plasticity proceeds through the destabilization and rearrangement of discrete units at given spatial positions (sometimes referred to as shear transformation zones). 19,20 Rearrangement of a unit at a given spatial position has an effect (through direct elastic interaction) on the stability of other units at different locations. It is thus clear that the form of the elastic interactions between different regions of the material play a central role in the problem. In this respect, it is crucial to realize that the elastic interaction kernel has alternating signs in different spatial regions. 21 The destabilization of a given site can thus have both a destabilizing or stabilizing effect on some other site.
The previous description shares many features with the well known case of depinning transitions of elastic manifolds on disordered energy landscapes. 22,23 However, a crucial difference steps from the fact that for depinning, destabilization of a given site always produces a destabilizing effect on any other site, namely the interaction kernel is positive definite. The mean field description of the depinning transition is well studied, and is known to produce a critical size distribution of avalanches of the form N (S) ∼ S −3/2 . It has to be mentioned that this value of τ stands for (at least) two quite different proto-cols used to trigger avalanches in the system. In one case (random loading) one site is chosen at random and the local stress is increased until the site destabilizes. In the second case (uniform loading) the stress on the system is increased uniformly over all sites, until the first site becomes destabilized. In both cases τ = 3/2 is found for depinning.
I will show that in the case of mean field plastic yielding models, the value of τ that is obtained is dependent on the loading mechanism applied. While for random destabilization the value τ = 3/2 is obtained, uniform loading produces a smaller value of τ , close to 1. The reason of this difference is the much larger stress that is added to the system in the uniform loading protocol compared to the random loading case. I present the model to be studied in the next Section. An analytical treatment of the problem is given in Section III, and results of numerical simulations are presented in Section IV. Section V contains the conclusions.

II. THE MODEL
Hebraud and Lequeux 24 considered a mean field model of the plastic yielding transition. The model considers a real valued variable y i , representing the stress at each spatial position of the system. 25 The coupling among y i variables is mean field like, i.e., any change in y i affects equally all other y's. I discuss here a small variation of the original HL model, in the limit of quasistatic loading. I consider N stresses y i , i = 1, ...N , which define the instantaneous configuration of the system. I will refer to the N variables i alternatively as "sites", or "particles". The configuration is stable as long as all y i are lower than a plastic threshold value, taken for simplicity to be uniform, of value 1. It will be more convenient to describe the model in terms of the variables x i ≡ 1 − y i . Each x i represents the amount to additional stress necessary to destabilize site i. The dynamics proceeds according to the following rules.
(a) On a stable configuration, an instability is produced (either by random or uniform loading), and an avalanche starts. Suppose i is the site that becomes unstable.
(b) x i is increased a quantity (1 + k)u, where u is a random quantity, of order 1 (below I take u to be exponentially distributed: p(u) = u −1 exp(−u/u), u > 0), and k is a non-conservation parameter (in such a way that criticality is expected as k → 0).
(c) All x j are decreased a quantity η/ √ N −u/N , where η is a Gaussian variable, with zero mean and fixed standard deviation η 0 . (d): Steps (b-c) are repeated over all destabilized sites, 26 until all sites become stable (this defines the end of the avalanche). Then I return to (a).
The instantaneous stress σ in the model is defined as The size S of an avalanche is conveniently defined as the sum of all u values generated during the avalanche. Note that according to the previous rules, an avalanche of size S generates a reduction of stress in the system of value kS/N . As k → 0 the average size of avalanches goes to infinity as 1/k. This corresponds to a broad distribution of avalanches sizes in this limit, which corresponds to criticality. For finite k, there is typically a maximum size S max of the avalanches that are generated.
The previous rules indicate that each site reaching the x = 0 border is reinserted at x = (1 + k)u, then x = 0 is an absorbing border. Each unstable site also produces random kicks (of intensity η/ √ N ) plus a systematic decrease (of value −u/N ) on all other sites. In terms of a fictitious time that counts the number of sites that reach the x = 0 border, these two terms can be described as diffusive, and convective, respectively.
In the thermodynamics (large N ) limit, the equilibrium distribution of x values is characterized by a function P (x), such that N P (x)∆x gives the number of sites with x within the small interval ∆x. The form of P (x) is determined as a balance between the sites that reach x = 0 and are thus redistributed, and the change in x that each site reaching x = 0 produces. In an equilibrium situation the total change ∆P produced by a single particle being destabilized and reinserted must vanish, and this means that P (x) satisfies the equation (for k → 0) (1) with D = η 2 0 /2 of step (c) above, and P (x) = 0 for x < 0. The three terms of the r.h.s. represent respectively the effect of diffusive and convective variations of x (the two terms in the change of x j at step (c)), and the effect of reinserted sites (step (b)). The amplitude of the reinserted term, namely DP ′ (0) + P (0) takes into account that sites can reach the x = 0 border due to the diffusive, or convective evolution.
There is an important difference in the equilibrium form of P (x) when D = 0, or D > 0. In the first case (in which the previous model describes depinning) P (x) is immediately obtained from (1) as In particular, P (x) goes to a non-zero value when x → 0. For D > 0, due to the diffusive term and the absorbing boundary condition, P (x) goes to zero at x = 0. Moreover, the particle rate across the x = 0 border due to the diffusive term is proportional to dP/dx| x=0 , and for this quantity to be finite, the form of P (x) must be P (x) ∼ Ax for small x. This is a mean field results that does not hold for short range interactions, where a form P (x) ∼ x θ with a non-trivial θ is expected. For k → 0, P (x) is independent of the form of the destabilizing mechanism (random or uniform loading) since at criticality the stress redistribution caused by direct external loading is negligible compared with reaccommodations during avalanches (this is not true for finite k, see the Appendix). The form of P (x) for different values of D are presented in Section IV.

III. RESULTS FOR AVALANCHE STATISTICS
In order to calculate the avalanche statistics, the continuous description provided by the P (x) function is not sufficient. In fact, what this distribution tells is that each single site reaching x = 0, produces changes in all other x such that on average an additional site reaches x = 0, i.e., a conservative critical situation is reached. However, in order to describe avalanche distribution we must consider fluctuations.

A. Random loading
It is convenient at this point to describe the dynamics of the system in term of a discrete "time unit" corresponding to one step of points (b) and (c) of the dynamical protocol described above. Namely, at each discrete time, one kick is given to every particle according to (c), and one destabilized particle is reinserted according to (b). Let us call n i the cumulative number of destabilized particles until time t ≡ i. The condition for an avalanche started at t = 0 to survive until time t = i is that n i ≥ i. The avalanche stops at the moment in which n i < i.
Since the sites carrying different values of x, once destabilized, are reinserted through a random process (implied by the random values of the u parameter) the exact values of x i are locally uncorrelated. Due to this uncorrelated behavior, the stochastic process n i corresponds to a Poisson process with a unitary rate per time step. Defining the stochastic process m i as m i = n i − i, we see that the avalanche lasts until the first time i 0 in which m i0 = 0. Note that as one single site is reinserted at each time step, i 0 measures also the avalanche size, i.e., S ≡ i 0 .
In the large i limit m i is a symmetric random walk, with probability distribution p(m i ) given by and the previous description corresponds to its survival probability in the presence of an absorbing boundary at the origin. This surviving probability goes as t −1/2 . The size distribution of avalanches is given by the probability density of the random walk being absorbed at time t, which gives N (S) ∼ t −3/2 , i.e., τ = 3/2. I emphasize that the previous method shows this result is valid for random loading, both in the case η 0 = 0 (depinning), and η 0 > 0 (plastic yielding) as it only depends on the fact that the number of sites reaching the instability threshold is an uncorrelated stochastic process with constant rate.

B. Uniform loading
In the uniform loading case, in order to start an avalanche we look for the smallest x value in the system, say x min , and add this additional load to the system. This is equivalent to say that all x's in the system are reduced in the quantity x min . The question is if this shifting of the P (x) distribution may bring some observable effect. We will see that the answer depends on whether we are considering the case η 0 = 0 or η 0 > 0.
For η 0 = 0 (depinning case) the form of P (x) tends to a constant as x → 0, i.e, P (x) ∼ P (0) for small x. When shifting the distribution by x min to start the avalanche, we modify P (x) to a new P (x) given by P (x) = P (x + x min ) ∼ P (0) (to leading order). This means that the form of P (x) is not greatly affected by the x min shift. In fact, results of numerical simulations described below show that both loading protocols produce the same avalanche distribution N (S) ∼ S −3/2 in the η 0 = 0 case. However, the situation is different in the case η 0 > 0. Now P (x) ∼ Ax for small x, and a shift of the distribution by x min transforms it to P (x) = P (x + x min ) ∼ Ax + Ax min , namely the probability distribution gets a constant correction near x = 0 which qualitatively modifies its form. This produces an important change in the avalanche size distribution as I will show now.
Suppose we have the distribution P (x) = Ax + A 0 . The linear part generates a flux of particles through the x = 0 border characterized by the Poissonian stochastic process n i . The average rate of this process is easily seen to be equal to AD, and since we have one particle arriving per unit time, we conclude that A = 1/D. The constant term A 0 produces the arrival at x = 0 of some additional particles, that contribute to the process. Let d i be the number of additional particles arrived at x = 0 at or before time i, originated in the additional constant term A 0 . Now the stochastic process to be considered in order to determine the size of the avalanches is m i = n i − i + d i , i.e, the avalanche survives as long as m i ≥ 0. We need to characterize the contribution d i from A 0 . Due again to the fact that sites that contribute to the constant A 0 are totally uncorrelated, d i is a (cumulative) Poissonian process, however now its average rate depends on time. The time dependence of the average rate can be calculated simply noticing that it gives the average number of particles absorbed at x = 0, starting from a distribution with constant value A 0 for all x > 0. The average number d i can be obtained by direct integration of the diffusion equation, and the result is The condition for an avalanche to survive until time i is then m i = m i + 2A 0 Di/π ≥ 0. In other words, the original condition m i ≥ 0 is now replaced by m i ≥ −2A 0 Di/π, i.e., the absorbing boundary that was originally at i = 0 now retracts in time, as ∼ i 1/2 . The effect of such a moving absorbing boundary condition on a random walk has been studied in the literature (see for instance Refs. [27][28][29] ). The important point that makes the problem analytically solvable is that the position of the absorbing boundary has the same dependence in time that the internal diffusive dynamics, and this allows to transform the problem to a single ordinary differential equation. 30 The survival probability is thus observed to remain a power law in time, but with a nonuniversal exponent. The value of τ that depends on the coefficient R ≡ 2A 0 D/π of the square root time behavior in Eq. 4 (assuming a normalized underlying random walk, as in Eq. 3). It is found that τ = 3/2 for R = 0, and τ → 1 when R → ∞. We see the full form of this dependence in Fig. 1.
In order to calculate the actual value of R, we must determine what value of A 0 is appropriate in our case. It turns out that we do not have a unique value of A 0 , but a continuous distribution. In fact, given a starting distribution P (x) = Ax, the lowest x value, which defines x min is itself a random variable, with distribution We can write The values of A 0 are obtained from A 0 = Ax min = √ Ax 0 , giving for R the value (remember that A = 1/D): We can estimate a typical value of τ as the one corresponding to x 0 equal to its average value (x 0 = π/2). This gives R = √ 2, that entered in Fig. 1 gives τ ≃ 1.09. However we must notice that different avalanches are produced with different x 0 values, and the final distribution of avalanche sizes has to be calculated as a composition from avalanches taken from distributions with different τ values according to where p(τ ) = p(x 0 )dx 0 /dτ (taken from data in Fig. 1) is the probability of different values of τ , and N (S(τ )) = C τ S −τ , with C τ a normalizing constant. Taking the minimum size of avalanches to be S min ≡ 1, it results C τ = (τ − 1), and a numerical evaluation of the expression (8) produces a non-perfect power law, that can be characterized by an effective, local exponent τ eff defined as The result of this evaluation can be seen in Fig. 2. As we move to larger values of S, τ eff is more dominated by the terms with the smallest τ in Eq. (8), and τ eff decreases. It has to be noticed however, that for this effect to be appreciable we must really move to extremely large values of S. As a rule of thumb, we can say that we expect a τ between 1.1 and 1.2 to be observed.

IV. NUMERICAL SIMULATIONS
I will present here results of numerical simulation of the model described in Section II, that will confirm the scenario described in the previous Section. In all the simulations presented the values of u are taken from an exponential distribution with mean value of 1, namely p(u) = exp (−u). In this case, in the thermodynamic limit the equilibrium form of P (x) can be readily found by integrating Eq. 1, and is given (for k → 0) by The form of P (x) for different values of η 0 at criticality (k → 0) is shown in Fig. 3. Note that it must be D < 1 for the solution to exist. As it was explained before, the main feature of these curves is their behavior near x = 0. We see that that P (x) ∼ x 0 for η 0 = 0, whereas P (x) ∼ x 1 for η 0 = 0.
Next, I present results for the avalanche size distributions. Results in Fig. 4 correspond to the depinning model (η 0 = 0). They were obtained in systems with N = 10 5 . It is seen that the size distribution N (S) is a power law N (S) ∼ S −3/2 that is cut-off at large avalanche size at some characteristic typical maximum size S max , that increases as k decreases. We can obtain the increase law of S max with k as follows. An avalanche of size S generates a reduction in the systems stress of value Sk. This quantity is (on average) compensated by the stress increase in the loading stage. For the η 0 = 0 case, this increase is given on average by u, which is equal to one in our case. So we get S = k −1 . Also, S is related to S max through S ∼ S 2−τ max , which gives S max ∼ k − 1 2−τ , and in the present case with τ = 3/2, we get S max ∼ k −2 , that is well satisfied by the results in Fig. 4. Note that the results for η 0 = 0 show no difference for the two different driving protocols, namely random or uniform loading, according to the arguments given in the previous Section. Now we move to the analysis of the plastic yielding case (η 0 = 0). In this case, the two loading protocols generate different results. For random loading we find a power law distribution with τ = 3/2 (Fig. 5), equivalent to that found in the depinning case.
In the case of uniform loading the results are different, as expected from the analysis of the previous section. In Fig. 6 we see the decay of number of avalanches with size. On the range of sizes displayed in Fig. 6, a decaying exponent τ ≃ 1.1 can be estimated as k is decreased. I notice the rather strong overshoot of the size distributions near S max , that I analyze in the Appendix. It is seen that as k decreases, a power law with a value smaller than the standard τ = 3/2 value is observed. It is also apparent the existence of a rather strong overshoot at the range on the largest avalanches observed.
To gain some additional insight on the results in Fig.  6, I present some additional analysis concerning the avalanche size and the values of x min at which they are triggered. In Fig. 7 we see the distribution of values of x min obtained along a simulation with N = 10 5 , k = 0.1. We see that the distribution follows very closely the distribution given by Eq. (5), which in turns originates in the linear distribution of x values near x = 0.
In addition, if avalanches are independent, and they follow the results from the theory of the previous section, we should expect a different decaying exponent for avalanches if they are classified according to the value of x min at which they were triggered. In order to check for this, in the inset to Fig. 8 I show partial size distributions of avalanches restricted to given intervals of values of x min that trigger them. If we concentrate in the low size part of this plots, the results follow a trend compatible with the theoretical analysis: avalanches started with smaller x min have a distribution with a decaying exponent close to 3/2, while those with larger x min produce a distribution with a slower decaying exponent. There seems to be however systematically lower values in the simulations compared with the theoretical results. There is an additional effect in the uniform loading plastic yielding case that deserves to be mentioned. For depinning, or random loading plastic yielding, I have argued that the mean size of the avalanches must satisfy the relation Sk ∼ 1, since the average external load between avalanches is order 1. In the case of uniform loading in plastic yielding, the average loading between avalanches is order N x min ∼ √ N . So in this case we have Sk ∼ √ N , or S ∼ √ N /k. Introducing the size of the largest avalanches in the system S max we get S max ∼ ( √ N /k) 1/(2−τ ) . This indicates that the size of the largest avalanches in the system are not completely determined by the value of k. There is also a direct dependence on N . In particular, at any fixed value of k, S max becomes arbitrarily large as N → ∞.

V. CONCLUSIONS
In this paper I have studied the avalanche size distributions in mean field models of the plastic depinning transition. The main conclusion is that contrary to the case of elastic depinning where the mean field exponent is well defined and given by τ = 3/2, here the result depends on details of the driving mechanism. I have shown that for plastic depinning, the τ = 3/2 exponent is recovered in the case of random loading, where each avalanche is started by making unstable a random site. However, the more standard loading mechanism corresponding to increase uniformly the stress in all the system until an instability is obtained, produces a significantly lower exponent τ ≃ 1.1 − 1.2. The origin of this effect was elucidated by mapping the present problem to the survival probability of a random walk in the presence of a moving absorbing boundary. In addition, this mapping provides insight into the size crossover for cases in which the system is not right at criticality.
The present results call the attention on the precise destabilization mechanism used in different models of plastic depinning. It is likely that the effects discussed here remain in more realistic cases (not mean field) where a more realistic symmetry and interaction range of the elastic kernel is considered. The kind of variation of the τ exponent in those cases remains to be elucidated.

VI. ACKNOWLEDGEMENTS
I thank Alberto Rosso for stimulating discussions, and P. Le Doussal for bringing to my attention Refs. 27,28 . This research was financially supported by Consejo Nacional de Investigaciones Científicas y Técnicas (CON-ICET), Argentina. Partial support from grant PICT-2012-3032 (ANPCyT, Argentina) is also acknowledged.
Appendix A: The approach to criticality, and the large size overshoot It has already been noticed (in connection with data in Fig. 6) the important size effect observed at large avalanche sizes, where an excess number of avalanches (compared with the results for random loading, Fig. 5) is observed. In order to study this effect in more detail, it is not enough to consider the critical (k = 0) case, but it is necessary to analyze the form in which the system evolves when a finite k is reduced towards 0.
I restrict here to the plastic depinning case (η 0 = 0). The modifications in the balance equation (1) produced by a finite k are different in the cases of uniform, or random loading. In the case of random loading, a finite k produces (according to the rules of Section II) that the reinsertion does not occur at the position u, but at u(1+k). This alters the argument of p and the normalization of the last term in (1). In addition, for finite k, the average size of the avalanche =1/k is finite, and we must take into account that 1 every 1/k particles is destabilized by the random mechanism that initiates the avalanche. As this random selection is made on the actual distribution P (x), this produces a negative term −kP/N per destabilized particle in the balance equation, that finally reads (A1) In the case of uniform loading, we still have the change u → u(1 + k) in the reinsertion term, but now the uniform loading mechanism produces a shift of P in a quantity x min every S destabilized particles. This gives a correction δP in the balance equation that is equal to δP = (dP/dx)dx = (dP/dx)(x min /S) = (dP/dx)(k/N ) where the last equality follows (on average) from the balance between the stress increase during loading (x min ), and the stress decrease during avalanche (kS/N ). Finally, the balance equation for uniform loading takes the form For our choice p(x) = exp(−x) the solutions to both equations are for random loading, and for uniform loading.
In the limit N → ∞, the distribution of avalanches is controlled by the form of these expressions near x = 0. Close to criticality, we must search for the leading terms in an expansion in powers of k. For random loading, linearizing Eq. (A3) near x = 0, and writing the coefficient to first order in k we obtain where we see that as k → 0, P (x) tends to the critical form P (x) = x/D. The finite value of k provides a cut off in the size distribution of avalanches. In fact, in the mapping to the random walk problem, the linear in k term of the previous expression corresponds to an absorbing moving wall at m i = ki/D (in the notation of Section III). This problem can be analytically solved in the large avalanche size limit, and provides an avalanche distribution given by This gives a size cutoff S max that increases a 1/k 2 as k → 0, and results that fit nicely the curves in Fig. 5. For uniform loading instead, it is easily verified that Eq. (A4) is independent of k to first order in x. 31 This means that in order to obtain an avalanche cutoff we must expand to second order in x. The result (to first order in k in the coefficient of x 2 ) is This form of P (x) allows to explain in detail the behavior observed in Fig. 6 for the uniform loading case. In fact, when k = 0 we are at criticality, and we can describe the avalanche distribution through the model of Section III, of the survival of a random walk, where the loading to x min corresponds to the presence of a retracting absorbing wall at ∼ −x min √ t. The linear in k term in Eq. (A7) corresponds to an additional current of particles through the absorbing wall that has to be subtracted (because it has a negative coefficient) from the total. The number of absorbed particles from this term until time t is given by −4k (1−D) √ Dπ t 3/2 . It represents an additional "forward movement" of the absorbing wall that is ultimately responsible for the decaying of the random walk in a finite time (i.e, avalanches have a maximum value). We end up with the problem of an unbiased random walk m(t) with an absorbing boundary b(t) moving with the law b(t) = −x min t 1/2 + kt 3/2 (A8) (some non-essential coefficients have been set to 1 in this expression). It is straightforward to numerically simulate this process, finding the size of the avalanches as the first time t 0 at which m(t 0 ) = b(t 0 ), namely when the walk is absorbed, and then making statistics on the values of t 0 . The results are shown in Fig. 9. Simulations are presented for different k values, whereas x min are taken from its known distribution (Eq. (5)). The results give a nice verification of the form of the cut off that was observed in the full simulations of Fig. 6. I notice that an analytical solution to the survival probability of a random walk in the presence of the absorbing wall as in Eq. (A8) is not known. However, the form of the cutoff can be heuristically determined by inserting Eq. (A8) in the asymptotic form expected for a solution of the diffusion equation of the random walk at large x, namely ∼ exp(−x 2 /2t), which gives for the cut off the form ∼ exp −k 2 S 2 /2 + x min kS (A9) We can thus expect the form of N (S) to be given by where τ is given as a function of x min in Section II. In Fig. 9 we see that this gives a remarkably good fitting of the numerical results. From Eq. (A10) it is clear that the intensity of the overshoot becomes larger as x min increases. 32 Note also that the scaling of S max as 1/k is immediate form expression (A10).