Criticality in elastoplastic models of amorphous solids with stress-dependent yielding rates

TL;DR

This paper investigates elastoplastic models of amorphous solids near the yielding transition, identifying universal static critical exponents and dynamic exponents sensitive to local yielding rules, with analytical and numerical support.

Contribution

It introduces two local yielding rules, establishes universal static critical exponents, and reveals how dynamic exponents depend on the specific yielding rate, supported by analytical calculations and simulations.

Findings

Universal static critical exponents are independent of yielding rules.

Dynamic exponents, including the Herschel-Bulkley exponent, depend on the local yielding rule.

Analytical support from the Hraud-Lequeux model confirms the numerical observations.

Abstract

We analyze the behavior of different elastoplastic models approaching the yielding transition. We propose two kind of rules for the local yielding events: yielding occurs above the local threshold either at a constant rate or with a rate that increases as the square root of the stress excess. We establish a family of "static" universal critical exponents which do not depend on this dynamic detail of the model rules: in particular, the exponents for the avalanche size distribution $P(S)\sim S^{-\tau_S}f(S/L^{d_f})$ and the exponents describing the density of sites at the verge of yielding, which we find to be of the form $P(x)\simeq P(0) + x^\theta$ with $P(0)\sim L^{-a}$ controlling the extremal statistics. On the other hand, we discuss "dynamical" exponents that are sensitive to the local yielding rule details. We find that, apart form the dynamical exponent $z$ controlling the duration of avalanches, also the flowcurve's (inverse) Herschel-Bulkley exponent $\beta$ ($\dot\gamma\sim(\sigma-\sigma_c)^\beta$) enters in this category, and is seen to differ in $\frac12$ between the two yielding rate cases. We give analytical support to this numerical observation by calculating the exponent variation in the H\'ebraud-Lequeux model and finding an identical shift. We further discuss an alternative mean-field approximation to yielding only based on the so-called Hurst exponent of the accumulated mechanical noise signal, which gives good predictions for the exponents extracted from simulations of fully spatial models.

Tables

Table 1. Table 1: Values of the β 𝛽 \beta exponent according to Eq. 25 for different models: Prandtl-Tomlinson (PT), Hébraud-Lequeux (HL) and the universe of 2 d 𝑑 d EPMs.

Rate type	PT	HL	2$d$-EPM
	(H=1)	(H=1/2)	(H=2/3)
Uniform	1	2	3/2
Progressive	3/2	5/2	2

Table 2. Table 2: Values of critical exponents determined from numerical simulations in two-dimensional systems. Exponents in general do not depend on the particular model analyzed. β 𝛽 \beta and z 𝑧 z differ for uniform and progressive rate rules.

	uniform	progressive
$\beta$	$1.5\pm 0.1$	$2.0\pm 0.1$
$\vartheta$ (from $\left(x_{\mathrm{𝚖𝚒𝚗}}\right)$)	$0.5\pm 0.05$	$0.5\pm 0.05$
$\theta$ (from $P\left(x\right)$)	$0.75\pm 0.07$	$0.75\pm 0.07$
${\tau }_S$	$1.33\pm 0.03$	$1.33\pm 0.03$
$d_f$	$1.08\pm 0.05$	$1.08\pm 0.05$
$H$	$0.67\pm 0.03$	$0.67\pm 0.03$
$z/d_f$	$0.54\pm 0.02$	$0.43\pm 0.04$
$z$	$\sim 0.583$	$\sim 0.464$

Equations

$$\frac{\partial {\sigma }_i(t)}{\partial t}=\mu {\dot{\gamma }}^{\mathtt{ext}}+\sum _j{G}_{ij}{n}_j(t)\frac{{\sigma }_j(t)}{\tau };$$

$$\frac{\partial {\sigma }_i(t)}{\partial t}=\mu {\dot{\gamma }}^{\mathtt{ext}}-g_0{n}_i(t)\frac{{\sigma }_i(t)}{\tau }+\sum _{j\ne i}{G}_{ij}{n}_j(t)\frac{{\sigma }_j(t)}{\tau };$$

$$G_{\mathbf{q}}=-\frac{4q_x^2{q}_y^2}{(q_x^2+q_y^2{)}^2}.$$

$$G_{\mathbf{q}=0}=-\kappa$$

$$q_{x,y}^2\equiv 2-2\cos \left(\frac{\pi m_{x,y}}{L}\right)$$

$$n_i:\{\begin{array}{ll}0\to 1& \text{mbox}atrate\text{leavevmode}{\tau }_{\mathtt{on}}^{-1}\text{mbox}if{\sigma }_i>{\sigma }_{\mathtt{Y}}\\ 0\leftarrow 1& \text{mbox}atrate\text{leavevmode}{\tau }_{\mathtt{off}}^{-1}\end{array}$$

$$n_i:\{\begin{array}{ll}0\to 1& \text{mbox}atrate\text{leavevmode}{\tau }_{\mathtt{on}}^{-1}\text{mbox}if{\sigma }_i>{\sigma }_{\mathtt{Y}}\\ 0\leftarrow 1& \text{mbox}instantaneously\end{array}$$

$$n_i:\{\begin{array}{ll}0\to 1& \text{mbox}instantaneouslyif{\sigma }_i>{\sigma }_{\mathtt{Y}}\\ 0\leftarrow 1& \text{mbox}when\int d{t}^{\prime }|{\partial }_t\sigma (t^{\prime })/\mu +{\dot{\gamma }}^{pl}(t^{\prime })|\ge {\gamma }_c\end{array}$$

$$\lambda \sim ({\sigma }_i-{\sigma }_i^{th}{)}^{1/2}.$$

$$\frac{\partial \overline{\sigma }(t)}{\partial t}=\mu {\dot{\gamma }}^{\mathtt{ext}}-\frac{\kappa }{\tau }\overline{n\sigma }(t);$$

$$⟨x_{\mathtt{min}}⟩\propto N^{-\varphi }\text{leavevmode},\text{leavevmode}\text{ with }\text{leavevmode}1<\varphi <2$$

$$⟨x_{\mathtt{min}}⟩\sim N^{-1/(1+\theta )}=L^{-d/(1+\theta )}.$$

$$\xi (kx)\sim k^H\xi (x)$$

$$P(\delta \xi )\sim \frac{1}{|\delta \xi {|}^{\mu +1}}$$

$$T\sim N^{\alpha }{S}^{(1-\alpha )/2}$$

$$\dot{\gamma }\propto (⟨\sigma (\dot{\gamma })⟩-{\sigma }_c{)}^{\beta },$$

$${\tau }_S=2-\frac{\vartheta }{\vartheta +1}\frac{d}{d_f}.$$

$$P(\delta \xi )=\frac{A}{N}{\left|\delta \xi \right|}^{-\mu -1}$$

$$H=\frac{1}{\mu }$$

$$\theta =\mu /2$$

$$\frac{1}{H}=\beta -1+\frac{1}{\omega }.$$

$$\dot{\gamma }\sim (\sigma -{\sigma }_c{)}^{\eta +2}$$

$${\dot{e}}_i=-d{V}_i/e_i+\sum _j{G}_{ij}{e}_j+\sigma .$$

$${\dot{e}}_i=-d{V}_i/e_i+{\eta }_i(t)+\sigma .$$

$$\beta =1+\frac{1}{H}-\frac{1}{\omega }.$$

$${\dot{e}}_i\sim -e_i^2-\delta \sigma$$

$$\begin{array}{rl}{\partial }_t\mathcal{P}(\sigma ,t)=& -G_0\dot{\gamma }(t){\partial }_{\sigma }\mathcal{P}+D_{\mathtt{pl}}(t){\partial }_{\sigma }^2\mathcal{P}\\ & -\nu (\sigma ,{\sigma }_c)\mathcal{P}+\Gamma (t)\delta (\sigma )\end{array}$$

$${\nu }_{\mathtt{HL}}(\sigma ,{\sigma }_c)\equiv \frac{1}{\tau }\Theta (\sigma -{\sigma }_c)$$

$$\Gamma (t)=\frac{1}{\tau }{\int }_{{\sigma }^{\prime }>{\sigma }_c}d{\sigma }^{\prime }\mathcal{P}({\sigma }^{\prime },t),$$

$$D_{\mathtt{pl}}(t)=\tilde{\alpha }\Gamma (t)$$

Full text

I Introduction

Last years have witnessed major advances in the understanding of the deformation of amorphous solids from the statistical physics point of view Karmakar et al. (2010); Lin et al. (2014). The so-called yielding transition, an out-of-equilibrium (driven) transition between a solid-like and a flowing phase has been analyzed making use of comparisons and analogies with other well known problems in the literature. In particular, the depinning transition of elastic manifolds in disordered media Fisher (1998) has largely influenced the interpretation of yielding Dahmen et al. (2009); Tyukodi et al. (2016); Lin et al. (2014). In analogy with the classical velocity-force characteristics of depinning $v\sim(f-f_{c})^{\beta}$ (see for example Ferrero et al. (2013)), yielding can be described by a critical behavior of the steady-state strain rate $\dot{\gamma}$, that is zero when the stress $\sigma$ is below a critical value $\sigma_{c}$ and becomes $\dot{\gamma}\sim(\sigma-\sigma_{c})^{\beta}$ when $\sigma>\sigma_{c}$. Additionally, in the transient regime, the onset of flow starting from a given glass has been shown to display two qualitatively different kinds of behavior according to the annealing level, separated by a critical point, akin to athermally driven random-field Ising models Ozawa et al. (2018).

In the seek of simplification and generalization, elastoplastic models (EPMs) built at a coarse-grained level constituted the workhorse in the development of a statistical mechanics description of yielding Nicolas et al. (2018). Nevertheless, the highly phenomenological approach on which these models have been built supplied a broad variety of rules and details, yielding different quantitative results and obfuscating the establishment of a set of universal exponents for the yielding transition. Despite a broad recent activity on the study of yielding by means of EPMs, universal properties remain elusive. Although some consensus has been built in the numerical community around the avalanche statistics displayed being different from mean-field depinning Nicolas et al. (2018), quantitatively the reported critical exponents still differ. In particular, exponents such as the ones governing the flowcurve and the relation between avalanche size and duration show a wayward behavior. On the other hand, a comparison of EPMs with well known mean-field constructions and biased-random-walk problems, lead to a clearer expectation on where one would find universality and where model details should matter Jagla (2017). More importantly, experiments of yield-stress materials show themselves a broad variation of exponent values, e.g., $n\simeq 0.2-0.8$ for the Herschel-Bulkley exponent of the flowcurve Bonn et al. (2017) ($n\equiv\beta^{-1}$). While this dispersion may be originated in a number of reasons, including a variety of experimental and measurement protocol details, detailed predictions from a theoretical perspective, though for an idealized case, could be illuminating. More than an academic exercise, the forge of consensus around the critical properties of a problem is crucial in the practical development of the field. It is with that spirit that we address this issue.

In this work we contrast and analyze outputs from three different elastoplastic models previously used in the literature and further modify them to illustrate other three model cases. In particular, we focus on the modification of the rule that governs the onset of the increase of local plastic deformation, i.e., the local yielding. In general, when the local stress overcomes a preset threshold ($\sigma_{i}>\sigma_{i}^{th}$) there is a probability $\lambda$ per unit time that a plastic deformation occurs. In all versions of EPMs presented so far, the value of $\lambda$ is taken as a constant. However, another alternative seems more natural. The plastic strain increase when the local stress threshold is overcome can be described as the passage between a local state that becomes unstable to a new stable state. Therefore, the typical time needed to move to the new minimum, and equivalently the transition rate $\lambda$, should be a function of the “degree of instability” $\sigma_{i}-\sigma_{i}^{th}$. We will refer to this situation as diplaying “progressive rates”, as opposed to the constant, or uniform case. When this effect is quantitatively taken into account, the value of the inverse Herschel-Bulkley exponent $\beta$ is seen to differ in $\frac{1}{2}$ between the two cases.

We make a somewhat ad hoc classification of exponents governing scaling laws as static and dynamical exponents, relegating to the latter class those exponents that happen to depend on the particular rules of the model (constant or progressive rates), as it is the case for $z$ and $\beta$. Furthermore, we support our numerical observation by estimating analytically the impact of the “progressive rates” modification in the paradigmatic Hébraud-Lequeux model. Our results highlight an important ingredient to be considered in the interpretation of the broad range of experimental values observed for the flowcurve exponent of the yielding phenomenon.

In the athermal and overdamped limit in which the yielding transition is defined, the quasistatic dynamics is governed by a priori uncorrelated collections of plastic events tagged avalanches. Events shaping an avalanche, in turn, are supposed to be correlated, giving rise to a non-trivial dynamics depending on interaction kernel and dimensionality. Yet, standing on a distant point of the system and assuming ergodicity, all the physics could be described by the mechanical noise felt by this point due to the avalanches taking place elsewhere. By analyzing the time series of the mechanical noise accumulated on time at a distant location, we can interpret our critical exponents from the problem of a biased random walk; e.g., giving an explicit expression for $\beta$ in terms of the Hurst exponent of the accumulated mechanical noise signal, which constitutes in itself an alternative mean-field proposal for the study of yielding.

II Models and simulation protocols

We consider amorphous materials at a coarse-grained-level description, laying in between the particle-based simulations and the continuum-level description provided by classical approaches such as Soft Glass Rheology Sollich (1998) or the Hèbraud-Lequeux mean-field model Hébraud and Lequeux (1998). Full background, context and historical development of these so-called elasto-plastic models (EPMs) can be found in a recent review article Nicolas et al. (2018). Briefly, the amorphous solid is represented by a coarse-grained scalar stress field $\sigma({\bm{r}},t)$, at spatial position ${\bm{r}}$ and time $t$ under an externaly applied shear strain. Space is discretized in blocks (e.g., square lattice). At a given time, each block can be “inactive” or “active” (i.e., yielding). This state is defined by the value of an additional variable: $n({\bm{r}},t)=0$ (inactive), or $n({\bm{r}},t)=1$ (active). An over-damped dynamics is imposed for the stress on each block, following some basic rules: (i) The stress loads locally in an elastic manner while the block is inactive. (ii) When the local stress overcomes a local yield stress, a plastic event occurs with a given probability, and the block becomes “active” ($n({\bm{r}})$ is set to one). Upon activation, dissipation occurs locally, and this is expressed as a progressive drop of the local stress, together with a redistribution of the stresses in the rest of the system in the form of a long-range elastic perturbation. A block ceases to be active when a prescribed criterion is met. The auxiliary binary field $n({\bm{r}},t)$ shows up in the equation of motion for the local stress $\sigma({\bm{r}},t)$, defining a dynamics that is typically non-Markovian. While the structure of the equation of motion for the local stresses is almost unique in the literature, both its parameters and the rules governing the transitions of $n({\bm{r}})$ ($0\rightleftharpoons 1$) show a variety of choices.

We define our EPMs as a $2$-dimensional scalar field $\sigma({\bm{r}},t)$, with ${\bm{r}}$ discretized on a square lattice and each block $\sigma_{i}$ subject to the following evolution in real space

[TABLE]

where $\dot{\gamma}^{\tt ext}$ is the externally applied strain rate, and the kernel $G_{ij}$ is the Eshelby stress propagator Picard et al. (2004).

It is sometimes convenient to explicitly separate the $i=j$ term in the previous sum, as

[TABLE]

where $g_{0}\equiv-G_{ii}>0$ (no sum) sets the local stress dissipation rate for an active site. The form of $G$ is $G({\bm{r}},{\bm{r}}^{\prime})\equiv G(r,\varphi)\sim\frac{1}{\pi r^{2}}\cos(4\varphi)$ in polar coordinates, where $\varphi\equiv\arccos(({\bm{r}}-{\bm{r}}^{\prime})\cdot{\bm{r}}_{\dot{\gamma}^{\tt(ext)}})$ and $r\equiv\left|{\bm{r}}-{\bm{r}}^{\prime}\right|$. For our simulations we obtain $G_{ij}$ from the values of the propagator in Fourier space $G_{\bf q}$, defined as

[TABLE]

for $\bf q\neq 0$ and

[TABLE]

with $\kappa$ a numerical constant (see below). Note that in our square numerical mesh of size $L\times L$, $q_{x}^{2}$, $q_{y}^{2}$ must be understood as

[TABLE]

with $m_{x,y}=0,...,L-1$.

The elastic (e.g. shear) modulus $\mu=1$ defines the stress unit, and the mechanical relaxation time $\tau=1$, the time unit of the problem. The last term of (2) constitutes a mechanical noise acting on $\sigma_{i}$ due to the instantaneous integrated plastic activity over all other blocks ($j\neq i$) in the system. The picture is completed by a dynamical law for the local state variable $n_{i}=\{0,1\}$. We define hereafter three different rules corresponding to three different models:

Picard’s model Picard et al. (2005)

[TABLE]

where $\tau_{\tt on}$ and $\tau_{\tt off}$ are parameters and $P({\sigma_{\tt Y}})=\delta({\sigma_{\tt Y}}-1)$.

• 2.

Lin’s model Lin et al. (2014)

[TABLE]

where $\tau_{\tt on}=1$ and $P({\sigma_{\tt Y}})=\delta({\sigma_{\tt Y}}-1)$. The plastic stress release is instantaneous and the block becomes inactive immediately after being activated, all $j\neq i$ neighbors receive a kick that is proportional to the value of the local stress drop and mediated by the Eshelby propagator. The local stress drop value is chosen to be $\sigma_{i}\pm\epsilon$, with $\sigma_{i}$ the local stress value just before yielding and $\epsilon$ a uniformly distributed random variable with amplitude $\epsilon_{0}=0.1$, to avoid periodic dynamics effects.

• 3.

Nicolas’ model Nicolas et al. (2014)

[TABLE]

here $P({\sigma_{\tt Y}})$ is exponentially distributed (as in Liu et al. (2016)) and the integral over $dt^{\prime}$ accounts for the accumulated plastic deformation of the block after the last local yielding.

Besides their differences, one thing that these EPMs as-found-in-the-literature have in common is that they have a fixed or constant transition rate for plastic activation $\lambda_{\tt fix}=\tau_{\tt on}^{-1}$, be it finite or arbitrarily large as in Nicolas’ model. As already mentioned, a more natural alternative would be to associate a stress-dependent typical time with the passage between stable states. This is seen more convincingly by considering models with continuous local disorder potentials. This analysis (postponed to Section IV.4) reveals that, for a smooth form of the effective local confining potential, $\tau_{\tt on}\sim(\sigma_{i}-\sigma_{i}^{th})^{-1/2}$. If we want to maintain an implementation in terms of transition rates, $\lambda\equiv\tau_{\tt on}^{-1}$ should be then expressed as

[TABLE]

Notice that a situation of constant rate is recovered if the local disorder potential is assumed to produce a jump in the force at the transition point, a situation that occurs, for instance, if the potential is formed by a concatenation of parabolas. In this case, the time it takes for the local stress to reach the new minimum when $\sigma_{i}>\sigma_{i}^{th}$ is roughly constant, independent of the degree of instability. We can then associate the prescription of a constant transition rate in previous EPMs to a local disorder potential formed by the concatenation of parabolic pieces.

Therefore, apart form the “classical” implementations of the three models described above, we further include in our present analysis the same three models but modified with stress-dependent transition rates of the kind $\tau_{\tt on}^{-1}=\lambda_{\tt prog}=(\sigma_{i}-{\sigma_{\tt Y}}_{i})^{1/2}$. Other choices of $\eta$ for $\lambda_{\tt prog}\propto(\sigma_{i}-{\sigma_{\tt Y}}_{i})^{\eta}$ could be possible. Yet, in the analogy, the disordered potential that would produce something different from $\eta=1/2$ or $\eta=0$ turns out to be very difficult to justify (see Section IV.4). Anyhow, for $\eta>0$ the intuitive expectation is that local regions that have exceeded their stress threshold by larger amounts will take precedence in their fluidization moment over others. As we will see, the change from $\eta=0$ to $\eta=1/2$ largely impacts the rheological properties of EPMs.

II.1 Finite strain-rate protocol

Starting from an initial configuration $\{\sigma_{i}\}$ that observes mechanical stability the system is evolved according to the equation of motion (2) with an elementary discretized time step $dt=10^{-2}$ (or smaller). After each time step integration of $\sigma_{i}$, the $n_{i}$ are updated according to the corresponding rules (6, 7 or 8, or one of their ‘progressive rate’ equivalents). The calculation of the convolution between $n_{j}\sigma_{j}$ and $G_{ij}$ is done in Fourier space. If the ${\bf q}=0$ mode of the propagator is set to zero, i.e. $\kappa=0$ (Eq. 4), then $\sum_{i}G_{ij}=0$, and the dynamics is stress conserving. However, in order to be able to control the strain rate $\dot{\gamma}^{\tt ext}$ we are interested to work with non-conserved stress protocols. This is accomplished by taking $\kappa>0$. By spatially averaging Eq. 1 we get

[TABLE]

that clearly indicates how $\kappa$ produces a stress decay in the system as long as there are active sites ($n\neq 0$). We use $\kappa=1$, as in previous strain-controlled EPMs implementations Martens et al. (2012); Nicolas et al. (2014); Liu et al. (2016), unless otherwise specified.

It has to be emphasized that initial conditions play no role in our study, as all measurements are done in the steady state, after a very long straining stage where the memory of the initial condition is totally lost. The transient where we don’t take any measurements is sufficiently long so that every block has yielded several times already. For the flowcurves construction each fix strain-rate has its own long transient strain period.

II.2 Quasistatic protocol

For the analysis of avalanche statistics, it is convenient to have a protocol that allows for the triggering and unperturbed evolution (no driving) of avalanches until they stop (what is guaranteed by a degree of stress non-conservation $\kappa>0$). This is the quasi-static protocol described here.

Starting from any stable configuration, i.e., no site is active and no site stress is above its local threshold ($n_{i}=0$ and $\sigma_{i}<{\sigma_{\tt Y}}_{i}$ for all sites), the next avalanche of plastic activity is triggered by globally increasing the stress by the minimum amount necessary for a site to reach its local threshold. That site (the weakest) is activated at threshold with no stochastic delays; it perturbs the stress values of other sites and the rest of the avalanche evolves without any external drive following the dynamics prescribed by Eq. (2) (and the corresponding activation rule) with ${\dot{\gamma}}=0$. The avalanche stops once there are no more active sites and all stresses are below their corresponding thresholds again. At this point the loading process is repeated. For each simulation run, data is collected only in the steady-state.

We have noticed that Nicolas’ model as originally proposed Nicolas et al. (2014) is ill-defined in the quasistatic limit. Eventually one arrives to the pathological case in which a site is active and the criterion to recover local elasticity (8) is never met (any finite ${\dot{\gamma}}$ guarantees it, but not the quasistatic protocol). To overcome this issue, we add a maximum duration bound to plastic events, yet large enough to guarantee the full relaxation of the local stress.

III Results

III.1 Quasistatic avalanche distributions

Starting from a system that has been deformed at a slow strain rate, we apply a quasistatic protocol. Size $S$ and duration $T$ of each triggered avalanche are measured. $S$ is simply calculated from the total stress drop $\Delta\sigma$ caused by the avalanche of plastic events as $S=\Delta\sigma L^{d}$, while $T$ is the time elapsed until the avalanche ceases its activity, measured in units of $dt$.

Figure 1 shows avalanche size distributions for all three models in their two transition rate variants, at various system sizes. All distributions have the form $P(S)\sim S^{-\tau_{S}}f(S/S_{\tt max})$, with $f(x)$ a rapidly decaying function (a compressed exponential). Taking a single EPM no difference between uniform and progressive rates is observed. Furthermore, the same value of $\tau_{S}\simeq 1.33$ characterizes all avalanche distributions, irrespective of the six model variants analyzed. Good collapse of the distributions for different system sizes is found when plotting $P(S)L^{\tau_{S}d_{f}}$ vs $S/L^{d_{f}}$ with the so-called “fractal dimension” $d_{f}=1.08$. However, explicit and acceptable fittings of $P(S)=AS^{-\tau_{S}}\exp(-(S/S_{\tt max})^{b})$ with $b\simeq 1.5\ldots 2$ may result in $S_{\tt max}\propto L^{d_{f}}$ with $d_{f}\simeq 1.1\ldots 1.5$ as well. In any case, we notice that $d_{f}$ is independent on the particular form of the rates used.

Apart from the uncertaintes in the determination of $d_{f}$, a word should be said about the accurate quantitative estimation of $\tau_{S}$. We have noticed that its estimation can be sensitive to some technicalities. The avalanche sizes distribution shape changes with the parameter of stress (non-)conservation $\kappa$ (the value chosen for $-G_{{\bf q}=0}$). When $\kappa=0$, stress is conserved by the dynamics. In that limit, we will soon or late observe a never ending avalanche once the loading takes the system above the critical threshold. For small but finite values of $\kappa$ any stress excess eventually decays, but large avalanches are difficult to stop and this results in a plateau or bump in the $P(S)$ distribution at large $S$. For relatively small systems sizes, as the ones one can usually simulate, an important part of the distribution is impacted by this protocol detail (the value of $\kappa$). Fig. 2 illustrates (for Lin’s model with uniform rates) how the choice of $\kappa$ affects the quantitative estimation of $\tau_{S}$ in a typical finite-size avalanche size distribution. If one chooses a power-law window in $P(S)$ neglecting very small and very large $S$ (usually associated with finite time-step $dt$ discretization and finite system size $N$, respectively), we would observe $\tau_{S}$ varying from $\tau_{S}\simeq 1.33$ for $\kappa=1$ to $\tau_{S}\simeq 1.25$ for $\kappa=0.01$. Still, $\tau_{S}$ should be independent of the precise value of $\kappa$ in the thermodynamic limit. The particular value of $\kappa$ has an effect on $P(S)$ in the range of $S$ for which the last term in Eq. (10) gives a total stress drop comparable with the present value of the stress at the yielding sites times $\kappa$. For an avalanche of size $S$ in a system of size $N\equiv L^{d}$ this occurs when $\Delta\sigma L^{d}\sim{\sigma_{\tt Y}}L^{d_{f}}\kappa$ or larger. This means that for any $\kappa$ we should have a consistent value of $\tau_{S}$ as long as we fit $P(S)$ for an avalanche range not reaching this limit. For the system in Fig. 2, this is $S/\kappa\sim 512^{d_{f}}$. On the left of such a value, all curves seems consistent with the power-law found in Fig.1 (red lines). This is to argue that, although different exponent values are found in the literature Talamali et al. (2011); Budrikis et al. (2017); Liu et al. (2016); Lin et al. (2014); Nicolas et al. (2018) and we cannot be conclusive about $\tau_{S}\simeq 1.33$ being the universal value, the sensitivity of $P(S)$ to protocols with effectively different values of $\kappa$ may justify the spread of reported values of $\tau$. Systematic finite-size analysis are needed to reach an agreement on a unique universal exponent. We have done a first step here by clearly demonstrating that, at least when fixing $\kappa$, $\tau_{S}$ and $d_{f}$ do not depend on the yielding rate rule, nor in other details defining the particular EP model.

III.2 Extremal statistics and density of shear transformations $P(x)$

Directly related to the avalanche size statistics during quasistatic simulations is the stress $x_{\tt min}$ needed to trigger the next avalanche each time ($x_{\tt min}={\sigma_{\tt Y}}_{w}-\sigma_{w}$, where $w$ stands for the “weakest” site). It was first observed by molecular-dynamics simulations Maloney and Lemaître (2004) and then contextualized in a framework of extreme statistics both in MD Karmakar et al. (2010) and EPMs Lin et al. (2014), that the finite size scaling of the average loading stress needed to trigger avalanches is sub-extensive, i.e.,

[TABLE]

with $N=L^{d}$. This is at odds with the intuition that we get from analogies with sand-pile models or other stick-slip dynamics as the one related to the depinning transition. There, the intensive local variables (pile height or force at site $i$, equivalent to the stresses $\sigma_{i}$) can only increase as the intra-avalanche dynamics proceeds and therefore a larger system will always show an equally weaker site on average ($\phi=1$).

This phenomenological sub-extensiveness was interpreted as a consequence of marginal stability in the driven amorphous solid by Lin et al.Lin et al. (2014). In average, starting at any given steady state configuration, sites far away from their thresholds have a much larger life time before yielding than sites that are close to it. Since plastic activity provides signed kicks to each site stress, the problem of local yielding is that of a first passage time, and the probability of finding the walker away from the boundary is much larger than finding it very close to it. Analyzing the distribution $P(x)$ of local stresses needed to reach local thresholds ($x_{i}={\sigma_{\tt Y}}_{i}-\sigma_{i}$)–tagged the “density of shear transformations”–and assuming it to be characterized by a power-law $P(x)\sim x^{\theta}$ at vanishing $x$, it can be shown that the distribution of the minimal values ($x_{\tt min}$) among randomly chosen sets of $N\gg 1$ independent samples from $P(x)$ is a Weibull distribution with mean value

[TABLE]

In this way, any $\theta>0$ provides a sub-extensive $\left<x_{\tt min}\right>$, and in fact, $\theta>0$ was estimated Karmakar et al. (2010); Lin et al. (2014); Nicolas et al. (2018) for the yielding of amorphous solids in contrast with the $\theta=0$ value of the depinning of elastic manifolds. Fig. 3 shows $\left<x_{\tt min}\right>$ as a function of system size for different elastoplastic models obtained in the steady state of quasistatic simulations. We can observe that results for the uniform and progressive rates variants of the models are totally consistent. Furthermore, all the simulated EPMs display the same scaling law $<x_{\tt min}>\sim L^{-1.33}$, compatible with $\theta\simeq 0.5$ in Eq. (12).

At this point, it is appropriate to clarify a dichotomy in the determination of exponent $\theta$. It turns out that $P(x)$ does not keep a form $P(x)\sim x^{\theta}$ when $x\to 0$ as it is commonly assumed. Instead, it deviates from such power-law. In Fig. 4 we show the results obtained for $P(x)$, as determined from quasistatic simulations collecting data right after avalanches have finished and before a new loading phase begins, in the steady state. The deviation from a common power law, and the establishment of a system-size dependent plateau at $x\to 0$ is clear, as it is also the fact that this plateau occurs systematically at smaller values of $x$ as $L$ increases. Note that this effect can be easily overseen if an arbitrary lower cutoff of the histograms is imposed.

Therefore, the form of the distribution of stress distances to local threshold is rather $P(x)\simeq P(0)+x^{\theta}$. The emergence of this plateau at small $x$ disentangles $\theta$ from the exponent $\phi$ that rules the finite size scaling of $\langle x_{\tt min}\rangle$ in Eq.11. So we give a new name to the ‘apparent’ $\theta$ holding Eq. 12: $\vartheta=d/\phi-1$, while keeping $\theta$ for a behavior $P(x)\sim x^{\theta}$ in the intermediate range of $x$ where it holds. The inset of Fig. 4 shows an empirical scaling of the form $P(x)=P(0)+x^{\theta}$ with $P(0)\sim L^{-a}$. We obtain $a\simeq 0.55$. $\theta$ is in the range $0.6-0.77$ depending on the criterion chosen to identify the power-law regime. Such power-law range is usually limited by finite size effects. In our case, only for the biggest $L$ simulated we can say that $P(x)$ displays a power-law in more than a decade, with exponent $\theta\simeq 0.75$.

In Fig. 5 we plot $P(x)$ for different EP models with $L=1024$. Again, the values of $\theta$ that one could obtain from such curves vary between $0.56-0.77$ according to the criterion used to identify the power-law regime. Yet, in any case, the extracted exponent value is independent on each activation rate used, as the correspondent distributions for the same model lay one on top of the other. Beyond the numerical uncertainty, the values of $\theta$ observed are clearly above the value of $\vartheta=0.5$ extracted from the scaling of $\langle x_{\tt min}\rangle$. Furthermore, when lines proportional to $x^{\theta}$ with $\theta=0.75$ are proposed, all model variants show a good compatibility with them in a given range of $x$ (see Fig.5). This value of $\theta$ happens to be consistent with the exponents found for the flowcurve, as we discuss in Sec.IV.2.

The same global form of $P(x)$, with an apparent plateau at $x\to 0$, has been independently and simultaneously reported by Tyukodi et al. Tyukodi et al. (2019) and later confirmed in MD simulations Ruscher and Rottler (2019) while this work was under review. We further discuss this in Sec. IV.1.

III.3 Accumulated mechanical noise

Standing on an arbitrary site $i$ in the system we can define $\xi_{i}(t)$ as the accumulated mechanical noise on $i$ at time $t$ due to avalanches occurring elsewhere. We would like to study the statistical properties of this noise and relate it to the observed values of the avalanche statistics and the rheological exponent $\beta$. In particular, we are interested in describing $\xi(t)$ as a correlated noise with statistical properties defined by the so-called Hurst exponent $H$, so that $\xi(x)$ is a (stochastic) homogeneous function of degree $H$,

[TABLE]

Note that a standard random walk has $H=1/2$. In general, in a mean-field description Lin and Wyart (2016), $H=1/\mu$ where each contribution $\delta\xi$ to the accumulated noise $\xi$ is assumed to come from a distribution with a long tail

[TABLE]

and with $1\leq\mu\leq 2$.

Taking into account the action of individual plastic events through the Eshelby propagator, Lin and Wyart (2016) arrived to the conclusion that $\mu=1$ is the only value with “physical” sense. Nevertheless, if we discretize time in such a way that we can see the occurrence of avalanches as the elementary event contributing to the mechanical noise, other values of $\mu$ may acquire a physical sense. In fact, a distant site at a given moment in the relevant time scale does not feel the mechanical kick produced by a single site, but by a full avalanche of sites yielding. A coarse-graining in time is done by the system itself, which does not necessary constraint the plastic events to occur one at a time and well separated but overlaps them in a sudden burst of activity. In a quasistatic loading protocol, the time scale separation between the duration of one avalanche and the loading time to trigger the next one allows us to make this approximation and set the time-scale in avalanche counts.

Figure 6 displays different realizations of $\xi$ measured for various of our EPMs. We can already notice with the bare eye the characteristic long steps or Lèvy flights that distinguish this kind of signal from standard random-walks. More importantly, we will show that in all cases these signals show values of $H$ that are very close to each other, pointing to a robust value of $H$ that defines an intrinsic property of the system. To compute $H$ we can analyze the signal in the following way: for each window of size $\epsilon$ in the avalanche sequence compute the average value of the absolute noise difference in that window $\left<\delta\right>(\epsilon)=\frac{1}{M}\sum_{i}^{M}|\xi(i)-\xi(i+\epsilon)|$. $H$ is then the exponent of the relation $\left<\delta\right>\sim\epsilon^{H}$. This is basically what the Detrended Fluctuation Analysis (DFA) (Nolds Schölzel (2018)) does after having subtracted from the signal a global trend. Therefore, we collect large series of the accumulated mechanical noise at different points in the system and analyze these signals with the DFA algorithm.

Results for the histograms of $H$ values observed are shown in Fig. 7 for the different models and rate rule variants. We observe that all systems show some spread around a mean value $\sim 0.67$, but already the concurrence of all model variants around such a particular value is non-trivial. The DFA analysis is sensitive to details as the total signal length and the minimum segment length for polynomial fitting, installing a non negligible uncertainty in its output. We believe that a good theoretical proxy for the Hurst exponent in 2d-EPMs is $H\simeq\frac{2}{3}$, signaled by a vertical dashed line in Fig. 7. This corresponds to $\mu=\frac{1}{H}\simeq 1.5$ in Eq. 14, subscribing the idea that a coarse-graining in time can give a physical meaning to values of $\mu$ different from $\mu=1$. The Hurst exponent is independent of the rate rule, what puts it at the level of other static critical exponents such as $\tau_{S}$, $\theta$ or $d_{f}$. In the Section IV.2 below we provide arguments to believe that the exponent $H$ characterizing the mechanical noise in a quasistatic measure is useful for the estimation of the $\beta$ exponent of the flowcurve as it departs from the critical yielding point.

III.4 Avalanche durations

We now consider a critical exponent for which a dependence on the rate law may be expected. This is the dynamical exponent $z$ that, combined with $d_{f}$, relates the avalanche size with the avalanche duration. The avalanche duration $T$ is expected to scale as a power law of the avalanche size $S$, namely $T\sim S^{1/\delta}$. If we assume that both $S$ and $T$ are controlled by a correlation length $\ell$ through $S\sim\ell^{d_{f}}$ and $T\sim\ell^{z}$, then $T\sim S^{z/d_{f}}$ and $\delta=d_{f}/z$.

Fig. 8(top) shows the $T$-$S$ relation for a large set of avalanches of quasistatic simulations at a fix system size ($L=1024$) from Lin’s model with both constant and progressive rates. We can observe for both cases a broad cloud of points, clearly different for the two different rate rules. Although scaling relations are difficult to guess from these clouds, it is clear that the progressive rate allows for a much larger spread of avalanche durations for the same avalanche size, in particular, it allows for avalanches much longer in time. Now, averaging the values of $T$ within small $S$ intervals, for different system sizes we obtain the curves in Fig. 8(bottom). For uniform rates, the results show a consistent power law with an exponent $z/d_{f}\simeq 0.54$. Data for all system sizes overlap on the same curve. Using the previously determined value $d_{f}\simeq 1.08$, we can obtain $z\simeq 0.58$. For progressive rates, the results are definitely different. At a fix system size, a grow of $\left<T\right>$ vs $S$ with exponent $z/d_{f}\simeq 0.43$ is established for large enough avalanches, yielding for the dynamical exponent $z\simeq 0.46$, i.e., lower than the uniform case. Moreover, an additional feature is observed. The results for increasingly large system sizes do not simply extend the region in which a power law is observed but also produce a shift of the average duration towards larger values.

This unexpected result can be rationalized in the following way. Each time an avalanche triggers a new plastic event, its life is expanded and some time will be added to its final duration. For uniform activation rates, this time is on average $dT=\lambda^{-1}=1$. On the other hand, for progressive rates, the time added to the avalanche duration depends on the stress excess of the new yielding site; on average $dT=\lambda^{-1}=(\sigma_{i}-{\sigma_{\tt Y}}_{i})^{-\frac{1}{2}}$. The events that most contribute to the avalanche duration then, are the ones that occur very close above their threshold. While their individual probability of yielding is low, the observation of these events increases with the number of sites susceptible of being in that situation; this is, it increases with increasing system size. This previously unnoticed phenomenon is also present at a mean-field level, where a more quantitative estimation of the value of $z$ can be given and an exact scaling law derived (see Appendix A.2), leading to the $N$ and $S$ dependence of the avalanche duration:

[TABLE]

We have tested this scaling in the inset of Fig. 8(bottom), and it works very well with $\alpha\simeq 0.16$. Notice that in the numerical determination of the durations $T$ reported in this section, the time needed to destabilize the first site in the avalanche was not considered. This was done in this way since, as the first site is destabilized by an infinitesimal quantity, it would add a diverging contribution to the total time when progressive rates are at play, completely spoiling the duration estimation. For uniform rates the first site only adds an additional time unit, but for consistency we do not considered it either in this case.

III.5 Flowcurves

The results of the previous section about the dependence of $z$ on the activation protocol is not surprising as $z$ is the prototypical dynamical exponent of the model. Much more unexpected is the fact that the flowcurve exponent $\beta$ is also dependent on the activation rate law, as we will discuss now. To build flowcurves we use the strain rate controlled simulation protocol. We first deform with a large strain rate (e.g. ${\dot{\gamma}}\simeq 0.5$) until a global steady stress value is reached, and then average the measured stress in time. After a measurement at a given strain rate, the strain rate value is reduced and the protocol is repeated, covering a proper range in strain rates to draw the flowcurve in log-scale.

When arriving from the “fluid” side the approach to the yielding point is continuous and the transition is critical. As the strain rate ${\dot{\gamma}}$ vanishes, the average stress reaches a finite value that we call the (dynamical) yield stress $\sigma_{c}$. In the vicinity of the limit ${\dot{\gamma}}\to 0$ the flowcurve obtained by a strain rate control protocol $\left<\sigma({\dot{\gamma}})\right>=\sigma_{c}+A{\dot{\gamma}}^{n}$ can be also written as

[TABLE]

where $\beta=1/n$.

Fig. 9 displays flowcurves for all three EPMs with both constant and progressive transition rates for local yielding. While the upper panel displays average stress vs. strain rate, in the lower panel we plot the strain rate as a function of the stress excess above the dynamical yield stress. One important detail in the construction of this curves is the appropriate estimation of $\sigma_{c}$, for which we have implemented separately quasistatic measurements and averaged the stress in this limit over long deformation windows. Although it is a delicate fitting procedure, not free from finite-size effects, once $\sigma_{c}$ has been properly determined, all EPMs with uniform rate show a flowcurve exponent compatible with $\beta=3/2$. as for the models with piece-wise parabolic disorder potential analysed in  Fernández Aguirre and Jagla (2018). Interestingly, when progressive rates are used instead, the exponent changes to $\beta=2$ in all the three models taken from the literature111 The progressive case of the Nicolas’s model in Fig. 9 shows a rather narrow range where $\beta=2$ seems to be valid. We believe this is due to the relatively large time the variable $n_{i}$ can remain in the ”active” (i.e. $n_{i}=1$) state in this case, masking the effect of progressive rates unless $\dot{\gamma}$ is very small, suggesting that a local fluidization rate of the form $\lambda_{\tt prog}=(\sigma-{\sigma_{\tt Y}})^{\frac{1}{2}}$ takes us to the case of ‘smooth’ disorder potentials Fernández Aguirre and Jagla (2018) (see Section IV.4). In any case, Fig.9 clearly shows that the particular dynamical rule that is used for local yielding has a direct impact in the value of the flowcurve exponent.

It can be noticed in passing that, by using progressive rates, the exponent obtained in our 2$d$-systems is $\beta=2$, as the one derived for the mean-field Hébraud-Lequeux (HL) model Hébraud and Lequeux (1998); Agoritsas et al. (2015); Agoritsas and Martens (2017). Nevertheless, in the HL model the activation rate is constant. The reason for this accidental coincidence rests in the fact that the flowcurve exponent is not only determined by the transition rates, but also by the mechanical noise statistics (which is dimension-dependent), as discussed in the next Section.

IV Discussion

IV.1 The residual density of sites at the verge of yielding

In Sec.III.2 we have shown data for the distribution $P(x)$ of $x={\sigma_{\tt Y}}-\sigma$, the local distances to threshold, showing that it develops a ‘plateau’ as $x\to 0$. This feature has been ignored or widely overlooked in the literature so far, but has been put forward recently not only by a draft version of this work but also in two other independent an very recent preprints by Tyukodi et al. Tyukodi et al. (2019) (see also Tyukodi (2016)) and Ruscher&Rottler Ruscher and Rottler (2019). Despite differences in the modeling approach, in particular the treatment of periodic boundary conditions, our results and the results of Ref. Tyukodi et al. (2019) are quite consistent. We both observe for quasistatic simulations in the steady state that a plateau develops for the distributions $P(x)$ for different system sizes at small $x$ and that the $\langle x_{\tt min}\rangle$ values lie between the plateau and power-law regimes (see crosses in Fig.4, an idea borrowed from Ref. Tyukodi et al. (2019)). They observe a scaling of the plateau level $p_{0}\equiv P(x=0)$ with system size of the form $p_{0}\sim L^{-a}$ with $a=0.6$, while we report $a\simeq 0.55$ in Sec.III.2. Also, the $\phi$ exponents fitted are quite close, 1.35 and 1.33 in Ref. Tyukodi et al. (2019) and in our case, respectively. The small discrepancy in the value proposed for $\theta$ (their 0.67 vs. our 0.75) could be explained by the range of system sizes explored. In MD simulations of 2D binary Lennard-Jones mixtures and using a ‘frozen matrix method’ the authors of Ref. Ruscher and Rottler (2019) find a scaling analogous to the one of Fig.4 with $a\simeq 0.65$ for the stationary regime (they also analyze the MD ‘as-quenched’ configurations prior to deformation), and a value of $\theta\simeq 1.05$ even larger than ours. They suggest the development of correlations in a preferential direction at the microscopic level to justify a ‘depinning-like’ form of $P(x\to 0)$. We have explored several spatial configurations of $P(x)$ as-left after an avalanche and we cannot hypothesize about the presence of correlations. Rather than thinking that way, we interpret the plateau effect at small $x$ as a consequence of working with a discrete time dynamics and a finite system size. For a discrete time random walk in the presence of an absorbing border we know that the probability amplitude right at the border is finite and its value decreases with the amplitude of the time step. Yet, we have performed simulations with time steps down to $dt=10^{-4}$ and the modification of the obtained $\langle x_{\tt min}\rangle$ values is very slow. The true origin of the finite size scaling for the $P(x)$ plateau at $x\to 0$ should have to do therefore not with the discretization in time of the individual events evolution, but with the discreteness of the smallest kick felt by any site in the system due to an avalanche happening elsewhere Ferrero and Jagla (2019).

It could be argued that being an effect that decreases with system size, the plateau in $P(x)$ for small $x$ must eventually become irrelevant: for sufficiently large system sizes the full $P(x)$ would acquire the expected pure $P(x)\sim x^{\theta}$ form, and the scaling in Eq. (12) (with $\theta$ instead of $\vartheta$) would hold. But this is not so. In Fig. 3 we see that $\langle x_{\tt min}\rangle$ is always in between the power-law and the plateau region of $P(x)$, and then it is mostly the scaling of this plateau with system size what determines the value of $\vartheta$ and the scaling in Eq.12. Thinking the other way around, we do not find any strong reason why the value of $\vartheta$ determined from Fig. 3 should be equal to that detrmined from the asymptotic form of $P(x)$ (the red dashed line in Fig. 4). In fact, these values do not coincide as far as we can tell: while from Fig. 3 we obtain $\vartheta\simeq 0.5$, the value from the asymptotic form in Fig. 4 is rather $\theta\gtrsim 0.7$. Actually, this situation is not new in the literature, where $\vartheta\sim 0.45-0.5$ has been determined from the scaling of $\langle x_{\tt min}\rangle$ but larger values ($\theta\simeq 0.6$) have been obtained from the full $P(x)$ in two-dimensional systems Karmakar et al. (2010); Lin et al. (2014); Nicolas et al. (2018). It should be always kept in mind that the exponent ruling the avalanche dynamics in the quasistatic limit is $\vartheta$, namely the one obtained from $\left<x_{\tt min}\right>$, since this is the average value of stress that is actually added to the system after each avalanche to maintain it in a stationary state. On the other hand, it is the value of $\theta$ (extracted from $P(x)$) the one that one can relate to the flowcurve exponent Lin and Wyart (2016), as we justify in the next section.

IV.2 Relations among exponents

In the light of our results, we recall and discuss in this section some scaling laws proposed in the literature. To start with, an interesting relation among exponents comes from a very simple argument originally presented by Lin et al. (2014). In the steady state, the stress has a well defined average and describes a jerky plateau as a function of strain for any finite system size. Therefore, on average, the stress increases in the loading phases must be balanced by the stress drops during avalanches. From the avalanches distributions $P(S)\sim S^{-\tau_{S}}f(S/L^{d_{f}})$, with $f(y)$ a fast decaying function for $y\gg 1$ and considering $S\equiv\Delta\sigma L^{d}$, one obtains $\langle\Delta\sigma\rangle\propto L^{d_{f}(2-\tau_{S})-d}$. On the other hand, the stress increment needed to trigger a new avalanche is what we have called $x_{\tt min}$ which scales as $\langle x_{\tt min}\rangle\propto L^{-\frac{d}{\vartheta+1}}$. Equating $\langle\Delta\sigma\rangle\propto\langle x_{\tt min}\rangle$, leads to

[TABLE]

By considering the values of $\vartheta\simeq 0.5$ (taken from Fig. 3), and $\tau_{S}\simeq 1.33$ and $d_{f}\simeq 1.08$ (from Fig. 1), we verify that this relation fairly holds for our data. This is not surprising since, as we explained, Eq. 17 is originated in a simple requirement of stress stationarity in the system. Note however that Eq. 17 is not satisfied using the $\theta$ value obtained from $P(x)\sim x^{\theta}$ in our data.

When the flowcurve exponent enter into the discussion, things are less clear. Different relations have been proposed in the literature, in particular Lin et al. (2014) $\beta=1+\frac{z}{d-d_{f}}$ and Lin and Wyart (2018) $\beta=1+\frac{1}{d-d_{f}}$, that show either a weak agreement with simulations or an asymptotic agreement in dimensions higher than $d=3$. We start our discussion by considering the exponent $\mu$ governing the fat tails of the mechanical ‘kicks’ $\delta\xi$ felt by a given site in the system, through the distribution Lin and Wyart (2016, 2018); Fernández Aguirre and Jagla (2018); Jagla (2018)

[TABLE]

with appropriate cutoffs. When integrated in time, this noise produces a time signal that can be characterized as a fractional Brownian motion with Hurst exponent $H$ that is given by

[TABLE]

In addition, it is also found analytically Lin and Wyart (2016, 2018) that such fractional Brownian motion generates a distribution $P(x)$ that behaves at low $x$ as $x^{\theta}$, with

[TABLE]

Lin et al. Lin and Wyart (2016, 2018) have suggested that $\mu$ is also the exponent $\beta$ of the flowcurve if it happens to be larger than one.

Our numerical results, together with those in Fernández Aguirre and Jagla (2018); Jagla (2018) (summarized in Section IV.4), strongly suggest that two dimensional yielding should be well described as a mean field transition characterized by a value of $H=2/3$. In fact, that leads to a consistent value of $\theta=3/4$ as the one we measured in Sec.III.2.

Furthermore, this generalized mean field analysis relates $H$ (and therefore $\theta$) to the flowcurve exponent $\beta$ through

[TABLE]

where, $\omega$ is related to the form of the transition rates: $\omega=1$ for uniform rates and $\omega=2$ progressive rates. Note that the prediction in Lin et al. Lin and Wyart (2016, 2018) (i.e. $\mu=\beta$) is consistent with this equation for the usual case of constant rates ($\omega=1$). In this case the equations above lead to $\theta=\beta/2$ 222Notice that in a previous publication by one of us Fernández Aguirre and Jagla (2018), the relation $\theta=\frac{1}{H}-1$ was proposed. However, the derivation of this relation was based on a wrong assumption, and therefore is not correct. , an extremely simple relation between the flowcurve exponent and the power exponent of the density of shear transformations, that we verify in our data (see Sec. III.5). It makes perfect sense: The steeper the $P(x)$, the sooner many sites will yield when we set the system in motion, the faster strain rate will increase and (equivalently) the slower the global stress will increase with $\dot{\gamma}$ due to the frequency of local fluidization. Notice that such a simple expression is also valid for the Hébraud-Lequeux model ($\beta=2$, $\theta=1$).

The values of $\beta$ according to formula (21) for different values of $H$ are summarized in Table 1. The case $H=1$ corresponds essentially to the absence of stochastic noise, and this is realized in the standard one-particle Prandtl-Tomlinson model. In this case the values of $\beta$ are well known. The value $H=2/3$ was argued before to correspond to the 2$d$-EPM, and we observe in fact that the values on Table 1 are exactly the flow exponents that were obtained from the numerical results in Fig. 9. The case of Gaussian noise $H=1/2$ with uniform rates is known to correspond to the Hébraud-Lequeux model. The value $\beta=2$ obtained from Eq. 21 in fact constitutes a well-known result for this model ($n=1/\beta=0.5$) Hébraud and Lequeux (1998); Agoritsas et al. (2015); Agoritsas and Martens (2017).

To complete the verification of the values contained in Table 1, which is to say, the accuracy of Eq. 21, it remains to consider the HL model for progressive rates. An analytical derivation of $\beta=5/2$ for that case is presented with some detail in the Appendix (A.1). In brief, by extending the standard rate of plastic events $\nu_{\tt HL}(\sigma,\sigma_{c})\equiv\frac{1}{\tau_{S}}\Theta(\sigma-\sigma_{c})$ to $\nu(\sigma,\sigma_{c})\equiv\frac{1}{\tau_{S}}(\sigma-\sigma_{c})^{\eta}\Theta(\sigma-\sigma_{c})$, so allowing for a progressive transition rate depending on the local stress excess, we follow the derivation of the flowcurve ending up with

[TABLE]

for general $\eta$. The case corresponding to smooth potentials has $\eta=\frac{1}{2}$, obtaining $\beta=5/2$. In Table 2 we summarize all the values of the exponents we have obtained in the numerical simulations of the present work.

IV.3 Brief digression about the dynamical critical exponent $z$

Since the instantaneous long-range Eshelby interaction is responsible for the dependence on microscopic properties and ultimately, for the different values of $\beta$ and $z$ found, one can argue (as Lin and Wyart (2018)) that spurious effects will disappear when this ‘unphysical’ instantaneous interaction is eliminated (by adding a mechanism for propagating waves, for example). While this might be the case, the use of an instantaneous interaction is not necessarily an unphysical assumption that should be improved, but rather makes sense in many different physical systems when restricted to the appropriate spatial and temporal scales. In the depinning of magnetic structures for instance, long-range interactions (typically dipolar) are considered, although these interactions cannot travel faster than the speed of light! Also, the problem of contact line depinning is known to have a long-range elastic interaction decaying as $1/r^{2}$. This interaction is generated through the elastic propagation of interactions in the bulk of the material, and in the end its propagation velocity must be limited to the speed of sound in the material. Nevertheless, that finiteness for the elastic propagation velocity has not been a major obstacle for the theory there. In our case, an Eshelby instantaneous interaction is an approximation that considers the speed of sound in the system to be infinite. This is clear for instance in the analysis of Cao et al. (2018), where in order to derive the Eshelby interaction an infinite bulk modulus –and thus an infinite sound velocity– is considered. This approximation produces some results that are certainly not correct in a more realistic situation. For instance a value of the dynamical exponent $z$ lower than one indicates a propagation of perturbations as $r=Ct^{1/z}$, and thus at velocities larger that any finite value if considering sufficiently large times or distances. The consideration of a finite sound velocity $c_{s}$ must provide ultimately a value of $z\geq 1$. Yet the values we find for $z$ and $\beta$ and, furthermore their different values for different dynamical rules, still make sense if we are analyzing cases (avalanche durations and sizes, for instance) for which $r/t=Ct^{1/z-1}=C^{z}r^{1-z}\ll c_{s}$.

IV.4 Relation to potential-energy-landscape-based models of yielding and justification of progressive rates

A different class of mesoscopic models to study the yielding transition , based on a continuous description of the strain in the system and the explicit inclusion of a quenched disorder landscape, was proposed Fernández Aguirre and Jagla (2018). A brief analysis of them will allow to clarify the meaning of uniform and progressive rates (Eq.9) we have used in EPMs, and the origin of Eq. 21. These models can be written through an equation giving the time evolution of the local strain in the system, $e_{i}$, in the form

[TABLE]

where $G_{ij}$ is the same elastic propagator used in Eq. 1, $\sigma$ is the applied stress, and $V_{i}$ indicate a set on random pinning potentials on which $e_{i}$ evolve.

A generalized mean-field treatment for this kind of models was proposed in Jagla (2017); Fernández Aguirre and Jagla (2018); Jagla (2018). It accounts to replace the mechanical noise term in 23 by a stochastic term $\eta_{i}(t)$ with correlation characterized by its Hurst exponent $H$, leading to decoupled equations of the form

[TABLE]

A detailed analysis of this stochastic one-particle modelJagla (2018) shows that the value of the flowcurve exponent $\beta$ can be calculated once the values of $H$ and the analytical form of the pinning potentials at the transition point between potential wells are know. The main result is:

[TABLE]

where $\omega=2$ for smooth potentials, and $\omega=1$ for potential with discontinuous forces. We now link this with EPMs.

Notice that all EPMs we have studied in this paper share the same kind of dynamics for the local stress $\sigma_{i}$: When $\sigma_{i}$ overpasses the local yield stress ${\sigma_{\tt Y}}$, it enters a stage in which there is a probability per unit time that the site becomes plastic. The average time for this stochastic activation process is note the “activation time” $\tau_{\tt on}$. The inverse of $\tau_{\tt on}$ time is the “activation rate” $\lambda$. Once in the plastic state ($n_{i}=1$) the site evolves towards a new relaxed position, decreasing its stress and accumulating plastic strain. The time related with this second stage of the local yielding event varies enormously with the particular EP model.

Yet, it is expected that the activation time will dominate the total duration of the plastic event.

The same kind of behavior for the local stress is realized in models like the one of Eq. 23.

Consider the limiting case of a ‘narrow wells potential’ $V$. In this case, the local strain $e_{i}$ evolves in a random disordered pinning potential as sketched in Fig. 10(a). The wells are assumed to be extremely narrow, in such a way that the value of $e_{i}$ is constant within the well 333In the analogous EPM, $n_{i}=0$ while the site is in the well and $n_{i}=1$ when it is out.. Each well is characterized by the force that it has to be applied in order for the site to escape from it. This threshold force is the analogous of the local yield stress ${\sigma_{\tt Y}}$ in EPMs. As long as the locally applied stress $\sigma$ (originated in the externally applied stress and the contribution from all other sites through $G_{ij}$) is lower than ${\sigma_{\tt Y}}$, the value $e_{i}$ remains constant. When $\sigma>{\sigma_{\tt Y}}$ the value of $e_{i}$ increases, and eventually reaches the value corresponding to the next potential well. Due to the increase of $e_{i}$ there is also a decrease of the local stress, similarly to what happens in EPMs. The important point to consider is the time $\tau_{e}$ it takes for this process to occur, i.e., for $e_{i}$ passing from a given potential well to the next one. To evaluate $\tau_{e}$, its is not enough to know that potential wells are narrow, we must consider also their shape. Referring to Fig. 10(b), let’s suppose that we have some smooth pinning well defined by some potential energy function.

The maximum of the derivative of this potential will define the local yield stress. If the overstress $\delta\sigma\equiv\sigma_{i}-{\sigma_{\tt Y}}$ is small, the time to escape from the well $\tau_{e}$ can be estimated by solving an equation of the form

[TABLE]

with initial condition $e_{i}\sim 0$ at $t=0$, and looking for the time at which $e_{i}\sim 1$. A direct solution of the equation, or a dimensional analysis Jagla (2018), provides $\tau_{e}\simeq\delta\sigma^{-1/2}$. This is the reason to consider –in the model analyzed in this paper– a “progressive” transition rate of the form $\lambda_{prog}\sim\delta\sigma^{1/2}$.

The previous argument can be generalized to potential wells with a more general form at the escape point, like $|e|^{\psi}$, instead of $e^{2}$ in Eq. 26. The analysis of this case yields $\lambda_{prog}\sim\sigma^{(\psi-1)/\psi}$. However, to have potential wells with this characteristic requires an extremely fine tuning of the potential, that is required to be non-analytical at the point of its maximum derivative. The only sensible possibility (in addition to the smooth case) is that of a potential that sharply vanishes at a finite value of $e$ (Fig. 10(c)). This case corresponds in fact to $\psi=1$ and provides $\lambda=cte$

; and so, justifies our assertion that constant activation rates can be assimilated to pinning potentials with successive wells separated by singular points in which there is a jump in the pinning force.

Given this equivalence, Eq. 25 can be considered to be valid also for EPMs, identifying the $\omega=1$ case with constant transition rates, and $\omega=2$ with progressive transition rates.

V Summary

In this paper we have investigated different versions of elasto-plastic models (EPMs) discussed in the literature addressing their critical properties at the yielding transition. Accurate numerical simulations have revealed that for the three cases analyzed all critical exponents are the same, thus suggesting that these values are universal, independent of model details. Nevertheless, we have identified a dynamical rule in EPMs on which “dynamical” critical exponents depend upon. This is the form of the rate law used to fluidize sites that have exceeded the local yielding threshold. Previous implementations of EPMs in the literature considered this rate $\lambda$ to be a constant, independent of the degree of overstress above the threshold value. We have gone beyond that, analyzing also the case of a progressive rate that depends on the degree of overstress, in particular $\lambda\sim(\sigma_{i}-{\sigma_{\tt Y}}_{i})^{1/2}$. This change of dynamical rule has a strong effect on the flowcurve exponent $\beta$ (which changes from $\beta=3/2$ to $\beta=2$) and the dynamical exponent $z$. The change of rule impacts in exactly the same manner on the three models analyzed. Furthermore, we have investigated also the effect of progressive rate in mean-field. For the well known Hébraud-Lequeux model the inclusion of progressive rates transforms the flowcurve exponent from the standard $\beta=2$ value to $\beta=5/2$ in the progressive case. Critical exponents that we have called “static” do not depend on the yielding rate rule: in particular, the exponents $\tau_{S}$ and $d_{f}$ for the avalanche size distribution and the exponents describing the density of sites at the verge of yielding. For the latter, we find a distribution of the form $P(x)\simeq P(0)+x^{\theta}$. Therefore, $\vartheta=d/\phi-1$ extracted from the finite-size scaling of $\langle x_{\tt min}\rangle\sim N^{-\phi}$, happens to be different from $\theta$. Both are independent on model and rate-rule details. Also the Hurst exponent $H$ derived from the accumulated mechanical noise signal is a static and model-independent exponent. It further allows to propose an alternative mean-field approximation to yielding which we find to give good predictions for the exponents extracted from simulations of fully spatial EP models in this work.

Conceptually, we believe that progressive yielding rates (i.e., smooth pinning potentials) should be more representative of real systems. It would be more difficult to think of a real system as having an effective sharp potential. Even though a direct comparison with experiments could be difficult, molecular dynamics simulations can be set to deal with interaction potentials of the smooth or sharp type, and investigate if the same variation of $\beta$ we find here is also found in those atomistic simulations.

The finding of the critical exponent $\beta$ depending on microscopic details of the disorder potential for the yielding transition has arisen somehow as a surprise, since it is commonly expected (based in the evidence provided by renormalization group arguments for the depinning transition) that these microscopic details should be irrelevant. A second thought, however, shows that this influence of dynamical rules that we observe in yielding is analogous to the one observed in depinning in the mean-field fully-connected case, where $\beta$ also depends on the particular shape of the disordered potential (see for instance Kolton and Jagla (2018)). Actually, the elastic coupling $G_{\bf q}$ for yielding (Eq. 3) is in fact independent of the absolute value of ${\bf q}$ (exactly as in the fully interacting depinning case, where $G_{\bf q}=-A$, with $A$ a constant for all $\bf q$), generating a long range interaction in real space $G(r)\sim 1/r^{d}$. This ‘sorority’ between the yielding transition of amorphous solids and the fully-connected mean-field depinning classical problem is subject of analysis of a separate work Ferrero and Jagla (2019).

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

EEF acknowledges support from grant PICT 2017-1202, ANPCyT (Argentina).

Appendix A Appendix

A.1 Progressive rates in the Hébraud-Lequeux model

The Hébraud-Lequeux model Hébraud and Lequeux (1998) is defined by the evolution of the probability distribution function ${\mathcal{P}(\sigma,t)}$ of local stress $\sigma$ at time $t$, under an external strain rate ${\dot{\gamma}(t)}$ as

[TABLE]

where the rate $\nu$ is given by

[TABLE]

(with $\Theta$ the Heaviside function and $\delta$ the Dirac distribution) i.e., the rate of plastic events is approximated by a fixed value $1/\tau$ in any overstressed region, and the plastic activity is defined from

[TABLE]

Finally, a key ingredient of the model is that the diffusion coefficient $D_{\tt pl}(t)$ is assumed to be proportional to the plastic activity:

[TABLE]

where $\tilde{\alpha}>0$ is an ad-hoc parameter of the model. In a stationary situation, as the plastic activity is proportional to $\dot{\gamma}$ for low shear rates Agoritsas and Martens (2017), the previous equation also implies that $\dot{\gamma}\sim D_{\tt pl}$. For any $\tilde{\alpha}$ smaller than a critical value $\tilde{\alpha}_{c}=\frac{\sigma_{c}^{2}}{2}$, a non-trivial solution exist for the (nonlinear) evolution in ${\mathcal{P}(\sigma,t)}$ that translates in a flowcurve ($\left<\sigma\right>=\sigma_{c}+A\dot{\gamma}^{n}$) with $n=1/2$ (i.e., $\beta=2$).

In the following we will show that modifying Eq.(28) to

[TABLE]

i.e., introducing a stress-dependent progressive rate in the HL model, the flowcurve results in

[TABLE]

Rather than attempting a full solution of Eq.(27) with the modified rate (31), we simply concentrate in tracking the impact of the functional form of the local plastic rate in the scaling properties of the stress close to the critical point. First, notice that the result $\dot{\gamma}\sim D_{\tt pl}$ at first order Agoritsas et al. (2015); Agoritsas and Martens (2017) is still valid in this case. At the same time, the stress excess $\Delta\sigma\equiv\sigma-\sigma_{c}$ when we impose a small strain rate is proportional to the probability density at the local threshold $\mathcal{P}(\sigma_{c})$.

Looking for a stationary solution of Eqs. (27,31) at zero order in $\dot{\gamma}$, we can write for $\sigma>\sigma_{c}$

[TABLE]

The solution for $\mathcal{P}$ should have the form

[TABLE]

with $f(0)=1$. Now, asking for the usual closure equation of the HL model (30) that provides self-consistency, we obtain

[TABLE]

The last integral is just a number, then we conclude that

[TABLE]

and

[TABLE]

In other words,

[TABLE]

for the HL model modified with progressive rates for local yielding of the form $\sim(\sigma-\sigma_{c})^{\eta}$. In the particular case corresponding to smooth potentials we must take $\eta=\frac{1}{2}$, obtaining $\beta=5/2$ which is precisely the value predicted by Eq. (25) for this case.

A.2 Avalanche duration vs. size in mean-field

A toy model with discrete pinning sites and stress dependent transition rates can be used to investigate analytically the dependence between duration and spatial extent of avalanches. Imagine a pinning potential consisting of narrow wells in which an elastic interface is pinned. Each well has a threshold force value that is necessary to apply to unlock the interface from it. Once a local threshold is exceeded the interface locally jumps to the next well with probability $\lambda$ (the inverse of the time it will take a particle to move from one well to the next). $\lambda$ may be a constant or a quantity depending on the force excess over the threshold: $\lambda=(f_{i}-f_{i}^{th})^{\alpha}$, where $f_{i}$ is the actual value of the force applied at site $i$, $f_{i}^{th}$ is the local threshold force, and $\alpha$ is a numerical exponent. A configuration of the system is characterized by the values of $f_{i}$, which are the elastic forces at every site. The configuration is stable if $f_{i}<f_{i}^{th}$. On this configuration an avalanche is triggered by increasing uniformly the force up to the point where $f_{i}$ becomes equal to $f_{i}^{th}$ at some particular $i$, which then becomes unstable. After some time the unstable site jumps to the next well, $f_{i}$ is reduced by some quantity $\delta$, and all $f_{j}$ move up a quantity $\delta/N$. This may produce some new sites to become unstable (i.e., to overpass their $f^{th}$). The process continues until there are no more unstable sites. At this point an avalanche of some duration $T$ and some spatial extent $S$ has occurred.

We can obtain the relation between $T$ and $S$ for this avalanche as follows. Suppose we plot the number of unstable sites $n$ as a function of the number of sites $k$ that have already jumped to their new positions within the avalanche. Such a plot is a random walk. The avalanche size $S$ is the number of sites that have jumped when there are no more unstable sites, namely $n(S)=0$.

In order to calculate $T$ for a given $S$ we have to sum all time intervals $dt_{k}$ between successive values of $k$, namely

[TABLE]

Let us consider first the simplest case in which $\lambda$ is constant. This means that any individual site jumps in a time that is always $\sim 1$. If there are $n(k)$ unstable sites at some moment, the typical time for the first one of them to jump is $dt_{k}\sim 1/n(k)$. In addition, the average form of $n(k)$ for a random walk between two zero crossings at 0 and $S$ is $n(k)\sim\sqrt{k(S-k)}/L^{1/2}$ then we obtain

[TABLE]

that for sufficiently large avalanches can be written (using $x\equiv k/S$) as

[TABLE]

i.e., we obtain $z=1/2$ which is the standard value of the dynamical exponent for depinning in mean-field.

The previous calculation can be extended to the case of an arbitrary value of $\alpha$ in $\lambda=(f_{i}-f_{i}^{th})^{\alpha}$. It requires the evaluation of $dt_{k}$ in this situation. The calculation is now non-trivial, because every unstable site has a different value of $y_{i}\equiv f_{i}-f^{th}$, and then a different transition rate. The value of $dt_{k}$ can be calculated as

[TABLE]

(note that if $\alpha=0$, $dt_{k}=1/n(k)$ and we return to the previous case). To calculate $dt_{k}$ the actual distribution of $y_{i}$ values if needed. We focus on its calculation now.

Consider a situation with $n$ unstable sites with values $y_{1}$ … $y_{n}$. We will work in a continuum limit, with $n$ large. We want to calculate the expected distribution of $y$, that will be characterized by a function $P(y)$, such that

[TABLE]

As the rate goes as $y^{\alpha}$, we can write

[TABLE]

The equilibrium condition for $P(y)$ is that when one of the $y_{i}$ jumps (and then disappears as an unstable site) and all the distribution is shifted by $1/N$ the configuration remains stable. This leads to

[TABLE]

The constant $C_{0}$ remains yet undetermined, and will be fixed by normalization below. Combined with the previous equation this gives

[TABLE]

and from here

[TABLE]

The two constants are determined from the normalization condition (Eq. 45), and from the fact that $P(0)=N$, since on average, a shift in $y$ of $1/N$ produces the appearance of one new unstable site. The final result is

[TABLE]

with $a_{0}$ a numerical constant. The continuous form of Eq. (44) allows to calculate the average time interval up to the next particle jump $dt_{k}$ as

[TABLE]

Doing the integration, the result is

[TABLE]

Using this expression to perform the same analysis that was done before in Eqs. (42), (43) leads to (assuming $\alpha<1$)

[TABLE]

[TABLE]

This is the final result. It provides the value of the dynamical exponent as $z=(1-\alpha)/2$ and, at the same time, a dependence on the system size $N$. The latter disappears in the uniform rate case ($\alpha=0$), corresponding to sharp pinning potentials.