# Critical nucleation length for accelerating frictional slip

## Abstract

The spontaneous nucleation of accelerating slip along slowly driven frictional interfaces is central to a broad range of geophysical, physical and engineering systems, with particularly far-reaching implications for earthquake physics. A common approach to this problem associates nucleation with an instability of an expanding creep patch upon surpassing a critical length . The critical nucleation length is conventionally obtained from a spring-block linear stability analysis extended to interfaces separating elastically-deformable bodies using model-dependent fracture mechanics estimates. We propose an alternative approach in which the critical nucleation length is obtained from a related linear stability analysis of homogeneous sliding along interfaces separating elastically-deformable bodies. For elastically identical half-spaces and rate-and-state friction, the two approaches are shown to yield that features the same scaling structure, but with substantially different numerical pre-factors, resulting in a significantly larger in our approach. The proposed approach is also shown to be naturally applicable to finite-size systems and bimaterial interfaces, for which various analytic results are derived. To quantitatively test the proposed approach, we performed inertial Finite-Element-Method calculations for a finite-size two-dimensional elastically-deformable body in rate-and-state frictional contact with a rigid body under sideway loading. We show that the theoretically predicted and its finite-size dependence are in reasonably good quantitative agreement with the full numerical solutions, lending support to the proposed approach. These results offer a theoretical framework for predicting rapid slip nucleation along frictional interfaces.

## I Introduction

The process of rupture nucleation in which slowly driven frictional interfaces (faults) spontaneously develop elastodynamically propagating fronts accompanied by rapid slip is of fundamental importance for various fields, with far-reaching implications for earthquake physics. Quantitatively understanding the nucleation process is essential for predicting the dynamics of frictional interfaces in general and for earthquake dynamics in particular. There exists some observational evidence, based on seismological records (Scholz, 1998; Ohnaka, 2000; Harris, 2017), and some experimental evidence, based on laboratory measurements (Dieterich, 1979; Ohnaka and Kuwahara, 1990; Kato et al., 1992; McLaskey and Kilgore, 2013; Latour et al., 2013), which suggest that rapid rupture propagation accompanied by a marked seismological signature is preceded by precursory aseismic slip. This precursory aseismic slip is commonly associated with a slowly expanding creep patch defined as a slipping segment of finite linear size , embedded within a non-slipping fault. Accelerating slip is expected to emerge once surpasses a critical nucleation length . We note that other nucleation scenarios have been considered in the literature, see for example Ben-Zion (2008), but are not discussed here.

Various theoretical and computational works have indicated that the nucleation of accelerating slip is related to a frictional instability (Ruina, 1983; Yamashita and Ohnaka, 1991; Ben-Zion and Rice, 1997; Scholz, 1998; Ben-Zion, 2001; Lapusta and Rice, 2003; Uenishi and Rice, 2003; Ben-Zion, 2008; Kaneko and Lapusta, 2008; Kaneko et al., 2016). From this perspective, the critical nucleation length corresponds to the critical conditions for the onset of instability that leads to accelerating slip and to the spontaneous propagation of elastodynamic rupture fronts. A major challenge is to understand the relations between the critical instability conditions and . In this Letter, we propose a theoretical approach for predicting which differs from the conventional approach.

The conventional approach, based on a single degree-of-freedom spring-block analysis extended to deformable bodies using various model-dependent fracture mechanics estimates, is discussed in the framework of rate-and-state constitutive laws in Sect. II. Our approach, based on the stability of homogeneous sliding of elastically-deformable bodies, is introduced in Sect. III and is shown to yield a significantly larger for elastically identical half-spaces and rate-and-state friction. In Sect. IV we show that the proposed approach is naturally applicable to bimaterial interfaces, which are of great interest in various contexts (Weertman, 1980; Andrews and Ben-Zion, 1997; Ben-Zion and Andrews, 1998; Cochard and Rice, 2000; Adams, 2001; Ben-Zion, 2001; Gerde and Marder, 2001; Ranjith and Rice, 2001; Rice et al., 2001; Shi and Ben-Zion, 2006; Rubin and Ampuero, 2007; Ampuero and Ben-Zion, 2008; Allam et al., 2014; Brener et al., 2016; Aldam et al., 2017), and derive analytic results for in this case, indicating that the bimaterial effect decreases compared to available predictions in the literature. Finally, in Sect. V we show that the proposed approach is applicable to finite-size systems and test our predictions against inertial Finite-Element-Method calculations for a finite-size two-dimensional elastically-deformable body in rate-and-state frictional contact with a rigid body under sideway loading. The theoretically predicted and its finite-size dependence are shown to be in reasonably good quantitative agreement with the full numerical solutions, lending support to the proposed approach. Section VI offers some concluding remarks and discusses some prospects.

## Ii A conventional approach to calculating the nucleation length

As stated, the most prevalent approach to the nucleation of rapid slip at frictional interfaces associates nucleation with an instability of a slowly expanding creep patch. The creep patch features a non-uniform spatial distribution of slip velocity, in the quasi-static regime (where inertia and acoustic radiation are negligible), due to some external loading. It is assumed to be stable as long as its length is smaller than a critical nucleation length . When , the patch becomes unstable and transforms into a rupture front, accompanied by accelerated slip and dynamic propagation (where inertia and significant acoustic radiation are involved). As creep patches are non-stationary objects that involve spatially varying fields, determining their stability — and hence — is a non-trivial challenge that typically requires invoking some approximations.

The most common approximation proceeds in two steps (Dieterich, 1986, 1992; Lapusta et al., 2000; Kaneko and Lapusta, 2008). First, the creep patch and the two elastically deformable bodies that form the frictional interface are replaced by a rigid block of mass in contact with a rigid substrate and attached to a Hookean spring of stiffness . That is, all of the spatial aspects of the problem are first neglected. The external loading and the typical slip velocity within the patch are mimicked by constantly pulling the Hookean spring at a velocity . The rigid block is pressed against the rigid substrate by a normal force , which gives rise to a frictional resistance force , where is described by the friction law, which may depend on the block’s slip , its time-derivatives and the state of the frictional interface.

This single degree-of-freedom spring-block system is described by the force balance equation where each superimposed dot denotes a time-derivative. We assume that can be described by the rate-and-state constitutive framework, where is a function of the slip velocity and of an internal state variable . The latter, which quantifies the typical age/maturity of contact asperities, evolves according to , where is a memory lengthscale and the function satisfies and . For example, two popular choices, i.e. (Ruina, 1983; Marone, 1998; Nakatani, 2001; Baumberger and Caroli, 2006; Bhattacharya and Rubin, 2014) and (Ruina, 1983; Gu et al., 1984; Bhattacharya and Rubin, 2014), feature .

Consider then a steady sliding state at a constant driving velocity such that . A standard linear stability analysis implies that this steady state becomes unstable if (Rice and Ruina, 1983; Ruina, 1983; Gu et al., 1984; Lapusta et al., 2000; Baumberger and Caroli, 2006; Bhattacharya and Rubin, 2014)

(1) |

where an inertial term proportional to has been neglected. That is, an instability is predicted when the spring stiffness is smaller than a critical stiffness . Note that since generically , a necessary condition for instability is , i.e. that the sliding velocity belongs to the velocity-weakening branch of the steady state friction curve (Ruina, 1983).

In the second step, the analysis is extended to spatially varying fields and elastically deformable bodies — relevant to realistic creep patches — by identifying the spring stiffness in the spring-block system with an -dependent effective stiffness in the spatially varying and elastically deformable system. This is typically done through some fracture mechanics estimates which yield (Dieterich, 1986; Rice, 1993)

(2) |

where is the shear modulus, is the nominal contact area and the dimensionless number is a model-dependent pre-factor. As expected physically, the effective stiffness of the overall system, , is a decreasing function of the length of the creep patch, . Using then of Eq. (1) as an instability criterion, one obtains

(3) |

where . The numerical pre-factor is model-dependent (e.g. it depends on the crack configuration, dimensionality and loading configuration) and its value varies between and in the available literature (Dieterich, 1992, see Table 1). The nucleation criterion in Eq. (3), with close to unity, is widely used in the literature, though we are not aware of computational or experimental studies that quantitatively and systematically tested it. Next, we present a different approach for calculating .

## Iii An approach based on the stability of homogeneous sliding of elastically-deformable bodies

Our goal here is to propose an alternative approach to calculating the critical nucleation length . In the proposed approach, nucleation is viewed as a spatiotemporal instability occurring along the creep patch which is assumed to be stable from the fracture mechanics perspective, i.e. to propagate under stable Griffith energy balance conditions (Freund, 1990). Since, in general, an elastic body can be thought of as a scale-dependent spring, one expects short wavelength (large wavenumber ) perturbations to be stable and instability — if it exists — to emerge beyond a critical (minimal) wavelength (i.e. below a critical wavenumber ). Consequently, when the size of the expanding creep patch is small, , we expect it to be stable. A loss of stability is expected when an unstable perturbation can first fit into the creep patch, i.e. when the patch size satisfies .

In this physical picture, the major goal is to calculate the critical wavenumber . There is, however, no unique and general procedure to study the stability of non-stationary (time-dependent) and spatially varying solutions such as those associated with an expanding creep patch. Consequently, we invoke an approximation in which the spatially varying slip velocity within the creep patch is replaced by a homogeneous (space-independent) characteristic slip velocity . With this approximation in mind, we need to study the stability of steady-state homogeneous sliding of an infinitely long system (in the sliding direction) in order to calculate . Applying the result to the actual creep patch, accelerating slip nucleation is predicted to occur when . This idea has been introduced, pursued and substantiated in the context of thin layers sliding on top of rigid substrates in Bar-Sinai et al. (2013). Our aim here is to significantly generalize the idea to any frictional system.

We consider a long elastic body in the -direction of height in the -direction steadily sliding with a relative slip velocity on top of a long elastic body of height . The bodies may be made of different elastic materials and are pressed one against the other by a normal stress , see Fig. 1 (left). As we are interested in the response of the system to spatiotemporal perturbations on top of the homogeneous sliding state at a velocity , we define the slip displacement and the slip velocity , where is the displacement field and is the fault plane (the superscript +/– means approaching the fault plane from the upper/lower body side, respectively). for each body satisfies the Navier-Lamé equation , with its own shear modulus , Poisson’s ratio and mass density Landau and Lifshitz (1986). The Cauchy stress tensor field was related to the displacement field through Hooke’s law and each superimposed dot represents a partial time derivative.

The fault at is assumed to be described by the rate-and-state constitutive relation . Fault opening or interpenetration are excluded, i.e. we assume , and and are continuous across the fault. The internal state field evolves according to , with and , as in Sect. II. We then introduce interfacial slip perturbations of the form , where is the complex growth rate and is the wavenumber. The shear and normal stress perturbations are related to using the solution of the quasi-static Navier-Lamé equation, and take the form , , is the shear modulus of the upper body. We focus on the quasi-static regime, i.e. excluding inertia, because nucleation generically takes place in this regime. The quasi-static elastic transfer functions and , see Supporting Information (Geubelle and Rice, 1995), contain all of the information about the system’s geometry, the elastic properties of the sliding bodies and loading conditions (e.g. velocity vs. stress boundary condition). The perturbation in the frictional resistance takes the form , where we used , and the definitions , and (the latter two are evaluated at ), see Supporting Information. Note that can be both positive (velocity-weakening friction) and negative (velocity-strengthening friction) depending on the materials, the sliding velocity and physical conditions (e.g. temperature) (Bar-Sinai et al., 2014). For the small slip velocities regime of interest here we assume that friction is velocity-weakening, hence we consider .

The linear perturbation spectrum is determined by the perturbation of the constitutive relation, which reads

(4) |

Substituting the results for , and , we obtain an equation for

(5) |

Once solutions are obtained, instability is implied whenever , corresponding to an exponential growth of perturbations. Consequently, is determined as the largest wavenumber (smallest wavelength) for which and the critical nucleation length is estimated as .

Solutions to Eq. (5) for some cases are available in the literature. Most notably, for two identical half-spaces we have and (see Supporting Information), where the latter represents the absence of a bimaterial effect for elastically identical materials of the same shape/geometry. Plugging these transfer functions into Eq. (5), one can readily obtain a known result for the critical wavenumber (Rice and Ruina, 1983), which reads . Using our proposed criterion , we obtain

(6) |

where was defined in Eq. (3). This prediction for the critical nucleation length is identical to the one in Eq. (3), which basically follows from dimensional considerations, once the pre-factor is identified as done above (and the definitions of and are recalled). This value of the pre-factor is times larger than the largest value we have been able to trace in the available literature based on the conventional approach, hence we conclude that the proposed approach predicts a significantly larger nucleation length for identical half-spaces as compared to the conventional approach. Indeed, some numerical simulations of earthquake nucleation indicated that the conventional prediction with quite significantly underestimates the observed (Lapusta and Rice, 2003).

The physical picture of nucleation developed in this section suggests that the origin of nucleation is a linear frictional instability, while the outcome of nucleation is typically strongly nonlinear. In particular, the critical nucleation conditions coincide with the onset of linear instability when the patch size reaches , then the slip velocity increases exponentially in the linear regime until nonlinearities set in when the slip velocity is large enough. Finally, the patch breaks up into propagating rupture fronts. The linear stage of the instability is expected to be rather generic, and in particular nearly independent of the exact functional form of (with and ) within the rate-and-state constitutive framework and of the background strength of the fault quantified by the initial age , while the nonlinear stages that follow may depend on the details of the constitutive relation and the background fault strength.

These generic properties of the onset of nucleation will be explicitly demonstrated in Sect. V below. Furthermore, we note that the works of Rubin and Ampuero (2005); Ampuero and Rubin (2008) apparently focus on the nonlinear stages of nucleation, which is consistent with the fact that they find differences between different friction laws and that their patches can shrink/expand during the nonlinear evolution of the instability. The nonlinear stages – on the route to rupture propagation – cannot take place, though, if the patch does not reach first the size determined by the linear instability. Hence, we believe that the above defined is the relevant nucleation length, and not any other length that might characterize the nonlinear evolution of the instability.

## Iv Application to bimaterial interfaces

The general framework laid down in the previous section, unlike the conventional approach, can be naturally applied to bimaterial interfaces. We consider then two half-spaces made of different elastic materials, the upper half-space is characterized by a shear modulus and Poisson’s ratio and the lower half-space by a shear modulus and Poisson’s ratio . It corresponds to Fig. 1 (left), once the limits and are taken. Defining and (i.e. the shear modulus of the upper body is denoted by , as before), the elastic transfer functions for this bimaterial system take the form (Rice et al., 2001) (see also Supporting Information)

(7) |

where

(8) |

plays the role of an effective bimaterial modulus, which approaches in the identical materials limit, and . , which appears in but not in , vanishes in the identical materials limit (and consequently vanishes in this limit as well) and hence it quantifies the bimaterial effect.

The presence of a bimaterial contrast, , introduces a new destabilization effect associated with a coupling between slip and normal stress perturbations, in addition to the the destabilizing effect associated with velocity-weakening friction, . Hence, on physical grounds one expects to decrease with increasing bimaterial contrast. To test this, we insert of Eq. (7) into Eq. (5) and calculate , obtaining the following expression for

(9) |

The first multiplicative contribution on the right-hand-side, , is obtained by replacing in our result in Eq. (6) by the effective modulus . A similar replacement has been proposed by Rubin and Ampuero (2007) in the context of a different heuristic estimate of the critical nucleation length for bimaterial interfaces. Consequently, we plot in Fig. 2 of Eq. (9), normalized by , as a function of for various values of . It is observed that for bimaterial interfaces is generically smaller than the conventional estimate , indicating that bimaterial interfaces may be more unstable than previously considered. We note in passing that Eq. (9) remains valid also in the presence of velocity-strengthening friction, , for which it predicts that for sufficiently strong bimaterial contrasts, , instability is implied even for velocity-strengthening friction (Rice et al., 2001).

## V Application to Finite-size systems and comparison to inertial Finite-Element-Method calculations

The general framework laid down in section III, unlike the conventional approach, can be naturally applied to finite-size systems. To demonstrate this, we consider here a system that features both finite dimensions and a bimaterial contrast. In particular, we consider a long deformable body of height , and of elastic constants and , in rate-and-state frictional contact with a rigid substrate under the application of a compressive stress and a shear stress . This configuration corresponds to Fig. 1 (left), once the limit is taken. In this case, the elastic transfer functions appearing in Eq. (5) take the form (see Supporting Information)

(10) |

When substituted in Eq. (5), we obtain a complex equation which is not analytically tractable, but rather is amenable to numerical analysis. Let us denote the solution by and the corresponding prediction for the critical nucleation length by .

Equation (5), with of Eq. (10), does admit an analytic solution in the limit , i.e. when the system height is small compared to field variations parallel to the interface characterized by a lengthscale . In this limit, we find and . Using these in Eq. (5), we obtain

(11) |

predicts the small behavior of and constrains any numerical calculation of to be quantitatively consistent with it in this limit. In addition, it is fully consistent with the results of Bar-Sinai et al. (2013). We numerically calculated for the following set of parameters: GPa, , , , , MPa, m, and m/s (the latter corresponds to an applied shear stress ). The result is plotted in the main panel of Fig. 3 (solid line). When of Eq. (11) is superimposed on it (dashed line), perfect agreement at small and significant deviations at larger are observed, as expected.

Our goal now is to quantitatively test the ability of the calculated to predict the critical nucleation length in a realistic situation in which a slowly expanding creep patch spontaneously nucleates accelerating slip. We would also like to test the theoretical prediction that is nearly independent of the specific friction law (in particular the aging vs. the slip evolution laws) and the background fault strength (the initial value of ). To these aims, we performed inertial Finite-Element-Method (FEM) calculations that are directly relevant for the geometrical configuration and material parameters discussed in the last two paragraphs. In particular, we consider a deformable body of height which is also of finite extent in the direction parallel to the interface and which is loaded (by an imposed velocity m/s that is initiated at ) at its lateral edge (defined as ), rather than at its top edge at , see Fig. 1 (right). The advantage of this sideway loading configuration is that it naturally generates a creep patch that slowly expands from along the interface, cf. the inset of Fig. 3. The interface is first described by the aging rate-and-state constitutive relation with and , where and the other parameters are as above. The initial conditions are and . The full constitutive relation used in the FEM calculations, which also allows a transition from stick () to slip (), can be found in the supporting information Hecht (2012).

Our theoretical approach predicts that the creep patch loses its stability and develops accelerating slip upon reaching a certain critical length. This is indeed observed in the inset of Fig. 3, where the slip velocity blows up when the creep patch reaches a certain length. We then measured the critical length in inertial FEM calculations for different system heights (in addition to the inset of Fig. 3, see also the supporting information for the details of the determination of in the numerical calculations) and superimposed the results for the aging law (red circles) on the main panel of Fig. 3. It is observed that the theoretical prediction for the critical nucleation length is in reasonably good quantitative agreement with the FEM results for the full range of system heights . This major result lends serious support to the approach developed in this Letter.

In order to test whether the theoretically predicted critical nucleation length is indeed nearly independent of the details of the friction law, we repeated the above described FEM calculations for m and m with the slip law instead of the aging law; that is, we used (the full constitutive relation can again be found in the supporting information). The resulting critical nucleation length (black triangles in main panel of Fig. 3) for both values exhibits only a small variation (less than ) compared to the results for the aging law. Furthermore, we repeated the above described FEM calculations for the aging law with m, except that we increased the initial age of the fault by three orders of magnitude, from s to s. The resulting critical nucleation length (brown square in main panel of Fig. 3) exhibits only a small variation (less than ) compared to the result for s. These results lend strong support to the idea that the critical nucleation length is determined by a linear instability that is reasonably predicted by the procedure developed in this Letter.

## Vi Concluding remarks

In this Letter we developed a theoretical approach for the calculation of the critical nucleation length for accelerating slip . The proposed approach builds on existing literature by adopting the view that nucleation is associated with a linear frictional instability of an expanding creep patch. It deviates from the conventional approach in the literature by replacing the problem of the stability of a spatiotemporally varying creep patch by an effective homogeneous sliding linear stability analysis for deformable bodies, rather than invoking a spring-block stability analysis supplemented with some fracture mechanics estimates for deformable bodies. The quality of the predictions emerging from the proposed approach therefore depend on the degree by which the creep patch can be approximated by spatially homogeneous fields. This approximation is expected to be reasonable in many cases in light of the weak/logarithmic velocity dependence of friction in many materials. The temporal aspects of the creep patch propagation are taken into account by the requirement that it becomes unstable upon attaining a length for which an unstable mode from the homogeneous linear stability analysis can be first fitted into.

The proposed approach is rather general and applies to a broad range of physical situations. For sliding along rate-and-state frictional interfaces separating identical elastic half-spaces, it has been shown to predict a significantly larger nucleation length compared to the conventional approach. For sliding along rate-and-state frictional interfaces separating different elastic half-spaces, the proposed approach has been shown to predict a bimaterial weakening effect which appears to be stronger than previously hypothesized, resulting in a smaller nucleation length. Finally, the proposed approach has been applied to finite-height systems. For this case, the scenario of a loss of stability of an expanding creep patch has been directly demonstrated using inertial FEM calculations and the predicted nucleation length has been shown to be in reasonably good quantitative agreement with direct FEM results for a range of system heights. The quality of the theoretical predictions has been shown to be nearly independent of the specific friction law used (aging vs. slip laws) and the background strength of the fault. These results offer a theoretical framework for predicting rapid slip nucleation along faults and hence may give rise to better short-term earthquake prediction capabilities. The proposed approach can and should be quantitatively tested in a wide variety of interfacial rupture nucleation problems, using both theoretical tools and extensive numerical simulations.

###### Acknowledgements.

E. B. acknowledges support from the Israel Science Foundation (Grant No. 295/16), the William Z. and Eda Bess Novick Young Scientist Fund, COST Action MP1303, and the Harold Perlman Family. R. S. acknowledges support from the DFG priority program 1713. We are grateful to Eric Dunham, one of the reviewers of the manuscript, for his valuable and constructive comments and suggestions. We also thank Robert Viesca for useful discussions in the context of nucleation on bimaterial interfaces. M. A. acknowledges Yohai Bar-Sinai for helpful guidance and assistance. The analytical formulae and numerical methods described in the main text and supporting information are sufficient to reproduce all the results and plots presented in the paper.Supplemental Materials for:

“Critical nucleation length for accelerating frictional slip”

## Introduction

The goal of this Supporting Information file is to provide additional details regarding the derivation of the theoretical results (Text S1) and Finite-Element-Method (FEM) calculations (Text S2) appearing in the main text.

## Text S1: Derivation of theoretical results appearing in the main text

### 1. The transfer functions and in plane-strain elasticity

We consider a two-dimensional (2D) elastic body that occupies the region and under plain-strain conditions. The bottom boundary of the body at is in contact with another body and hence some boundary tractions are generated at the interface. The latter are described by the interfacial stress components (constituting a vector where ). Our first goal would be to explicitly calculate the quasi-static (i.e. excluding inertia) relation between and the interfacial displacement in the form , i.e. to calculate the matrix . The body satisfies quasi-static equilibrium (momentum balance) and linear elasticity (Hooke’s law) (Landau and Lifshitz, 1986)

(S1) |

where is the infinitesimal strain tensor (not to be confused with the slip displacement discontinuity vector ), is Cauchy’s stress tensor, and and are the shear modulus and Poisson’s ratio, respectively.

At the top boundary the body is loaded by imposing a horizontal velocity and a compressive normal stress (with ). The homogeneous solution consistent with these boundary conditions reads

(S2) |

Since the equations of motion are linear, one can decompose a general solution to a sum of the steady solution of homogeneous sliding and a deviation from it, and write . The boundary conditions (BC) at are

(S3) |

Consider now a single Fourier mode, i.e. assume that all fields depend on and as , for which Eq. (S1) admits a solution of the form

(S4) |

where are (yet) unknown amplitudes which are determined by employing boundary conditions, and . Note that Eq. (S4) remains valid for plain-stress conditions, but with . Applying two boundary conditions at , cf. Eq. (S3), the solution takes the form

(S5) |

where , and the two amplitudes and remain unspecified (they depend on the contact interactions at the interface, which remain unspecified). Using Hooke’s law and momentum balance, cf. Eq. (S1), one can express the relation between the interfacial displacements and the interfacial stresses , for any and , in the form (Geubelle and Rice, 1995), where

(S6) |

Note that in case that the body under consideration occupies the region in space (with ) the analysis remains valid, but should be replaced by . This simply amounts to changing the sign of the diagonal entries of in Eq. (S6). In the limit , the matrix takes the form

(S7) |

The matrix characterizes a single body. Next, we aim at calculating the response of a general composite system, composed of two bodies made of different linear elastic materials of different heights in frictional contact (both bodies are assumed to be infinite in the -direction). The upper body, denoted by the superscript (1), is assumed to occupy the region and the lower one, denoted by the superscript (2), the region (with positive ). Since and may be different for the different bodies, they are also labeled with a superscript. Following the previous derivation, the displacements at the frictional interface are given as

(S8) |

Since both and are continuous at , the displacement discontinuity can be written as

(S9) |

Since we exclude opening gaps , but allow a slip discontinuity (cf. the main text), we obtain the following relation between and (Geubelle and Rice, 1995)

(S10) |

We thus write

(S11) |

as in Sect. 3 in the main text, with the definitions and . Note that in systems with reflection symmetry with respect to the interface is obtained from simply by taking . Therefore, the diagonal terms of and have opposite signs and the off-diagonal terms are identical, or in other words, is diagonal. Thus, is also diagonal, i.e. , and no coupling exists between tangential motion and normal traction in this case (i.e. for symmetric systems) (Ranjith and Rice, 2001). For example, in the case addressed in Sec. 3 of the main text we have , , and for these in Eqs. (S7) and (S10), one obtains (Rice et al., 2001)

(S12) |

### 2. Application to bimaterial interfaces

In Sec. 4 in the main text we apply our approach to the case of two half-spaces made of different materials. In this case, one can use Eq. (S7) for . For the ease of notation we define and , leading to (Rice et al., 2001)

(S13) |

This leads to

(S14) |

where

(S15) |

which identify with Eqs. (7)-(8) in the main text.

### 3. Application to finite height systems

In Sec. 5 in the main text a finite height body in frictional contact with an infinitely rigid substrate was considered. The requirement that the bottom layer is infinitely rigid leads to the boundary condition at the interface. In addition, since we consider here a constant shear stress at the top boundary, rather than a constant tangential velocity, the BC’s in Eq. (S3) are replaced by

(S16) |

### 4. The linear stability spectrum and finding the critical wavenumber

The transfer functions and are plugged in Eq. (4) in the main text. The latter transforms into a linear stability spectrum equation once is calculated, which is done as follows (Ruina, 1983; Rice and Ruina, 1983). As is a function of and , we have

(S19) |

where all the derivatives are evaluated at and . Next we insert and into the evolution equation , with and , and expand to leading order in to obtain

(S20) |

Using (evaluated at steady state, as described above), together with and Eq. (S20), we find

(S21) |

where we have defined , and (the latter two are evaluated at ). Substituting of Eq. (S21) into Eq. (4) in the main text, we obtain the general linear stability spectrum in Eq. (5) in the main text, which is reproduced here

(S22) |

with the extra notation which is omitted in the main text. Equation (S22) is then used to find the critical wavenumber. To that aim, we look for solutions to Eq. (S22) of the form , where is a real frequency (Rice and Ruina, 1983). Inserting this into Eq. (S22) and using the fact that all the quantities in the equation are real now, we split the equation into its real and imaginary part as follows