Spatiotemporal dynamics of frictional systems: The interplay of interfacial friction and bulk elasticity

Frictional interfaces are abundant in natural and engineering systems, and predicting their behavior still poses challenges of prime scientific and technological importance. At the heart of these challenges lies the inherent coupling between the interfacial constitutive relation -- the macroscopic friction law -- and the bulk elasticity of the bodies that form the frictional interface. In this feature paper, we discuss the generic properties of the macroscopic friction law and the many ways in which its coupling to bulk elasticity gives rise to rich spatiotemporal frictional dynamics. We first present the widely used rate-and-state friction constitutive framework, discuss its power and limitations, and propose extensions that are supported by experimental data. We then discuss how bulk elasticity couples different parts of the interface, and how the range and nature of this interaction are affected by the system's geometry. Finally, in light of the coupling between interfacial and bulk physics, we discuss basic phenomena in spatially-extended frictional systems, including the stability of homogeneous sliding, the onset of sliding motion and a wide variety of propagating frictional modes (e.g. rupture fronts, healing fronts and slip pulses). Overall, the results presented and discussed in this feature paper highlight the inseparable roles played by interfacial and bulk physics in spatially-extended frictional systems.


I. INTRODUCTION
The dynamics of frictional systems have been the focus of intense scientific research in the last few decades [1][2][3], with applications ranging from atomic and nano-scale systems [4][5][6], through engineering and control systems [7,8], to geophysical systems [9][10][11][12]. Despite the progress made in understanding and controlling friction at different scales [13][14][15][16], various basic aspects of frictional systems are not yet well-understood. As frictional systems and processes span a wide range of spatiotemporal scales, the problem can be approached from different perspectives. The focus of the present feature paper is on spatially-extended frictional systems approached from a macroscopic (continuum) perspective. Our major goal is to discuss the generic properties of macroscopic interfacial constitutive relations, within the rate-and-state friction framework [9,14,[17][18][19][20][21], and to demonstrate that the overall dynamics of frictional systems emerge from the inherent coupling between interfacial tribology and the bulk elasticity of the bodies forming the interface. As such, the physics and dynamics of frictional systems are highlighted to result from the inherent interplay between the local, micro-scale tribological behaviour of frictional interfaces and the non-local, long-range interaction between different parts of the interfaces, which are mediated by macro-scale bulk elasticity (cf. Fig. 1).
To see the nature and origin of this inherent coupling let us consider two bodies in frictional contact and the associated slip displacement δ, which is the difference between the displacement of the upper and lower bodies u ± , respectively, at the interface (cf. Fig. 1). In spatially-extended frictional systems, δ generically varies along the interface, and hence should be in fact represented by a field δ(x, t), where x is the coordinate along the interface that is located at y = 0 and t is the time (cf. Fig. 1 and note that we do not consider the z direction, which is nevertheless marked in the figure). A coarse-grained constitutive law relates the frictional stress τ (x, t) (sometimes termed the frictional strength/resistance) to quantities such as the slip displacement field δ(x, t), the slip velocity field v(x, t) = ∂ t δ(x, t), the interfacial normal (compressive) stress σ(x, t) and possibly other interfacial fields, and can be expressed as Here the dots represent yet unspecified internal state fields, which satisfy their own evolution equations. The latter, together with the functional f [·] that appears as a proportionality factor ("friction coefficient") between σ(x, t) and τ (x, t), define the coarse-grained constitutive relation. The purely local nature of f [·] is manifested by the absence of spatial derivatives of δ(x, t) in it. When bulk elasticity is not taken into account, i.e. when the bulks are considered infinitely rigid as in spring-block models [14,21,22], Eq. (1) takes the form τ 0 = f σ 0 , where τ 0 and σ 0 are the far-field applied shear stress and normal (compressive) stress, respectively (cf. Fig. 1). In this case, the frictional problem has a single degree of freedom, namely δ(t). A major goal of this feature paper is to understand how spatial degrees of freedom, emerging from bulk elasticity, change/extend the relation τ 0 = f σ 0 and to demonstrate their far-reaching implications for frictional dynamics. When the elasticity of the bodies forming the interface is taken into account, Eq. (1) serves as a boundary condition that relates two components of the stress tensor field σ ij (x, y, t) as the interface is approached, y → 0. In particular, τ (x, t) equals the bulk shear stress, σ xy (x, y, t), and σ(x, t) equals the bulk normal (compressive) stress, −σ yy (x, y, t), in the limit y → 0. The interfacial shear stress and normal (compressive) stress, σ xy (x, y = 0, t) and −σ yy (x, y = 0, t), depend on the far-field applied shear stress and normal (compressive) stress, τ 0 and σ 0 respectively, and on stresses that emerge from all processes that are mediated by the presence of the bulk of the bodies forming the interface. The latter depend on gradients of the slip displacement field, ∂ x δ(x, t), and on the slip velocity field, v(x, t) = ∂ t δ(x, t), and can be expressed as F τ [∂ x δ(x, t), ∂ t δ(x, t)] and F σ [∂ x δ(x, t), ∂ t δ(x, t)], respectively. In case the relevant bulk physics is faithfully represented by linear elastodynamics, as is widely assumed elsewhere and as is also adopted here, the interfacial stresses take the form The linear functionals F τ [∂ x δ(x, t), ∂ t δ(x, t)] and F σ [∂ x δ(x, t), ∂ t δ(x, t)] generically describe long-range spatiotemporal interactions (i.e. note the explicit appearance of spatial gradient in ∂ x δ(x, t), which schematically represents also higher order spatial derivatives and even a spatial integral) and interface-bulk radiation effects. They are affected by dimensionality, by the geometry of the bodies forming the interface and by their constitutive parameters (cf. Fig. 1). Using the equalities τ (x, t) = σ xy (x, y = 0, t) and σ(x, t) = −σ yy (x, y = 0, t), Eq. (1) takes the form which sets the framework for the discussion in this paper. This basic equation clearly demonstrates that frictional dynamics, which correspond to solutions for δ(x, t), inherently depend on both the interfacial constitutive relation f [δ(x, t), ∂ t δ(x, t), ...] (along with the internal state fields and their own evolution equations, to be discussed below), and Here we plot the general case in which the bodies are made of different materials (hence the different colors used) and have different geometries. The system is loaded normally and tangentially (σ 0 and τ 0 , respectively) far from the interface. Stress conditions at the interface, which give rise to frictional motion, result from the interplay between local frictional rheology (see zoom in) and forces mediated by the elastic bulks, induced by the loading and by differential sliding in distant parts of the interface (see curved double-arrow). The height H and length L of one of the bodies are marked and a coordinate system is introduced. The local slip displacement δ, the difference between the local displacement of the upper/lower bodies (u ± ), is also schematically denoted.
on the bulk constitutive contributions F τ [∂ x δ(x, t), ∂ t δ(x, t)] and F σ [∂ x δ(x, t), ∂ t δ(x, t)]. Consequently, the dynamics of frictional systems inseparably emerge from the coupling between interfacial and bulk physics, as encapsulated in Eq. (3). In the following parts of the paper we first discuss the generic properties of the interfacial constitutive relation f [δ(x, t), ∂ t δ(x, t), ...], including internal state fields and their evolution equations, and then consider the bulk constitutive contributions F τ [∂ x δ(x, t), ∂ t δ(x, t)] and F σ [∂ x δ(x, t), ∂ t δ(x, t)], for various dimensionalities, physical situations and geometric configurations relevant for a broad range of frictional systems. While in so doing we build on extensive effort and progress accumulated in the literature over the last few decades, we are also strongly biased toward our own recent work. As such, this paper inevitably presents our personal perspective on the issues at stake. Yet, we hope that the emerging physical picture is sufficiently broad and general to be useful for scientists interested in a wide variety of frictional systems. Whenever possible, we make an effort to highlight basic concepts and generic physical effects, while minimizing the technical complexity of the presentation and referring the readers to a more technically involved literature.

II. GENERIC PROPERTIES OF INTERFACIAL CONSTITUTIVE LAWS
A widely-used class of interfacial constitutive laws, termed rate-and-state friction (RSF) models, originated from the seminal works of Dieterich [17,18,23] and Rice and Ruina [20,24] in the late 1970's, which elaborated on ideas formulated by Bowden and Tabor [1] and Rabinowicz [25] in the 1950's. These models introduce one or more coarse-grained internal state variables to account for the evolution of the interfacial structural state (and hence they carry memory of the interfacial history), and rate dependence to account for interfacial rheology. In this section, we will briefly describe the conventional RSF model, discuss what we view as its power and possible shortcomings, and consequently propose modifications. Our ultimate goal is to provide an overview of the generic properties of interfacial constitutive laws.

A. Conventional RSF models
It is now well-established that a crucial physical quantity affecting the dynamics of spatially-extended dry frictional interfaces is the real contact area A r , which is typically orders of magnitude smaller than the nominal (apparent) contact area A n [14,[26][27][28][29]. That is, a spatially-extended interface separating two macroscopic bodies in frictional contact is composed of a collection of discrete contact asperities (i.e. it is a multi-contact interface, cf. Fig. 1) whose total area is typically much smaller than the nominal overlap of the bodies. A major implication of this physical picture is that the frictional stress τ can be written, following Bowden and Tabor [1], as τ = (A r /A n )τ s (v), where τ is the frictional stress of Eq. (1), τ s (v) is a rate-dependent asperity-level shear stress and v is the interfacial slip velocity. The normalized real contact area, A ≡ A r /A n , is commonly expressed in terms of another state variable φ, of time dimensions, as Here σ is the interfacial normal stress of Eq. (1), σ H is the material hardness [14,27,30] and b is a dimensionless material parameter, typically of the order of 10 −2 . In this formulation, φ * is a rather arbitrary time scale required for dimensional consistency. Below we discuss this point in detail. Equation (4) properly describes the outcome of a broad range of experimental works that have demonstrated that the real contact area (or the static frictional resistance, which is proportional to A according to the Bowden and Tabor relation τ = A(φ)τ s (v)) grows logarithmically with the time of contact in the absence of sliding [27,[31][32][33][34][35]. This is the well-known logarithmic aging of the static frictional resistance. Thus, Eq. (4) allows one to write the simple relation ∂ t φ = 1 in the absence of sliding, v = 0. Within this physical interpretation, φ is the typical lifetime of a contact, and therefore, if the typical slip distance required to destroy an ensemble of contacts is D, we expect φ(v) = D/|v| under persistent (steady-state) sliding, v = 0. This implies that A decreases logarithmically with increasing sliding velocity v, in agreement with the prevalent observation that the steady-state frictional resistance (to be discussed below), for a broad range of velocities, decreases logarithmically with v [17,[36][37][38][39]. To actually show this, one should consider τ s (v), which we do next.
The rate-dependent contribution to the frictional resistance, τ s (v), is usually associated with a stress-biased thermal activation process [14,38,[40][41][42][43]. Within a simplified model of a single energy barrier E, linearly biased by the applied stress according to E = E 0 − Ω τ , where Ω is an activation volume and E 0 is the bare energy barrier, τ s (v) takes the form [40,43] Here k B is Bolzmann's constant, T is the temperature and v c is a velocity scale setting the upper limit for the validity of the thermal activation model. Combining Eqs. (4) and (5), we obtain In the limit where E 0 k B T , which is quite typical, it can be shown that this equation takes the well-known form of RSF friction [40,43] where f (φ, v) = τ (φ, v)/σ as in Eq. (1), higher order logarithmic terms have been omitted, and the dimensionless parameters f 0 , α and β can be expressed in terms of physical parameters defined above (see Appendix A). The second logarithmic contribution on the right hand side of Eq. (7) emerges from the evolution of the real contact area, i.e. the state-dependent contribution of Eq. (4). The first logarithmic contribution, on the other hand, emerges from the interfacial rheology, i.e. the rate-dependent contribution of Eq. (5). This logarithmic v-dependent behavior is directly supported by experiments in which the slip velocity v changes abruptly, such that φ does not change appreciably, and the change in the frictional stress is measured on a short time scale. Under steady-state conditions, φ is slaved to the slip velocity according to φ = D/|v|, and the overall rate-dependence of f (φ = D/|v|, v) = f ss (v) is proportional to log(v), as mentioned above. This logarithmic dependence is velocity-weakening for α < β and velocitystrengthening for α > β. This issue will be further discussed in Sect. II B 2. Finally, the RSF model is fully defined once a dynamical evolution equation for φ, which connects the no-sliding φ = t and sliding φ = D/|v| fixed-points, is provided. This is commonly done through the popular "aging law" [9,23] which together with Eq. (7), constitute the conventional RSF model. Other φ evolution equations,φ = Φ(φ|v|/D) with other functions Φ(·), are considered and studied in the literature [20,[44][45][46], but are not discussed here.
B. Limitations of conventional RSF models and proposed modifications RSF models have been quite extensively used in macroscopic modeling of frictional phenomena. Yet, we believe that while these models are quite successful in capturing many physical aspects of frictional phenomenology, they are also deficient in various important ways. Our goal here is to highlight these limitations and to propose modifications to be incorporated into an extended friction model.

Linear reversible response and an extended RSF model
Conventional RSF models predict the frictional stress, cf. Eq. (6), which should provide physically meaningful information both in the presence of sliding (v = 0) and in its absence (v = 0). Using conventional terminology, it should describe both static and dynamic friction. However, Eq. (6) predicts τ (φ, v = 0) = 0 in the limit v → 0, i.e. it predicts vanishing static friction. This happens because τ s (v) in Eq. (5) corresponds to a Newtonian fluid behavior at small sliding velocities, i.e. τ s (v → 0) ∝ v. One way in which this problem is circumvented in the literature is to change the v-dependent term in Eq. (7) to be α log(1 + v/v c ), such that it vanishes for v = 0, together with interpreting the remaining part, f 0 + β log(t/φ * ) (note that φ = t for v = 0), as an upper threshold. That is, in this case f (φ, v = 0) does not correspond to the frictional stress, but rather to the threshold for the onset of sliding. Mathematically, it means that f (φ, v = 0) is used to define an inequality imposed on a frictional stress that is calculated from other considerations (in fact, from the bulk shear stress σ xy of Eq. (2)). At the same time, f (φ, v) of the very same Eq. (7), is interpreted as the actual frictional stress for v = 0, implying that f (φ, v) has a different physical meaning for v = 0 and for v = 0.
Another common way to circumvent the vanishing static friction problem is by forcing the frictional interface to slide with a small background velocity v = 0 and use Eq. (7) as is [45,47]. While these ad hoc solutions to the stated problem might provide a practical way out under some circumstances, from a basic physics perspective we expect τ (φ, v), or f (φ, v), to continuously describe the frictional resistance stress under both static and dynamic conditions. In particular, it should properly account for the transition between the two behaviors, e.g. for the onset of sliding and for sliding direction reversal (which is relevant for cyclic frictional motion, appearing in many engineering systems).
In our view, these problems originate from the fact that the conventional RSF model does not feature a linear reversible response at low stresses, corresponding to the elastic (reversible) deformation of contact asperities (note that, as explained above, the model does feature a Newtonian viscous behavior, i.e. a linear irreversible response). The existence of a reversible response has been experimentally demonstrated [22,[48][49][50][51], see also Fig. 2a here, and is expected on general grounds [52]. That is, we expect the asperity-level stress τ s to be finite and reversible under small static shear loading conditions due to the elastic deformation of contact asperities. This should result in τ (φ, v) that is valid for both v = 0 and v = 0, as well as for v that changes sign dynamically, which should resolve the problems with conventional RSF models stated above. Indeed, some authors have incorporated an interfacial elastic response in their modeling efforts [51,[53][54][55][56][57][58][59][60][61] and noted its potential importance for various frictional phenomena, such as the onset of sliding [22] and rupture mode selection [55]. Here we describe an extended RSF-type model that incorporates a linear reversible response at low stresses. The model has been briefly introduced and studied in previous publications [58,60,61].
The basic idea is to incorporate into the RSF framework an elastic interfacial response at small stresses, within a Kelvin-Voigt-like viscoealstic model [62,63] together with threshold dynamics. To that aim, we reinterpret the asperity-level stress in Eq. (5) as a viscous-like stress τ vis s (φ, v), which is a partial stress constituting only one contribution to the total asperity-level stress τ s . The complementary contribution is the asperity-level elastic stress τ el s , such that τ s = τ vis s + τ el s . While this way of introducing asperity-level elasticity into the friction model is not unique, it suffices for our present purposes, as will be shown below. Using the Bowden and Tabor relation τ = A(φ)τ s (v) to translate the asperity-level stress to a coarse-grained stress, we obtain where, as explained above, τ vis (φ, v) = A(φ)τ s (v) with τ s (v) of Eq. (5). In Eq. (9), τ el is a dynamic variable by itself, in addition to φ and v. Therefore, we need to derive its evolution equation. As explained above, at small slip displacements we expect the interface to respond linear elastically. To capture this basic physical behavior, we treat the interface as a boundary-layer of height h, much smaller than the macroscopic size H of the system, which is characterized by an effective shear modulus µ 0 . This effective modulus may be different from the shear modulus µ of the bulk due to the existence of stress-free surfaces along the interface (i.e. surface roughness). Under deformation, the total shear accumulated in the boundary layer is δ/h, producing a shear stress µ 0 δ/h. Such an elastic response of the interface has been directly measured [48,50,64], as shown in Fig. 2a, and has been predicted theoretically [52]. Defining a coarse-grained elastic stress τ el through A n τ el = A r µ 0 δ/h and differentiating it with respect to time, we obtainτ el = µ 0 A v/h (where we assumed that A varies much slower than τ el such that its time-derivative can be neglected). The term µ 0 A v/h constitutes only one contribution toτ el ; the other contribution is discussed next.
Irreversible sliding drives contact annihilation and renewal, and hence is accompanied by a reduction in the average elastic stress stored in the contact asperities. For simplicity, we assume a linear form for this reduction, which is similar to that of φ in Eq. (8). Thus, the evolution equation of τ el is taken to bė where the function g(τ, v) controls the transition between the reversible and irreversible regimes. In fact, the very same function should appear in a modified version of Eq. (8), which takes the formφ = 1 − φ|v|g(τ, v)/D and is adopted here. When g 1, the system responds elastically, i.e. τ el = µ 0 Aδ/h. When g(·) ∼ O(1), irreversible stress and area reduction take place. Therefore, g(·) plays the role of an effective threshold for the onset of irreversible slip, as in conventional elasto-plastic models that feature a yield stress. Explicit expressions for g(·) will be discussed below in relation to various applications.
To directly demonstrate the physical relevance and importance of interfacial elasticity, we focus on the experimental results of [50], which are reproduced here in Fig. 2a. In these experiments, the multi-contact interface between two PMMA blocks has been subjected to two time-dependent shear loading-unloading protocols below the 'static' threshold for gross slip (which for this experimental system corresponded to a shear force of F static = 10N). In the first protocol (the squares in Fig. 2a, equally spaced in time), the shear force F S (t) is externally increased at a constant (and small) rate up to 4N (i.e. about 40% of the nominal 'static' friction threshold) and then decreased back to zero at the same small rate. The observed response, i.e. the slip displacement δ(t), is nearly linear in the driving force F S (t), featuring very little dissipation (i.e. hysteresis) and almost full recovery of the initial state. These observations explicitly demonstrate the existence of a linear reversible (elastic) response. In the second protocol (the circles in Fig. 2a, equally spaced in time), the shear force F S (t) is externally increased at a constant (and small) rate up to 9N (i.e. about 90% of the nominal 'static' friction threshold), it is then maintained at this level for sufficiently long time and eventually it is decreased back to zero at the same small rate. Under this load-hold-unload protocol, the slip displacement δ(t) is no longer linear in the driving force F S (t) for values larger than 7N. Moreover, when F S (t) is kept at its maximal value (here of 9N), δ(t) exhibits creep dynamics that may indicate arrest after a finite time (as suggested by the increasing density of the circles). Finally, as the shear force is reduced back to zero, a finite residual (irreversible) slip remains. This behavior corresponds to a viscoealstic response with threshold dynamics.
The friction model discussed above, which incorporates an interfacial elastic response at small driving forces, naturally and semi-quantitatively reproduces these experimental observations. To see this, we neglect all bulk degrees of freedom, whose role has been deliberately minimized in the experiments of [50] in order to isolate the interfacial response (this is the only place in this paper where bulk degrees of freedom are neglected), and focus on the equation Here F S (t) is the applied shear force (i.e. corresponding to each one of the two protocols discussed above), F N is the constant applied normal force and A n is the nominal (macroscopic) contact area (such that σ 0 = F N /A n is the normal stress, not used here). τ (φ,δ, τ el ) of Eq. (11) is given in Eq.

Short time cutoff and the steady-state friction curve
The second problem we identify with conventional RSF models is related to Eq. (4). Imagine that a fresh interface is formed by bringing two bodies into frictional contact. In the absence of irreversible slip, for whichφ = 1, the age φ can be replaced by the time t elapsed from the formation of the interface. Equation (4) then predicts that the contact area will grow logarithmically with time, as observed in a broad range of experiments [35,37,65,66]. However, when t → 0 it becomes singular, which is clearly unphysical. Evidently, a short time cutoff exists and must be incorporated. A simple way to do so is to replace Eq. (4) by With this modification, φ * is no longer an arbitrary time scale -it is an intrinsic interfacial/material parameter which can be interpreted as a typical cutoff time for the onset of logarithmic aging. The existence of this cutoff has been directly verified experimentally [35,65,67]; it has been suggested by Dieterich [17,18], mentioned by Baumberger, Caroli and coworkers [14,34,38] and extensively discussed recently in [43,60]. While this modification might appear at first glance rather innocent, it has major implications for persistent sliding because the short time cutoff is translated, through the steady-state relation φ = D/|v|, to a high-velocity cutoff. To see this, note that Eq. (12) predicts that A decreases logarithmically with the steady sliding velocity v. The presence of a short time cutoff implies that this decrease saturates at velocities D/φ * . At higher velocities the real contact area approaches a constant, and the velocity dependence of friction is dominated by τ s (v). Since the latter is intrinsically an increasing function of v (cf. Eq. (5)), friction becomes velocity-strengthening in this range of velocities. Put differently, the short time cutoff makes the steady sliding frictional resistance a non-monotonic function of the sliding velocity, reaching a minimum at v min D/φ * . As will be thoroughly discussed below, the crossover from velocity weakening to velocity-strengthening has important implications on many aspects of the spatiotemporal dynamics of frictional systems.

FIG. 3:
The generic form of steady-state frictional resistance f ss (v) (see text and Appendix A for details). At extremely low sliding velocities (yellow-shaded), steady-state friction is velocity-strengthening. This regime, associated with v ≈ 0 in the figure, corresponds to stick conditions (effectively no-sliding) either because the slip velocity is minute or because steady-state is dynamically inaccessible due to extremely long relaxation timescales (see text for additional discussion). At intermediate velocities (purple-shaded), friction is velocity-weakening. Finally, at higher velocities (green-shaded), above a local minimum occurring at v min (corresponding to a frictional resistance τ min /σ 0 ), friction becomes velocity-strengthening again. When loaded at stresses higher than this minimal value, τ 0 > τ min (horizontal line), three homogeneous solutions (fixed-points) exist; two stable (on the velocity-strengthening branches, marked by squares) and one unstable (on the velocity-weakening branch, marked by a circle), see text for additional discussion.
We conclude this discussion of conventional RSF models and their extensions (the complete set of equations defining the constitutive law is collected in Appendix A) by stating what we perceive to be a generic picture of the steadysliding friction curve. The latter, denoted by f ss (v) and derived in Appendix A, is schematically plotted in Fig. 3. The function f ss (v) has a generic "N-like" shape, composed of three distinguishable regimes. First, at very small slip velocities, the frictional stress increases linearly with v, corresponding to a viscous-like creep motion [30]. We note that at these low velocities, dynamically reaching the true steady-state might be very difficult/long because the relaxation time scales as D/|v| (though there exist experimental measurements of this regime [68]). This velocity-strengthening branch (i.e. an increase in the frictional resistance with increasing steady-state sliding velocity) corresponds in practical terms to 'stick conditions', i.e. to situations in which the frictional interface is essentially locked-in; the reason for this is either the extremely low slip velocities in this regime or that it takes so long to reach this low velocity steady-state that the interfacial response is in fact predominantly elastic (representing an extremely long transient). This regime is terminated by a maximum of the curve, followed by a logarithmic velocity-weakening regime (corresponding to the α < β case discussed in Sect. II A in relation to Eq. (7)). Finally, the curve reaches a minimum around v min D/φ * and then increases logarithmically (i.e. follows a velocity-strengthening behavior, as discussed above).
The validity of the N-shaped steady-state friction curve of Fig. 3 is supported by experimental data for a broad range of materials [43] and its implications for the dynamics of frictional systems will be reviewed below. We note that while we consider here logarithmic velocity-strengthening above the minimum of the steady-state curve, strengthening can in fact become stronger than logarithmic, as discussed in [14,43,69]. We also note that at yet higher slip velocities (compared to those shown in Fig. 3), the heat generated by frictional dissipation might not have enough time to diffuse away from the interface, leading to thermal softening accompanied by a substantial drop in the frictional resistance [70,71]. This regime is not shown in Fig. 3 and is not further discussed in this paper. Finally, note that the velocity-weakening and post-minimum velocity-strengthening regimes of the steady-state friction curve are still approximately described by Eq. (7) once β log(φ/φ * ) is replaced by β log(1 + φ/φ * ) (following Eq. (12)), but with redefined f 0 , α and β (see Appendix A for details). For α > β, the steady-state friction curve is velocity-strengthening for all slip velocities (again, if extremely high-velocity thermal softening is excluded from the discussion). With this, we conclude our discussion of the generic properties of interfacial constitutive laws.

III. BULK ELASTICITY IN FRICTIONAL SYSTEMS
In the previous section we focused on the generic properties of the interfacial constitutive law, i.e. on f [δ(x, t), ∂ t δ(x, t), ...] in Eq. (3). Our goal here is to discuss the properties of the bulk-mediated interaction contributions to this basic equa- , represents the coupling between interfacial slip and normal stress variations along the interface. It is now well-established that this coupling exists only when reflection symmetry with respect to the frictional interface is broken [40,[72][73][74]. This symmetry can be broken by material contrast (i.e. having two blocks made of different materials, cf. Fig. 1), by geometric asymmetry (i.e. having two blocks of identical materials, but with different shapes, cf. Fig. 1) or by the loading configuration. In the first part of the discussion in this section, we assume perfect material and geometric symmetry, and perfectly anti-symmetric loading, such that The shear stress contribution generally features complex spatiotemporal properties. The spatial properties, in particular the range of bulk-mediated interaction between different parts of the interface, is determined by the system dimensions. We assume that the system length L is much larger than its height H (cf. Fig. 1 for a visual definition of L and H), and focus on the effect of the latter on where µ is the shear modulus (already defined above) and c s is the shear wave-speed (s(x, t) will be defined and briefly discussed below). The right-hand-side of Eq. (13) represents the infinite half-space linear elastodynamic Green's function [75][76][77], corresponding to the momentum balance equation ρü = ∇·σ (where ρ is the mass density and the stress tensor σ is related to displacement gradients ∇u through Hooke's law), and contains two contributions. The contribution µ∂ t δ(x, t)/2c s is 'instantaneous', i.e. it depends on the slip velocity ∂ t δ(x, t) alone, and corresponds to wave radiation from the interface into the bulks forming it [45,47,78,79]. It is known as the 'radiation-damping' term because mathematically it appears as an effective linear drag/viscous stress that the bulks apply on the interface. The contribution s(x, t) in Eq. (13) is infinitely long-ranged in space, i.e. it involves integration over the whole interface with power-law kernels, and involves time integration that is limited by causality (formally it is a functional of δ(x, t), not just a function of x and t). In general, s(x, t) does not admit an analytic real-space representation and even its spectral (Fourier) representation may be highly complex [75][76][77]. In order to nevertheless gain some physical intuition into s(x, t) and to support the statements just made about it, we consider situations in which the frictional dynamics are significantly slower compared to elastodynamic effects. That is, we consider the quasi-static limit in which the inertia of the bulks plays no role. In this quasi-static limit the radiation-damping term is also negligibly small and we have (see, for example, [80]) The infinitely long-range bulk-mediated interaction between different parts of the interfaces and the associated powerlaw kernel are evident, as is also schematically illustrated in Fig. 1. Mathematically handling the infinitely long-range interaction in the H → ∞ limit generally raises significant difficulties and in many case a variety of numerical methods are invoked, e.g. spectral techniques [75][76][77]. The spatiotemporal range of interaction encapsulated in F τ [∂ x δ(x, t), ∂ t δ(x, t)] diminishes as the system height H is reduced. A particularly useful and analytically tractable expression is obtained in the limit of small H. That is, if one assumes that H is sufficiently smaller than any other relevant lengthscale in the problem, a systematic expansion (reminiscent of the commonly used shallow water approximation in fluid mechanics [81]) yields [60,61,82] Hereμ is proportional to µ, with a prefactor that depends on Poisson's ratio [61]. In this thin-systems approximation, sometimes termed the quasi-1D limit, the spatiotemporal interaction becomes short-ranged (i.e. captured by second order differential operators, with no integral terms) and the bulk -described by a scalar wave equation with Hdependent coefficients -is effectively inseparable from the interface. Precisely for this reason, the radiation-damping term of Eq. (13) also disappears (i.e. there is no distinct bulk to radiate energy into) and F σ = 0. While the short-range spatiotemporal interaction is a limitation, the relative simplicity of the small-H formulation allows a comprehensive, quantitative and predictive analytical treatment of complex spatiotemporal frictional dynamics, in addition to extremely efficient numerical calculations [60,61,69], as will be demonstrated below.
When H is neither infinitely large nor small, interesting new effects related to the finite system height emerge, as will be discussed below in relation to various physical phenomena. Moreover, when the perfect symmetry assumption is relaxed, and reflection symmetry with respect to the frictional interface is broken by material contrast, geometric asymmetry or the loading configuration, F σ [∂ x δ(x, t), ∂ t δ(x, t)] does not vanish. In such situations, which are quite prevalent, one needs to consider both F τ [∂ x δ(x, t), ∂ t δ(x, t)] and F σ [∂ x δ(x, t), ∂ t δ(x, t)], and rich dynamical effects emerge from Eq. (3). While the analysis of finite systems can be pursued analytically in some cases, numerical methods such as the Finite Element Method (FEM) are required in many other cases. Finally, we note that finite-H systems, i.e. realistic systems, are well described by the infinite-H formulation discussed above on time scales shorter than the typical wave reflection time from outer boundaries (which is determined by H/c s ). In the next section we will see how the various forms of bulk elasticity discussed in this section couple to the interfacial constitutive law to give rise to rich spatiotemporal frictional dynamics.

IV. BASIC PHENOMENA IN SPATIALLY-EXTENDED FRICTIONAL SYSTEMS
Our goal in this section is to demonstrate, through a broad range of examples, how basic phenomena in spatiallyextended frictional systems emerge from the inherent coupling between the interfacial constitutive law (discussed in Sect. II) and bulk elasticity (discussed in Sect. III).

A. The stability of homogeneous sliding
A central question in relation to spatially-extended frictional systems concerns the stability of homogeneous sliding [19,40]. That is, suppose that two bodies in frictional contact slide at a constant homogeneous slip velocity v 0 . When this sliding motion is unstable, interesting phenomena such as squeaky door hinges [83], squealing car brakes [84] or earthquakes along geological faults [85,86] might emerge. What are the conditions for instability? How are they influenced by the interplay between interfacial friction and bulk elasticity? In this section, we scratch the surface of these important questions through two examples that highlight the interplay between interfacial friction and bulk elasticity in relation to the stability of frictional systems. Mathematically, one considers small spatiotemporal perturbations δv(x, t) on top of the homogeneous sliding state v 0 , i.e. v(x, t) = v 0 + δv(x, t), and expand all relevant equations to leading (linear) order in δv(x, t). Due to linearization, one can consider each Fourier mode separately, i.e. express the perturbation as δv(x, t) ∼ exp(Λt − ikx), where Λ is the complex frequency and k is the wavenumber. An instability emerges whenever the real part of Λ becomes positive, i.e. [Λ] > 0, which implies that perturbations exponentially grow in time.
The complex frequency Λ and the wavenumber k can be related by perturbing Eq. (3) to linear order in δv(x, t). The result takes the schematic form where δ{·} (not to be confused with the slip displacement δ(x, t)) represents a variation with respect to v.
The first example we consider in this context concerns the stability of homogeneous sliding in the small H limit in the quasi-static regime, i.e. when Eq. (15) reads F τ −μH∂ xx δ(x, t). In the small H limit there exists no coupling between slip and normal stress variations, hence F σ = 0 in Eq. (16). Finally, v 0 is assumed here to be on the velocityweakening branch of the steady-state friction curve in Fig. 3 In this case, δf in Eq. (16) is destabilizing, i.e. from the perspective of the friction law alone, perturbations are expected to be amplified. However, as will be shown next, the bulk contribution δF τ is stabilizing and consequently it generates a stability regime that depends on an emerging elasto-frictional lengthscale. Performing the analysis [19,87], one finds that homogeneous sliding under the stated conditions is unstable for 0 < k < k c , where k c = 2π/L c is determined by a critical length of the form [61] This result shows that indeed stability is enhanced (i.e. the range of instability shrinks) as the bulk combinationμH increases compared to the interfacial friction combination dfss(v) d log v /D, and that the competition between the two (along with the applied normal stress σ 0 ) gives rise to a characteristic lengthscale. These properties are demonstrated in Fig. 4a. Finally, note that a similar physical picture emerges in the H → ∞ limit, where the critical length of Eq. (17) is replaced by an H-independent length of the form L c ∼ The second example we consider concerns the stability of homogeneous sliding of two finite height bodies made of the same material, where fully inertial dynamics are taken into account (unlike the quasi-static dynamics discussed in the previous example). This time, the homogeneous slip velocity v 0 is assumed to be on the velocity-strengthening branch of the steady-state friction curve in Fig. 3, i.e. dfss(v) d log v > 0 at v 0 . In this case, δf in Eq. (16) is stabilizing, i.e. from the perspective of the friction law alone perturbations are expected to be attenuated, and δF τ is still stabilizing. Consequently, any possible instability can arise only from a destabilizing normal stress contribution δF σ . Let us denote the height of the upper body by H and that of the lower one by ηH [73]. For η = 1, i.e. for a perfectly symmetric system, we have δF σ = 0, and homogeneous sliding under the stated conditions is categorically stable [88]. This is demonstrated in Fig. 4b. Consider then the case in which η 1, i.e. there exists significant geometric asymmetry in the system, giving rise to a destabilizing normal stress contribution δF σ that competes with the velocity-strengthening stabilizing friction contribution. Depending on the values of the problem's parameters the former can even overcome the latter, implying an instability. This scenario is demonstrated in Fig. 4b, where an instability emerges at a critical geometric asymmetry level η c > 1 and a finite range of unstable modes with wavenumbers k ∼ H −1 exists for larger geometric asymmetry levels, showing that bulk geometry and elasticity can destabilize frictional sliding that is otherwise stable [40,73]. as a function of wavenumber k (both properly non-dimensionalized). When [Λ] > 0, sliding is unstable. (a) In the small-H limit, sliding on top of a rigid substrate becomes unstable when k is smaller than a critical value, determined by the critical nucleation length L c of Eq. (17). Curves are shown for two parameter sets, demonstrating that stability is promoted (i.e. the range of unstable modes shrinks) when the body is stiffer (largeμ) and when friction is less velocity-weakening (smaller σ 0 df ss /d log(v)). Details of the derivation of the linear stability spectrum can be found in the Supporting Information of [87], though similar calculations appeared earlier [19,40]. Note that the kink in the spectrum is related to its complex-plane structure. (b) For a system of two elastic bodies of different height (quantified by the height contrast η), velocity-strengthening sliding can be destabilized by geometric contrast/asymmetry alone (adapted from Fig. 5 of [73]). In the absence of geometric contrast (η = 1), sliding is categorically stable. At a critical contrast (η c ≈ 2), an instability emerges (first zero crossing). Finally, for large contrast (η 1), a range of unstable modes with k ∼ H −1 exists (see text for additional details). The insets show the geometric configuration in each panel.
There are other physical ways in which the interplay between interfacial friction and bulk elasticity affects the stability of frictional systems. For example, the combined effect of bimaterial contrast, finite geometry and inertial dynamics, which strongly affect δF σ , has been shown to play crucial roles in the development of instabilities of velocitystrengthening interfaces (i.e. when the frictional contribution δf is stabilizing [89]). In particular, it has been shown that there exists a universal instability with a wavenumber k ∼ H −1 and a maximal growth rate [Λ] ∼ f c s /H, in the presence of strong bimaterial contrast (e.g. a deformable body sliding on top of a rigid one). This universal instability has been shown to be mediated by waveguide-like modes that are intrinsically coupled to the system height H. Furthermore, in has been shown that in the limit of large systems, H → ∞, there exist -in addition to quasistatic instability modes that are relevant at small sliding velocities [40] -also fully dynamic instabilities. These are mediated by propagating plane-wave-like modes, featuring interesting directionality effects with respect to the sliding direction [89]. All in all, these examples show how the interplay between interfacial friction and bulk elasticity determines the stability of homogeneous frictional sliding.

B. The onset of sliding motion: Creep patches and their stability
The discussion in the previous subsection focused on homogeneous sliding that results from homogeneous loading conditions. However, in many engineering applications and laboratory experiments the loading on a frictional system is inhomogeneous, e.g. it is applied to the trailing edge of the bodies involved. To address this situation, we consider a body of height H in frictional contact with an infinitely rigid substrate that is being pushed sideways, say at x = 0, with a velocity v d [61,87]. The latter is assumed to be on the velocity-weakening branch of the steady-state friction curve in Fig. 3 [87]. Reasonable and non-trivial agreement between the theory, which is based on homogeneous linear stability of steady sliding, and the numerics, based on transient inhomogeneous simulations, is demonstrated (see [87] for additional discussion).
Do the creep patches propagate indefinitely as long as v d is maintained finite? The stability analysis of homogeneous sliding presented above suggests that this is not the case. While creep patches are characterized by inhomogeneous slip velocity distributions (cf. Fig. 5a), varying from v(x = 0, t) = v d to v(x > l(t), t) = 0 at any time t, the values themselves are predominantly on the velocity-weakening branch of the steady-state friction curve; overall, the physical situation is not expected to dramatically differ from that of a homogeneous patch characterized by an average slip velocitȳ v(t) = [l(t)] −1 l(t) 0 v(x, t)dx ∝ v d that belongs to the velocity-weakening branch. Consequently, the stability analysis discussed above predicts that creep patches undergo a linear instability when their size surpasses the critical (minimal) length for instability, i.e. when l(t = t c ) = L c (H). That is, when they become large enough to contain the smallest unstable mode. This prediction is tested in Fig. 5a and is fully supported by the FEM calculation, i.e. when l(t) surpasses L c (H) (for the H employed in this calculation), an instability characterized by an exponential growth of slip velocities is observed.
To test this prediction over a broad range of H values, L c (H) has been theoretically calculated using the homogeneous linear stability analysis discussed in Sect. IV A [87]. Note, though, that the present problem is more complicated than the case discussed in Sect. IV A due to the existence of bimaterial contrast, i.e. the present problem involves a deformable body sliding on top of a rigid one. The scaling properties of the solutions reported on in Sect. IV A are nevertheless preserved, i.e. L c (H) ∼ √ H for small H and L c (H) → const. for large H, as shown in Fig. 5b. The comparison between the size in which creep pathes lose their stability, measured in dynamical simulations of propagating creep patches for various H values, and the theoretical prediction for L c (H) is also shown in Fig. 5b. The dynamics of creep patches and their accompanying instabilities determine the onset of global frictional motion, occurring when slip first reaches the leading edge of the body. These rich and realistic spatiotemporal dynamics crucially depend on the interplay between interfacial friction and bulk elasticity (as well as other related analyses in the literature, e.g. [82,[90][91][92][93]), and are qualitatively different from the description of the onset of frictional motion using spring-block models that lack spatially-extended degrees of freedom. In fact, the transfer of slip to the leading edge of the body also involves dynamically propagating modes that are excited by the instabilities of creep patches (e.g. the one shown in Fig. 5a), not discussed up until now. Such propagating frictional modes are discussed next.

C. Propagating frictional modes
Propagating frictional modes play crucial roles in the dynamics of frictional systems [3,10,11,27,28,93,94]. For example, following the previous subsection, propagating frictional modes that are generated by instabilities of creep patches play important roles in the onset of global frictional motion, relevant to a wide variety of engineering systems. Our goal in this subsection is to understand how such propagating frictional modes emerge from the interplay between interfacial friction and bulk elasticity and what forms can they take. To this aim, we consider a frictional system subjected to homogeneous shear stress and normal (compressive) stress loading, τ 0 and σ 0 respectively (cf. Fig. 1).
Suppose then that τ 0 is larger than τ min ≡ σ 0 f min , where f min corresponds to the minimum of the steady-state friction curve in Fig. 3. Under these conditions, there exist 3 spatially homogeneous fixed-points (already mentioned above) that correspond to the intersection points of the τ 0 /σ 0 line with f ss (v), cf. Fig. 3. The first and last intersection points, on the velocity-strengthening branches, are stable fixed-points from the perspective of the friction law, while the intermediate one, on the velocity-weakening branch, is an unstable fixed-point (it has been discussed in the context of creep patches in Sect. IV B). Generally speaking, frictional modes that propagate at a steady-state speed c spatially connect the two stable fixed-points. The steady nature of these propagating objects imply that they can be described in a co-moving frame defined through ξ ≡ x − ct, and consequently that ∂ t = −c∂ ξ and ∂ x = ∂ ξ such that all of the relevant equations can be expressed in terms of the co-moving coordinate ξ.
We discuss below 3 types of steady-state propagating frictional modes. First, we consider modes that connect the very low-v stable fixed-point (v ≈ 0, termed the 'stick state') with the high-v stable fixed-point (v > 0, the 'sliding/slipping state'). If the latter is realized at ξ → −∞ and the former at ξ → ∞, the modes propagate from left to right with c > 0, and are characterized by a lengthscale around ξ = 0 over which the 'stick state' changes to the 'sliding/slipping state'. These propagating frictional modes are sometimes termed frictional rupture fronts. Second, when the stick state invades the sliding state, we obtain healing ('inverse') frictional fronts. Finally, when the stick state exists both at ξ → ±∞, with a spatially-compact region of v > 0 around ξ = 0, we obtain slip pulses. As will be explained below, slip pulses can be constructed from the interaction of rupture and healing frictional fronts (we note in passing that there exists a related frictional mode, composed of a periodic train of pulses [55,74,95], which is not discussed here). In all cases, propagating frictional modes involve the self-selection of both the characteristic transition length and the propagation speed c. Our goal below is to understand how they are determined by the interplay between interfacial friction and bulk elasticity.
As a first step, we show that and c satisfy a generic scaling relation. To see this, note that the time it takes for a propagating frictional mode to travel over a given point on the interface is /c. During this time, the slip accumulated at this point scales as v( ) /c ∼ v p /c, where we estimated the slip velocity in the transition region v( ) by the peak (maximal) slip velocity v p . Finally, since the accumulated slip scales with the typical slip distance D (while the accumulated slip itself is quite significantly larger than D [96]), we obtain This relation can be somewhat more formally rationalized using Eq. (8) (the same reasoning applies to Eq. (10)). Note that the relation between and c in Eq. (18) is independent of the loading conditions and dimensionality as a scaling relation, though an implicit dependence emerges through v p = v p (τ 0 /σ 0 , H), as will be discussed below. To see how these physical quantities emerge in realistic frictional systems, we first discuss steady-state propagating rupture modes. To that aim, we consider a long deformable body of height H in frictional contact with a rigid substrate, subjected at y = H to homogeneous shear and normal stress, τ 0 and σ 0 respectively (cf. Fig. 1). We then obtain steady-state propagating rupture modes using a treadmill FEM routine detailed in Appendix C. The results, demonstrating steady-state rupture modes propagating from left to right for a broad range of H values, are presented in Fig. 6. In panel (a) we present τ (ξ)/σ 0 , focusing on the transition region of size that clearly emerges near ξ = 0. It is observed that far ahead of the transition region τ (ξ) simply equals to the applied shear stress τ 0 , as expected. Far behind the transition region, τ (ξ) approaches a residual stress τ r , which identifies here with the loading shear stress τ r → τ 0 ; this is expected since the two stable fixed-points that the rupture mode is connecting feature the same stress τ 0 (cf. the intersection points marked by black squares in Fig. 3). The fact that the residual and applied shear stresses are identical for steady-state rupture modes will be further discussed below. Finally, in the transition region, τ (ξ) goes through a maximum, which we denote by τ p .
The corresponding slip velocity v(ξ)/v min is shown in Fig. 6b. The residual slip velocity left behind the propagating rupture mode, v r , corresponds to the high-v stable fixed-point in Fig. 3, τ 0 /σ 0 = f ss (v r ), while far ahead of it the slip velocity is vanishingly small, v ≈ 0 (formally also a solution of τ 0 /σ 0 = f ss (v)). In between, the two are smoothly connected over a lenthscale , where the connection is monotonic for small H and nonmonotonic for larger H, going through a maximum denoted by v p . When the transition is monotonic, we simply have for the maximal slip velocity v p = v r . In panel (c) the normalized real contact area A(ξ) is presented, directly justifying the term 'rupture front' as a low A value state ('partially ruptured contact') is observed to invade a high A value state (reference real contact area). Finally, in panel (d) we present σ(ξ)/σ 0 , which is simply unity for small H (since δF σ = 0 in this limit) and deviates from unity as H increases due to the bimaterial effect manifested in δF σ = 0 (recall that we are considering a deformable body sliding on top of a rigid one).
To gain deeper insight into and c, going beyond the scaling relation in Eq. (18), we follow the scaling analysis of [60] in the small-H limit (H ), to obtain where the so-called dynamic stress drop is defined as ∆τ p−r ≡ τ p − τ r . Note that strongly inertial effects, which emerge when c approaches the elastic wave speeds, have been neglected, that we used the fact that v p (τ 0 /σ 0 , H) = v r (τ 0 /σ 0 ) in the small-H limit, and that and c indeed satisfy the scaling relation in Eq. (18). The relations in Eq. (19) clearly demonstrate how and c emerge from the interplay between interfacial friction and bulk elasticity. The prediction ∼ √ H is tested in Fig. 6 through rescaling ξ by √ H. It is observed that for small H, but over a large range (here H varies 16-fold, from 625µm to 1cm), the fields nearly collapse on H-independent curves, verifying the ∼ √ H prediction. For yet larger H values, some deviations are observed, as will be discussed next. The scaling prediction for c in Eq. (18) will also be tested below.
With increasing H, we expect the small-H predictions in Eq. (19) to break down gradually, eventually being replaced by their large-H counterparts (H ). The procedure to obtain the latter from Eq. (19) is in principle straightforward; one should just replace H in the first relation in Eq. (19) by and solve for the latter to obtain ∼ µ D/∆τ p−r (the small-H shear modulusμ is replaced here by µ). Then the latter is used to replace H in the second relation in Eq. (19), leading to c ∼ v p µ/∆τ p−r . The only remaining issue is the H-dependence of v p (τ 0 /σ 0 , H), which is clearly observed for the largest H values used in Fig. 6b. This pronounced increase of v p with H is nothing but a signature of the build-up of the famous square root singularity of crack-like objects [97], which is absent in the quasi-one-dimensional limit of small H. Consequently, using conventional fracture mechanics in long systems (recall that we assumed L → ∞, i.e. our system features a 'strip geometry' H L) [97], the residual slip velocity v r is amplified in the transition region according to v p (τ 0 /σ 0 , H) ∼ v r (τ 0 /σ 0 ) H/ . Taken together, we obtain for the large-H limit the following predictions This analysis seems to suggest that the propagation speed c follows the scaling relation c ∼ √ H in both the small and large H limits, cf. Eqs. (19)- (20). While at present we cannot comprehensively test this prediction, simply because very large H values are computationally inaccessible, we can consider the available results of Fig. 6, where deviations from the small-H limit are observed (e.g. the lack of collapse of the curves shown in panel (a) for the largest H values used and the appearance of a significant overshoot in the slip velocity profiles shown in panel (b) there). Consequently, we plot in Fig. 7a c(H) for the results presented in Fig. 6. It is observed that the scaling prediction c ∼ √ H is accurately followed up to the largest H value considered here. Upon further increasing H, c cannot increase without bound and we expect inertial effects to intervene, leading to the saturation of c at the elastic wave speed. As just explained, probing this limit using a treadmill FEM routine is practically impossible due to the computational costs associated with large H values. Yet, we do demonstrate the inertial saturation effect by considering the dependence of c on the applied shear stress τ 0 , a dependence that has not been addressed up to now, using numerical calculations in the small-H limit. In fact, note that results of such calculations (using F τ of Eq. (15)) and efficient numerics [60,61,69]) have already been included in Fig. 6 and Fig. 7a, demonstrating perfect agreement with the small-H FEM results. According to Eqs. (19)- (20), the τ 0 dependence of c is fully contained in v r (τ 0 /σ 0 ), which is a property of the high-v velocity-strengthening branch of the steady-state friction curve (determined through τ 0 /σ 0 = f ss (v r )).
The 'spectrum' of propagation speeds c(τ 0 ) is shown in Fig. 7b. Note that c is well-defined for τ 0 > τ min , as only in this regime steady-state rupture fronts exist [60]. Interestingly, as discussed extensively in [60], in the limit τ 0 → τ min the propagation speed c attains a finite value that is much smaller than the elastic wave speed c s (cf. inset), a result that might be related to the long-debated problem of slow earthquakes/slip phenomena [27,[98][99][100]. As τ 0 is increased away from τ min , the propagation speed c increases according to the theoretical prediction of Eq. (19), which corresponds to the dashed line. However, inasmuch as c cannot increase indefinitely with H, it cannot also increase indefinitely with τ 0 . In particular, we expect that once c becomes comparable to c s , inertial effects become dominant and c → c s . This expectation is fully supported by the results of Fig. 7b. In the inertia-dominated regime other effects that are not discussed here emerge. For example, the characteristic length undergoes relativistic (Lorentz) contraction that follows a 1 − c 2 /c 2 s (in the large H limit, c s is replaced by the Rayleigh wave speed [94,97]). As mentioned above, one can consider in addition to frictional rupture fronts also frictional healing fronts. In principle, the discussion and procedures presented above in relation to rupture fronts can be applied to healing fronts as well, where the stick state invades the sliding state instead of vice versa. Here we focus on just one aspect of the results concerning healing fronts, the one that addresses the spectrum of propagation speeds c(τ 0 ). The major difference compared to the results presented in Fig. 7b, where c(τ 0 ) for rupture fronts is shown to be an increasing function of τ 0 , is that for healing fronts c(τ 0 ) is a decreasing function of τ 0 [101,102]. This implies that the two propagation speed spectra intersect at some loading level, say τ * . The latter is of physical significance because it is related to the third type of steady-state propagating frictional modes mentioned above, i.e. to slip pulses (sometimes also termed self-healing pulses) [101,103,104]. The point is that when the system is driven at τ 0 = τ * , rupture and healing fronts propagate at the same speed and hence can be superimposed without interaction to form a steady-state pulse of infinite width. As τ 0 is increased above τ * , the two fronts interact, leading to pulses of finite width L(τ 0 ), with L → ∞ for τ 0 → τ * . An example of such a finite-width propagating slip pulse in the small-H limit (i.e. using F τ of Eq. (15)) is presented in Fig. 8a [102].
Up until now, we focused on steady-state propagating frictional modes, i.e. on solutions that depend only on the comoving coordinate ξ = x−ct. These modes can be also studied in fully dynamic numerical simulations in which steadystate conditions are not enforced. A snapshot from such a simulation in the infinite-H limit (i.e. using F τ of Eq. (13)) is shown in Fig. 8b [102]. The figure presents a slip pulse that appears to propagate over long distances without appreciably changing its shape and propagation speed. Such pulses are sometimes termed sustained pulses [102], though determining whether they are truly steady-state objects or not is difficult based on dynamic simulations. The quite pronounced differences in the shape of the slip pulses presented in Figs. 8a-b are related to dimensionality (i.e. the former is obtained in the quasi-1D limit of small H and the latter in the 2D infinite-H limit).
Looking for steady-state solutions by using the co-moving coordinate ξ = x − ct has other implications. Most notably, by so doing one excludes from the discussion the possibility that some of the obtained steady-state solutions are in fact unstable. It turns out that small-H steady-state pulse solutions, an example of which is presented in Fig. 8a, are actually unstable [102]. In particular, it has been shown that small-H steady-state pulse solutions can be used to sort out perturbations based on their spatial extent; suppose that a homogeneous state corresponding to the stick fixed-point of the steady-state friction curve at a given τ 0 > τ * is introduced with a slip velocity perturbation. If this perturbation is narrower than the characteristic width of the corresponding steady-state slip pulse (i.e. its spatial extent is smaller than L(τ 0 )), then it decays in the subsequent dynamics (where steady-state is not enforced). On the other hand, perturbations characterized by a spatial extent larger than L(τ 0 ) are amplified and develop into propagating rupture fronts [102].
The unstable nature of small-H steady-state slip pulses has led to the idea that they might in fact serve as critical nuclei for rapid slip propagation mediated by rupture fronts, in analogy to first order phase transitions. That is, it has been suggested that slip pulses, and in particular their characteristic size L(τ 0 ), determine the minimal size of perturbations needed for a homogeneous stick state to develop rupture fronts that bring the system to the high-v velocity-strengthening fixed-point for the same τ 0 . This idea has also been tested in the 2D infinite-H limit, where the critical nucleus size L(τ 0 ) has been identified (obviously this L(τ 0 ) is larger than the width of the sustained pulses observed at the same τ 0 , cf. Fig. 8b). An example of a rapid rupture front emerging from introducing a perturbation larger than L(τ 0 ) into a homogeneous stick state, in the infinite-H limit, is shown in Fig. 8c. This phase transition-like nucleation process of rupture fronts at locked-in (i.e. stick state) frictional interfaces is qualitatively different from the linear frictional instability nucleation process at velocity-weakening interfaces, discussed in Sect. IV B (also compare to the nucleation scenario proposed in [105]).  Fig. 2 in [102]). Shown is the slip velocity profile as a function of the co-moving coordinate ξ used in steady-state calculations (normalized by L c of Eq. (17)). (b) A sustained slip pulse obtained through a spectral boundary integral method calculation in the 2D infinite-H limit using F τ of Eq. (13) (extracted from supplementary movie S4 in [102]). The curve is truncated in the vertical axis for visual clarify. (c) A non-steady 2D infinite-H rapid rupture front generated through a spectral boundary integral method calculation (adapted from Fig. 3 in [106]). The curve is truncated in the vertical axis for visual clarify. For all panels, see text for additional details.
The 2D infinite-H rapid rupture front shown in Fig. 8c highlights another interesting bulk elasticity effect in frictional systems, which concludes our discussion. Comparing the rupture front in Fig. 6a to the one in Fig. 8c reveals a qualitative difference; the former leaves behind it a residual stress τ r that equals the applied shear stress τ 0 (achieved ahead of the front), while in the latter we observe τ r < τ 0 . Put differently, the former features a vanishing stress drop ∆τ 0−r ≡ τ 0 − τ r → 0 (∆τ 0−r should be distinguished from the dynamic stress drop ∆τ p−r discussed above), while the latter features a finite stress drop ∆τ 0−r . The existence of stress drops has important implications for frictional rupture [3,[106][107][108]. The origin of this qualitative difference is that F τ of Eq. (13), whose value behind the rupture fronts is simply ∆τ 0−r (cf. Eq. (3)), is finite for transient, out of steady-state rapid rupture (due to both the radiation-damping term and the long-ranged interaction term s(x, t) in Eq. (13)), while it vanishes under steady-state conditions in finite-H systems [106]. These findings demonstrate that even a quantity such as the residual stress τ r , which is traditionally regarded as a property of interfacial friction law [109,110], is in fact strongly affected by bulk elasticity and bulk dynamics.

V. CONCLUDING REMARKS
Frictional systems pose great conceptual, physical and mathematical challenges, with far-reaching scientific and technological implications. A frictional system is generically composed of a low-dimensional part -the contact interface -and a higher-dimensional part -the bulks that surround the interface. The major theme of this feature paper is that in order to understand the behavior of frictional systems one has to take into account both the interfacial physics that characterizes the frictional contact and the bulk-mediated interactions between different parts of the interface, and that this inherent interface-bulk coupling is indispensable.
In the context of macroscopic interfacial constitutive laws, we discussed the widely-used rate-and-state friction modeling framework. We highlighted the power of this coarse-grained approach and also some of its limitations, and then proposed physically-motivated modifications to overcome these limitations. This part of the paper emphasized the generic properties of interfacial constitutive relations, where the discussion focused on spatially-extended, multicontact dry frictional interfaces. We did not explicitly take into account third-body effects, wear or various lubrication processes that may give rise to additional important and interesting interfacial physics. It is worthwhile mentioning that the steady-state friction curve of Fig. 3 features all of the salient properties of the Stribeck curve [111], which generically describes fluid-lubricated contact interfaces. Consequently, we expect some aspects of our results to be relevant also to lubricated interfaces, offering some unifying concepts.
In the context of bulk-mediated interactions, we focused on a linear elastodynamic bulk constitutive relation, highlighting how macroscopic geometry and material contrast affect the range and nature of the emerging spatiotemporal interactions. While we believe that linear elastodynamics is generally relevant and applicable to frictional systems, other bulk physics might be relevant and important in various situations. The latter may include bulk viscoelasticity (e.g. [112][113][114]), poroelasticity (i.e. the interaction between fluid flow and solid deformation in porous materials, e.g. [74,115]), plasticity and distributed/off-fault damage (e.g. [57]), which can affect the overall behavior of some frictional systems.
Finally, we have shown how basic phenomena in frictional systems -including instabilities of homogeneous sliding, the onset of sliding motion and a wide variety of propagating frictional modes (e.g. rupture fronts, healing fronts and slip pulses) -emerge from the inherent interplay of interfacial and bulk physics. We believe that these phenomena are generic and essential for understanding the behavior of a broad range of frictional systems. We therefore hope that the approach discussed and advocated in this feature paper will be useful for solving problems in fields ranging from engineering tribology to geophysics. This model contains the following parameters (in the order they appear above): , D, f 0 ≡ Dµ0 hσH and whatever additional parameters that appear in the threshold function g(·). Note that α defined here is the same one used in Eq. (7). The values of these parameters are discussed below in Appendix B and Appendix C.
In most of the studies discussed in this work (exceptions are mentioned in Sect. II B 1 and in Appendix B), g(v) = 1 + (v * /v) 2 is used. This choice, which introduces a single additional parameter v * , ensures that |v|g(v) → v * in Eq. (A2) is extremely small for |v| v * , once an extremely small v * is chosen, and that |v|g(v) → |v| for v v * . For this choice of g(·), the emerging predictions are insensitive to the precise value of v * , as long as it is significantly smaller than all other velocity scales in a given problem. In some applications of the model, most notably when large slip displacements are involved, one assumes that τ el quickly relaxes to its steady-state. Under these conditions, we setτ el = 0 inside Eq. (A2) and substitute the resulting τ el in Eqs. (A1)-(A2) to obtain The steady-state friction curve f ss (v), which is sketched in Fig. 3, is obtained once we setφ = 0 and substitute the solution for φ in f (φ, v) of Eq. (A3). Finally, note that the resulting f ss (v) reduces to the conventional RSF expression in Eq. (7) once the two constitutive extensions introduced in this paper are removed, i.e. log(1 + φ/φ * ) → log(φ/φ * ) and µ 0 → 0 (the latter impliesf 0 = 0) are substituted in Eq. (A3), and then the steady-state of φ is used. To see this, set g(·) = 1 ( g(·) accounts for the low-v velocity-strengthening branch in Fig. 3, which is not described by Eq. (7)), expand sinh −1 (x) ≈ log(2x) for large argument x and drop the resulting O(log 2 ) term. The outcome identifies with Eq. (7), with the parameters f 0 ≡ E0 Ω σH and β ≡ bf 0 , where α is the same as above (and note thatv is not used here, i.e.v ≡ v c exp[− E0 k B T ] has been substituted in Eq. (A3) before sinh −1 (·) has been expanded).

Appendix B: Theoretical predictions for frictional dynamics under small stresses
The set of equations used to obtain the results reported on in Fig. 2b is listed and discussed in detail in Sect. II B 1. The only additional input needed in order to reproduce Fig 2b is the function g(·) and the determination of the material parameters. The experimental data in Fig. 2a clearly indicate that the frictional response up to shear forces of 5−6N is predominantly linear elastic, i.e. that τ vis is negligible in Eq. (9) (and consequently also in Eq. (A1)) in this range of loading forces. While the onset of irreversibility is not sharp, it is clear that significant deviations from a linear elastic behavior emerge only for shear forces larger than 7N. To capture these observations in a phenomenological manner, we set g(τ ) = Θ(τ − A(φ)τ c ), where Θ is the Heaviside step function. Here τ c is a characteristic interfacial yield stress, whose value is chosen such that the onset of irreversibility takes place around 7N (in particular, given A(φ(t = 0)) of Eq. (12) as detailed below, we set τ c = 70MPa). With this choice of g(·), which has also been used in [60], the v ≈ 0 regime in the steady-state friction curve of Fig. 3 becomes a strictly vertical (elastic) v = 0 line (as shown in Fig. 1a of [60]).
The friction parameters for the PMMA plates used in the experiments of [50], and which enter Eqs. (A1)-(A2), are the same as those reported on in Table I below (note that v * there is irrelevant for the present analysis as g(·) here is characterized by τ c , as just discussed). The external parameters are as reported explicitly in [50] (specifically we have A n = 0.01m 2 , F N = 24N and |dF S /dt|= 0.5N/s) and the initial conditions are φ(t = 0) = 100s (which together with b of Table I allows to determine A(φ(t = 0)) mentioned above),δ(t = 0) = 0 and τ el (t = 0) = 0. The only missing parameter is the interfacial elasticity ratio µ 0 /h, which is directly extracted from the initial linear slope in the F S (δ) experimental data (dashed line in Fig. 2a) according to dF S /dδ = µ0F N hσH [1 + b log(1 + φ(t = 0)/φ * )], resulting inf 0 = Dµ0 hσH = 0.209 (which is not the same as the one listed in Table I, though the two are quite close to each other). Then Eqs. (11), (A1) and (A2) are solved for φ(t),δ(t) and τ el (t) (for the two applied loading protocols F S (t)). Finally, δ(t) of Fig. 2b is obtained (for the two protocols) by integratingδ(t) with δ(t = 0) = 0.
Appendix C: A treadmill FEM routine for obtaining steady-state frictional modes The FEM results reported above have been obtained using the FEM software package FreeFem++ [116]. The partial differential equations are solved in a monolithic scheme on a triangular mesh using a real space implementation. The fully inertial evolution equations, which include the linear elastodynamic equations (momentum balance with Hooke's law), are expressed in weak form in the spirit of the finite element method. The time derivatives are expressed via finite differences. The frictional boundary condition is implemented via a semi-implicit integration scheme. Details about the results reported on in Fig. 5 can be found in [87], while here we focus on the results reported on in Fig. 6 and Fig. 7a.
In these calculations, a deformable body of height H and length L H sliding on top of a rigid substrate has been considered. A shear stress of magnitude τ 0 , which is larger than the minimum of the steady-state friction curve in Fig. 3, and a normal stress σ yy = −σ 0 are applied at y = H (cf. Fig. 1). The lateral boundaries are traction free and the rigid substrate implies that we have u y = 0 at the interface (which ensures no interpenetration into and detachment from the substrate). In order to efficiently reach steady-state conditions, a "treadmill" routine in which the system is shifted such that the propagating mode always remains in the lateral center of the sliding body has been employed. Consequently, we effectively obtain a description in a co-moving reference frame. This requires imposing additional boundary conditions on the fields φ and τ el (cf. Eq. (A2)) at the leading edge, which are chosen according to the stable low-velocity fixed-point corresponding to the intersection of τ 0 the steady-state friction curve in Fig. 3. The treadmill scheme also allows to locally increase the spatial resolution near the front position, a main region of interest, without the need for remeshing during propagation, as the front always remains at the same position.
For the largest simulated systems with H = 16cm we used a slider of length up to L = 9.6m. Altogether, we typically used non-uniform discretizations with up to 800×20 grid points, and checked that the results are insensitive against further mesh refinement. To avoid numerical instabilities, the dynamics is suppressed near the trailing and leading edge of the slider, and its length L is chosen to be sufficiently large to avoid their influence on the front dynamics. We used an initially large timestep, ∆t = 10 −4 s, which allows for rapid progress in the initial stage of the simulations. This large initial timestep also stabilizes the dynamics, but leads to reduced accuracy. Consequently, we varied the timestep such that it relaxes to a final value of ∆t = 10 −7 s in the course of the simulations, and we checked that the results are then insensitive against further reduction of the timestep. The simulations typically run until the front has traversed the system length at least once to get rid of sound waves and perturbations which result from the initialization of the system, and which interfere with the front dynamics. After this time, steady-state conditions are reached, and the spatial profiles of all relevant fields are extracted. The results in Fig. 6 and Fig. 7a) involve the friction law of Eqs. (A1)-(A2) (together with g(v) = 1 + (v * /v) 2 ), with the bulk and interfacial parameters listed in Table I. The same interfacial parameters are relevant also for Fig. 2b, see the caption of Table I Fig. 6 and Fig. 7a. Note that ν is Poisson's ratio (not defined in the manuscript), which is needed for Hooke's law, and that an unrealistically small mass density ρ is used (leading to unrealistically large wave speeds), as strongly inertial effects have not been of interest in these particular calculations. The applied normal stress in these calculations has been set to σ 0 = 1MPa. Finally, note that σ H is not needed for solving the equations, just for calculating A of Eq. (12), the result of which is presented in Fig. 6c. The very same interfacial parameter, except for v * andf 0 , have been used to generate the results in Fig. 2b (see Appendix B for details).