How collective asperity detachments nucleate slip at frictional interfaces

TL;DR

This paper presents a microscopic model of frictional interfaces showing how asperity detachments lead to slip nucleation, revealing a pseudo-gap distribution and finite size effects that influence stick-slip behavior.

Contribution

It introduces a novel asperity-level disorder model with elastic interactions and inertia, linking frictional slip nucleation to amorphous plasticity concepts and finite size effects.

Findings

Slip nucleation governed by a Griffith criterion for asperity avalanches.

Presence of a pseudo-gap in asperity yield distances with a non-universal exponent.

Stick-slip behavior is a slow finite size effect with a diverging nucleation radius.

Abstract

Sliding at a quasi-statically loaded frictional interface can occur via macroscopic slip events, which nucleate locally before propagating as rupture fronts very similar to fracture. We introduce a novel microscopic model of a frictional interface that includes asperity-level disorder, elastic interaction between local slip events, and inertia. For a perfectly flat and homogeneously loaded interface, we find that slip is nucleated by avalanches of asperity detachments of extension larger than a critical radius $A_c$ governed by a Griffith criterion. We find that after slip, the density of asperities at a local distance to yielding $x_\sigma$ presents a pseudo-gap $P(x_\sigma) \sim (x_\sigma)^\theta$, where $\theta$ is a non-universal exponent that depends on the statistics of the disorder. This result makes a link between friction and the plasticity of amorphous materials where a pseudo-gap is also present. For friction, we find that a consequence is that stick-slip is an extremely slowly decaying finite size effect, while the slip nucleation radius $A_c$ diverges as a $\theta$-dependent power law of the system size. We discuss how these predictions can be tested experimentally.

Equations

$$P(S)\sim S^{-\tau }$$

$$S\sim A^{d_f}$$

$$P(A)\sim A^{-d_f(\tau -1)-1}$$

$$A_c\sim (\sigma -{\sigma }_c{)}^{-2}$$

$$n_{a}P(A>A_c)\sim 1$$

$$P(A>A_c)\sim A_c^{d_f(1-\tau )}\sim (\sigma -{\sigma }_c{)}^{-2d_f(1-\tau )}$$

$$P(x_{\sigma })\sim (x_{\sigma }{)}^{\theta }$$

$${\Phi }_a\sim {\int }_0^{{\Delta }_{\sigma }}(x_{\sigma }{)}^{\theta }d{x}_{\sigma }\sim (\sigma -{\sigma }_c{)}^{\theta +1}$$

$$n_a=N{\Phi }_a\sim N(\sigma -{\sigma }_c{)}^{\theta +1}$$

$${\sigma }_n-{\sigma }_c\sim N^{\frac{-1}{2d_f(\tau -1)+\theta +1}}\sim N^{-0.16}$$

$$A_c(\sigma ={\sigma }_n)\sim N^{\frac{2}{2d_f(\tau -1)+\theta +1}}\sim N^{0.32}$$

$$\rho\;\vec{a}(\vec{r})=\vec{\nabla}\cdot\bm{\sigma}\big{(}\bm{\varepsilon}(\vec{r})\big{)}-\alpha\vec{v}(\vec{r})$$

$$\mathbf{\sigma }\equiv \frac{K}{d}\mathrm{tr}\left(\mathbf{\epsilon }\right)\mathbf{I}+G{\mathbf{\epsilon }}_{\mathrm{d}}$$

$$\mathbf{\epsilon }\equiv \frac{1}{2}\left(\vec{\nabla }\vec{u}+(\vec{\nabla }\vec{u}{)}^T\right)$$

$${\mathbf{\epsilon }}_{\mathrm{d}}\equiv \mathbf{\epsilon }-\frac{1}{d}\mathrm{tr}(\mathbf{\epsilon })\mathbf{I}$$

$$\epsilon \equiv \sqrt{\frac{1}{2}{\mathbf{\epsilon }}_{\mathrm{d}}:{\mathbf{\epsilon }}_{\mathrm{d}}}$$

$$V(\epsilon )=G{\epsilon }^2$$

$$V\Big{(}\varepsilon_{\mathrm{y}}^{(i)}\leq\varepsilon<\varepsilon_{\mathrm{y}}^{(i+1)}\Big{)}=G\,\bigg{[}\,\Big{[}\,\varepsilon-\varepsilon_{\mathrm{min}}^{(i)}\,\Big{]}^{2}-\Big{[}\,\Delta\varepsilon_{\mathrm{y}}^{(i)}\,\Big{]}^{2}\,\bigg{]}$$

$$\bm{\sigma}(\bm{\varepsilon})=\frac{K}{d}\mathrm{tr}\left(\bm{\varepsilon}\right)\bm{I}+G\,\Big{[}\,\varepsilon-\varepsilon_{\mathrm{min}}^{(i)}\,\Big{]}\;\bm{N}_{\mathrm{d}}\qquad\mathrm{for}\;\varepsilon_{\mathrm{y}}^{(i)}\leq\varepsilon<\varepsilon_{\mathrm{y}}^{(i+1)}$$

$${\mathbf{N}}_{\mathrm{d}}\equiv \frac{{\mathbf{\epsilon }}_{\mathrm{d}}}{\epsilon }$$

$$P(\Delta {\epsilon }_{\mathrm{y}})=k(\Delta {\epsilon }_{\mathrm{y}}{)}^{k-1}\exp \left[-(\Delta {\epsilon }_{\mathrm{y}}{)}^k\right]$$

$$\int\limits_{\Omega}\rho(\vec{r})\;\delta\vec{u}(\vec{r})\cdot\vec{a}(\vec{r})\;\mathrm{d}\Omega=\int\limits_{\Omega}\delta\vec{u}(\vec{r})\cdot\Big{[}\,\vec{\nabla}\cdot\bm{\sigma}(\vec{r})-\alpha\vec{v}(\vec{r})\,\Big{]}\;\mathrm{d}\Omega$$

$$\int\limits_{\Omega}\rho(\vec{r})\;\delta\vec{u}(\vec{r})\cdot\vec{a}(\vec{r})\;\mathrm{d}\Omega=-\int\limits_{\Omega}\big{[}\,\vec{\nabla}\delta\vec{u}(\vec{r})\,\big{]}:\bm{\sigma}(\vec{r})\;\mathrm{d}\Omega-\int\limits_{\Omega}\alpha\;\delta\vec{u}(\vec{r})\cdot\vec{v}(\vec{r})\;\mathrm{d}\Omega\qquad\forall\;\delta\vec{u}(\vec{r})\in\mathbb{R}^{d}$$

$$\vec{u}(\vec{r})$$

$$\delta \vec{u}(\vec{r})$$

$$\underbrace{\int\limits_{\Omega^{h}}\rho(\vec{r})\;\varphi_{k}(\vec{r})\;\varphi_{l}(\vec{r})\;\mathrm{d}\Omega^{h}}_{\displaystyle M_{kl}(\vec{r})}\;\vec{a}_{l}=-\underbrace{\int\limits_{\Omega^{h}}\big{[}\,\vec{\nabla}\varphi_{k}(\vec{r})\,\big{]}:\bm{\sigma}(\vec{r})\;\mathrm{d}\Omega^{h}}_{\displaystyle\vec{f}_{k}(\vec{r})}-\underbrace{\int\limits_{\Omega^{h}}\alpha(\vec{r})\;\varphi_{k}(\vec{r})\;\varphi_{l}(\vec{r})\;\mathrm{d}\Omega^{h}}_{\displaystyle D_{kl}(\vec{r})}\;\vec{v}_{l}$$

$$\underline{\underline{M}}\underline{\vec{a}}=-{\underline{\vec{f}}}_{\mathrm{int}}=-\underline{\vec{f}}-\underline{\underline{D}}\underline{\vec{v}}$$

$${\underline{\vec{u}}}^{(n+1)}={\underline{\vec{u}}}^{(n)}+{\Delta }_t{\underline{\vec{v}}}^{(n)}+\frac{1}{2}{\Delta }_t^2{\underline{\vec{a}}}^{(n)}$$

$$\hat{\underline{\vec{v}}}^{(n+1)}=\underline{\vec{v}}^{(n)}+\tfrac{1}{2}\Delta_{t}\Big{[}\,\underline{\vec{a}}^{(n)}+\underline{\vec{a}}\big{(}\underline{\vec{u}}^{(n+1)},\underline{\vec{v}}^{(n)}+\Delta_{t}\underline{\vec{a}}^{(n)},t^{(n+1)}\big{)}\,\Big{]}$$

$$\underline{\vec{v}}^{(n+1)}=\underline{\vec{v}}^{(n)}+\tfrac{1}{2}\Delta_{t}\Big{[}\,\underline{\vec{a}}^{(n)}+\vec{a}\big{(}\underline{\vec{u}}^{(n+1)},\underline{\hat{\vec{v}}}^{(n+1)},t^{(n+1)}\big{)}\,\Big{]}$$

$$\rho {\partial }_t^{2}u=\frac{1}{2}G{\partial }_x^{2}u-\alpha {\partial }_{t}u$$

Full text

Significance statement

Understanding how slip at a frictional interface initiates is important for a range of problems including earthquake prediction and precision engineering. The force needed to start sliding a solid object over a flat surface is classically described by a ‘static friction coefficient’: a constant established by measurements. It was recently questioned if such constant exists, as it was shown to be poorly reproducible. We provide a model supporting that it is stochastic even for very large system sizes: sliding is nucleated when, by chance, an avalanche of microscopic detachments reaches a critical radius, beyond which slip becomes unstable and propagates along the interface. It leads to testable predictions on key observables characterising the stability of the interface.

Keywords: Friction; Inertia; Avalanches; Fracture; Stick-slip

Introduction

The sliding of a block that rests on a flat surface starts when the applied tangential force passes some threshold $F_{S}$, which is proportional to the normal force $F_{N}$. Their ratio defines the friction coefficient $\mu\equiv F_{S}/F_{N}$, which typically decreases with increasing sliding velocity when the latter is small (Scholz and Engelder, 1976; Baumberger and Caroli, 2006; Rabinowicz, 1956; Marone, 1998; Heslot et al., 1994). This phenomenology leads to stick-slip, whereby driving a system quasi-statically results in periods of loading that are punctuated by sudden macroscopic slip events. Experimental observations support that these events proceed by ‘fracture’ (Xia et al., 2004; Rubinstein et al., 2004; Ben-David et al., 2010; Passelègue et al., 2013): after a nucleation phase in which slip appears locally and evolves slowly (Ben-David and Fineberg, 2011; Ohnaka and Kuwahara, 1990), a well-defined rupture front appears, that travels ballistically across the frictional interface, unzipping it. This front is accompanied by a stress field in the elastic bulk that is well described by that of a propagating crack (Svetlizky and Fineberg, 2014; Svetlizky et al., 2016). By contrast, the nucleation phase is much less understood. It is observed that (i) its spatial extension $A_{c}$ decreases with increasing shear stress (Ben-David and Fineberg, 2011), (ii) there is a considerable variability in the tangential force magnitude at which macroscopic slip nucleates (Ben-David and Fineberg, 2011; Popov, 2010; Rabinowicz, 1992) and (iii) acoustic emission (McLaskey and Glaser, 2011; Johnson et al., 2013) supports that nucleation occurs by bursts of spatially resolvable McLaskey and Glaser (2011) detachments of micrometer-sized asperities (Bowden and Tabor, 1954; Dieterich and Kilgore, 1994; Hyun et al., 2004). Explaining these facts is relevant for earthquake predictions (Brace and Byerlee, 1966) as well as to forecast the variability of the measured friction coefficient (Ben-David and Fineberg, 2011; Popov, 2010; Rabinowicz, 1992), of importance for precision engineering (Armstrong-Hélouvry et al., 1994).

At a continuum level, rate-and-state models (Dieterich, 1979; Rice and Ruina, 1983; Ruina, 1983; Scholz, 1998) are powerful phenomenological descriptions of frictional interfaces, in which friction depends on the sliding velocity as well as some history-dependent state $\phi$ of the interface. Several length scales appear in these approaches Ruina (1983); Ohnaka and Kuwahara (1990), including a Griffith length beyond which sustained slip-pulses can propagate, as well as a larger length at which these pulses can nucleate fracture, reminiscent of a first-order phase transition111In this dynamical phase transition, the order parameter is the strain rate while the control parameter is the stress. A first-order transition therefore refers to a discontinuous strain rate vs stress behaviour. (Brener et al., 2018). Yet these descriptions are coarse-grained and phenomenological, and connecting them with asperity-level phenomena where disorder, that arises from surface roughness, is preponderant remains a challenge. Likewise, the precise meaning of the state variable $\phi$, often thought as capturing the ageing of contacts, remains to be clarified. One microscopic view is that a sliding frictional interface shares similarities with the plastic flow of amorphous materials Baumberger et al. (1999). Interestingly, it was recently observed that the ‘state’ of bulk amorphous materials (Lemaître and Caroli, 2007; Karmakar et al., 2010) can be quantified by the density of soft spots about to yield locally, which scales as $P(x_{\sigma})\sim(x_{\sigma})^{\theta}$, where $x_{\sigma}$ is the stress increment required for a given spot to yield222The scaling $P(x_{\sigma})\sim(x_{\sigma})^{\theta}$ holds only for small $x_{\sigma}$, namely for the soft spots. The term “pseudo-gap” refers an exponent $\theta>0$ that corresponds to a singular depletion for small $x_{\sigma}$ Müller and Wyart (2015).. The existence of a non-trivial exponent $\theta>0$ was shown to be a necessary consequence of the long-range elastic interactions and of the non-monotonic (varying in sign) stress redistribution triggered by plastic events Lin et al. (2014a, b, 2015); Lin and Wyart (2016). This parallel raises the intriguing possibility that frictional interfaces are characterised by some exponent $\theta$ as well.

Our goal is to propose a description of frictional interfaces that captures disorder at the asperity level, long-range elastic interactions between local slip events, and inertia. When the latter is absent, the physics is well understood and falls into the universality class of the depinning transition333The depinning transition occurs for example when an elastic manifold is pulled through a disordered medium (Fisher, 1998). with monotonic interactions (Fisher, 1998; Kardar, 1998; Ferrero et al., 2013), for which a continuous transition at a unique, well-defined, macroscopic critical force $F_{c}$ (Middleton, 1992) separates a flowing and an arrested phase. There exists no macroscopic stick-slip, and at $F_{c}$ motion corresponds to power law avalanches in which many local slip events act in concert. However, what happens to this scenario when inertia matters (as it does at a frictional interface) is a matter of debate. The popular view is that inertia destroys criticality: the transition becomes first order, stick-slip appears, and for significant inertia macroscopic slip events are nucleated by a few asperities acting together (Fisher et al., 1997; Dahmen et al., 1998). However, another scenario has been proposed by Schwartz et al. (Schwarz and Fisher, 2003) based on a simplified cellular automaton model describing short-range elasticity, in which stick-slip is a slowly decaying finite size effect that vanishes in the thermodynamic limit, for which the transition is continuous. Nevertheless Schwartz et al. later argued (Maimon and Schwarz, 2004) that for physical systems the scenario developed in (Fisher et al., 1997) was presumably correct, and that the conclusions of (Schwarz and Fisher, 2003) on the absence of hysteresis in infinite systems were non-generic and only valid for a finely tuned model.

In this work we introduce a novel numerical model of flat, homogeneously loaded, frictional interfaces where inertia is properly treated by discretising the bodies in contact by finite elements. Asperities at the interface are described by elements endowed with a random potential that represents the presence of surface roughness, allowing for sudden local slip events when a local (random) threshold stress is reached444This abstraction presents similarities to the treatment of the local rearrangement of particles in amorphous solids as shear-transformation zones (Argon, 1979), and allows similar numerical treatment (Homer and Schuh, 2009; Jagla, 2007, 2017).. Our model thereby differs from existing ones as no microscopic constitutive friction model (such as slip-weakening, velocity-weakening, or rate-and-state; see e.g. the spring–block models of (Andrews, 1976; Trømborg et al., 2014, 2015; Amon et al., 2017) and references therein) is presumed. In contrast, slip-weakening emerges as a consequence of mechanical noise (elastic waves generated by inertia) emitted when asperities detach, which can cause the nucleation of macroscopic slip.

Our main findings are that: (a) surprisingly, the scenario developed in (Schwarz and Fisher, 2003) is correct: stick-slip is a finite-size effect, although we find its power law decay with system size to be so slow that it is significant even in very large systems. In that regime, the friction coefficient is intrinsically stochastic, consistent with (ii) above. (b) Despite the interaction being monotonic, due to inertia the interface presents a non-trivial exponent $\theta$ characterising a pseudo-gap in the density of asperities about to yield. We argue that this conclusion will hold more generally to depinning problems with inertia and long-range elasticity. (c) Nucleation is triggered when avalanches get bigger than a critical radius, governed by a Griffith criterion, that diverges as a power law of the system size. Experimentally, the presence of avalanches is consistent with the measured distribution of acoustic emission (iii), while Griffith’s criterion is supported by a decreasing nucleation length with increasing stress (i). We relate the exponents associated with properties (a) and (c) to $\theta$ and those characterising avalanches, and confirm our predictions numerically. Finally, we propose experiments to test our results and measure $\theta$.

Model

The geometry of our set-up is illustrated in Fig. 1. It comprises a frictional interface (in red) embedded between two identical isotropic linear elastic materials (in blue), all discretised using finite elements and interacting in the same manner. The elements along the frictional interface (referred to as “blocks”) are plastic: they respond elastically (with the same elastic constants as the elastic bodies) up to a local yield strain (see below). To mimic energy leakage at the boundaries (by the transmission of elastic waves), we consider a viscous damping in the bulk whose magnitude is such that waves travel on the order of the system size before decaying (see Methods). The system is periodic in the horizontal direction, while the top and bottom boundaries are used to impose an event-driven quasi-static simple shear. In this protocol, the strain is increased up to the next plastic event (the response to this increase is purely elastic and in mechanical equilibrium), after which an infinitesimal strain increment is applied, triggering (an avalanche of) plasticity. Once motion stops, this sequence is repeated.

The frictional interface consists of $N$ ‘elasto-plastic’ blocks (finite elements) of linear size $h$, each representing one or a few asperities555More specifically, each block corresponds to the so-called Larkin length (Cao et al., 2018) below which asperities always collectively rearrange. Our predictions below apply if the Larkin length is much smaller than the whole system size. In the experiments of (Ben-David and Fineberg, 2011; Svetlizky et al., 2016; Svetlizky and Fineberg, 2014; Passelègue et al., 2013; Ben-David et al., 2010; Rubinstein et al., 2004), the nucleation length is found to be quite smaller than the system size, consistent with this assumption. See Appendix F for quantitative statements.. Similar blocks are used in models of plasticity of amorphous materials (Jagla, 2007, 2017). Each block is characterised by a random potential $V(\varepsilon)$ function of the equivalent shear strain $\varepsilon$ (the norm of the deviatoric part of the strain tensor), see Fig. 1(bottom). $V(\varepsilon)$ is constructed from a sequence of quadratic potentials of identical curvature, whose intersections define the yield strains. Disorder is introduced by randomly drawing the yield strains $\Delta\varepsilon_{\mathrm{y}}$ from some distribution, independently for each block. We chose a Weibull distribution $P(\Delta\varepsilon_{\mathrm{y}})=k\;(\Delta\varepsilon_{\mathrm{y}})^{k-1}\exp\left[-(\Delta\varepsilon_{\mathrm{y}})^{k}\right]$ with $k=2$. To acquire statistics we consider an ensemble of independent realisations and focus on the two-dimensional case where larger systems can be reached (see Methods and Appendix B for details).

A single plastic event

Under shear loading, a block responds linearly up to reaching the local yield strain, corresponding to a cusp in $V(\varepsilon)$. Passed that point, the block releases some of its elastic energy, and settles in a new equilibrium position determined by the potential energy of the element and the interaction with its surroundings. Such plastic shear strain leads to a permanent redistribution of shear stress in the system, that decays as a force dipole $1/r^{d}$ (Eshelby, 1956, and Appendix C), with $r$ being the distance from the block and $d$ the dimension of the space (here $d=2$). Along the weak layer the kick in shear stress is strictly positive (and decays in space as $1/r^{d}$), corresponding to a monotonic interaction. This effect alone can destabilise other blocks, leading an avalanche of yielding events.

In addition, each yielding event emits elastic waves, causing a transient stress, whose amplitudes decay in space as a force monopole $1/r^{d-1}$. This effect can trigger yielding of blocks that would have remained stable otherwise, causing a dynamical weakening effect: plastic activity leads to more inertial mechanical noise, which in turn creates more plastic activity.

Avalanches as precursors of macroscopic slips

A typical stress–strain response is shown in Fig. 2(a). The system first responds elastically, followed by a steady state stick-slip behaviour (highlighted in grey). The stick-slip phase consists of loading intervals punctuated by macroscopic slip during which all blocks yield many times, on average causing the stress to drop from $\sigma_{n}$ to $\sigma_{c}$ (see Fig. 2(a)). Such macroscopic slips are fracture-like, as supported by the time evolution that presents a ballistic propagation front, that travels at a super-shear velocity, consistent with the recent experiments (Rubinstein et al., 2004; Ben-David et al., 2010; Svetlizky and Fineberg, 2014; Svetlizky et al., 2016), see Appendix C.

As in experiments (Johnson et al., 2013), we observe microscopic activity during the loading phases. It corresponds to events that failed to nucleate macroscopic slip, and as such are important to analyse. The distribution of slip sizes $\tilde{S}$, defined as the total number of times that blocks yield during an event shows a clear separation in two types of events: macroscopic slips at $\tilde{S}\gg N$ (indicating that blocks have yielded many times), and avalanches that occur during the loading phase666Macroscopic slips are events in which all blocks yield at least once, and avalanches are all other, localised, events.; see Fig. 2(c) and the sketch in Fig. 2(b). However, the occurrence of avalanches is too rare to be insightful.

To gain more information about the avalanches, we manually trigger events at different stresses $\Delta_{\sigma}\equiv\sigma-\sigma_{c}$, by locally applying a shear displacement perturbation to a randomly selected block along the weak layer (see Appendix B). If all blocks were elastic, the displacement would simply snap back to the original equilibrium configuration. But for the elasto-plastic blocks an avalanche can be triggered, leading to a new equilibrium state.

We first focus on $\Delta_{\sigma}=0$. The distribution of avalanches sizes, $P(S)$, is obtained by eliminating events that result in macroscopic slip (defined as an event in which all $N$ blocks yielded at least once). Strikingly, we find a power law distribution of avalanches at $\sigma=\sigma_{c}$:

[TABLE]

with the exponent $\tau\simeq 1.5$ (Fig. 3(a)) and a fractal dimension $d_{f}\simeq 1.7$. The latter relates spatial extension $A$ (the number of sites that yielded at least once) of an avalanche to its size:

[TABLE]

(Fig. 3(c)). These results imply

[TABLE]

as confirmed in Fig. 3(b). We conclude that the stress $\sigma_{c}$, after macroscopic slip, is a critical point at which the distribution of avalanche sizes is scale free.

Mechanism for ‘fracture’ nucleation

Our central observation in Figs. 3(a,b) is that increasing $\sigma$ above the critical point $\sigma_{c}$ leads to a smaller and smaller cutoff $A_{c}$ and $S_{c}$ for the distribution $P(A)$ and $P(S)$. At first glance this is surprising, since at large stresses one may expect avalanches to be bigger. In fact, this cutoff signifies that large avalanches run away, and lead to macroscopic slip (not included in these distributions). Thus the cutoff $A_{c}$ and $S_{c}$ characterise the size of the avalanches required to nucleate a macroscopic slip event.

We now propose a scaling relationship for $A_{c}$ as a function of $\sigma-\sigma_{c}$. We posit that $\sigma_{c}$ is the maximum stress that the frictional layer can locally carry in the presence of endogenous inertial mechanical noise. This noise is generated by the ballistic pulses of stress emitted by failing blocks when the interface is in the process of plastically rearranging locally. Now consider triggering an avalanche at $\sigma>\sigma_{c}$. Avalanches are compact objects (as $d_{f}>1$), implying that each block yields many times, inducing a large inertial mechanical noise. On average, this will reduce the stress inside the avalanche to $\sigma_{c}$ (see sketch of Fig. 4(a)), while at large distances from the avalanche the stress remains $\sigma>\sigma_{c}$. This mismatch leads to stress concentrations at avalanche’s edges proportional to a stress intensity factor $(\sigma-\sigma_{c})\sqrt{A}$. As postulated by Griffith (Anderson, 2005; Griffith, 1921), a fracture instability777Note that in contrast to an opening crack, which cannot carry any stress, a stress $\sigma_{c}$ can still be carried during macroscopic slip. will take place when the intensity factor reaches a threshold, implying:

[TABLE]

(for any $d$). We confirm this result in Fig. 4(b), supporting our hypothesis that the stress inside the avalanche on average drops to $\sigma_{c}$. The departure from scaling in Fig. 4(b) at small $\Delta_{\sigma}$ is due to the value of $A_{c}$ being so large that our measurements suffer from finite size effects, see Appendix E. Note that we measure $A_{c}$ using the ratio of successive moments to extract $A_{c}\equiv\langle A^{p+1}\rangle/\langle A^{p}\rangle$. In practice we use $p=4$ to be more sensitive to the biggest avalanches while still having good statistics, but our results are robust to different choices, see Appendix E.

Macroscopic slip

Macroscopic slip nucleates when, by chance, an avalanche exceeds the nucleation radius $A_{c}$. As stress increases, more and more avalanches are triggered, while concurrently the nucleation radius shrinks. Nucleation of macroscopic slip thus becomes more and more likely with increasing stress. Typically, macroscopic slip will have happened when the stress is sufficiently large such that

[TABLE]

where $n_{a}$ is the number of triggered avalanches that have occurred as a result of a stress increment $\Delta_{\sigma}=\sigma-\sigma_{c}$, and ${P(A>A_{c})}$ is the fraction of those avalanches that exceed the radius at which macroscopic slip is nucleated. Eq. (5) thus sets the typical value of stress, $\sigma_{n}$, at which macroscopic slip occurs.

The probability that an avalanche has a radius larger than $A_{c}$ follows from Eqs. (3,4):

[TABLE]

as verified in Fig. 4(c).

The number of avalanches $n_{a}$ follows from the distribution $P(x_{\sigma})$ of the stress increment $x_{\sigma}$ required for a given block to yield for the first time after a big slip event and trigger an avalanche. Let us assert for the moment (and confirm below) that this distribution follows a power law:

[TABLE]

In that case, the fraction of blocks that triggers an avalanche upon increasing the stress by $\Delta_{\sigma}=\sigma-\sigma_{c}$ scales like

[TABLE]

This allows us to measure $\theta$ by counting the number of avalanches during the loading periods. We find a non-trivial exponent $\theta\simeq 3.7$, as shown in Fig. 5. For the number of avalanches, we thus get:

[TABLE]

Inserting Eqs. (6,9) into Eq. (5) leads to:

[TABLE]

It follows from this argument that: (i) The stress $\sigma_{n}$ at which macroscopic slip nucleates is stochastic, as embodied by Eq. (6). (ii) The stick-slip amplitude $\sigma_{n}-\sigma_{c}$ eventually vanishes as the number of asperities $N\rightarrow\infty$. Stick-slip is thus a finite size effect, yet the decay is so slow that it is expected to persist in realistic systems. In a truly infinite system avalanches should be power law distributed. (iii) The fracture nucleation radius diverges as:

[TABLE]

Argument for pseudo-gap $P(x_{\sigma})\sim(x_{\sigma})^{\theta}$

The stability distribution, $P(x_{\sigma})$, can in general obey one of three scenarios at small $x_{\sigma}$: (i) Depinning: a finite number of blocks can yield after a small increase of stress, characterised by an exponent $\theta=0$ (Lin et al., 2014a; Fisher, 1998). (ii) A pseudo-gap: the number of blocks that can yield vanishes only at $x_{\sigma}=0$, i.e. $\theta>0$. (iii) Gap: a small depleted region at small $x_{\sigma}$, such that $P(x_{\sigma}<D)=0$ for some small but finite $D$, thus requiring a finite increase of stress to destabilise any block. This scenario appears to be required to get true stick-slip as ${N\rightarrow\infty}$.

Our data in Fig. 5 and other measurements below support scenario (ii). We now exclude the depinning scenario based on a stability argument. In the presence of inertia, the temporary stress overshoot can destabilise blocks that would otherwise stop at small $x_{\sigma}$. Stability of the system requires that the number of blocks that are destabilised by one event does not diverge when the system size goes to infinity. This leads to the condition $\theta>0$ as follows: When a block fails, it emits a temporary stress overshoot $\sigma_{I}\sim 1/r$ (in 2D). The probability that this will destabilise other blocks is $P(x_{\sigma}<\sigma_{I})\sim r^{-(\theta+1)}$. Consequently, in a system of size $R$ the number of destabilised blocks $n_{f}\sim\int_{h}^{R}r^{-(\theta+1)}dr$ diverges as $R\to\infty$, unless $\theta>0$ (see Appendix D for a more general argument).

We currently do not have a theory for the value of exponent $\theta$, but preliminary observations indicate that $\theta$ is non-universal. Building a theory to understand $\theta$ should explain the following observations (presented in detail in Appendix D): (a) The blocks for which $x_{\sigma}$ is very small following a macroscopic slip event typically lie in a shallow well followed by another shallow well in the block (de Geus et al., 2015). (b) As a consequence, when triggered they tend to lead to small slips, and are less likely to trigger slip in other sites. As a result, there exists another exponent $\theta^{\prime}\simeq 2.5\leq\theta$ characterising the density of sites at a distance $x_{\sigma}$ to yield, unconditioned to subsequently triggering an avalanche (our argument and measure in Fig. 5 is conditioned to sites triggering an avalanche). (c) The exponent $\theta$ is not universal and depends on the specific choice of disorder, in particular on the parameter $k$ entering the Weibull distribution, and characterising the probability to find narrow wells. Using $k=1.2$ instead of $k=2$, we find $\theta^{\prime}\simeq 1.4$.

The presence of the strong depletion of the number of the almost unstable blocks, induced by the macroscopic slips, ($\theta>0$ for Eq. (7)) can be a possible explanation of the observed exponent $\tau\simeq 1.5$ that characterises the distribution of the avalanche sizes in Eq. (1). For the depinning transition, in the overdamped limit, we know that $\tau_{\text{dep}}=1.28$ (Bonamy and Bouchaud, 2011; Moulinet et al., 2004). However, we also know that each such avalanche is a collection of spatially disconnected slipping regions, called clusters. When treated as separate events, the distribution of avalanche sizes of individual clusters is also scale free, but with a larger exponent, $\tau_{\text{clus}}\simeq 1.56$ Laurson et al. (2010), close to our measured $\tau\simeq 1.5$. The existence of these disconnected clusters is a consequence of the long-range nature of the elastic interactions Joanny and de Gennes (1984): a slipping block is a source of instability for the neighbourhood, but also for blocks far away that are very close to their yield stress (have a small $x_{\sigma}$). The presently observed strong depletion of $P(x_{\sigma})$ for small $x_{\sigma}$ implies that there are very few blocks close to yielding, thus reducing the likeliness of triggering a ‘secondary’, disconnected, avalanche.

Discussion

We have introduced a model of a frictional interface that includes microscopic disorder at the asperity scale, long-range elastic coupling between local slip events, and the propagation of inertial waves. Our results support a description unifying collective avalanches of asperity detachments and fracture-like macroscopic slip events, in which the former nucleates the latter once a critical avalanche size is reached. These predictions are compatible with existing observations: the presence of avalanches is consistent with the measured distribution of acoustic emission (McLaskey and Glaser, 2011; Johnson et al., 2013), while Griffith’s criterion is supported by a decreasing nucleation length with increasing stress (Ben-David and Fineberg, 2011). Two surprises emerge from our predictions. First, a key aspect of the interface is the distribution of asperities about to yield, which is very much depleted and characterised by a non-trivial exponent $\theta$ after a macroscopic slip event. Second, we find that the transition to sliding is a continuous transition in the thermodynamic limit, but that finite size effects decay extremely slowly: the stress drop is a stochastic quantity whose typical scale decays as $N^{-0.16}$ and will thus persist in very large system, leading to a slowly diverging nucleation radius $A_{c}\sim N^{0.32}$.

Our predictions could be quantitatively tested in nearly flat and homogeneously loaded samples, which may be achievable experimentally using the apparatus of (Sahli et al., 2018; Bureau et al., 2000). In particular, microscopic slip events could be measured using (an array of) mechanical or acoustic sensors like in (McLaskey and Glaser, 2011; Rubinstein et al., 2004). Their cumulative number while quasi-statically loading the sample by a stress increment $\Delta_{\sigma}$ after a macroscopic slip event is proportional to $(\Delta_{\sigma})^{\theta+1}$, thus allowing one to access empirically the pseudo-gap exponent $\theta$. Moreover, well-separated avalanches could be acquired using our trick of triggering avalanches at different stress levels after macroscopic slip, for instance by supplying a focused acoustic signal to the system and measuring the magnitude of the mechanical or acoustic response. We expect the distribution of the magnitude to display a power law $P(S)\sim S^{-\tau}$ with a cutoff $S_{c}$ decreasing as $S_{c}\sim A_{c}^{d_{f}}\sim(\Delta_{\sigma})^{-2d_{f}}$. Beyond $\tau$, such a measurement would thus also yield an estimate of the fractal dimension of the avalanches $d_{f}$, without the need to spatially resolve the avalanches. As a reference, we document the statistics needed to extract these exponents reliably using our model in Appendix F.

There is an apparent opposition between the description presented here, and rate-and-state models where velocity-weakening is assumed to hold in the continuous limit, and nucleation stems from a first order transition. It would be very interesting to study how these two scenarios evolve when disorder is present at all scales (including the fact that the surface can have a roughness exponent, and the loading can be very heterogeneous). It is possible that in our approach as well, the transition becomes first order for certain statistics of the disorder. We view it as an important extension of the present work. Another important extension is the inclusion of creep. It may be readily achievable by putting our model in contact with a thermal bath, since in that case individual asperities will age to find a deeper nearby well.

Finally, it is interesting to ask which class of dynamical transitions can become first order due to inertia, and which cannot. The role of inertia has been studied recently in amorphous materials (Karimi et al., 2017; Nicolas et al., 2016; Salerno and Robbins, 2013; DeGiuli and Wyart, 2017; Vasisht et al., 2018), where it leads to a large pseudo-gap exponent $\theta$ comparable to ours (Karimi et al., 2017) (and much larger than the one present in the absence of inertia in these materials). It has been proposed that depending on the amount of damping, different universality classes could exist Nicolas et al. (2016); Salerno and Robbins (2013), but that for strongly underdamped systems the transition appears to become first order (Karimi et al., 2017). If confirmed, we speculate that the cause of the difference between amorphous solids and frictional interfaces is that avalanches are compact objects (having a fractal dimension $d_{f}>1$) only in the latter case. If that would not be true, our assumption that the inertial noise within an avalanche is comparable to that occurring in a macroscopic slip event may not hold, possibly leading to different physics.

Acknowledgement

T.G. was partly financially supported by The Netherlands Organisation for Scientific Research (NWO) by a NWO Rubicon grant number 680-50-1520. M.W. thanks the Swiss National Science Foundation for support under Grant No. 200021-165509 and the Simons Foundation Grant ($\#$454953 Matthieu Wyart). We acknowledge an anonymous referee for useful comments on experimental validation.

Methods

We consider two ensembles, each consisting of independent realisations comprising an approximately square box characterised by $N$ blocks along the weak layer, with $N=3^{6}$ ($300$ realisations) and $N=3^{6}\times 2$ ($1000$ realisations). The mechanical response is approximately incompressible, which allows us to focus on the shear response. The box is assumed periodic in horizontal direction. Quasi-static shear is applied by fixing the displacement of the bottom boundary to zero, while incrementing the displacement of the top boundary in very small steps (though we efficiently skip periods in which no yielding takes place, by homogeneously distributing the shear strain). Loading is stopped when the local strain exceeds a maximum.

After each step the energy is minimised according to the following equation of motion:

[TABLE]

(see Appendix A for nomenclature). From left to right, this equation comprises (i) an inertial term, in which $\rho$ is the mass density and $\vec{a}=\partial_{t}^{2}\vec{u}$ is the acceleration (where $\vec{u}$ is displacement and $t$ is time); (ii) the divergence of the stress tensor $\bm{\sigma}$; and (iii) a non-Rayleigh damping term, where $\alpha$ is the damping coefficient and $\vec{v}=\partial_{t}\vec{u}$ is the velocity. The stress $\bm{\sigma}$ follows from strain $\bm{\varepsilon}$ (which is the symmetric gradient of the displacement $\vec{u}$) using the constitutive model outlined in the main text. We set $\alpha$ such that kinetic energy is effectively leaked at the (periodic) boundaries.

Eq. (12) is solved in the weak form by discretising in space and time. In space, we discretise using finite elements. These elements coincide with the elasto-plastic blocks along the weak layer, while the elastic domain is discretised using elements that are conveniently chosen to increase in size with increasing distance to the weak layer to save computational costs. We discretise in time using the velocity-Verlet protocol.

Note that we formulate our model under the small strain assumption. To respect this assumption but still acquire a decently long steady state response, we choose the yield strains to be very small. This fixes the absolute strain and stress values to be small, which we rescale for visualisation to be of order one. See Appendix B for details. Furthermore, note that the numerical implementation is open-source de Geus (2018a, b) and that we have made all data underlying this manuscript freely available de Geus (2019).

Appendix A Nomenclature

[TABLE]

Appendix B Model

B.1 Constitutive model

Linear elasticity

The constitutive model, as illustrated in Fig. 1, is based on linear elasticity. Such behaviour is provided by Hooke’s law

[TABLE]

where $K$ and $G$ respectively are the bulk and shear modulus, and $d$ is the number of dimensions ($d=2$ in our case). Here $\bm{\varepsilon}$ is the (linear) strain tensor, defined as the symmetric gradient of the displacement $\vec{u}$, i.e.

[TABLE]

Finally, $\bm{\varepsilon}_{\mathrm{d}}$ is the strain deviator which contains all shear components of the strain, i.e. all strain components that do not lead to change of volume:

[TABLE]

We assume that yielding is a shear transformation. The volumetric response (through $K$) is therefore assumed to be purely elastic. Plasticity is defined through the potential energy. To obtain the shear part of Eq. (13) in an energetic picture, we need to introduce an equivalent shear strain

[TABLE]

This quantity thus characterises the magnitude of the shear strains encompassed in the strain deviator $\bm{\varepsilon}_{\mathrm{d}}$. With this, we retrieve Eq. (13) with a quadratic potential energy for the deviatoric part of the strain:

[TABLE]

(where $\bm{\sigma}_{\mathrm{d}}\equiv\partial V/\partial\bm{\varepsilon}_{\mathrm{d}}$). For completeness we introduce the work-conjugate equivalent shear stress $\sigma^{2}\equiv 2\bm{\sigma}_{\mathrm{d}}:\bm{\sigma}_{\mathrm{d}}$.

Elasto-plasticity

Following Jagla (2017), a shear transformation is introduced by using a manifold of quadratic potentials, given by:

[TABLE]

see Fig. 1. The elastic response is always governed by the shear modulus $G$. The intersections of the potentials are set by the sequence of yield strains $\varepsilon_{\mathrm{y}}^{(i)}$ (where $\varepsilon_{\mathrm{y}}^{(i+1)}>\varepsilon_{\mathrm{y}}^{(i)}$). Furthermore, $\Delta\varepsilon_{\mathrm{y}}^{(i)}\equiv(\varepsilon_{\mathrm{y}}^{(i+1)}-\varepsilon_{\mathrm{y}}^{(i)})/2$ and $\varepsilon_{\mathrm{min}}^{(i)}\equiv\varepsilon_{\mathrm{y}}^{(i)}+\Delta\varepsilon_{\mathrm{y}}^{(i)}$. Compared to Ref. (Jagla, 2017) our model remains isotropic, which corresponds to a yield strain that bounds the elastic domain using a sphere in principal deviatoric strain space. Finally, we obtain the following expression for the stress tensor:

[TABLE]

whereby the direction of shear is contained in

[TABLE]

Note finally that $\varepsilon_{\mathrm{p}}=\varepsilon_{\mathrm{min}}^{(i)}$ may be interpreted as a plastic strain, and thus that $\varepsilon_{\mathrm{e}}=\varepsilon-\varepsilon_{\mathrm{min}}^{(i)}$ corresponds to an elastic strain.

Parameters

The yield strains are drawn from a Weibull distribution

[TABLE]

for which we use $k=2$, see Fig. 6. To avoid difficulties in our algorithm (which is freely available (de Geus, 2018b)) we add a small offset $\delta_{\mathrm{y}}$ to each drawn yield strain, such that $P(\Delta\varepsilon_{\mathrm{y}}<\delta_{\mathrm{y}})=0$. To minimise the duration (in strain) of the transient initial loading preceding the steady state (highlighted in Fig. 2(a)), the first yield strain of each block is taken from a uniform distribution $\Delta\varepsilon_{\mathrm{y}}=[\delta_{\mathrm{y}},1)$. Incompressibility is approximated by $K/G=10$, whereby the shear modulus $G=1$. Furthermore, the shear wave speed $c_{s}=\sqrt{G/(2\rho)}=1/\sqrt{2}$ (note that the factor $1/\sqrt{2}$ appears because of our definition of the shear modulus). A crucial final point is that the absolute strains (and stresses) are fixed by the yield strains. To stay, as much as possible, within the small strain limit we scale the yield strains by $\varepsilon_{0}=10^{-3}/2$ (and use $\delta_{\mathrm{y}}=10^{-5}/2$). All strains and stresses that appear in diagrams have been rescaled by this typical strain and corresponding typical stress making them $\mathcal{O}(1)$.

B.2 Finite element discretisation

Overview

We solve Eq. (12) by discretising space according to the Finite Element Method (FEM). Besides our elasto-plastic blocks, we discretise also the elastic region using quadrilateral elements. This is illustrated in Fig. 1, where it is observed that we systematically coarsen the regions that present less interesting physics to reduce computational costs. FEM treats Eq. (12) in a weak sense. The resulting volume integral is evaluated element-by-element using numerical quadrature, in our case using four Gauss points. It is at these points that the stress and strain are evaluated, and thus where the potentials that are illustrated in Fig. 1 are defined. To fix the disorder to the scale of the blocks, all four Gauss points in one element get the same local potential energy landscape. The discretised weak form of Eq. (12) is, furthermore, discretised in time using the velocity-Verlet protocol. This results in incremental updates for the nodal velocity $\underline{\vec{v}}$ and displacement $\underline{\vec{u}}$, based the nodal acceleration $\underline{\vec{a}}$ that results from solving the discrete equation of motion taking into account the fixed displacement and periodic boundary conditions (as illustrated in Fig. 1). Note that our FEM code is also freely available (de Geus, 2018a).

Weak form

The crux of the Finite Element Method is to solve Eq. (12) in its weak form. For this, one has to satisfy

[TABLE]

for any possible test function $\delta\vec{u}(\vec{r})$. Note that $\Omega$ is the volume of the box. We note once more that the stress $\bm{\sigma}(\vec{r})$ depends in some specific way on the strain $\bm{\varepsilon}(\vec{r})$: the symmetric gradient of the displacement $\vec{u}(\vec{r})$ (see Appendix B.1 for details). In order to be able to evaluate the integral element-by-element, partial integration is employed next. This results in

[TABLE]

whereby the boundary integral, that incorporates the external forces that appear from partial integration has been omitted as its contribution is irrelevant because we fix all displacements along the top and bottom boundaries, and assume periodicity along the rest of the boundary. The discretised external forces needed to sustain the prescribed displacements can be easily retrieved, as discussed in Appendix B.3.

Discretisation in space

The problem is now discretised in space using a set of nodes that are connected through elements. Shape functions $\varphi_{k}(\vec{r})$ are used to interpolate the nodal displacement and test functions throughout the discretised domain $\Omega^{h}$. Note thereby that $\varphi_{k}(\vec{r})$ is locally supported, being non-zero only in the elements that contain the node $k$ (see Fig. 7 for a one-dimensional example). Furthermore, the shape functions constitute to a partition of unity. Their expression is standard, and for our four-noded quadrilateral elements they are bilinear.

For the displacement field and test functions, we thus have that

[TABLE]

where $k$ loops over all nodes. When applied to Eq. (23) we get

[TABLE]

Which is automatically satisfied for all nodal test functions. The integrals are finally evaluated at a discrete set of quadrature points. For our quadrilateral elements this step is exact when using four so-called Gauss points. We make one approximation here for the integrals that result in $M_{kl}$ and $D_{kl}$, by choosing equally weight quadrature points that coincide with the nodes. By making this approximation both matrices become diagonal, and their numerical treatment, including inversion, very cheap. Physically this corresponds to concentrating the mass as point masses at the nodes (whereby the mass still depends on the corresponding volume). We finally note that the strains at the Gauss points, that are needed to compute the stresses there, are obtained from the interpolation of the nodal displacements using Eq. (24). Finally, we introduce a short hand notation that reads

[TABLE]

Discretisation in time

To solve the second order differential equation in time, we proceed by discretising time. For this we employ the velocity-Verlet protocol, which:

Computes the position at time $t^{(n+1)}=t^{(n)}+\Delta_{t}$:

[TABLE] 2. 2.

Estimates the velocity at time $t^{(n+1)}=t^{(n)}+\Delta_{t}$ (by solving Eq. (27)):

[TABLE] 3. 3.

Corrects $\hat{\underline{\vec{v}}}^{(n+1)}$ (by solving Eq. (27)):

[TABLE] 4. 4.

Computes $\underline{\vec{a}}^{(n+1)}$ by solving Eq. (27) (using $\underline{\vec{u}}^{(n+1)}$ and $\underline{\vec{v}}^{(n+1)}$).

Parameters

To set the background damping such that waves whose wavelength is longer than the system’s size are critically damped we use the following one-dimensional wave equation

[TABLE]

Critical damping is found for wave numbers

[TABLE]

substituting $q=2\pi/(Nh)$ gives us the value of $\alpha$ that we seek.

The time step is taken much smaller than the time needed for a shear wave to travel to the shortest length scale in our problem: the block’s size. Consequently we take $\Delta_{t}=(1/(\sqrt{2}c_{s}q_{h}))/10$, with where $q_{h}=2\pi/h$.

B.3 Event-driven protocol

The box (e.g. in Fig. 1) is sheared, in simple shear, by prescribing the displacement of the top boundary. We stop the simulation when the local equivalent shear strain reaches $\varepsilon=0.5$ anywhere in the system (note that this is the absolute strain, before rescaling by the typical strain $\varepsilon_{0}$). As we seek to carefully measure the avalanches triggered by a local yielding event, we prescribe a very small equivalent shear strain change ($\delta_{\varepsilon}=10^{-7}$, also in terms of the absolute strain) during each loading increment. This coincides with the quasi-static protocol. The protocol is thereby that we affinely displace the entire box, such that the strain increment is homogeneous (see below). The top (and bottom) boundaries are then kept fixed, while the rest of the system evolves until energy has been minimised.

To enhance efficiency we make use of the relatively large strain intervals in which the entire box responds elastically. Since we drive very slowly and we know exactly how to distribute the strains to reach equilibrium (simply homogeneously in this case). We can thus transverse the entire elastic regime in one step, allowing us to run an event-driven quasi-static loading protocol. This protocol consists of two steps. In the first step, an affine displacement is added to the entire box such that the point that was closest to yield, is brought to the verge of yielding. In particular, the affine displacement is applied such that the increment in equivalent shear strain, $\Delta\varepsilon$, satisfies:

[TABLE]

(below we describe how this translates to an affine displacement increment). Since this step is purely elastic (and the displacement is affine), energy is instantaneously minimised. Then, in the second step, a small ‘kick’ $\Delta\varepsilon=\delta_{\varepsilon}$ is given to the system (again by applying an affine displacement), that causes yielding in at least one point.

An affine shear displacement field $\Delta u_{x}(y)=y\,\Delta\gamma$ is applied (where $x$ is the horizontal coordinate and $y$ is the vertical coordinate). This leads to the following local strain deviator:

[TABLE]

where $\varepsilon_{\mathrm{ss}}$ and $\varepsilon_{\mathrm{ps}}$ indicate the local contributions in simple shear (ss) and pure shear (ps) – the two principle deviatoric strains. This gives us an expression for the equivalent shear strain

[TABLE]

which can be solved exactly for $\Delta\gamma$, thereby taking the strain components from the point that is closest to yielding (see Eq. (33)).

The kick in strain is followed by an energy minimisation, using the equation of motion in Eq. (27) until all residual forces are sufficiently small. In particular we satisfy

[TABLE]

where the reaction forces, $\vec{F}_{k}$, measured at the nodes whose displacement is fixed (i.e. those at the top and bottom boundaries), are used for normalisation. They are

[TABLE]

B.4 Triggering

We measure the response at different stresses above $\sigma_{c}$ by manually triggering events after macroscopic slip (during which all blocks yielded at least once). In particular, we trigger at different $\sigma$ above $\sigma_{c}$ after macroscopic slip, as is illustrated in Fig. 2(a,b). This state can be reached exactly by applying an affine deformation to the last equilibrium state following each macroscopic slip (and an arbitrary number of avalanches) for which the stress is still smaller than $\sigma$, as sketched in Fig. 2(b). This protocol provides us with one equilibrium state following each macroscopic slip (provided that it ended at a sufficiently low stress, and that the next macroscopic slip nucleated at a sufficiently high stress). To acquire statistics we obtain more than one measurement from each relevant equilibrium state by triggering different blocks per realisation: each $73^{\mathrm{th}}$ block along the weak layer.

We trigger the event by temporarily applying a displacement fluctuation to the selected block, see Fig. 8. Note that the boundary conditions are not changed in any way, the system is free to possibly reach a new minimum governed by the same boundary conditions. The amplitude of the displacement fluctuation is such that the strain in the block is just above the yield strain; the same strain ‘kick’ is used as in Appendix B.3.

Appendix C (Typical) Response

C.1 A single plastic event

When a single block yields, it releases its built-up potential energy, accompanied by an increase in shear strain. This triggers dipolar force field on the surrounding blocks, which leads a change of stress whose amplitude decays in space as $1/r^{d}$, with $r$ being the distance to the yielding block and $d$ the number of dimensions ($d=2$ in our case). Along different directions this change is positive or negative, but along our weak layer the stress is strictly increased. This corresponds to the classical results by Eshelby (1956), accessibly summarised in a recent review of elasticity (Clouet et al., 2018).

Inertia introduces a second, transient, effect. Upon the sudden release of elastic energy, the surrounding blocks transiently experience the effect of a series of monopolar forces. This causes a series of temporary stress overshoots, followed by temporary stress undershoots in the surrounding blocks, whereby the amplitude of the stress overshoot decays in space as $1/r^{d-1}$.

We can observe these classical results in our system by triggering yielding of a single block embedded in an otherwise homogeneous elastic system (using the same simple-shear drive as is being used in the main text, see Fig. 1). Fig. 9 shows the response in the neighbouring blocks along the horizontal direction. The scaling of the permanent stress increase measured at long time (denoted by $\Delta\sigma$, in blue) and the temporary maximum stress overshoot (denoted by $\sigma_{I}$, in red) are consistent with our prediction.

C.2 Ballistic rupture front

To support our picture of an avalanche that nucleates a fracture, it is instructive to consider the time evolution during a single system-spanning event (corresponding to macroscopic slip characterised by a macroscopic stress drop, see Fig. 2(a)). In particular we expect to first see the avalanche as a fractal object (that is rather compact as $d_{f}>1$). Beyond a critical radius $A_{c}$ this object transitions into a clear rupture front. Since the system is finite, $N=3^{6}$ for the results below, nucleation happens at a relatively high stress, or a small nucleation radius $A_{c}$. We now consider a typical time evolution, in Fig. 10, whereby each marker corresponds to a local yield event. Indeed, after the first yielding event, an avalanche is observed in which the blocks yield over and over. Then, after some time, the avalanche succeeds in nucleating a fracture that is characterised by a well-defined rupture front. Although not studied here, we remark that the rupture front is supersonic, its velocity $\simeq 2c_{s}$ (which is, in accordance with the laws of elasticity, below the compressive wave speed, that for the used elastic constants is $\sqrt{10}\,c_{s}$, see Appendix B.1). After this front has crossed the entire box, yielding continues, corresponding to macroscopic slip.

Appendix D Stability & number of avalanches

D.1 Argument for pseudo-gap $P(x_{\sigma})\sim(x_{\sigma})^{\theta}$

The probability that a block will become unstable by a stress increase $\Delta_{\sigma}$ is

[TABLE]

The stress increase caused by the failure of a single block decays in space as

[TABLE]

with $q=d-1$ for the inertial stress overshoot and $q=d$ for the permanent stress increase, see Fig. 9. The number of blocks that will be destabilised by the failure of a block therefore reads

[TABLE]

The only way that $n_{f}$ does not diverge for $R\rightarrow\infty$, for neither the inertial stress overshoot with $q=d-1$ nor the permanent stress increase with $q=d$, is when $\theta>0$. Note that the number of failing blocks remains finite because the microscopic length scale, $h$, is finite. This proves that $P(x_{\sigma})$ displays a pseudo-gap, with $\theta>0$, or a gap (though the scenario of a gap is excluded by our data).

D.2 Distance to yielding

We distinguish: (I) the distance to yielding of individual blocks from (II) the distance to triggering an avalanche. Namely, blocks can yield without triggering an avalanche, which can alter the relevant value of exponent $\theta$. We measure both distributions independently. For the first distribution (I) we measure, for each block $i$, the additional amount of stress needed to yield: $x_{\varepsilon}^{(i)}\equiv\varepsilon_{\mathrm{y}}^{(i)}-\varepsilon^{(i)}$ (without loss of generality we directly use strain, which is uniquely related to the stress using the shear modulus: $x_{\sigma}\sim x_{\varepsilon}$). The distribution displays a pseudo-gap

[TABLE]

with an exponent $\theta^{\prime}\simeq 2.5$, as shown in Fig. 11(b) and 12(a,b). Note that this measurement is consistent with the scaling of the cumulative probability of yield events upon a stress increase $\Delta_{\sigma}=\sigma-\sigma_{c}$, that scales as $(\Delta_{\sigma})^{\theta^{\prime}+1}$ as shown in Fig. 11(a).

For the second distribution (II), we measure the probability that an avalanche occurs (that yielding occurs more than once) when manually triggering yielding of a block at a certain $x_{\varepsilon}$. We find that this probability scales like $(x_{\varepsilon})^{1.2}$, see Fig. 11(c). These measurements are thus consistent with $\theta\simeq 2.5+1.2=3.7$ (cf. Fig. 5).

D.3 Microscopic details matter

To test if the value of $\theta$ is universal we inspect the energy landscape around the configuration that is in mechanical equilibrium after macroscopic slip. In our case, this local energy landscape is fully defined by the yield strains in each block. On average, we find that the yield strains are large (40% larger than the average yield strain drawn from the yield strain distribution, which we set to $1$ after rescaling the data, see Appendix B.1). Since $\theta$ characterises the blocks close to yielding, we focus on the energy landscape around the blocks displaying a small $x_{\varepsilon}^{(i)}$. To this end, we compute the average yield strain bounding the local energy minimum of those blocks characterised by a small $x_{\varepsilon}^{(i)}$. We also compute the average yield strains bounding the next and previous local minima in the same block, as well as its left and right neighbours. Following (de Geus et al., 2015), we can compute:

[TABLE]

To bias this weighted average to those blocks that are closest to yield, we use

[TABLE]

where

[TABLE]

takes only the positive part of $\bullet$. As a reference, we include the average yield strain without biasing for small distance to yielding ($\tilde{x}_{\varepsilon}=x_{\varepsilon}$, denoted without any subscript). All results are normalised by the ensemble average yield strain $\langle\Delta\varepsilon_{\mathrm{y}}\rangle_{0}\approx 1$.

We find that those blocks that are close to yielding after macroscopic slip, have strong neighbours while they themselves are weak (their average local yield strain is 46% lower than the typical one), and also their next yield strain is equally low, see Fig 12. This indicates that the block can fail due to mechanical noise, but can also move back because of the low local yield strain. In the dying activity, such block can behave as if there was no inertia, thereby avoiding a gap in $P(x_{\sigma})$. This suggests that the presence of pseudo-gap and/or its exponent may depend on microscopic details (as hinted on in (Schwarz and Fisher, 2003)). We test this by considering a different distribution of yield strains, representing a different surface roughness, for which it is more likely to find two sequential low yield strains in the local energy landscape. To this end we use the Weibull distribution of Eq. (21) with $k=1.2$ (instead of $k=2$). Fig. 13 shows the stability distribution for an ensemble of $40$ realisations of $N=3^{6}$. As expected $\theta^{\prime}$ is decreased with respect to our other results.

Appendix E Verification of robustness

E.1 Cutoff

The proposed critical radius at which an avalanche nucleates macroscopic slip, in Eq. (4), can also be expressed in terms of avalanche size using the fractal dimension $d_{f}$. This corresponds to a critical avalanche size

[TABLE]

This scaling is verified in Fig. 14(D.1). Like for $A_{c}$, for our measurement we use $S_{c}\equiv\langle S^{p+1}\rangle/\langle S^{p}\rangle$ with $p=4$ to be mostly sensitive of the biggest avalanches that did not grow unstable, but the scaling is robust also for different choices of $p$ and for two other protocols, see below. With the measured $S_{c}$ we can count the fraction of events (avalanches or macroscopic slip) whose size $\tilde{S}>S_{c}$. Based the power law distribution of avalanche sizes in Eq. (1) we expect:

[TABLE]

as verified in Fig. 14(D.3).

Because the measurements of the scaling of the cutoff radius $A_{c}$ and of the cutoff size $S_{c}$ are crucial to test our theory, we test the robustness of these measurements. Above we have used

[TABLE]

with $p=4$. Here we compare the scaling for different values of $p$, and additionally consider

[TABLE]

and

[TABLE]

The latter is the average maximum when all independent measurements are separated in 20 ensembles of measurements. The results are shown in Fig. 14 for all crucial scaling measurements. As observed, all measurements are consistent with our theory, though in particular for $p=1$ the measurement is polluted by avalanches of small $A$ and $S$.

E.2 Critical stress $\sigma_{c}$ and system’s size $N$

In the main text, we measured the scaling of $A_{c}$, $S_{c}$, and the number of avalanches (characterised by the exponent and $\theta$) as a function of an increase of stress $\Delta_{\sigma}=\sigma-\sigma_{c}$ compared to the ensemble averaged critical stress $\sigma_{c}\equiv\langle\sigma_{c}(s)\rangle$ (where $s$ is an index that loops over all macroscopic slips in the ensemble). To verify the robustness of this protocol, we measure the scaling of these quantities as a function of an increase of stress $\Delta_{\sigma}^{\prime}\equiv\sigma(s)-\sigma_{c}(s)$ compared to the local critical stress: the value of stress directly after the last macroscopic slip. The results are fully consistent with the measured scaling relationships in the main text, as shown in Fig. 15. Note that, for consistency, we denoted $\sigma-\sigma_{c}^{\prime}\equiv\Delta_{\sigma}^{\prime}$.

Furthermore, in Fig. 15 we show all results for two system sizes $N=3^{6}\times 2$ (used in the main text) and $N=3^{6}$. Note that even smaller systems do not offer a perspective of validation as those systems do not display a clear separation between avalanches and macroscopic slip, see Appendix F. The results of $A_{c}$ and $S_{c}$, in Fig. 15(e,f), furthermore, emphasise that departure from the predicted scaling of $A_{c}$ and $S_{c}$, at small $\sigma-\sigma_{c}^{\prime}$, is a finite size effect. In particular, large $A_{c}$ and $S_{c}$ are better approximated for the largest system. The departure from scaling happens when $A_{c}$ is larger than $N/2$. For such a radius, the periodic repetitions are certainly felt by a propagating avalanche. This causes macroscopic slip to be nucleated sooner, which leads to the observed smaller than predicted $A_{c}$ and $S_{c}$.

Appendix F Required statistics

To aid setting up (experimental) validations of our theoretical predictions we list the (statistical) details of our measurements. We emphasise that because each (experimental) protocol comes with its own sources of uncertainties these numbers should be used merely as guideline. We, furthermore, emphasise that the number of avalanches scales with the system size $N$ through the stability exponent $\theta$ as in Eq. (9), a fact that can be used to optimise the employed (experimental) protocol.

• •

System size: $N$.

Throughout the text we use a system size of $N=3^{6}\times 2\approx 1500$ blocks, whose size is equal to the Larkin length. We find that all our results are robust for a system size of $N=3^{6}\approx 750$ (see Appendix E). We find that smaller systems do not display a clear distinction between avalanches and runaway slip events. Such finite size effects also appear in the distribution of the local distance to yielding, $P(x_{\varepsilon})$, as shown in Fig. 16.

• •

Stability distribution: $P(x_{\sigma})\sim(x_{\sigma})^{\theta}$, see Eq. (7).

The result in Fig. 5 is based on 8115 avalanches obtained from 2279 steady state stick-slip cycles. Note that using less than 2000 avalanches (or 570 stick-slip cycles) the power law scaling was not obvious.

• •

Avalanche exponent: $P(S)\sim S^{-\tau}$, see Eq. (1), or $P(A)\sim A^{-d_{f}(\tau-1)-1}$, see Eq. (3).

Fig. 3(b) is based on 4513 (triggered) avalanches. Note that we could not extract the correct exponent using less than 500 (triggered) avalanches.

Without triggering, the number of avalanches $n_{a}\sim N(\sigma-\sigma_{c})^{\theta+1}$ (Eq. (9)). In our model $\theta+1\simeq 4.7$ so that avalanches at stresses $\sigma\approx\sigma_{c}$ are rare. For our model, only 0.2% of the naturally formed avalanches are in the first bin of Fig. 3(b) (i.e. at a stress $\sigma\leq\sigma_{c}+\frac{1}{6}(\sigma_{n}-\sigma_{c})$), one thus needs roughly $6\cdot 10^{4}$ steady state stick-slip cycles if triggering is not used.