Mean-field description of plastic flow in amorphous solids

Failure and flow of amorphous materials are central to various phenomena including earthquakes and landslides. There is accumulating evidence that the yielding transition between a flowing and an arrested phase is a critical phenomenon, but the associated exponents are not understood, even at a mean-field level where the validity of popular models is debated. Here we solve a mean-field model that captures the broad distribution of the mechanical noise generated by plasticity, whose behavior is related to biased L\'evy flights near an absorbing boundary. We compute the exponent $\theta$ characterising the density of shear transformation $P(x)\sim x^{\theta}$, where $x$ is the stress increment beyond which they yield. We find that after an isotropic thermal quench, $\theta=1/2$. However, $\theta$ depends continuously on the applied shear stress, this dependence is not monotonic, and its value at the yield stress is not universal. The model rationalizes previously unexplained observations, and captures reasonably well the value of exponents in three dimensions. Values of exponents in four dimensions are accurately predicted. These results support that it is the true mean-field model that applies in large dimension, and raise fundamental questions on the nature of the yielding transition.

Amorphous solids such as emulsions, glasses or sand are yield stress materials, which fail and flow if a sufficient shear stress is applied. In the solid phase, plasticity can be conceived as consisting of elementary rearrangements, the so-called shear transformations [1][2][3][4][5]. Shear transformations are localized but display long-range elastic interactions [6], and organize dynamically into elongated highly plastic regions [4,[7][8][9][10]. Above some threshold stress, failure occurs and one enters a fluid phase where a stationary flow can be maintained. In various materials, rheological properties appear to be controlled by a critical point at the yield stress Σ c where the flow stops: at that point, flow curves relating shear stress and strain rate are singular [11], and the dynamics displays long-range spatial correlations [9,12,13]. Despite the importance of these properties in a variety of phenomena including earthquakes and landslides, a quantitative microscopic description is lacking. As is generally the case in condensed matter systems, one expects that the density of elementary excitations strongly affects such properties. For amorphous solids this corresponds to the density P (x) of shear transformations, characterized by the additional shear stress x required to trigger them [14,15]. One finds empirically a pseudo-gap, i.e. P (x) ∼ x θ , for small x with θ > 0. The value of θ was argued to control the singular rheological properties and diverging length scale of flow just above the yield stress Σ c [16], and implies crackling (system spanning avalanche-type response) in the entire solid phase Σ < Σ c [17], where Σ is the applied shear stress. However, little progress has been made to predict the value of θ theoretically. Most recent data indicate that (i) for quasi-static flows (i.e. at Σ c ) θ ≈ 0.5 and θ ≈ 0.35 in two and three dimensions respectively [16,18], (ii) right after a fast quench at zero stress, θ ∈ [0.4, 0.6] [15,19] and (iii) as the stress increases within the solid phase, θ rapidly drops initially, and then slowly rises again as Σ c is approached from below [17].
Pseudo-gaps appear in a variety of glassy systems [20]. The associated exponent is constrained by stability, as occurs in electron glasses [21], mean field spin glasses [22][23][24][25][26] and hard sphere packings [27][28][29][30][31]. However the associated stability bound can be proven to be saturated (a scenario referred to as marginal stability) only if the interaction is sufficiently long-range [20]. For amorphous solids, stability toward extensive avalanches requires θ > 0 [19]. This inequality proves the existence of a pseudo-gap but is not saturated, and thus does not yield the exponent value. Another route seeks progress by considering mean-field models, that would allow one to compute θ in large spatial dimension d. Hebraud and Lequeux (HL) [32] introduced a popular model where all shear transformations interact with each other, with a similar magnitude. A pseudo-gap is predicted, but one finds θ = 1 which is far from the values observed in two and three dimensions, and does not depend on the applied stress. Lemaitre and Caroli (LC) pointed out that since elastic interactions decays as an inverse power of distance, the magnitude between two randomly chosen shear transformations is broadly distributed [14]. Including this effect led to a mean-field model which was numerically shown to display a pseudo-gap, but which was not solved analytically.
In this letter we introduce a class of mean-field models that interpolate continuously between these two cases, and solve them using a combination of probabilistic arguments and analysis. In our models, the variables x describing the stability of shear transformations undergo biased one-dimensional Lévy-flights of index µ with absorbing conditions outside a compact interval, and reinsertion within this interval. HL model corresponds to µ ≥ 2, whereas the more physical LC model corresponds to µ = 1. Our findings are that (a) for µ > 1, θ is independent of system preparation and follows θ = µ/2 for µ ∈ (1, 2]. (b) For µ < 1, θ = µ/2 after an isotropic (Σ = 0) quench but θ = 0 if Σ > 0, in particular at the yield stress Σ c . (c) For the physical case µ = 1, θ = 1/2 after a quench, but θ is not universal for Σ > 0, and is shown to drop immediately as Σ increases from zero, and then increases continuously with the applied stress. These predictions are confirmed numerically. They are remarkably consistent with the finite-dimension observations (ii, iii) listed above, and support that the nonuniversal value θ(Σ c ) tends to decrease with d, in agreement with (i), but will never reach a well-defined value above some critical dimension.
Mean field model: an unstable site (red line) returns to xi = 1 after a time τr. Concomitantly, other sites implement a drift toward negative x to conserve the mean stress, and a random jump of symmetric distribution w(ξ). The dashed line represents P (x), assuming Σ > 0. A pseudogap is represented at x = 0.

MEAN FIELD MODELS
Following [12,19,[33][34][35], we describe amorphous materials as consisting of N blocks, each characterized by a local shear stress σ i and a local yield stress σ th i , which we choose to be unity. The shear stress applied on the material is Σ = i σ i /N , and we assume that Σ is constant in time (we relax this hypothesis below). A block becomes unstable if |σ i | > σ th ≡ 1. This condition is easily expressed by defining x i ≡ 1 − σ i and corresponds to x i ∈ [0, 2]. If the variable x i exits this interval at some time t, the block i is declared "unstable". After some constant time interval τ r (describing the time for a local plastic event to occur) σ i relaxes to zero, i.e. x i → 1 at time t + τ r . For this choice of dynamics, x i / ∈ [0, 2] are always absorbing conditions for the variable x i (another popular choice of dynamics assumes that if x i / ∈ [0, 2], there is a finite evaporation rate at which sites become unstable. Such a dynamic leads to an absorbing condition outside [0, 2] only in the quasi-static limit, where our results below must apply. Indeed in this limit the pseudo-gap exponent does not depend on the choice of dynamics [16]).
Such a relaxation event corresponds to a plastic strain in site i of magnitude ∆ i = σ i /E, where E is the shear modulus. This relaxation causes a global plastic strain ∆ = ∆ i /N . Moreover, the stress in i is redistributed on other sites, leading to: where G ji ( r i − r j ) is the interaction kernel which a priori depends on the position r j of the sites, and the integer l numbers plastic events in chronological order. If the stress is fixed one must have j =i G ij = −G ii = −1. In the following, we set E = τ r = 1. For amorphous solids, G(r) is well-approximated by the Eshelby kernel, of magnitude G(r) ∼ 1/r d and whose sign depends on the relative directions between ( r i − r j ) and the imposed shear [6,36,37]. This property implies that the distribution ρ(∆x) of kicks ∆x j = x j (l + 1) − x j (l) at each plastic event is broadly distributed, since site i can be close or far from site j. Assuming σ i (l) = 1 in Eq.(1), one readily finds [14] that ρ(∆x) ∼ |∆x| −1−µ with µ = 1. It is straightforward to extend this result to the general case G(r) ∼ 1/r α , where we find µ = d/α. Extending [14], mean field models can now be constructed where the distribution of kicks amplitudes preserves that of the original problem, but where all sites interact statistically in an identical way. This is achieved by replacing the prescription Eq.(1) for the relaxation of site i by the new rule: x i (l + 1) = 1 In the last equation, the first term on the right side (referred to as drift term below) ensures conservation of stress, i.e. i x i (l) does not depend on l. The random variable ξ has zero mean ξ i = 0, and its distribution w(ξ) mimics that of the finite dimensional model: with a lower cut-off ξ c = ( 2A µ ) 1/µ N −1/µ fixed by normalization and an upper cut-off ξ m = ( 2A µ ) 1/µ (such a cut-off is present in the finite-dimensional model, and corresponds to the amplitude of the kick given by an adjacent site). The dynamical rule of the mean field models are illustrated in Fig.1.
Such mean field models behave qualitatively like standard elastoplastic models: there exists Σ c > 0 such that for Σ < Σ c , the dynamics eventually stops, corresponding to the solid phase. For Σ > Σ c , the dynamics never stops in the thermodynamic limit, and is characterized by a rate of plastic strain˙ = N u σ u /N , where N u is the instantaneous number of unstable blocks, and σ u = 1 − x u their mean stress value when they relax. σ u has both a positive contribution from x < 0 and a negative one from x > 2, and the symmetry between these is broken as soon as Σ > 0. In our convention, for Σ > 0, x < 1 and more sites become unstable at the boundary x < 0, leading to˙ > 0. Ultimately, our model describes Lévy flights with absorbing conditions for x / ∈ [0, 2]. Due to the drift term in Eq.(2) these flights are biased, which tends to bring them toward the unstable region x < 0 if Σ > 0. For Σ > Σ c where a stationary state is reached, computing the pseudo-gap in this mean-field approximation requires to obtain the stationary distribution of the stable sites P (x) ≡ stable δ(x − x i )/N of biased Lévy flights near an absorbing boundary. Table I. Summary of results: mean field values (MF) of θ at the yield stress Σc and after a quench at Σ = 0 as a function of the Lévy index µ, the random kick amplitude A, and the bias v. For comparison we also report θ(marginality) corresponding to the saturation of the stability bound derived in [19].

CONTINUOUS DESCRIPTION
We consider the limit N → ∞ while keeping the variable γ ≡ l/N fixed, with γ 1 (γ is essentially a measure of the accumulated plastic strain , as = γ σ u where σ u is of order one for Σ ∼ Σ c ). In this limit Eq.(2) becomes: where the random kicks ξ(γ) satisfy the probability dis- We assume that v > 0 and will relax this hypothesis when discussing thermal quenches at Σ = 0. Eq.(4), together with our rule that unstable sites are reinserted in x = 1 leads to the master equation: with the condition that P (x) = 0 if x / ∈ [0, 2], and where w (ξ) = lim γ→0 w γ (ξ)/γ.
Case µ ≥ 2: w(ξ) has then a finite variance, and w(k) = 1 − ξ 2 k 2 + O(k 4 ). Forw γ (k) to converge in the large N limit, one must choose . w (ξ) then converges to a Gaussian, and Eq.(4) leads to the standard Fokker-Planck equation: Solutions of such a diffusion equation vanish linearly near the absorbing condition, i.e. P (x) ∼ x at small x, as found in HL [32]. This corresponds to θ = 1 Case µ < 2: In that case, one recovers the well-known results for Lévy distribution [38]: with an upper cut-off at ξ m , see Appendix.A for details, Eq.(4) then leads to:

ASYMPTOTIC SOLUTIONS
We seek stationary solutions of Eq.(8) of the form P (x) = P 0 + C 0 x n for x 1, with n > 0. Defining s = y/x, Eq.(8) reduces to: We denote the last term T 3 . In the limit x 1, T 3 converges to: Case µ < 1: T 3 is always negligible in Eq. (9). Keeping only the dominant terms, we obtain: implying n − 1 = −µ. Thus we find: corresponding to θ = 0. Case 1 < µ < 2: In that case, Eq.(9) cannot be satisfied if P 0 > 0, and T 3 is therefore not negligible. Equating the dominant terms leads to: If n > µ, the first terms tends to 0, while the second term remains O(1). Thus n ≤ µ, and: We checked that the unique solution of this equation is n = µ/2, a result which has a simple probabilistic interpretation, as discussed below. Thus P (x) ∼ x µ/2 and θ = µ/2. Case µ = 1: The most important case physically is also the richest. Solution can exist only if P 0 = 0 and n < 1, and Eq.(9) implies asymptotically: The last integral yields I 1 = ∞ 0 s n −1 |1−s 2 | ds = 1 − πn cot(πn). Thus: implying that θ continuously depends on the drift v and the magnitude of the noise A.

INTERPRETATION
A detailed probabilistic derivation of these results is given in Appendix B. In a nutshell, P (x) for x 1 is proportional to the number of walks starting from x = 1 which ends up in position x after a time γ = O(1), without having crossed the absorbing condition x < 0. For a Brownian motion (corresponding to µ ≥ 2 in our model) it is well known that this number vanishes linearly in x. This result is independent of the bias for a Brownian motion, because on small length scales of order x (or equivalently on small time scale of order γ ∼ x 2 ), fluctuations always dominate the bias, which is therefore irrelevant. Fluctuations also dominate the bias for Levy flights if µ > 1, and thus P (x) at small x must again be independent of the bias, in agreement with our result that θ = µ/2 if µ ∈ (1, 2]. The case µ = 1 however is marginal: bias and fluctuations are always comparable on all time scales, and both affects the value of θ as shown in Eq. (16). Finally, for Levy flights with µ < 1, the bias dominates on small time scales: typically walkers are essentially convected toward the origin, leading to a non-vanishing P (x) at x = 0.

NUMERICAL TESTS
We now test our predictions numerically, considering first the case Σ = Σ c > 0, implying v > 0. To compute numerically P (x) at Σ c , we use the extremal dynamics method, see e.g. [35]. Starting from a small value of the stress, it is increased each time the dynamics stops, so as to trigger new avalanches of plasticity. During avalanches, the stress is lowered. Using such dynamics the system spontaneously reaches the stationary state where Σ = Σ c . Fluctuations of stress vanish in the thermodynamic limit, and we expect our fixed stress predictions to hold, as we confirm numerically. Fig.3 shows P (x) for µ = 0.8 and µ = 4/3 for two different choices of kick amplitude A. For µ = 4/3, we measure θ by fitting the part of the curves that overlap for different system sizes, and find θ = 0.61±0.04 (A = 0.15) and θ = 0.63 ± 0.03 (A = 0.35). These results are slightly smaller but close to the predicted value θ = 0.667. For µ = 0.8, we fit P (x) by the functional form P 0 + C 0 x 1−µ predicted in Eq. (12). The fit is very good as shown in Fig.3(b).
For µ = 1, θ continuously depends on the kick amplitude, A and the bias v. We plot the measured value of θ and the theoretical prediction in Fig.4, and find once again a very good agreement.

TRANSIENT BEHAVIOR
Consider a liquid state with Σ = 0. We model it as a configuration where many blocks are unstable, due to thermal fluctuations. The total initial distribution (including stable and unstable sites) P t (x) must be symmetric around x = 1, and display tails in the unstable regions x < 0 and x > 2. Next, we suddenly quenched the system by setting the temperature T to zero. Importantly, the symmetry Σ = 0 imposes that in the dynamics that follows, the same number of sites become unstable at x < 0 and x > 2, implying v = 1− x u = 0. According to Eq.(16) we thus expect θ = 1/2 in our mean-field approximation. This prediction is consistent with the molecular dynamics simulations of [15] which find θ ≈ 0.6 after a quench both for d = 2 and d = 3. It is tested numerically in our model in Fig.5 where we find θ = 0.53±0.03 for an initial condition where P t (x) is uniform in [−1, 3]. Numerically we found consistent results as long as enough unstable sites are initially present.
This situation dramatically changes however as soon as Σ increases from 0. Avalanches are then triggered, and the stress v.s. plastic strain curves (experiments generally report the stress vs the total strain, which is a plastic strain plus an elastic contribution Σ/E), although smooth in the thermodynamic limit, consists of steps as shown in inset(a) of Fig.5. Inside these avalanches (horizontal segment in inset), the stress is fixed and one can measure a drift v, as shown in inset(b) of Fig.5. However, when the stress goes up in between avalanches (vertical segment in inset), all sites are shifted toward negative x, leading to an additional contribution to the drift. Its magnitude in the thermodynamic limit follows dΣ/dγ = vdΣ/d . This contribution is large (and dominant) initially and vanishes at Σ c due to the shape of the stress-strain curves displayed in inset (a) of Fig.5. Using Eq.(16), we obtain the mean-field prediction for θ in a transient: This prediction is tested in Fig.5 and works remarkably well. Most importantly, Fig.5 is very similar to what is found in finite-dimensional elasto-plastic model [17]. It would be very interesting to obtain such data experimentally (in particular via finite-size effects such as in [15]) or in molecular simulations.

ROLE OF SPATIAL DIMENSIONS
In finite dimension the interaction kernel G is welldescribed by the Eshelby kernel [6,36,37]. For d = 2 for example it follows G ji ( r i − r j ) ∼ cos(4φ)/| r i − r j | 2 where φ is the angle between the shear direction and r i − r j . To quantify finite-dimensional effects, we consider the mean-field model obtained by shuffling G ji randomly at each event, i.e: where k l (j) is a random permutation of all indices j = i. We then measure θ at Σ c , considering both the twodimensional and three-dimensional Eshelby kernel, and show our results in Fig.6.a. We find that θ SF 2D = 0.39 ± 0.02 and θ SF 3D = 0.29±0.02. These values are in very good agreement with the prediction Eq.(16), θ M F 2D = 0.37 ± 0.05, θ M F 3D = 0.29 ± 0.04, using direct measurements of v (which is very close to 1) and extracting A directly from the spatial kernel, as done in Fig.6.b.
These mean-field values for θ are systematically smaller than observations in finite dimensions where θ 2D ≈ 0.55 and θ 3D ≈ 0.35 [16,18]. This result indicates that finitedimensional effects still matter, in particular for d = 2. For d = 3, the difference is small (and perhaps not significant), and we expect it to vanish as d increases. It is presently unclear if this difference goes to zero above some critical dimension. Our analysis implies that even if this is the case, the values of θ will not be a simple number even for larger dimensions, as it is not universal. CONCLUSION We have analytically solved a mean-field model of the plasticity in amorphous solids, focusing on the exponent θ characterizing the density of shear transformations and the stability toward avalanches. The behavior of this model is very rich: θ is found to be universal after an isotropic quench, but to be otherwise history-dependent and non-universal. The remarkable similarity between these predictions and previous observations in finite dimensions, that improves with increasing dimension, supports that the model investigated is the "good" meanfield model, i.e. that it will capture correctly exponents in high dimensions. Although we have focused on amorphous solids, it is very plausible that this model applies to disordered crystals as well, where plasticity is mediated by dislocations whose motions interact with the same Eshelby kernel studied here. It will thus be very interesting to test our predictions in both classes of materials.
Our work is consistent with the notion that the yielding transition at Σ c is a dynamical phase transition, but supports that it is a transition of a curious kind, where exponents can depend continuously on parameters. It is still unclear if exponents in finite dimensions can be computed via a perturbation around some critical dimension, where the mean-field solution becomes exact but nonuniversal. Finally, to our knowledge the first example where violation of marginal stability (i.e. the fact that the pseudo-gap exponent is strictly larger than what required by stability in a glassy system) can be proven. The concept of marginality has been very influential in electron glasses [21] but its validity is still debated in that context [20,39]. Introducing mean-field models of the type discussed here may shed light on these questions.
It is a pleasure to thank J. The probability distribution of ξ is with the lower cut-off at |ξ| c = ( 2A where B 0 = 1 2(2−µ)Iµ . In the large N limit, the coarsegrained distribution w γ (k) becomes where l µ = (2AI µ ) 1/µ . We are interested in w γ (ξ) in the limit γ → 0, so it is convenient to definẽ Making use of Eq.(A.5), we can decompose w γ (ξ) as Here w 1 = 1 lµγ 1/µ L µ (ξ), and in the limit we are interested in γ → 0, it reduces to We discuss the w 2 contribution in the two regimes: Case I: 1 ξ ξ m . In this regime, we can expand the exponential terms to the first order, and obtain 10) which turns out to be much smaller than Eq.(A.9) and therefore negligible in this regime! Case II: 1 ξ m ξ . In this case, we can extends the integral to infinity, because of the oscillating factor, where the second term cancels w 1 , and the resulting distribution w γ (ξ) in this limit is which is much smaller than Eq.(A.9) at ξ ≈ ξ m , and behaves as an upper cut-off in w γ (ξ). Our numerical calculation of w γ (ξ) is consistent with the above theoretical analysis, shown in Fig.A.1.

B: PROBABILISTIC INTERPRETATION
We now present a probabilistic interpretation of these results. The dynamics of a single block is equivalent to a random walker with Lévy index µ in x space with a bias −v towards the absorbing boundary at x = 0. We introduce the Hurst exponent for the mean absolute displacements without bias for Lévy flights H = 1/µ if µ < 2, and H = 1/2 if µ ≥ 2. We define P + (x, x 0 , −v, γ) as the probability to find the walker starting from x 0 with a bias −v at position x after a time γ, without having crossed the x < 0 half-line (in our argument the presence of another wall at positive x is irrelevant for the small x behavior of interest here).
A block that has just relaxed starts from x 0 = 1, and reaches x after some time γ. In the stationary state, we obtain For small x, this integral is dominated by γ 1. Indeed v 1, thus for γ 1 the walker is still far from the origin and has negligible chance to hit the origin due to the random noise. Instead, for γ 1 the walker is very likely to have entered the unstable region due to the bias. Thus we have for small x: For Lévy flight with no bias, it is known that the absorbing boundary leads to P + (x, x 0 , 0, γ) ∼ x θ at small x [40,41]. We define the persistence probability S(x 0 , γ) as the probability for the walker starting at x 0 > 0 to have remained in the positive axis at time γ. In the limit x 0 γ H , S(x 0 , γ) ∼ γ −κ , and κ is called the persistence exponent. There is a general scaling relation in the absence of bias, θ = κ/H [42]. In addition, for Lévy-flights according to the Sparre Andersen theorem [43] κ = 1 2 , so for a random walker without bias (v = 0), θ = µ/2 for µ < 2 and θ = 1 for µ ≥ 2.
We now extend these results to the case where there is a bias. Time-reversal symmetry implies: Thus we transfer the problem to the backward process, where x is now the starting position, with a drift away from the absorbing boundary. We expect P + (x 0 , x, v, γ) to have the scaling form: Integrating over x 0 , we obtain the persistence probability S(x, v, γ) = γ H−H G( x γ H , v γ H−1 ) where G is some function. Because S(x, v, 0) = 1, we must have H = H, therefore We define the normalized surviving probability density as p x,v (x 0 , γ) = P + (x 0 , x, v, γ) where we use the scaled variable y = x/γ H , ω = v/γ H−1 , and p y,ω (y 0 ) = p x,v (x 0 , γ)γ H . It is clear from its definition that p y,ω (y 0 ) must converge to a constant p 0,ω (y 0 ) as y → 0. Thus the asymptotic behavior of F in the small y limit is that of the surviving probability G(y, ω). S(x, v, γ) ∼ G(y, 0), i.e. the drift term v is irrelevant at small times. In the absence of bias it is known that G(y, 0) ∼ y 1/2H [40]. As γ increases toward one, the effect of the bias becomes of the order of the noise, we thus expect it to only affect the survival probability by a numerical prefactor, implying that the result: holds even with a finite bias, as proven above. Case II: µ < 1. In this case, H > 1, and the bias is relevant at small times according to Eq.(A.18). In that limit we thus get S(x, v, γ) ∼ G(y, ∞) = 1. Once again when γ increases and becomes of order one, the effect of the noise becomes of order of the bias, and will affect the survival probability by some numerical prefactor. We thus get G(y, 1) ∼ 1 for small y, leading to θ = 0 as derived above.