Frictional sliding without geometrical reflection symmetry

The dynamics of frictional interfaces play an important role in many physical systems spanning a broad range of scales. It is well-known that frictional interfaces separating two dissimilar materials couple interfacial slip and normal stress variations, a coupling that has major implications on their stability, failure mechanism and rupture directionality. In contrast, interfaces separating identical materials are traditionally assumed not to feature such a coupling due to symmetry considerations. We show, combining theory and experiments, that interfaces which separate bodies made of macroscopically identical materials, but lack geometrical reflection symmetry, generically feature such a coupling. We discuss two applications of this novel feature. First, we show that it accounts for a distinct, and previously unexplained, experimentally observed weakening effect in frictional cracks. Second, we demonstrate that it can destabilize frictional sliding which is otherwise stable. The emerging framework is expected to find applications in a broad range of systems.


I. INTRODUCTION
Understanding frictional sliding is a long-standing challenge with important practical and theoretical implications. It is relevant in diverse physical systems spanning a broad range of scales, from the nano-scale to the planetary-scale. A complete analytic treatment of sliding frictional interfaces is generally a formidable task. Two major factors are responsible for the complexity of the problem. First, the friction law, i.e. the constitutive relation that describes the shear traction at the frictional interface, poses experimental challenges and depends on the slip rate and slip history in a highly nonlinear fashion [1][2][3][4][5][6][7][8][9][10]. The second factor is the elastodynamics of the sliding bodies, i.e. the time-dependent long-range stress transfer mechanisms between different points along the interface. It is particularly challenging when the two bodies are made of different materials and in the generic case in which spontaneously-generated interfacial rupture fronts dynamically propagate along the interface [11][12][13][14][15][16][17][18][19][20][21][22][23].
A significant simplification in relation to the second factor is obtained when the system possesses reflection symmetry across the interface, i.e. when the two materials are identical, the geometry is symmetric, and the loading configuration is antisymmetric (here and elsewhere we consider macroscopic geometry. Differences in small-scale roughness typically exist and are effectively incorporated into the interfacial constitutive relation). A prototypical example of such a situation is that of two semi-infinite half-spaces made of identical elastic materials, a situation that was extensively studied in the literature (see, for example, [24][25][26][27][28][29][30][31][32][33][34][35][36]). The main simplification comes from the fact that such a symmetry precludes a coupling between tangential slip and variations in the normal stress. The lack of such symmetry has important implications on the stability of slid- * M. Aldam and Y. Bar-Sinai contributed equally to this work. ing [11,15,16,20,21,23,37,38], the failure mechanism and rupture directionality [17,21,23,[37][38][39][40][41][42][43][44][45]. Physically, this happens because sliding can enhance (reduce) the normal stress, which in turn can inhibit (facilitate) frictional sliding. The origin of the absence of reflection symmetry is traditionally assumed to be constitutive in nature, i.e. sliding of dissimilar materials is usually considered. This is known as the bi-material effect. Sliding along such bi-material interfaces has been quite extensively studied in the literature and this material contrast is thought to have important implications for frictional dynamics [11, 15-17, 20-23, 37-45]. The purpose of this paper is to explore the possibility of asymmetry of a geometric origin, i.e. sliding of two bodies made of the same mate-rial without geometrical reflection symmetry. Examples of such geometries are depicted in Fig. 1: sliding of two blocks with different thickness in the direction orthogonal to sliding, an experimental setup that was used in various recent works [46][47][48][49][50] and will be theoretically addressed below (panel a); sliding of a block of finite height H over a semi-infinite bulk, a simple example to be analyzed in depth in this work (panel b); finally, an idealized sketch of tectonic subduction motion is shown (panel c), a situation in which one lithospheric plate is subducted beneath another one and is responsible for most of the large magnitude earthquakes ("megathrust") occurring on the Earth's crust [51-56, for example]. Obviously, many other sliding geometries which lack reflection symmetry can be conceived. Generally speaking, this situation is expected to be the rule rather than the exception, since no physical system features perfect reflection symmetry.
In this paper we lay out a rather general theoretical framework to address frictional sliding in the absence of geometrical reflection symmetry and support it by extensive experiments. A major outcome is that the effect of geometric asymmetry resembles, sometimes qualitatively and sometimes semi-quantitatively, that of material asymmetry. Consequently, many results obtained for bi-material interfaces are also relevant to interfaces separating bodies made of identical materials with different geometry. As first applications, two main results are obtained within the newly developed framework: • A novel explanation of a sizable weakening effect observed in recent experiments on rupture fronts propagation along frictional interfaces [49]. The weakening effect is directly linked to geometric asymmetry and is shown experimentally to disappear in its absence. This result has important implications for the failure dynamics of frictional interfaces.
• We demonstrate that geometric asymmetry can destabilize frictional sliding which is otherwise stable. We consequently expect geometric asymmetry to play an important role in frictional instabilities.
The emerging framework should find additional applications in a broad range of frictional systems.

II. GENERAL FRAMEWORK
Consider two blocks in frictional contact. At this stage, the discussion remains completely general, allowing the two blocks to be made of different materials, to feature different geometries and to experience general external loadings. We denote the displacement vector fields in the two blocks as u (1) (x, t) and u (2) (x, t), where the superscripts correspond to the upper and lower blocks, respectively. Each of these satisfies the momentum balance equation ∇·σ = ρü, where ρ is the mass density of each block. Cauchy's stress tensor σ is related to the displacement gradient tensor ∇u according to the isotropic Hooke's law (1 + ν)µ [∇u+(∇u) T ] = σ − ν(I tr σ − σ).
Here I is the identity tensor, ν is Poisson's ratio and µ is the shear modulus of each block. The coordinates are chosen such that the interface lies along the x-axis, which is also the direction of sliding, see Fig. 1. The direction normal to the interface is the y-axis and the interface is the surface y = 0. The z-axis is in the thickness direction, where z = 0 is the center line. While the formulation below and the analysis in Sect. IV are twodimensional (2D), we shall see that three-dimensional (3D) effects involving the z-coordinate play an important role in Sect. III.
Since the bulk equations are linear, one can analyze separately each interfacial Fourier mode, i.e. write u (n) (x, y = 0, t) = u (n) e ik(x−ct) , where k > 0, c is the complex phase/propagation velocity and n = 1, 2. The relation between the interfacial displacements and stresses is also linear and can be written as u where the matrix M (n) can be obtained from the Green's function of the corresponding medium and σ (n) yi are the interfacial stresses, i.e. at y = 0. For example, under quasi-static conditions for semi-infinite blocks in 2D, this relation for the lower block (i.e. the block at y < 0) takes the form [57] The essence of frictional motion is that the displacement field is discontinuous across the interface. We denote the slip discontinuity at the frictional interface by i (x) ≡ u (1) i (x, y = 0 + ) − u (2) i (x, y = 0 − ) .
On the interface, y = 0, no separation or inter-penetration between the bulks implies y = 0 and continuity of σ yi . Together with the known dynamic response matrices M (n) , these requirements can be used to calculate the relation between the slip discontinuity and the interfacial stresses of the composite system that consists of both bulks. Following [58], this is done by noting that for y = 0 we have σ yi = (M (1) −M (2) ) −1 ij j . Thus, the response of the composite system in 2D reads (in Fourier space) where we defined the elastic response functions yx . That is, the G i 's can be expressed as functions of the response coefficients of both bulks. Note that the imaginary unit i is included in Eq. (3) for convenience. Note also that we use the same notation for a function and its Fourier transform, as they are easily distinguishable by the context or the stated arguments (e.g. k or x).
The central player in the analysis to follow is G y , which represents the elastodynamic coupling between tangential slip and normal traction along the interface. In systems with complete reflection symmetry along y = 0, this coupling is precluded by symmetry. To see this, note that in this case the off-diagonal elements of M (1) and M (2) are identical [57], and thus M (1) − M (2) , as well as its inverse, is diagonal. This immediately implies G y = 0. In what follows, we study two important frictional problems in which the sliding bodies are made of identical materials, i.e. ρ (1) = ρ (2) ≡ ρ, µ (1) = µ (2) ≡ µ and ν (1) = ν (2) ≡ ν, yet reflection symmetry relative to the interface is absent due to asymmetry in the geometry of the bodies, leading to G y = 0. These problems highlight the importance of geometrical asymmetry to frictional sliding and its relation to the conventional bi-material effect.

III. "THIN-ON-THICK" SYSTEMS AND THE PROPAGATION OF FRICTIONAL CRACKS
The first problem that we examine, which is directly motivated by recent experimental observations [44,49], is depicted in Fig. 1a. In this system, a thin block of width W = 5.5 mm is pushed along its length (the x-axis) on top of a significantly thicker block (here 30 mm). This "thin-on-thick" experimental setup was used in various studies [46][47][48][49][50], where a transparent glassy polymer (poly(methyl-methacrylate), PMMA) was used. The transparent material allows a direct real-time visualization and quantification of a fundamental interfacial quantity: the real contact area, A r . The latter is the sum over isolated micro-contacts formed due to the small scale roughness of macroscopic surfaces.
A r is typically orders of magnitude smaller than the nominal contact area, A n . Their ratio, A ≡ A r /A n 1, plays a critical role in interfacial dynamics [2,3,[46][47][48][49][59][60][61][62] because the frictional resistance/stress is proportional to A, σ xy ∝ A, i.e. the larger the real contact area the larger the frictional resistance. A itself depends on the normal stress and also on the slip history of the interface according to where ψ is an internal variable characterizing the state of the interface. Frictional sliding leads to reduction of ψ, i.e. to a reduction of the contact area [1-3, 6, 61, 63]. In the absence of sliding, ψ (and hence A) grows logarithmically with time, a process known as frictional aging [1,3,48,64,65]. In [49] it was found that sliding is mediated by a succession of crack-like rupture fronts propagating along the frictional interface and that these fronts are surprisingly well described by the classical theory of shear cracks propagating along an interface separating identical materials, Linear Elastic Fracture Mechanics (LEFM). The variation of A along a few of these fronts is shown here in Fig. 2a. It is seen that the rupture fronts involve a significant overall reduction of the contact area, which weakens the interface (i.e. reduces ψ) and facilitates sliding. We would like to focus our attention on a distinct feature of these curves: As observed in Fig. 2a, fronts which travel at 90% of the Rayleigh wave-speed c R , here c R 1237 m/s (for plane-stress conditions [57]), or slower (not shown), feature a monotonic decrease of A. However, in fronts propagating even closer to c R , A features a non-monotonic behavior, i.e. A undershoots the asymptotic value A ∞ (i.e. A as x → −∞) and then rapidly increases, at a rate way too high to be explained by slow frictional aging. This non-monotonic behavior remained unexplained in [49], where it was stated that "the nonmonotonic behavior of A ... suggests interesting dynamics as c → c R ...".
In Fig. 2b we show the spatial profiles of the slip velocity v, corresponding to the contact area profiles shown in Fig. 2a. These profiles were calculated from the experimental data using the simplest cohesive zone model [66][67][68] which is consistent with the measurements of the fracture energy and cohesive zone size (see [57] for more details). This model, while generally used to describe identical materials, is motivated by the empirical observation [57] that the strain fields are, to first order, quite similar in the thin-on-thin and thin-on-thick setups. This approximation would, of course, have to be modified in cases of strong material contrast, where the fields on both sides of the interface differ strongly [44].
We denote the maximum slip velocity in these profiles by v m . Next, in order to quantify the non-monotonic effect, we define the magnitude of the undershoot ∆A as the difference between the asymptotic value A ∞ and the minimum of the profile over the range −5 mm < x < 0, which is the typical spatial range for which ∆A > 0 is observed in the thin-on-thick setup, see Fig. 2a. ∆A/A ∞ is plotted vs. v m in Fig. 2c (red symbols), demonstrating that the former is quasi-linear (i.e. predominantly linear) in the latter. Note that the spread in the data does not allow to identify any systematic deviations from linearity. We stress that the effect is not only qualitatively novel, i.e. the existence of a non-monotonic contact area behavior ∆A/A ∞ > 0, but it is also quantitatively important. As Fig. 2c shows, the local reduction in the real contact area ∆A/A ∞ can reach nearly 25%. This is a large quantitative effect, compared to other documented frictional effects, implying the existence of significant local frictional weakening which can significantly influence interfacial dynamics.
What is the source of this non-monotonicity, why does it scale quasi-linearly with the slip velocity and why does it appear only at sonic propagation velocities? Such behavior has recently been observed in [44] when investigating the frictional motion of bi-material interfaces in a geometrically symmetric system. There, a very large local reduction of A was observed at sonic propagation velocities, but it entirely disappeared when the upper and lower blocks were made of the same material. We propose that the same happens in our case, only here it is due to geometric asymmetry. That is, we suggest that the non-monotonicity of A stems from the absence of geometrical reflection symmetry of the two blocks, i.e. from the difference in their thickness. If true, then the fast  (a) Snapshots of the spatial profile of the contact area A of rupture fronts in the "thin-onthick" setup (see Fig. 1a and [49] for additional details). These fronts propagate to the right at velocities c indicated in the legend of panel b (0.900cR < c < 0.993cR, see [57]), where x = 0 corresponds to the tip of each rupture front. The contact area is normalized by its value A0 before the passage of the front. (b) The slip velocity profiles corresponding to the snapshots in panel a (see [57] for details). (c) ∆A/A∞, where ∆A is the magnitude of the real contact area undershoot and A∞ is the asymptotic value (see inset), vs. the maximal slip velocity vm (see panel b) for both the "thin-on-thick" setup (red symbols) and the geometrically symmetric "thin-on-thin" setup (blue symbols). Different symbols correspond to different experiments and their size roughly corresponds to the measurement error. The red line is the best linear fit for the red symbols. (inset) The contact area profile for c = 0.993cR in the "thin-on-thick" setup (red line, already appearing in panel a) and in the geometrically symmetric "thin-on-thin" setup (blue line).
non-monotonic variation of A is not an intrinsically frictional phenomenon, i.e. a result of the dynamics of the state of the interface ψ, but rather an elastodynamic effect emerging from the coupling between slip and normal stress variations, solely induced by geometrical effects. In terms of Eq. (4), we propose that ψ is monotonic and that the non-monotonicity of A results from a non-monotonic behavior of σ yy .
Our strategy in testing and exploring this idea is twofold. First, our idea can be directly tested by a definitive experiment. That is, we expect that when the width of the lower block equals that of the upper one, i.e. in a "thin-on-thin" setup, the non-monotonicity in A disappears altogether even in the limit c → c R . We performed this experiment, as in [44], and present a representative example (for c = 0.993c R ) in the inset of Fig. 2c (blue line). The curve is indeed monotonic. Moreover, note that the asymptotic value A ∞ is the same as that in the "thin-on-thick" setup (cf. Fig. 1a, for the same propagation velocity), even though the latter exhibits a large undershoot. In Fig. 2c (main panel) we added ∆A/A ∞ of many rupture fronts in the "thin-on-thin" setup (blue symbols). ∆A/A ∞ is indeed very close to zero (small negative values simply correspond to monotonic behavior), i.e. all of the A profiles in the "thin-on-thin" setup are monotonic. This direct experimental evidence provides unquestionable support of our basic idea that the non-monotonic behavior corresponding to the "thin-onthick" setup data in Fig. 2a emerge from the absence of geometrical reflection symmetry.
Next, our aim is to develop a theoretical understanding of the origin of non-monotonicity. The challenge is to explain both the fact that it emerges at asymptotic propagation velocities (c → c R ) and the quasi-linear relation between ∆A/A ∞ and v m . The complete problem, involving a thin block sliding atop a thicker one, is a very complicated 3D elastodynamic problem. We approach the problem by breaking it into two steps. First, we perform a simplified analysis, invoking physically-motivated approximations, which allow us to reduce the mathematical complexity of the problem and gain analytic insight into it. The major simplification is to consider the corresponding quasi-static problem instead of the full elastodynamic one. The physical rationale for this is clear: the absence of geometrical reflection symmetry should manifest its generic implications also in the framework of static elasticity and hence the simplified analysis is expected to reveal the origin of the non-monotonicity of the real contact area. Then, in the second step, we use the static results in an effective dynamic calculation, to be explained below.
The main outcome of the first step is that the static 3D problem can be approximately mapped onto a 2D problem involving two elastically dissimilar materials. That is, we show that the geometric asymmetry can be approximately mapped onto an effective constitutive asymmetry, i.e. an effective material contrast. To see how this emerges, we assume that both blocks are infinite in the ydirection and that the thicker (lower) block is also infinite in the z direction. That is, the lower block is replaced by a semi-infinite 3D half-space, which allows us to use the well-known interfacial Green's function [69]. More specifically, the 3D real-space Green's function matrix M 3D (r − r ) [69] allows us to express the interfacial displacements at a point r = (x, y = 0, z = 0) on the symmetry line, u(r) = (u x , u y ), induced by a point force applied by the upper block at r = (x , y = 0, z ), F (r ) = (F x , F y ). Note that the latter is assumed not to contain an outof-plane component, i.e. F z = 0, which in principle could emerge from frustrated Poisson expansion at the interface. It is reasonable, though, to neglect it to leading order.
We physically expect shear tractions to be uniform across the thickness W , hence they are taken to be constant for |z| ≤ W 2 (and of course to vanish for |z| > W 2 ). Thus, we obtain where the effective 2D response matrix M eff (k) of the thicker (lower) block is given by the Fourier transform of M 3D over the strip |z| ≤ W 2 , The integration can be carried out, resulting in where q ≡ kW and B(q) = π −1 q 0 K 0 (q /2)dq (K 0 (z) is the modified Bessel function of the second kind of order 0). The outcome of the analysis, which is presented in full detail in [57], is that M eff (k) appears to identify with the 2D response matrix of Eq. (1), if one defines the effective elastic moduli of the lower (thicker) block as These are plotted in Fig. 3a.
The mapping of the 3D problem onto an effective 2D problem is formally valid as long as the interfacial stresses (and hence displacements) in Eq. (5) are approximately localized in Fourier k-space. Otherwise, Eq. (5) will not identify with Eq. (1) due to the extra k-dependence of µ eff (kW ) and ν eff (kW ), which is a result of the 3D nature of the original problem. We note in passing that in the limit q = kW 1, µ eff → µ and ν eff → ν, which corresponds to 2D plane-strain conditions [71]. This is expected for small wavelengths, for which the thinner block also appears infinitely thick, and hence is a consistency check on our calculation. The important observation, though, as is clearly seen in Fig. 3a, is that for the thicker block µ eff (k) > µ for all experimentally relevant k's [72]. This suggests that the thicker block is effectively stiffer than the thinner one, as hypothesized in [73,74] where the thicker block was assumed to correspond to plane-strain conditions in numerical simulations. That is, the main physical insight gained from the performed analysis is that geometric asymmetry gives rise to an effective material contrast.
With this physical insight in hand, we aim now at addressing the non-monotonicity of A discussed in panels a and c of Fig. 2. The 3D static analysis presented above may not yield quantitatively accurate predictions when strongly elastodynamic 2D interfacial rupture fronts are considered. Yet, we believe that the insight embodied in the relations µ eff (k) > µ and ν eff (k) > ν is physically robust and hence try to explore their quantitative implications in relation to the experimental observations in the dynamic regime.
To accomplish this, we consider the 2D dynamic transfer function G y (c, k) in Eq. (3) and take it to approximately describe the experimental system when the effective moduli µ eff (k) > µ and ν eff (k) > ν are used for the thicker (lower) block and plane-stress conditions [71] are assumed for the thinner (upper) block. Note that it is justified to treat the heights of the two blocks as infinite since the experimental rupture fronts are so fast that they do not interact with the upper and lower boundaries before traversing the whole system. Therefore, Eq. (3) can be rewritten as where we used v =˙ x = −i c k x for a constant propagation velocity c.  7). (inset) The variation of the effective Poisson's ratio ν eff (q) with q = kW . In both we used ν = 0.33, which is relevant for PMMA [49,70]. (b) The response function Gy, quantifying the effective bi-material contrast according to µ eff (q) and ν eff (q) (for the thicker block, the thinner one is represented by plane-stress conditions), corresponding to selected values of q = kW . The corresponding values of the elastic moduli µ eff (q) > µ and ν eff (q) > ν are marked in panel a and its inset using the same color code. (c) ∆σyy given in Eq. The 2D infinite-system dynamic transfer function G y (·) in Eq. (8) was calculated by Weertman for sliding of dissimilar materials quite some time ago [12]. We reiterate that the basic idea here is to use a known result for dissimilar materials to represent a system composed of identical materials with geometric asymmetry, utilizing the effective moduli derived in Eq. (7), µ eff (k) > µ and ν eff (k) > ν. In the presence of any contrast between the shear moduli of the materials, G y (c) is finite and increases significantly at elastodynamic velocities (in fact, it diverges when c approaches the shear wave-speed c s of the more compliant material), as shown in Fig. 3b. Thus, we expect rupture fronts that propagate at near-sonic velocities to be accompanied by a significant reduction in the local normal stress as implied by Eq. (8), reducing locally the real contact area. In turn, this reduces the interfacial strength, which facilitates sliding. This is consistent with the experimental observations of Fig. 2a, where the non-monotonicity of A becomes substantial at asymptotic propagation velocities (c → c R ). This normal stress reduction is also remarkably similar to the recent observations of [44] in bi-material systems, a similarity that further strengthens the analogy between geometric asymmetry and material asymmetry.
The connection between geometric and material asymmetries is yet further strengthened when the directionality of rupture is considered. The sub-Rayleigh (c < c R ) rupture fronts, shown in Fig. 2b, propagate from left to right, in the direction of sliding of the thinner (upper) block (see also Fig. 1a). Sub-Rayleigh rupture fronts that are accompanied by normal stress reduction are known to propagate in the direction of sliding of the more compliant material in a bi-material setup, the so-called "preferred direction" [12,17,44]. This is fully consistent with our result that the thinner (upper) block is effectively softer than the thicker (lower) block (or alternatively, that the thicker block is effectively stiffer than the thinner one).
The quasi-linearity of ∆A with the (maximal) slip velocity, observed in Fig. 2c, naturally emerges from Eq. (8). To see this, note that ∆A ∝ ∆σ yy according to Eq. (4) (recall that ψ in that equation is expected to be monotonic) and that c remains close to c R to within a few percent. In this regime (c c R ), G y does not change appreciably as a function of c, while the maximal v varies quite substantially (cf. Fig. 2a). That means that while c c R is required for the existence of the weakening effect, its variability is mainly determined by v. Put together, we obtain ∆A ∝ v. To obtain some estimate of the proportionality factor between ∆A and v along this line of reasoning, we interpret ∆σ yy in Eq. (8) to be a function Fig. 3a).
The results for ∆σ yy (v) with kW = 3, 4, 5, normalized by the experimentally applied normal stress σ 0 , are shown in Fig. 3c. The slope of the kW = 5 line is very close to the slope of the linear fit in Fig. 2c, which was added to Fig. 3c for comparison (gray dashed line). Note that the experimental line features a finite v intercept, which is absent in the theoretical one. This is expected since the undershoot, ∆A, is generally susceptible to variations both of σ yy and the fracture of contacts (variations of ψ in Eq. (4)). For low values of v, variations of σ yy should be small, and the spatial profile of A is therefore dominated by variations of ψ. A(x) should therefore be monotonic in space, similar to the spatial profile in the "thin-on-thin" setup, thus rendering any undershoots (i.e. ∆A) to be unmeasurable. This quantitative agreement should be taken with some caution in light of the various approximations invoked above. Yet, the existence of a characteristic wavenumber kW = 5 is not unreasonable as the typical scale of the velocity peaks (see Fig. 2b), the spatial scale of the undershoot in the contact area (see Fig. 2a) and W are all in the mm-scale. Furthermore, the relative magnitudes of the slopes in Fig. 3c provide a testable prediction for how the slope decreases with increasing W . This should be experimentally tested in the future. Finally, as W increases and approaches the width of the lower block, the non-monotonicity is predicted to disappear, as demonstrated experimentally in Fig. 2c (blue symbols).
The results presented in this section demonstrate that global geometric features of the sliding bodies in a frictional problem, here a difference in their thickness, affect the frictional resistance to sliding and in fact makes it easier for interfacial rupture fronts that mediate sliding to propagate. In fact, the effect of geometric asymmetry is maximal at the extreme rupture velocities that are the norm in frictional sliding. This reduction in frictional dissipation applies to any engineering or tribological system involving identical materials and geometric asymmetry. As such, it implies that the design and friction control of any real-life tribological application must take into account not only the interfacial properties, but also the relative size of the sliding bodies. In the next section we show that the same concept applies to another class of important sliding friction problems, where a different form of geometric asymmetry controls the dynamic response of the system.

IV. STABILITY OF FRICTIONAL SLIDING
We now focus on a different, yet conceptually related, physical situation in which geometric asymmetry plays a crucial role as well. While in Sect. III geometric asymmetry was associated with a difference in the thickness of the sliding bodies, here its origin is a difference in their height. Moreover, while in Sect. III we addressed the propagation of spatially-localized interfacial cracks, here the focus will be on the stability of homogeneous sliding. Yet, in both cases a geometry-induced coupling between interfacial slip and normal stress variations, encapsulated in the function G y in Eq. (3), is the dominant physical player.
We consider an elastic block of a height H (1) sliding atop a block of height H (2) = ηH (1) (with a dimensionless positive η, 0 < η < ∞), both made of the same material under plane-strain conditions [71], as depicted in Fig. 1b. Note that η = 1 corresponds to a symmetric system. The blocks initially slide at a fixed velocity and all of the fields are assumed to reach steady state. A homogeneous compressive normal stress σ (1) yy = −σ 0 is imposed at both y = H (1) and y = −H (2) . In addition, a constant velocitẏ u (1) x = v in the positive x-direction is imposed at y = H (1) andu (2) x = 0 at y = −H (2) . In this problem, unlike the problem considered in Sect. III, the interfacial dynamics are coupled to the boundaries at y = H (1) , −H (2) , and hence the heights H 1,2 are expected to play a central role here.
To fully define the problem, one needs to specify the frictional boundary condition at the interface. Friction is commonly modeled as a linear relation between the interfacial normal stress and the interfacial shear stress (frictional resistance/stress), i.e.
where f (·) represents the friction law. Our major goal here is to understand the destabilizing effect associated with geometric asymmetry, i.e. η = 1, which to the best of knowledge has not been studied before. Consequently, in order to isolate the geometric effect, we will focus below on situations in which friction is intrinsically stabilizing such that any instability, if exists, is associated with the absence of geometrical reflection symmetry.
To achieve this, we proceed in two steps. First, in Sec. IV A, we present a simplified analysis involving a simple velocity-dependent friction law and strong geometrical asymmetry. This will allow us to gain much insight into the role of geometric asymmetry in frictional sliding and to clearly identify the physical origin of instability. Then, in Sec. IV B, we present a significantly generalized analysis for a realistic friction law, including an internal state variable and an interfacial memory length, and for any level of geometric asymmetry. The emerging results strengthen the findings of Sec. IV A and extend them.
A. Simplified analysis: Velocity-dependent friction and large geometric asymmetry As a primer, we use here a simple friction law where f (·) in Eq. (9) depends only on the interfacial slip velocity v ≡˙ x , i.e. f (v). We focus on velocity-strengthening interfaces, f (v) > 0, because in this case sliding is unconditionally stable for symmetric systems [57,75,76] and thus the origin of any emerging instability must be associated with the absence of geometrical reflection symmetry. Moreover, steady state velocity-strengthening friction has been recently shown to be a generic feature of dry interfaces over some velocity range [6]. Finally, to simplify the analysis further we consider the case in which the lower block is much higher than the upper one, η 1. That is, we take the limit Under what conditions is homogeneous sliding stable? This question, which is of fundamental importance in a broad range of frictional problems (see, for example [3,20,23,63,75,[77][78][79][80][81][82][83][84]), is first investigated in the context of the simplified problem defined above. As the interface is characterized by velocity-strengthening  friction, f (v) > 0, friction itself tends to stabilize sliding. Consequently, the only possible destabilizing piece of physics can be the geometric-asymmetry-induced coupling between interfacial slip and normal stress variations, encapsulated in the function G y (cf. Eq. (3)), which also played a crucial role in Sect. III. Can geometric asymmetry destabilize velocity-strengthening frictional interfaces in much the same way as material asymmetry (the bi-material effect) can [23]?
To address the stability question, we perturb Eq. (9) to linear order, obtaining [57] µ G x (c, k) + i µ f G y (c, k) + i c σ 0 δf /δv = 0 , (10) which is an implicit equation defining the linear stability spectrum c(k). In the simple velocity-dependent friction case considered here, we have δf /δv = f (v) (more general interfacial constitutive laws are considered in Sect. IV B). Perturbations with [c] > 0 are unstable and will grow exponentially, while perturbations with [c] < 0 are stable (remember that k > 0). An explicit calculation shows that G y reads [57] where c s and c d are respectively the shear and dilatational wave-speeds and α 2 s,d ≡ 1 − c 2 /c 2 s,d was introduced. The limit H → ∞ amounts to a symmetric system, in which case η → 1, and indeed G y vanishes in this limit. We can thus expect the system to be unconditionally stable for H → ∞. G y also vanishes in the limit H → 0. Similarly, G x takes the form [57] .
Equipped with the results for the dynamic response functions G i (c, k), the implicit equation for the spectrum, Eq. (10), can be in principle solved, at least numerically. The equation admits a few solution branches, and in general its analysis is far from trivial. However, since the purpose of the present discussion is not a complete analysis of Eq. (10), but rather a demonstration of the qualitative effect of the absence of geometrical reflection symmetry, we focus here on a particular branch of solutions which is shown in Fig. 4a. It is observed that for a range of parameters, and for a finite range of wavenumbers, the solutions are unstable ( [c] > 0). This is a direct numerical evidence that geometric asymmetry can destabilize systems which are otherwise stable (remember that f (v) > 0).
It seems natural at this point to ask under what conditions this instability is observed. What are the conditions on the various system parameters such that there will be a range of k's for which [c(k)] > 0? As a prelude, we perform a dimensional analysis. Clearly, the only lengthscale in the problem is H and indeed the wavenumber k only appears in the dimensionless combination kH. Thus, large (small) k is equivalent to large (small) H and since G y vanishes in both limits H → 0 and H → ∞, we expect to find unstable modes only in a finite range k min < k < k max , if any.
Another dimensionless combination is γ ≡ µ/(σ 0 c s f (v)), which is the ratio of the elastodynamic quantity µ/c s -proportional to the so-called radiation damping factor for sliding [23,26,75,85] and the response of the frictional stress to variations in the sliding velocity. As such, γ quantifies the importance of elastodynamics, which tends to destablize sliding when geometrical asymmetry is present, relative to velocity-strengthening friction, which generically stabilizes sliding. We thus expect large γ to promote instability, if G y = 0. In addition, as G y is the only possible source of instability in the problem, the appearance of f G y is associated with destabilization (because f and G y enter the spectrum in Eq. (10) only through the combination f G y ). Finally, the ratio of the two wave-speeds β ≡ c s /c d = (1 − 2ν)/(2 − 2ν) is also a dimensionless parameter of the system which depends only on the bulk Poisson's ratio.
To obtain analytic insight into the instability presented in Fig. 4a, note that solutions in this instability branch are located near the Rayleigh wave-speed, as shown in Fig. 4b (note that here c R 0.95c s ). Consequently, we expand Eq. (10) to linear order around c = c R + δc, obtaining an explicit expression for δc(kH) [57]. A consequence of this expansion is that the transition between stable and unstable modes occurs for k's which approximately satisfy [57] γf G y (c R , k) ≈ −c R /c s .
This approximate stability criterion explains the existence of an instability and in fact gives reasonable quantitative estimates for its onset.
To see this, note that since G y (c R , k) (which is negative, cf. Eq. (11) and [57]) vanishes for both k = 0 and k = ∞, and attains a global minimum for k of order H −1 , Eq. (13) admits solutions only for certain values of the product χ ≡ γf . When χ is smaller than a critical value χ c , no solutions exist and this branch of solutions is stable for all wavenumbers. Note that this criterion has exactly the expected structure: the instability is indeed governed by G y , and large γ or f promote instability, which only happens at a finite range of wavenumbers. These predictions are quantitatively verified in Fig. 4c. In addition, the real and imaginary parts of the approximate solution for δc(kH) [57] are added to Figs. 4a-b (dashed lines), demonstrating reasonable quantitative agreement with the full numerical solution for various parameters.
The results presented in this section demonstrate the destabilizing role that the absence of geometrical reflection symmetry may play in frictional dynamics. In the next section, we significantly extend the analysis to include more realistic friction laws and any geometric contrast.
B. Generalized analysis: State dependence, memory length and arbitrary geometric asymmetry The analysis presented in the previous section adopted two simplifying assumptions, i.e. that the frictional response depends only on the instantaneous slip velocity v and that the lower block is much higher than the upper one, η → ∞. Frictional interfaces, however, are known to depend also on the state of the interface, not just on the slip velocity, and obviously the sliding bodies can feature any geometric asymmetry, i.e. the system can attain any value of η. Consequently, our goal here is to relax these simplifying assumptions and to present a significantly generalized analysis applicable to a broad range of realistic frictional systems.
It is experimentally well-established that the response of frictional interfaces depends, in addition to the slip velocity v, on the state of the interface through the (normalized) real contact area A(φ) ∝ σ yy (1 + ψ(φ)) [1][2][3], as discussed in relation to Eq. (4). The auxiliary internal state variable φ, which represents the age/maturity of the contact and is of time dimensions, carries memory of the history of the interface. This implies that irrespective of the exact functional form of ψ(φ) (with dψ/dφ > 0) the frictional response f (·) in Eq. (9) depends on both v and φ, i.e. we have f (v, φ). Since f does not depend solely on the instantaneous sliding velocity, but also on φ, one should distinguish between where φ 0 (v) is the steady state value of φ. It is also well-established that after a rapid variation in v, accompanied by an instantaneous frictional response characterized by ∂ v f , a new steady state is established over a characteristic slip distance D, which can be regarded as an interfacial memory length. This generic behaviour is described by the following evolution equation with g(1) = 0 and g (1) < 0. While several functions g(·) were proposed and extensively studied in the literature [1,3], the only property that affects the linear stability is g (1). Note that if g(0) > 0 (corresponding to v = 0), the equation describes frictional aging (φ increases linearly with time under quiescent conditions) and that g(1) = 0 corresponds to steady state,φ = 0, implying φ 0 (v) = D/v. The latter describes contact rejuvenation, where the typical contact lifetime is inversely proportional to v. The physics incorporated in the distinction between ∂ v f and d v f , and in the memory length D -within the so-called rate-and-state friction constitutive framework -imply the existence of two dimensionless parameters that are absent in the simplified analysis of Sect. IV A Frictional interfaces generically feature ∂ v f > 0 [1-3], which is termed the "direct effect" (associated with thermally activated rheology [3,75]). As in Sect. IV A, we are interested in d v f > 0 (i.e. in steady state velocitystrengthening friction), which implies a positive ∆. In fact, ∆ varies in the range 0 < ∆ < 1 [23], while ξ can attain any positive value. Within this generalized framework, δf /δv of Eq. (10) takes the form [57] δf δv = In the limit ∆ → 1, i.e. when there is no distinction be-  To understand the effect of ∆ and ξ on frictional stability, we need to solve Eq. (10) using Eq. (17). As we also want to consider arbitrary values of the height ra-tio η, we should first derive expressions for the interfacial elastodynamic transfer function G x,y for any η. The generalized result takes the form [57] Note that Eqs. (11)-(12) are obtained from Eq. (18) by taking the η → ∞ limit, which amounts to setting tanh(ηkHα i ) to unity (since both k and [α i ] are positive). In addition, as expected, G y vanishes for symmetric systems, i.e. for η = 1.
We are now ready to study the effect of the geometric dimensionless parameter η, and of the constitutive dimensionless parameters ∆ and ξ, on the linear stability of frictional interfaces. That is, we aim at solving the implicit linear stability spectrum in Eq. (10), with Eqs. (17)- (18). The ultimate goal of such a generalized linear stability analysis is to derive the stability phasediagram in the γ (here ∂ v f replaces f (v) in the definition of γ in Sect. IV A), f , β, η, ∆ and ξ parameter space, where the stability boundary is a complex hypersurface in this multi-dimensional space.
As it is obviously impossible to visualize this highdimensional stability boundary and in order to gain clear physical insight, we analyze this hypersurface by studying its sections along various parameter directions. A first step was done in Sect. IV A, where the analysis was performed for fixed values of geometric asymmetry η, frictional resistance f and wave-speed ratio β, while γ varied. As a simple velocity-dependent friction model was adopted there, we also had ∆ = 1. As observed in Fig. 4a and analyzed theoretically in relation to Eq. (13), an instability emerges when γ becomes sufficiently large (here somewhere between γ = 2 and γ = 3). As γ = µ/(σ 0 c s ∂ v f ) quantifies the importance of elastodynamics relative to instantaneous velocity-strengthening friction, the instability emerges when elastodynamics becomes more dominant in the presence of large geometric asymmetry, η = ∞. Our next step is to isolate the geometric asymmetry effect embodied in η. We therefore use the parameters of Fig. 4a-b, together with γ = 3, and vary η over a very broad range, essentially from η = 1 (corresponding to a symmetric system) to η = ∞. [c(kH)/c s ], obtained by numerically solving Eqs. (10), (17) and (18), is shown in Fig. 5a. It is observed that for symmetric systems, η = 1, sliding is stable for all wave-numbers. As η is increased, [c(kH)/c s ] approaches the x-axis until they first intersect when η 3.3 at kH ∼ O(1), signaling the onset of instability. This result provides direct evidence for the destabilizing role played by geometric asymmetry in frictional sliding. As η is further increased, the system becomes more unstable in the sense of an increased range of unstable wave-numbers and a larger growth rate. Obviously, the result in the η = ∞ limit identifies with that of Fig. 4a. In fact, the η = ∞ analysis well captures the salient features of the instability spectrum for η values moderately above the critical value η 3.3.
Next, we would like to understand the effect of ∆, i.e. of a difference between the instantaneous response ∂ v f and the steady state response d v f , on the sliding stability in the presence of geometric asymmetry. For that aim, we plot in Fig. 5b [c(kH)/c s ] for various values of ∆, spanning the whole range 0 < ∆ < 1, and fixed ξ = 1 and η = ∞. It is observed that as d v f decreases relative to ∂ v f , i.e. as ∆ decreases, sliding becomes less stable, resulting in a broader range of unstable wavenumbers and a larger instability growth rate. This result demonstrates the stabilizing role played by steady state velocity-strengthening friction in frictional sliding. We note, though, that the qualitative properties of the instability spectrum are rather well captured by the ∆ = 1 analysis (i.e. for velocity-dependent friction, where no distinction is made between d v f and ∂ v f ). We stress that while ∆ affects the properties of instability, the origin of instability is still geometric asymmetry (i.e. sufficiently large η).
Finally, we explore the effect of varying the interfacial memory length D, corresponding to varying ξ, on frictional stability in the presence of geometric asymmetry. We plot in Fig. 5c [c(kH)/c s ] for a broad range of ξ values, and fixed η = ∞ and ∆ = 0.5. It is observed that increasing D (i.e. ξ) tends to stabilize sliding (i.e. shrink the instability range and growth rate) as it makes the real contact area less sensitive to slip velocity perturbations. We also stress here that while ξ affects the range and growth rate of instability, its origin is geometric asymmetry (i.e. sufficiently large η).
The results presented in this section provide a rather comprehensive physical picture of the implications of geometric asymmetry on the stability of frictional sliding, and of the interplay between geometric asymmetry and generic constitutive properties of frictional interfaces, most notably the effect of the state of the interface and of an interfacial memory length. The results significantly extend those presented in Sect. IV A, yet they show that the simplified analysis properly captured the destabilizing geometric asymmetry effect. We stress again that additional solutions to Eq. (10) (with Eqs. (17)-(18)) exist. These additional solution branches, along with a more detailed analysis of the multi-dimensional stability phase-diagram, will be presented in a follow-up report.
The results presented in this section regarding the stability of homogeneous sliding in the presence of geometric asymmetry may have far reaching implications for the dynamics of frictional interfaces in a variety of frictional systems. Under homogeneous loading applied to the top of long enough sliding bodies, as assumed in the analysis, we predict that no homogeneous steady state will be established experimentally under certain conditions that were carefully quantified. Instead, the interface separating geometrically asymmetric bodies will experience inhomogeneous slip related to the most unstable mode identified in the analysis. This will lead to spatiotemporal stick-slip-like motion, accompanied by distinct acoustic signature as in squeaking door hinges.
In frictional systems where the loading configuration promotes inhomogeneous slip, the obtained results may still be relevant. Inhomogeneous slip in slowly driven frictional interfaces typically takes the form of an expanding creep patch. The conditions under which an expanding creep patch spontaneously generates rapid/unstable slip, an important process known as nucleation, may be related to the minimal unstable wavelength in the stability analysis presented in this section for geometrically asymmetric systems. In particular, the minimal unstable wavelength may determine the size at which the expanding creep patch loses stability.
Finally, when rapid slip develops, it is typically mediated by the propagation of rupture modes. Which mode is actually realized in a given experimental system may be affected by the stability analysis presented here. In particular, extended crack-like rupture modes leave behind them a homogeneous sliding state, which may be precluded under certain conditions predicted by our analysis. Instead, localized pulse-like rupture modes may develop. Consequently, the results presented in this section may affect rupture modes selection, a basic open problem in the field of friction. Additional theoretical and experimental research should be carried out in order to fully explore these potential implications.

V. CONCLUDING REMARKS
In this paper, combining experiments and theory, we showed that frictional interfaces which separate bodies made of identical materials, but lack geometric reflection symmetry about the interface, generically feature coupling between interfacial slip and normal stress variations. This geometric asymmetry effect is shown to account for a sizable, and previously unexplained, normalstress-induced weakening effect in frictional cracks. New experiments support the theoretical predictions. We then showed that geometric asymmetry can destabilize homogeneous sliding with velocity-strengthening friction which is otherwise stable. These analyses demonstrate that the effect of geometric asymmetry resembles, sometimes qualitatively and sometimes semi-quantitatively, that of material asymmetry (the bi-material effect).
Since no system is perfectly symmetric, we expect the geometrically-induced coupling between interfacial slip and normal stress variations to generically exist in a broad range of man-made and natural frictional systems. Consequently, it should be incorporated into various theoretical approaches, into engineering models and employed in interpreting experimental observations. The implications in geophysical contexts, such as in subduction zone sliding (cf. Fig. 1c), call for further investigation.
Supplemental Materials for: "On the spatial distribution of thermal energy in equilibrium"

S-I. THE RESPONSE FUNCTIONS Mij IN PLANE-STRAIN ELASTICITY
The purpose of this section is to explicitly calculate the relation between the interfacial stress σ yi (which is a vector) and the interfacial displacement u for a two-dimensional (2D) elastic body that occupies the region −∞ < x < ∞ and 0 ≤ y ≤ H. The bottom boundary at y = 0 is a frictional interface and plane-strain conditions are assumed. Thus, the equations of motion are those of linear elasticity [S1], i.e.
where ε ij ≡ 1 2 (∂ i u j + ∂ j u i ) is the infinitesimal strain tensor (not to be confused with the slip displacement discontinuity vector i ), σ is Cauchy's stress tensor and µ, ν and ρ are the shear modulus, Poisson's ratio and mass density, respectively.
At the top boundary y = H the material is loaded by imposing a horizontal velocity v and a compressive normal stress σ yy = −σ 0 (with σ 0 > 0). The homogeneous solution u h consistent with these boundary conditions reads 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 u(x, y, t) = u h (y, t) + δu(x, y, t). The boundary conditions (BC) at y = H are ∂ t (δu x ) = 0 and δσ yy = 0 .
In what follows, we calculate the response function of the field δu. For easier readability we omit henceforth the notation δu and denote it simply by u.
Consider now a single Fourier mode, i.e. assume that all fields depend on x and t as ∝ e ik(x−ct) , for which Eq. (S1) admits a solution of the form where c s = µ ρ and c d = 2−2ν 1−2ν c s are the shear and dilatational wavespeeds and we defined α i ≡ 1 − c 2 /c 2 i where i ∈ {s, d}. A i are 4 unknown amplitudes which are determined by employing 4 boundary conditions. These are the 2 conditions at y = H, which are given in Eq. (S3), and 2 conditions at y = 0, given by the interfacial stresses σ yi e ik(x−ct) (which may arise from the frictional contact with another body or any other force-generating loading conditions). After calculating the amplitudes, one can express the relation between the interfacial displacements u i and the interfacial stresses σ yi in the form and T i ≡ tanh(kHα i ). Note that in case that the body under consideration occupies the region in space −H ≤ y ≤ 0 (with H > 0) the analysis remains valid, but H should be replaced by −H. This simply amounts to changing the sign of the diagonal entries of M in Eq. (S5). In the main text, a specific example of Eq. (S5) for a semi-infinite half-space y < 0 under quasi-static (QS) conditions was considered in Eq. (1). This is achieved from Eq. (S5) in two steps. First, a semi-infinite half-space y < 0 corresponds to H → −∞, which implies T i → −1, for which we obtain

S2
Second, the QS limit is obtained by taking c → 0. Note that since α i → 1 in this limit, all entries of the matrix in Eq. (S6) vanish, but the prefactor diverges. Their product approaches a finite limit, yielding This expression identifies with Eq. (1) in the manuscript.

S-II. THE RESPONSE OF A COMPOSITE SYSTEM
Next we aim at calculating the response of a general composite system, composed of two bodies made of different materials with different H's in frictional contact. Both bodies are assumed to be infinite in the x-direction. The upper material, denoted by the superscript (1), is assumed to occupy the region 0 < y < H (1) and the lower one, denoted by the superscript (2), the region −H (2) < y < 0 (with positive H (i) ). Since α s,d and T s,d (and µ) may be different for the different materials, they are labeled with a superscript, Following the previous section, the displacements at the frictional interface y = 0 are given as Since both σ xy and σ yy are continuous at y = 0, the displacement discontinuity i (cf. Eq. (2) in the main text) can be written as Since we demand y = 0, the relation between σ yi and x takes the form (S11) We thus write with the definitions G x ≡ G xx /(µ (1) k) and G y ≡ G yx /(iµ (1) k). Note that in symmetric systems M (2) is obtained from M (1) simply by taking H → −H. Therefore, the diagonal terms of M (1) and M (2) have opposite signs and the off-diagonal terms are identical, or in other words, M (1) − M (2) is diagonal. Thus, G is also diagonal, i.e. G y = 0, and no coupling exists between tangential motion and normal traction in this case (i.e. for symmetric systems).

S-III. LINEAR STABILITY ANALYSIS OF HOMOGENEOUS SLIDING
In this section we provide more details regarding the linear stability analysis discussed in Sec. IV of the main text. The basic constitutive relation for the frictional stress reads (S14)

S3
Taking the variation of this relation with respect to sliding velocity perturbations δv relative to steady state sliding at v, one readily obtains to linear order where we used the fact that to zeroth order we have σ yy = −σ 0 . Using Eqs. (S12) and the relation δv = −ickδ x we obtain where f is evaluated at the steady state sliding velocity v. This is Eq. (10) of the main text.
In the simple case of velocity-dependent friction, where f = f (v) (i.e. no state dependence), we simply have δf δv = f (v). In the general case, we have f = f (v, φ) withφ = g( vφ D ). In steady state we have φ = D/v such that g(1) = 0, and in addition we expect g (1) < 0. The perturbation of f takes the form Avoiding direct reference to ∂ φ f , we use the fact that in steady state φ = D/v and thus We then rewrite Eq. (S17) as In order to obtain δφ/δv, we perturbφ = g vφ D by setting v = v + δv and φ = D/v + δφ, which to leading order yields This is a linear equation that can be solved for δφ/δv. Plugging the solution into Eq. (S18), we finally obtain This is Eq. (17) in the manuscript.
A. Simplified analysis: η = ∞ and ∆ = 1 In Sect. IV A of the main text we examine the stability of the steady state sliding at a velocity v for which the frictional stress takes the form σ xy = f (v)σ 0 (i.e. no state-dependence, ∆ = 1). As detailed above, in this case the equation that defines the stability spectrum c(k) is We also assume the bottom layer is infinitely deep, which means η → ∞. As a result, T (2) s,d of Eq. (S13) should be replaced by −1.
As noted in the main text, Eq. (S21) admits multiple solution branches, a few of them might be stable or unstable, depending on the system parameters and on k. Here we focus on a solution branch which is located near the Rayleigh wavespeed c = c R in the complex c-plane, cf. Fig. 4b in the manuscript. First, we non-dimensionalize the equations by defining z = c/c s , q = kH and using the definitions introduced in the main text γ ≡ µ/(σ 0 c s f (v)) and β = c s /c d . With these, we obtain where α s = √ 1 − z 2 and α d = 1 − β 2 z 2 , and the implicit spectrum equation reads We now expand Eq. (S24) to linear order in δz, where z = z R + δz and z R ≡ c R /c s is the dimensionless Rayleigh wavespeed. The solution for δz reads where G j denotes ∂G j /∂z evaluated at z = z R . This solution is plotted in Fig. 4a-b of the main text. Since z R < ∼ 1, the functions G x (z R , q) and G y (z R , q) are real, as well as their derivatives. Thus, the real and imaginary parts of δz in Eq. (S25) are The assumption underlying the expansion to leading order in δz is |δz| z R , where z R is smaller than, but close to, unity.
[δz] is small by construction as we are interested in understanding the instability threshold, determined by a zero crossing of [δz]. To assess the smallness of [δz], note that α s (z R ) = 1 − z 2 R 1 due to the proximity of z R to unity (while α d (z R ) remains of order unity because of the factor β). Furthermore, note that G x (z R , q) contains a term proportional to α −1 s (z R ), cf. Eq. (S22). Since ∂ z (α s (z) −1 ) = z/α s (z) 3 , we have G x (z R , q)/G x (z R , q) ∼ α 2 s 1 and therefore the second term in the numerator of [δz] is negligible with respect to the first. This is demonstrated explicitly in Fig. S1. Thus, the criterion for the critical wavelnumber at which [δz] changes sign is approximately γf G y (z R , q) ≈ −z R , which is Eq. (13) of the main text. Finally, we note that the fact that G x (z R , q)/G x (z R , q) 1 is self-consistent with our working assumption that [δz] is small near the threshold. Near threshold, i.e. when γf G y + z R ≈ 0, we have

S-IV. RESPONSE FUNCTION OF THE "THIN-ON-THICK" GEOMETRY
In this section we calculate the response function of the "thin-on-thick" geometry discussed in Sec. III of the main text and depicted in Fig. 1a. The thinner (upper) block is assumed to be under plane-stress conditions, for which the elastic response is known to be identical to that of plane-strain, cf. Eq. (S7), but with renormalized elastic constants ν plane-stress = ν/(1 + ν) and µ plane-stress = µ [S2].

S5
To calculate the response of the thicker (lower) block, we model it as a semi-infinite half-space. Although in the main text (and in the experimental setup) the thicker block is the lower one, in order to conform to the notations of Sec. S-I we calculate here the response matrix of a semi-infinite half-space that occupies the upper half-space y > 0 and at the end transform the result to be valid for the lower half-space y < 0. As discussed in Sec. S-I, the response matrix of a material that occupies the lower half-space y < 0 is obtained by inverting the signs of the diagonal elements of the matrix.
The surface Green function of a half-space, which relates the displacement field u(x, y = 0, z) to a point force at the origin F δ(r), reads [S1,Eq. (8.19)] where r ≡ √ x 2 + z 2 . As explained in the main text, the following properties are used: (a) We examine the response along the symmetry line of the interface, i.e. y = z = 0.
(b) We set F z = 0, i.e. no out-of-plane forces emerge, consistent with the plane-stress assumption of the thinner block.
(c) F x , F y are constant along the width − W 2 ≤ z ≤ W 2 and consequently take the form where H is Heaviside's step function.
Since u z vanishes on the symmetry line z = 0 due to reflection symmetry through the x − y plane, we will only be interested in the x, y components of the displacement field. Points (a)+(b) imply Using point (c), we obtain the interfacial displacements by integrating over the contact region where we defined X ≡ x −x. The result of this integration is, by definition, the interfacial response matrix M (note that the e ikx prefactor is not included in M ). That is, if we define then M is given by The integration is easier in polar coordinates and it is more convenient to measure the polar angle from the z-axis, cf. Fig. S2. Then, the integrals of Eq. (S31) take the form The integrand is symmetric with respect to the reflection z → −z and thus the integral over θ can be performed on the domain − π 2 ≤ θ ≤ π 2 . Over this domain the cosine function does not change sign and we have where q ≡ kW was introduced. Employing the change of variables u = tan θ we obtain Using some straightforward manipulations, the explicit integration can be performed using [S3, Eq. (3.771. 2)], yielding finally where K 0 (z) is the modified Bessel function of the second kind of order 0. Using the expressions for I 1 , I 2 and I 3 in Eq. (S32), we obtain finally This analysis was performed for the upper half-space y > 0. As discussed above, the response matrix of the lower half-space, y < 0, is obtained by inverting the sings of the diagonal elements, yielding A. Mapping to an effective 2D material How does the response matrix of Eq. (S38) relate to that of an infinite 2D plane-strain material, i.e. Eq. (S7)? We want to write Eq. (S38) as One can show that B(q) approaches unity and C(q) vanishes in the limit q → ∞. Thus, in the large q limit (i.e. for wavelengths much smaller than W ), Eq. (S38) coincides exactly with Eq. (S39), which means in this limit the response of the thicker (lower) block is described by 2D plane-strain elasticity (as expected physically for wavelengths much smaller than W ). Such a mapping does not emerge as cleanly for finite q's. Clearly, for M of Eq. (S38) to have the same structure as M eff of Eq. (S39), the C(q) term in Eq. (S38) should be negligible with respect to the B(q) term. As shown in   S3, this is actually the case except at very small q. After neglecting the C(q) term, a mapping between Eq. (S38) and Eq. (S39) is obtained by equating the two independent terms in each matrix, i.e. by solving the two equations The q-dependent solution to these equations is which is identical to Eq. (7) in the main text. As stated above, in the limit q → ∞ we have B(q) → 1 and clearly µ eff → µ and ν eff → ν. In the opposite limit, q → 0, C(q) is no longer negligible compared to B(q) and a clean mapping to 2D does not emerge, i.e. the problem is fully 3D. The effective constants for intermediate values of q are plotted in Fig. 3 of the main text. It is seen that for the chosen value of ν = 0.33 we have µ eff > µ for all experimentally relevant values of q. For completeness, we note here that at very large q, when µ eff (q) approaches µ, µ eff minutely deviates from µ and approaches it from below in the limit. That is, the effective material contrast µ eff /µ is practically unity, but slightly smaller. For realistic values of ν, this effect is negligible and occurs at large q's: For ν = 0.3 the minimal value of µ eff /µ is 1 − 9.7·10 −7 and is obtained for q ≈ 19.3. The corresponding numbers for ν = 0.2 and ν = 0.4 are, respectively, min{µ eff /µ} = 1 − 7·10 −4 , 1 − 8·10 −8 which are obtained at q = 8.6, 45.

A. Sample construction
The experiments reported on in the main text study the frictional motion of two poly(methyl methacrylate) (PMMA) blocks and compare two geometrically different experimental setups. Both experimental setups were conducted using same upper block of dimensions 200 mm×100 mm×5.5 mm in the x, y and z direction, respectively (see Fig. 1a in the main text) while the lower block was of different geometry in the two setups. In the "thin-on-thin" (symmetric) experiment, a lower block of 250 mm ×100 mm×5.5 mm dimensions was used [S4]. The "thin-on-thick" experimental setup used a thicker lower block of 290 mm ×28 mm ×30 mm dimensions [S5]. The two blocks were pressed together by an external normal force (∼ 4.5 MPa nominal pressure).
The shear and longitudinal wavespeeds, c s and c d respectively, were obtained by measuring the time of flight of ultrasonic pulses, yielding c s = 1345 ± 10 m/s and c d = 2700 ± 10 m/s. Due to the high frequency (5 MHz) of the ultrasonic pulses used, the measured c d corresponds to plane-strain conditions (ε zz = 0). Using these measured values, c d for plane stress (σ zz = 0) was then calculated to be c d = 2333 ± 10 m/s. The corresponding Rayleigh wave speed is c R ≈ 1237 m/s. This velocity is indeed consistent with the maximal measured front velocities for the "thin-on-thin" setup. The measured maximal velocities for the "thin-on-thick" setup are systematically larger, by about 2%, quite S8 close to the values of c R for plain-strain conditions (1255 m/s). This value of c R , as well as the assumption of planestrain conditions, were used in previous work [S5] where the "thin-on-thick" setup was utilized. The experimental loading system, strain and contact area measurements are described in detail in [S5]. We specify here only the main differences of the current study.

B. Strain measurements
We used miniature Vishay 015RJ rosette strain gages for local strain measurements that were mounted 3.5mm above the frictional interface (top block only). Each rosette strain gage is composed of three active regions (each 0.34mm ×0.38mm size). Each active region provides a strain component, ε i , along the directions denoted by the yellow arrows in Fig. S4. S4. Geometry and dimensions (in mm) of a single rosette strain gauge. The black rectangles represent the active area of the measuring components, ε1, ε2 and ε3. Yellow arrows represent the direction of the measured strains.
Electrical resistance strain gages can be calibrated to a high precision when are used on very stiff materials such as various metals. However, when these strain gages are embedded on less stiff materials such as plastics, their presence might locally alter the strain field in their surroundings (see [S6] and references within). Analytical models and numerical efforts exist in the literature to estimate this effect and properly calibrate the measurement of strain.
For purposes of calibration, a rosette strain gage was glued at the center of 100 mm diameter PMMA disk (7.5 mm width). The disk was subjected to radial compression at various angles with respect to the rosette axis (y axis in Fig. S4). We assumed that a transformation that relates the altered strain field due to the rosette presence (here denoted by ε i ) to the "actual" strain field in its absence (ε i ) could be found. We indeed found that the calibration measurements can be described by a phenomenological transformation of the following form where a i are corrections for the gage factors, k i represent the transverse sensitivity of the strain gages and g i represent shear sensitivities. ( x, y) is the coordinate system rotated by 45 • relative to that of ε 1 (see Fig. S4). a 2 = 1 was chosen, as only the relative calibration of the components was of interest. Due to reflection symmetry with respect to the y axis, coefficients of ε 2 and ε 3 are identical and g 1 = 0. This reflection symmetry does not exist with respect to ε 2 and ε 3 , and hence shear sensitivity can not be excluded. As the effects of the elastic mismatch of the rosette configuration have not been previously considered, we note that shear sensitivity has not been discussed in the literature. Here, we find that shear sensitivity exists and is crucial for proper gage calibration. Our calibrations revealed that a 1 ≈ 0.95, k 1 ≈ −0.08, k 2 ≈ 0 and g 2 ≈ 0.1 (details of the calibration procedure will be published elsewhere). Once ε i are measured, ε i can be calculated using an inverse transformation.

C. Experimental results
Typical examples of strain measurements, for both experimental setups, are presented in Fig. S5. In previous work [S4, S5], it was found that the strains in vicinity of a rupture tip are well described by the singular Linear Elastic Fracture Mechanics (LEFM) solutions for ideal shear cracks with a single fitting parameter, the fracture energy Γ [S7]. It was found that for a wide range of rupture velocities, c, Γ is approximately constant (Γ ≈ 1.1 J/m 2 ). Some systemic discrepancies (at most 30%) are observed at extreme rupture velocities (cf. Fig. S5b) between the measured strains and LEFM predictions. These discrepancies may result from either errors involved in the strain gage calibration (see previous section), or, possibly, violations of our implicit 2D assumption. Comparison of strain and contact area measurements for "thin-on-thin" (blue symbols) and "thin-on-thick" (red symbols) geometries. a. Strain tensor variations, ∆εij, after subtracting the initial strains from εxx and εyy and the residual strain from εxy. Strains were measured 3.5 mm above the frictional interface and plotted with respect to the location of the rupture tip, xtip. The singular term of the LEFM solution, that is plotted in black (Γ = 1.1 J/m 2 , c is noted in the panels), describes rather well both geometrical setups. The apparent discrepancy in the shear component, ∆εxy, for x−xtip > 0 was shown to be related to nonsingular LEFM terms, as discussed in [S4]. These strain profiles correspond to two of the examples presented in Fig. 2b of the main text. b. Measured ∆εxx and ∆εyy were characterized by their peak values, ε m xx and ε m yy , respectively, as denoted in a(left). The prediction based on the singular LEFM solution, which corresponds to the black line, successfully captures the measurements with some systematic discrepancies at c > 0.98cR c. The dependence of the undershoot ∆A/A∞ on ε m xx for both geometries. Panel c here eventually transforms into Fig. 2c in the manuscript. This is done, as explained below (see text), in two steps. First, the measurements of ∆εxx at y = 3.5 mm are extrapolated to the interface, i.e. ∆εxx at y = 0. Then, ∆εxx(y = 0) is related to the slip velocity,˙ x, according to˙ x = −2c·∆εxx(y = 0).
In this work we are especially interested in the rupture dynamics at high rupture velocities 0.9c R < c < c R . While direct measurement of c can be performed in our system to ∼ 2% precision, we can significantly decrease this experimental uncertainty by exploiting the significant growth of the strain amplitudes as c → c R [S5]. Using this, we improve our measurements of c by fitting the measured ∆ε xx amplitudes to the singular solution, while assuming that Γ does not significantly change in the vicinity of c R and that the system obeys plane-stress boundary conditions. Using this method, c is the only fitting parameter. Results of this procedure are demonstrated in Fig. S5a and are employed to determine c in Fig. 2 of the main text.
Note that the assumption of plane-stress conditions should be violated for the "thin-on-thick" setup. As mentioned above, the measured asymptotic velocities for the "thin-on-thick" setup are about 2% above c R for plane-stress. Nevertheless, for simplicity, we have used the plane-stress assumption in the above analysis. This, therefore, may lead to systematic errors in our estimated values of c (for example, the directly measured velocities for the 3 highest velocities in Fig. 2a are all 1255 m/s). The use of ∆ε xx , however, enables us to quantitatively differentiate between the different high c measurements, despite possible systematic errors in determining the absolute values of c. Consequently, the rupture propagation velocities stated in the legend of Fig. 2b in the manuscript should be understood in relative terms when normalized with respect to the relevant c R . Figure S5c demonstrates that ∆A/A ∞ is correlated with the amplitude of ε xx , ε m xx , directly measured at 3.5mm above the frictional interface. Relating ∆A/A ∞ to the slip velocity, defined at y = 0, is of great interest. The LEFM singular solution, which describes our measurements well at y = 3.5mm, predicts that the slip velocity (and actually all strain and stress components) should be singular at y = 0 and x = x tip . These singularities are naturally regularized at the crack tip. In this section we will explain the underlying assumptions that enable us to estimate the slip velocity and relate the direct measurements of ∆ε xx presented in Fig. S5c to the extrapolated slip velocities in Fig. 2b of the main text. We first note that even in the extreme case (c = 0.993c R ) presented in Fig. S5, where ∆A/A ∞ ≈ 0.25, strain measurements obtained in both experiments with different geometrical setups are surprisingly similar, where only some differences are observed at x − x tip < 0. These strain differences are only minor when compared to the large qualitative difference in strain measurements presented in [S8], where strong material contrast is considered. This observation enables us to adapt a perturbative approach in which we invoke the simplest cohesive zone model valid for the "thin-on-thin" case to estimate the slip velocity at the interface for both geometrical setups. At this stage, however, we are not able to estimate the accuracy of this assumption for the "thin-on-thick" setup.

Slip velocity estimation
In the non-singular cohesive zone model [S9, S10], weakening initiates once the shear stress has reached a finite peak strength, τ p , above the residual value, τ r of the shear stress. The shear stress gradually decreases according to a prescribed shear stress profile, τ (x/X c ) = τ p · τ (x/X c ). X c is defined to be the cohesive zone size. Far ahead of the rupture tip the solution matches the square root singular form, i.e, σ xy (x X c , y = 0) = K II / √ 2πx. Therefore, τ p , X c and Γ = K 2 II /E, are related through [S10] In previous work [S5] it was argued that the length scale over which A is reduced provides an estimate of X c . It was shown that X c contracts as c → c R . As would be expected from elastodynamic theory, these measurements were