Elastic interfaces on disordered substrates: From mean-field depinning to yielding

TL;DR

This paper investigates a model of elastic interfaces on disordered substrates, revealing how critical exponents vary from mean-field depinning to finite-dimensional yielding, and highlighting the influence of microscopic pinning potentials.

Contribution

It introduces a generalized elastic model that interpolates between mean-field depinning and finite-dimensional yielding, showing smooth variation of critical exponents and dependence on microscopic details.

Findings

Critical exponents vary smoothly with elasticity modifications.

Herschel-Buckley exponent depends on the microscopic pinning potential.

Yielding in finite dimensions appears as a mean-field transition.

Abstract

We consider a model of an elastic manifold driven on a disordered energy landscape, with generalized long range elasticity. Varying the form of the elastic kernel by progressively allowing for the existence of zero-modes, the model interpolates smoothly between mean field depinning and finite dimensional yielding. We find that the critical exponents of the model change smoothly in this process. Also, we show that in all cases the Herschel-Buckley exponent of the flowcurve depends on the analytical form of the microscopic pinning potential. This is a compelling indication that within the present elastoplastic description yielding in finite dimension $d\geq 2$ is a mean-field transition.

Figures

Equations

$$\eta \frac{d{e}_i}{dt}=f_i(e_i)+\sum _j{G}_{ij}{e}_j+\sigma$$

$$\frac{d{e}_{\mathbf{q}}}{dt}=f(e){|}_{\mathbf{q}}+G_{\mathbf{q}}{e}_{\mathbf{q}}$$

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

$$\dot{\gamma }\equiv \frac{d\overline{e_i}}{dt}=\overline{f_i(e_i)}+\sigma$$

$$G_{\mathbf{q}}\equiv (1-\epsilon )G_{\mathbf{q}}^{\mathtt{Y}}+\epsilon G_{\mathbf{q}}^{\mathtt{MFD}}$$

$$G_{\mathbf{q}}^{3D}=\frac{2q_x^2(q_y^2+q_z^2)}{(q_x^2+q_y^2+q_x^2{)}^2}-1$$

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

$$e_L^{new}$$

$$e_R^{new}$$

$$f_i=-\sin (e_i-(e_L+e_R)/2).$$

$$\dot{\gamma }\equiv \frac{d\overline{e_i}}{dt}=\overline{f_i(e_i)}-\overline{e_i}\kappa +\sigma$$

$$f_i(e_i)=-(e_i-e_{0i})$$

$$e_{\mathbf{q}}=\frac{e_{0\mathbf{q}}}{1-G_{\mathbf{q}}}.$$

Full text

Statistical physics is built on analogies. The comparison of typically complex problems with a small number of simpler ones for which an exact solution is known is the first step in almost every argumentative construction. For instance, the so-called out-of-equilibrium phase transitions have been discussed in the mirror of equilibrium phenomena. The problem of depinning of an elastic manifold moving on a disordered landscape has been rationalized by analogy with the theory of critical phenomena Fisher (1998); Kardar (1998) and studied for over 30 years already. Once this problem has been reasonably understood, it serves in turn as the base model for a new analogy step. In this sense, depinning has shaped the theoretical endeavors in the understanding of the yielding transition of amorphous solids under deformation, that received full attention of the statistical physics community only recently. The problem with analogies is that, some times, they may prevent to see the big picture.

Members of the family of sand-pile problems, both depinning and yielding are paradigmatic examples of driven transitions and are intuitively very alike. Depinning is related to the movement of an elastic manifold in the presence of a quenched disordered potential, under the action of an external driving force. Yielding pertains to the flow of an amorphous solid upon the application of an external driving stress or deformation. In both cases, if the driving force is weak (and the possibility of thermal activation is excluded), the system remains in a frozen configuration; however, if a critical threshold is exceeded, the system reaches a dynamical state with a non-zero average velocity (depinning picture) or strain rate (yielding picture). The critical threshold defines the transition. In depinning the velocity-force characteristics of the system shows singular behavior at a critical force $f_{c}$. While $v=0$ for $f<f_{c}$, it behaves as $v\sim(f-f_{c})^{\beta}$ when $f$ increases above $f_{c}$, with $\beta$ a well defined number known as the flow exponent. In yielding the transition is characterized by the critical behavior of the strain rate $\dot{\gamma}$, which 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}$. The value of $\beta$ is referred again as the flow exponent. Its inverse $n\equiv 1/\beta$ is known as the Herschel-Buckley exponent.

The depinning transition finds a continuous model approach in the Quenched Edwards-Wilkinson equation, which allows for analytical treatment using Functional Renormalization Group (FRG) analysis. For yielding, elasto-plastic models (EPMs) built at a coarse grained scale Nicolas et al. (2018) provide a similar description Talamali et al. (2011); Lin et al. (2014a); Budrikis et al. (2017). Nevertheless, despite the analogous construction, an FRG treatment of EPMs has found limitations and analytical support for a theory of yielding is only provided so-far by mean-field variants Hébraud and Lequeux (1998); Agoritsas et al. (2015); Jagla (2015); Agoritsas and Martens (2017); Lin and Wyart (2016, 2018); Fernández Aguirre and Jagla (2018); Jagla (2018). In EPMs the instantaneous values of stress and plastic strain are evolved consistently. Under a condition of uniform load, the stress increases uniformly. When the stress exceeds locally a threshold value, the local plastic strain increases at that patch, causing a reduction of the local stress and also a perturbation of the stress in every other point in the system, following the action of elastic interactions. The form of this interactions is the one prescribed by the Eshelby propagator Eshelby (1957); Picard et al. (2004) of continuum mechanics; which has a $\sim 1/r^{d}$ spatial decay and thus it is a long-range interaction. Also, it has alternating signs depending of the direction, with a quadrupolar symmetry. This anisotropy is a curse for the FRG approach and the responsible for special avalanche correlations in the form of slip lines (or planes), that greatly determine the differences among the plastic yielding transition and its elastic depinning counterpart.

In the present paper we show that mean-field depinning and yielding transition in finite dimensions can be considered to be special cases of a generalized mean-field problem and, therefore, described within the same framework. Simply considering an elastic kernel as sum of two contributions $G^{\tt MFD}$ and $G^{\tt Y}$, corresponding respectively to a constant value propagator (in Fourier space) and the Eshelby propagator, we are able to smoothly interpolate ($\varepsilon G^{\tt MFD}+(1-\varepsilon)G^{\tt Y}$ with $0\geq\varepsilon\geq 1$) between mean-field depinning and yielding. In particular, we observe a smooth transition in the values of the critical exponents between the two limiting cases. Thus our work suggests an alternative view for the theoretical tackling of the yielding transition, interpreting it as a particular case of a general mean field problem that includes also the fully-connected mean-field depinning.

A general model for MF-depinning and yielding

The general model that allows to describe MF-depinning and yielding on the same footing is constructed in the following way. The variable of the model is a scalar field $e_{i}$ defined on the sites $i$ of a $d$ dimensional ordered lattice. For depinning $e_{i}$ represents the interface position at site $i$, whereas for yielding $e_{i}$ is the strain of an elemental volume of the system at site $i$. The dynamics is described by overdamped equations of motion of the form

[TABLE]

In the case of depinning, the terms $f_{i}(e_{i})\equiv-dV_{i}/de_{i}$ represent the force exerted by the external pinning potential $V_{i}$, whereas for yielding they describe the local thresholding behavior of a small piece of the amorphous material under deformation. In both cases the form of the functions $V_{i}$ are of the same kind: they have minima at different $e_{i}$ positions representing local equilibrium states. In Eq.(1) $G_{ij}$ represents the elastic interaction between $e$ values at different points. We restrict to cases in which this interaction preserves the homogeneity of the system, then $G_{ij}$ depends only on the difference between the (vector) positions $i$ and $j$. Also, $G_{ij}=G_{ji}$ is assumed. Elastic forces should be balanced in the system, therefore $\sum_{i}G_{ij}=0$ must be also satisfied. This still leaves us with a lot of freedom in the choice of a general form for $G_{ij}$. Nevertheless, an important additional constraint must be fulfilled: in the absence of local forces ($f_{i}\equiv 0$) the flat configuration of the interface $e_{i}={\tt cst}$ must be stable. This condition becomes more transparent in Fourier space, where Eq.(1) reads for ${\bf q}\neq 0$

[TABLE]

The stability condition is then $G_{\bf q}\leq 0$.

In the following we will mainly discuss the interaction kernel in Fourier space. One can consider “generalized mean field models” defined as cases in which the $G_{\bf q}$ is zeroth order homogeneous in $|q|$. These kernels produce a function $G_{ij}$ that is either independent of distance or decaying with $r_{ij}$ as $r_{ij}^{-d}$. In any of the two cases, the effect of a single site onto another site is negligible compared to the combined effect of all other sites in the lattice. Therefore, the dynamics of a given site can be solved by considering the existence of a (fluctuating) prescribed external field (see Fernández Aguirre and Jagla (2018)). In particular, the forms of $G_{\bf q}$ for mean-field depinning and yielding satisfy the prescription just mentioned. For mean field depinning $G^{\tt MFD}_{\bf q}=-1$ for ${\bf q}\neq 0$, whereas for yielding $G_{\bf q}$ is the Eshelby propagator that in $2D$ can be written as (${\bf q}\neq 0$)

[TABLE]

In both cases $G_{{\bf q}=0}$ is taken as zero in a stress conserved dynamics, as it follows from the condition $\sum_{i}G_{ij}=0$. The uniform mode in Eq.(1) is thus directly found from (we set $\eta\equiv 1$ for the rest of the paper)

[TABLE]

that defines the global strain rate $\dot{\gamma}$. The fact that both $G^{\tt MFD}$ and $G^{\tt Y}$ share the property of being $\mathcal{O}(q^{0})$, allows us to believe that mean-field depinning and yielding may share many common features.

With respect to the properties of the disorder term $f_{i}(e_{i})$, we restrict to the case of locally correlated potentials, where $\langle f_{i}(z)f_{i}(z+\Delta)\rangle$ decays to zero sufficiently fast with $\Delta$. Also, we consider the disorder site by site to be totally uncorrelated, namely $\langle f_{i}(z)f_{j}(z)\rangle=0$ for $i\neq j$. With these correlation properties, Renormalization Group theory teaches that the detailed form of $f_{i}(z)$ should be irrelevant when determining the critical properties of the transition, as long as the elastic interaction $G_{ij}$ decays sufficiently fast in space as a function of $r\equiv|r_{i}-r_{j}|$, concretely, if $G_{ij}\sim r^{-\alpha}$ with $\alpha>d$. In cases in which $\alpha\leq d$ (there included mean-field depinning and yielding) this result does not apply, and different values are obtained for the dynamical exponents when considering “cuspy” or “smooth” potentials Fisher (1998); Kolton and Jagla (2018). In the simulations presented below, we mainly focus in the case of a cuspy form of the pinning potential, taking $V_{i}$ as composed by a concatenation of parabolic pieces. A brief consideration of the smooth potential case is included at the end.

Results

We present simulations in $2D$ using a kernel

[TABLE]

that interpolates between the mean field depinning case (for $\varepsilon=1$) and the yielding case (for $\varepsilon=0$). The stress-controlled and quasistatic strain-controlled protocols to determine the flowcurves and the avalanche statistics are described in Appendix .1 in detail.

Both for an elastic interface undergoing depinning and for an amorphous solid at the onset of yielding, a singular behavior of $\dot{\gamma}$ (strain-rate or velocity) is expected at $\sigma_{c}$. Fig. 1 shows the flow curves for different values of $\varepsilon$. By plotting the data in logarithmic scale close to $(\sigma-\sigma_{c})$ 111$\sigma_{c}$ depends on the value of $\varepsilon$ and must be fitted for each curve, a clear power-law behavior allows to determine the flow exponent $\beta$. Going from $\varepsilon=1$ ($\tt MFD$) to $\varepsilon=0$ ($\tt Y$) we observe that the $\beta$ exponent moves from $\beta=1$ to $\beta\simeq 1.5$. Very importantly, this variation is smooth as the inset in Fig.1(b) shows, indicating the continuous evolution that exists between mean-field depinning and yielding.

We now discuss the avalanche size distribution $P(S)$ associated to the transition. When $P(S)$ is taken from large collections of avalanches obtained in a quasi-static simulation, it is expected to be power-law distributed, namely $P(S)\sim S^{-\tau}f(S/S_{\tt max})$ with the cutoff function $f(x)$ behaving as $f_{x\to 0}\to 1$ and $f_{x\to\infty}\to 0$, and $S_{\tt max}$ depending on the system size $L$ and the stress non-conserving parameter $\kappa$ used to define the value of $G_{{\bf q}=0}$ in the quasi-static protocol.

Avalanche size distributions are shown in Fig. 2 for different $\varepsilon$. 222For a better comparison, in Fig. 2 the values of $\kappa$ were adjusted to obtain curves with a similar cutoff. In depinning mean-field models $\tau=\frac{3}{2}$, and in fact, we obtain $\tau=1.5$ when $\varepsilon=1$. But as we decrease $\varepsilon$ moving towards yielding, $\tau$ diminishes, becoming $\tau\simeq 1.35$ at $\varepsilon=0$. Surprisingly, this change is continuous; meaning, without scale crossovers on the observable. The avalanche size distribution critical exponent is a smooth function of the parameter $\varepsilon$.

Directly related to the avalanche mean size is the loading stress needed to trigger avalanches $x_{\tt min}$. It is known for yielding Maloney and Lemaître (2004); Karmakar et al. (2010); Lin et al. (2014a), that its mean value scales sub-extensively with system volume $N$: $\overline{x_{\tt min}}\propto N^{-\alpha}$, with $1<\alpha<2$. This phenomenological sub-extensiveness in the plastic flow of amorphous solids under deformation was interpreted Lin et al. (2014b, a) as a consequence of a peculiar shape for the steady state distribution $P(x)$ of local distances to threshold $x$ 333In a driving protocol towards the right we define $x_{i}=e_{Ri}-e_{i}$ according to the notation in Appendix .1 . If this quantity has the form $P(x)\sim x^{\theta}$ as $x\to 0$, one can deduce Karmakar et al. (2010) that $\overline{x_{\tt min}}\propto N^{-1/(1+\theta)}$. Then, $\theta=0$ is expected for depinning-like models (where the kernel $G_{ij}$ is non-negative) and $\theta>0$ for yielding models (where the kernel $G_{ij}$ alternates in sign).

Figure 3 shows results for $\overline{x_{\tt min}}$ vs $L$ for systems with different $\varepsilon$. Power-law fits allow for a precise determination of the exponent values. Consistently with the expectation, $\theta=0$ for MF-depinning ($\varepsilon=1$) and a strictly positive value for the yielding case ($\varepsilon=0$). What is surprising again is that $\theta$ turns out to be a continuous function of the crossover parameter $\varepsilon$, going from [math] to $\simeq 0.5$ as we move from the MF-depinning limit to the 2$d$-yielding limit, as displayed in Fig. 3 inset. This tells us that we are working with a family of similar problems, each of them characterized by a given degree of sub-extensiveness of the load needed to trigger new avalanches.

Smooth pinning potentials– All results presented so far were obtained using a local disorder potential that has cusps in the transition from one potential well to the next one. Usually, according to renormalization arguments, this kind of details on the microscopic potential should not influence the critical properties of a system. In particular, the critical exponents of the depinning transition are expected to be independent on the potential being of the cuspy or smooth type. Nevertheless, the fully connected mean-field case is an exception (see discussion in Kolton and Jagla (2018)). There, we know that depinning displays a value $\beta=1$ for cuspy pinning potentials and a (different) value $\beta=3/2$ for smooth pinning potentials.

The smooth crossover of exponents that we observe between mean-field depinning and yielding in Fig.1, suggests that we will also find the above described dichotomy in the flow curve exponent value for the yielding case. Even more, we can expect to find larger values of $\beta$ using smooth potentials for any value of the crossover parameter $\varepsilon$.

Results of simulations contained in Fig. 4 confirm this result. Note that the $\beta$ value for smooth potentials always (i.e., for each $\varepsilon$) exceeds in $\frac{1}{2}$ the one for cuspy potentials, in full agreement with recent theoretical expectations derived from the Prandtl-Tomlinson model under stochastic driving Jagla (2018).

Three and larger dimensional cases– A scalar Eshelby kernel for $d=3$ in Fourier space (taking a diagonal component of the non-deviatoric stress tensor) can be written as

[TABLE]

One can notice that again $G_{\bf q}\sim{\bf q}^{0}$, and this is true in general for $d>1$. Therefore, we expect all our numerical observations and conclusions obtained in $d=2$ to be valid also in $d=3$ and higher dimensions. Preliminary simulations in $d=3$ for the cuspy potential give us, for instance, a smooth crossover between $\beta^{\tt Y}\simeq 1.3$ (not far from other estimations Lin et al. (2014a)) and $\beta^{\tt MFD}=1$ continuously moving with $\varepsilon$. The value of the pseudo-gap exponent $\theta$ also changes continuously with $\varepsilon$. The reduction of $\beta$ in passing from $2$ to $3$ dimensions can be rationalized as a consequence of the reduced density of zero modes in the elastic propagator in $3d$ compared to $2d$.

Why a smooth exponent crossover is surprising?– A remarkable point is that the critical exponents (values of $\beta$, $\theta$, and $\tau$ in particular) vary smoothly between those of $\tt MFD$ and $\tt Y$. This situation is not expected in general when studying crossovers between different asymptotic behaviors. Consider for instance the case of long range depinning. Choosing a kernel decaying in space as $G_{1}\sim 1/r^{\alpha_{1}}$ ($d<\alpha_{1}<d+2$) a set of critical exponents is obtained. For other decaying form of the kernel $G_{2}\sim 1/r^{\alpha_{2}}$ the exponents are different. However, if we combine the two kernels in the form $G=(1-\varepsilon)G_{1}+\varepsilon G_{2}$ the system will display the critical behavior corresponding to the lowest value of $\alpha$. In other words, if two different criticalities are mixed together the system will display at long enough scales the critical exponents corresponding to the longest range interactions. In order to have a variation of the critical exponents with $\varepsilon$, the long range weight of $G_{1}$ and $G_{2}$ must be similar. This non-common case is, in fact, what takes place in our kernel combination; clearly manifested in the ${\bf q}$ space form of the propagator: both the constant ($\tt MFD$) and the Eshelby ($\tt Y$) kernels scale as ${\bf q}^{0}$. The $\tt MFD$ kernel in Fourier space has a finite constant value, then the value of the interaction with any other site in real space is of order $\sim 1/L^{d}$. In real space, the Eshelby kernel decays as $1/r^{d}$ and, then, its typical value at distances of the order of the system size $L$ is also $\sim 1/L^{d}$.

Conclusions

We have studied a mesoscopic implementation of a scalar model for a generalized elastic manifold on a disordered landscape that is able to describe both depinning in mean-field ($\tt MFD$) or finite dimensional yielding ($\tt Y$) by changing the form of an elastic kernel interaction. The most important result is the observation of a smooth transition between $\tt MFD$ and $\tt Y$, as the kernel interpolates linearly between the two limiting cases. The identification of a common scenario for both transitions is assisted by recent reports in numerical simulations of yielding Fernández Aguirre and Jagla (2018); Jagla (2017); Ferrero and Jagla about phenomenological properties akin to mean-field depinning. In particular, dynamical critical exponents are seen to depend on the details of the local disorder potential, a result that has long been known for mean-field depinning Fisher (1998); Kolton and Jagla (2018).

In recent years, different attempts have been made to address’s the yielding transition of amorphous solids from an analytical perspective. Nevertheless, the noncompliance with the hypothesis needed for the Functional Renormalization Group analysis has largely confined these treatments to mean-field Hébraud-Lequeux-like approaches Agoritsas et al. (2015); Jagla (2015); Agoritsas and Martens (2017) and heavy-tail noise variants Lin and Wyart (2016, 2018); Fernández Aguirre and Jagla (2018); Jagla (2018). Those studies provided a common general picture but did not forge a consensus about critical exponents and scaling laws. One of the main conclusions of the present work is that the yielding transition, as described by a scalar model with an Eshelby interaction, can be treated as a special case of a generalized mean-field problem which has the very well known fully-connected MF-depinning problem as one limiting case. Our work suggest therefore that, instead of focussing on non-trivial correlations depending on the propagator properties and the dimension, an strategic angle of attack for theoretical studies of the yielding transition could be to start from a fully-connected depinning system and explore perturbations of the Eshelby type to the constant propagator.

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

.1 Simulation protocols

We do straightforward simulations of two-dimensional systems described by Eq.(1) integrated with a first order Euler method. In each integration step, the term $\sum_{j}G_{ij}e_{j}$ is treated in Fourier space, computed as $\sum_{\bf q}G_{\bf q}e_{\bf q}$. In this respect, note that in a square numerical mesh of size $L\times L$, quantities such as $q_{x}^{2}$ and $q_{y}^{2}$ must be understood as

[TABLE]

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

Concerning the form of the local forces $f_{i}=-dV_{i}/e_{i}$, we generate them ‘on the flight’ in the following way. We start with the value of $e_{i}$ such that $e_{L}<e_{i}<e_{R}$, where a parabolic potential is defined in term of $e_{L}$, $e_{R}$ and having an unitary curvature. In concrete, the force on $e_{i}$ coming from the disordered potential is taken as $f_{i}=-(e_{i}-(e_{L}+e_{R})/2)$. As soon as the dynamics makes $e_{i}$ larger than $e_{R}$, we set a new parabolic potential well for $e_{i}$ by choosing

[TABLE]

where the $\Delta$ is randomly chosen from a flat distribution between 0.5 and 1.5. This is what we call the ‘cuspy’ potential, since its composed by a concatenation of parabolic pieces and the transition from one to another produces a discontinuity in the force. In the case of ‘smooth potentials’ instead, the potential wells are defined and updated in the same form, but the force on each potential well is given by

[TABLE]

Note that in this form, the value of $f_{i}$ and also its derivative $df_{i}/e_{i}$ are continuous functions of $e_{i}$.

In constant stress simulations, the value of $\sigma$ in Eq.(1) is kept fixed, and the main output of the simulation is the value of $\dot{\gamma}$, Eq.(4). Such a protocol is used primarily to obtain the flow curves. A second important outcome of these simulations, is the distribution of local distances to instability, namely the distribution of the quantities $x_{i}\equiv e_{Ri}-e_{i}$. The average minimum value of $x_{i}$ is used to calculate the $\theta$ exponent.

On the other hand, we want to be able to access individual avalanches very close to the critical stress, and collect statistics of size, duration, etc. This can be better accomplished by using of a quasi-static protocol. To do so, we move away from the fix-stress modelling and modify $G_{\bf q}$ by defining $G_{{\bf q}=0}=-\kappa$, with $\kappa$ a constant parameter of order one. The equation for the evolution of the average strain is transformed to

[TABLE]

that can be interpreted as a progressive reduction of $\sigma$ (due to the term $\overline{e_{i}}\kappa$) as the average position of the interface moves forward. This stress reduction guarantees that any activity in the system will eventually stop, reaching a metastable static configuration. At this point the stress has to be increased again to trigger a new avalanche, and the process can be repeated. The evolution of stress along the simulation is sketched in Fig. 5. From such a simulation we collect the statistics of avalanche size ($S$) and stress increments needed to trigger new avalanches $x_{\tt min}$.

It should be mentioned, nevertheless, that quasistatic simulations in the form just described are rather inefficient: The dynamic evolution is continuous and we need to wait until the activity falls below a very low threshold to safely decide that the avalanche has stopped. In the same way, to trigger a new avalanche, the stress has to be increased very slowly to be sure to detect the precise beginning of the new avalanche. In the case of piece-wise parabolic potentials, an accelerated numerical scheme can be implemented as follows. Since we use a potential where all parabolic wells have the same curvature, the form of $f_{i}(e_{i})$ is simply given by

[TABLE]

with $e_{0i}=(e_{Li}+e_{Ri})/2$. Then Eq.(1) can be solved in a single step to obtain the new equilibrium position of the interface. The solution in Fourier space is given by

[TABLE]

If all the $e_{i}$ obtained by Fourier-inverting (13) lay within their potential well, i.e., $e_{Li}<e_{i}<e_{Ri}$, then the configuration found is a static solution to the problem. However, if some resulting $e_{i}$ happen to be outside the range $\left[e_{Li},e_{Ri}\right]$, it means that the corresponding $e_{0i}$ have to be adjusted and Eq.(13) solved again to find a new set of $e_{i}$. This process is repeated until all $e_{i}$ are within $\left[e_{Li},e_{Ri}\right]$. At this point the avalanche stops. In this scheme one actually loose the true continuous time evolution of the real dynamics, but it results to be computationally much more efficient and–as verified in test cases–it does not show any noticeable differences in the avalanche statistics with respect to the case in which the true dynamics is used.