Unconventional singularities and energy balance in frictional rupture

A widespread framework for understanding frictional rupture, such as earthquakes along geological faults, invokes an analogy to ordinary cracks. A distinct feature of ordinary cracks is that their near edge fields are characterized by a square root singularity, which is intimately related to the existence of strict dissipation-related lengthscale separation and edge-localized energy balance. Yet, the interrelations between the singularity order, lengthscale separation and edge-localized energy balance in frictional rupture are not fully understood, even in physical situations in which the conventional square root singularity remains approximately valid. Here we develop a macroscopic theory that shows that the generic rate-dependent nature of friction leads to deviations from the conventional singularity, and that even if this deviation is small, significant non-edge-localized rupture-related dissipation emerges. The physical origin of the latter, which is predicted to vanish identically in the crack analogy, is the breakdown of scale separation that leads an accumulated spatially-extended dissipation, involving macroscopic scales. The non-edge-localized rupture-related dissipation is also predicted to be position dependent. The theoretical predictions are quantitatively supported by available numerical results, and their possible implications for earthquake physics are discussed.

T he failure of frictional systems, composed of bodies interacting along contact interfaces, is mediated by the propagation of interfacial frictional rupture [1][2][3] . A prominent example of such spatiotemporal frictional rupture processes is earthquakes along geological faults [3][4][5][6][7] . A widespread framework for understanding frictional rupture invokes a close analogy to ordinary cracks  , despite notable differences in the underlying physics. Most importantly, an ordinary tensile (opening) crack that propagates inside a material loaded externally leaves behind it fully broken, stress-free surfaces, while the interface left behind a frictional rupture front remains in contact and hence features a finite frictional stress (strength) τ. Moreover, not only τ does not vanish as in ordinary tensile cracks, but in fact, is a dynamical field that varies in space and time, and that is self-selected by the failure dynamics; it depends on the local slip rate/velocity v and on the instantaneous structural state of the frictional interface.
The rather well-developed theory of ordinary cracks, the so-called linear elastic fracture mechanics (LEFM) 29,30 , offers powerful tools that would be very useful for understanding, interpreting, and quantifying frictional rupture if the analogy holds. LEFM is based on scale separation between edge-localized dissipation, which takes place on a short lengthscale ℓ, and linear elastic driving energy, which is stored on significantly larger scales, larger than the crack length L. In particular, cracks in the LEFM framework are characterized by edge-localized energy dissipation per unit area G c (the so-called fracture energy), which is balanced by an elastic energy flux G into the edge region. The latter is transported from large to small scales by singular fields that are characterized by a universal À 1 2 exponent 29,30 , valid at intermediate scales between ℓ and L. These relations between the singularity order, lengthscale separation, and edge-localized energy balance are illustrated in Fig. 1.
In relating frictional rupture to LEFM, one should consider the residual stress τ res , which is finite for frictional rupture, but vanishes for ordinary tensile cracks propagating inside bulk materials under external loading. τ res is not an intrinsic interfacial quantity but rather it is an emergent property that is self-selected by the dynamics of the system, through a coupling between the interfacial constitutive relation and bulk elastodynamics 22 . Once τ res is known, frictional rupture dynamics are quantified relative to a sliding state characterized by τ res . In particular, frictional rupture is then described by the difference between the frictional stress τ and τ res , i.e., by τ À τ res , as will become evident below.
While it is known that LEFM cannot be strictly valid for frictional rupture, where the frictional strength τ is self-selected and generally depends on the structural state of the frictional interface and on the slip velocity v, the conventional LEFM À 1 2 singularity, lengthscale separation, and edge-localized energy balance are extensively used in the context of modeling efforts, laboratory experiments, and field observations  . Yet, to the best of our knowledge, the range of validity of the approximated LEFM picture for frictional rupture, and the interrelations between the singularity order, lengthscale separation, and edgelocalized energy balance in frictional systems are still not fully understood. Our goal in this paper is to shed basic light on these fundamental issues by developing a comprehensive theory of rupture-related dissipation, lengthscale separation, and the singularity order of near rupture edge fields.
We show that the generic rate-dependent nature of friction leads to deviations from the conventional LEFM singularity, and that these deviations can be small if a properly identified dimensionless group of physical parameters is small. We also show that the emergence of unconventional singularities in frictional rupture is accompanied by the breakdown of scale separation, which leads to spatially extended dissipation that involves macroscopic scales. We show that when the deviation of the unconventional singularity order from the conventional LEFM À 1 2 one is small, edge-localized dissipation G c can be identified on a length ℓ, but G c can be significantly smaller than rupture-related dissipation (while they are identical in LEFM, cf. Fig. 1). Furthermore, the latter is shown to be positiondependent. The theory is quantitatively supported by extensive numerical simulations of rate-and-state-dependent frictional dynamics, including explaining recent puzzling observations 23 . Finally, some possible implications for earthquake physics are discussed.

Results
In order to study rupture-related dissipation and the associated scales involved, we consider the breakdown energy at an observation point x i along the rupture plane, defined as Here x i is a fixed position away from the hypocenter (the nucleation site of frictional rupture, whose instantaneous size is L(t), and nucleation occurred at t = 0), τ(δ; x i ) the frictional stress at that position and the slip displacement is δðt; where t x i is defined such that Lðt x i Þ ¼ x i . The term "breakdown" refers here to the fact that E BD involves stresses surpassing τ res , i.e., it does not account for the background frictional dissipation (heat) associated with sliding against the residual stress τ res . Consequently, E BD is the rupture-related dissipation. For ordinary cracks, E BD (t; x i ) is predicted to be independent of x i and to increase over a short timescale '=c r ðt x i Þ (for t > t x i , where c r ðtÞ _ LðtÞ is the instantaneous crack propagation velocity), until it A schematic illustration of dissipation-related lengthscale separation and singularity order in conventional cracks. The breakdown energy E BD is plotted as a function of the distance X behind a propagating crack (see inset), in a logarithmic scale. The breakdown energy E BD quantifies the energy dissipation associated with crack propagation and is given by E BD ðXÞ ¼ G c ðΔE BD ðXÞ þ 1Þ, where the spatial integral ΔE BD ðXÞ is given in Eq. (1) and G c is the fracture energy (defined in the text). E BD (X) increases over the spatial range 0 < X < ℓ (i.e., inside the so-called cohesive/process zone 29,30 ), but saturates at E BD (X) = G c for X ≥ ℓ (see dashed lines). That is, conventional cracks feature strict scale separation, where dissipation occurs only on a localization length ℓ. (inset) A schematic illustration of the slip velocity field v(X) behind a crack propagating at an instantaneous velocity _ L ¼ c r from left to right (L is the crack length and the dot stands for a time derivative. In addition, note that X here is increasing from right to left, unlike the main panel). v(X) features the conventional linear elastic fracture mechanics (LEFM) À 1 2 singularity (see triangle and note the logarithmic scale) on intermediate scales, ℓ ≪ X ≪ L. This square root singularity is directly related to the strict dissipation-related lengthscale separation illustrated in the main panel.
saturates at G c , as illustrated in Fig. 1. Our first goal is to develop a theory of the breakdown energy E BD for generic rate-and-state frictional interfaces.
Theory of the breakdown energy for rate-and-state frictional interfaces. We start by parameterizing the breakdown energy E BD (t; x i ) according to the distance X(t) ≡ L(t) − x i between the observation point x i and the rupture edge, instead of using t itself (see Fig. 2). Moreover, as ordinary cracks feature localized dissipation quantified by G c , we focus on the dimensionless excess breakdown energy, defined as ΔE BD ðX; x i Þ ðE BD ðX; x i Þ À G c Þ=G c , where the term "excess" refers here to the dissipation on top of the effective fracture energy G c . To calculate ΔE BD , consider a frictional rupture front steadily propagating at a constant velocity c r , for which the slip displacement increment at any point on the fault/ interface takes the form dδ = v(X; c r , L)dX/c r (unsteady rupture propagation will be discussed below). With this relation at hand, one can use the definition of E BD to define ΔE BD through the following spatial integral for ℓ ≤ X ≤ L (cf. Fig. 2, where ℓ, X, and L are illustrated), where we used the fact that the integral over 0 ≤ X < ℓ equals G c . Consider then a frictional interface that is described by a generic rate-and-state-dependent constitutive relation [32][33][34][35][36][37][38][39] , characterized by an N-shaped steady-state-friction curve τ ss (v) 40,41 and a single-structural-state field ϕ(x, t) 32-44 , as detailed in refs. 22,23 and in the "Methods". For a broad range of materials, τ ss (v) is characterized by a nonmonotonic and rather weak logarithmic rate dependence 41 . It is well-established that generic rate-and-state frictional interfaces host propagating rupture once the condition for rupture nucleation are met 2 . Suppose then that rupture nucleates at x = 0 (the hypocenter, cf. Fig. 2) at time t = 0, giving rise to two symmetrically propagating frictional rupture fronts. In Fig. 2, we present the frictional stress τ(X; c r , L) and slip velocity v(X; c r , L) fields of the right-propagating front at a later time t, as obtained by recent 2D anti-plane simulations of rate-and-state frictional interfaces 23 . In principle, the fields τ(X; c r , L) and v(X; c r , L) can be extracted from such a simulation and plugged into Eq. (1). Then the integral can be evaluated numerically to yield ΔE BD . Our goal, though, is to calculate ΔE BD analytically in order to gain insight into the underlying physics and then to test the resulting predictions against the simulational data.
The starting point for our development is the idea that for rateand-state frictional interfaces we have ½τðX; c r ; LÞ À τ res =τ res ( 1 for X > ℓ, as is indeed observed in Fig. 2. A quantitative criterion for this condition to hold is derived below. Had it been τ(X > ℓ; c r , L) = τ res , we would have ΔE BD ¼ 0 and the conventional slip velocity singularity vðX; c r ; LÞ $ 1= ffiffiffi ffi X p would have been exact for ℓ ≪ X ≪ L (as illustrated in Fig. 1 for ordinary cracks). Therefore, we treat the latter as a leading order solution and aim at expressing τ(X; c r , L) − τ res in terms of v(X; c r , L). We then assume that the evolution of the internal state field ϕ(X, t) is "fast", i.e., that it quickly equilibrates with v(X; c r , L). Under these conditions, we are left with τ(X; c r , L) = τ ss [v(X; c r , L)], where the latter is a nonlinear relation. To allow for an analytic treatment, we further assume that the smallness of ðτ ss ðvÞ À τ res Þ=τ res also implies the smallness of ðv À v res Þ=v res , presumably justifying a linearization of τ ss ðvÞ À τ res around v ¼ v res for the entire range X > ℓ.
With these ideas and assumptions in mind, we obtain the following expansion τ ss ðvÞ À τ res ' dτ ss ðv res Þ=dv where η dτ ss ðv res Þ=dv is an effective viscous-friction coefficient and τ res ) v res dτ ss ðv res Þ=dv, which is typically satisfied, has been assumed. As will be shown next, this effective linear viscousfriction relation allows to gain deep analytical insight into the physics of the problem at hand 45 Using then the conventional singular slip velocity field vðX; c r ; LÞ ' 2c r K=½μ α s ðc r Þ ffiffiffiffiffiffiffiffi ffi 2πX p for anti-plane conditions 29,30 , where α s ðc r Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 À c 2 r =c 2 s p is the relativistic Lorentz factor and K is the stress intensity factor, we can perform the integration to obtain which is expected to hold for ℓ ≪ X ≪ L, and where In deriving Eq. (3), we used the edge-localized energy balance 30 , which is associated with dissipation on the scale X~ℓ.
The effective viscous-friction coefficient η is positive for the Nshaped steady-state-friction curve τ ss (v) because v res typically Fig. 2 The near-edge fields of a rupture front propagating along with a rate-and-state frictional interface. A frictional rupture front that nucleated at x = 0 and propagates to the right in a 2D anti-plane rate-and-statefriction simulation (its symmetric counterpart, propagating to the left, is shown in Fig. 3a of ref. 23 , see details about the computer simulation therein). Its instantaneous half-length is L(t) and the instantaneous propagation velocity is c r (t) ≃ 0.94c s . Shown are the frictional stress (strength) field τ(X; c r , L) (green) and slip velocity field v(X; c r , L) (orange) left behind the propagating edge. Here, X(t) is a coordinate moving with the edge and pointing backward (cf. Fig. 1), whose origin (X(t) = 0) is defined according to v = 0. The two fields approach finite residual values, τ res and v res , respectively, far behind the propagating edge. (inset) A zoom in on the edge region, revealing a localization lengthscale ℓ associated with edgelocalized dissipation, resulting in an effective fracture energy G c (see Fig. 3a, for a precise definition of ℓ). v(x, t) follows, to a very good approximation, the conventional square root singularity of linear elastic fracture mechanics (LEFM) at an intermediate region, i.e., X > ℓ and prior to approaching v res (fit not shown here, see ref. 23 ). The same conventional singularity is featured by τ(X, t) ahead of the edge, X < 0 and |X| > ℓ, but it is not discussed here 23 .
resides on the velocity-strengthening branch of the friction law above its minimum, dτ ss ðv res Þ=dv > 0. While there is ample evidence that the N-shaped steady-state curve is a generic property of frictional interfaces 41 , hence η > 0, it is important to note that having η ¼ dτ ss ðv res Þ=dv < 0 does not violate any law of nature. The point is that ΔE BD G c is not the total dissipation, which includes also G c , the background frictional dissipation associated with sliding against the residual stress τ res and radiated energy. Together, these ensure positive total dissipation and in principle, one can have η < 0, which implies ΔE BD < 0. This would be the case if v res resides on a velocityweakening branch of the friction curve, dτ ss ðv res Þ=dv < 0.
The excess breakdown energy ΔE BD in Eq. (3), which constitutes one of our major results, shows that whenever friction is rate-dependent, i.e., η ∝ dτ ss /dv ≠ 0, the breakdown energy E BD deviates from the fracture energy G c (ΔE BD ≠ 0), dissipationrelated scale separation breaks down and ΔE BD explicitly depends on a macroscopic scale X. As X is generally orders of magnitude larger than ℓ, ΔE BD can, in general, be significantly larger than G c (the limiting/saturation level of ΔE BD will be discussed below). Our next goal is to understand the range of validity of the result in Eq. (3) and its relation to deviations from the conventional singularity of LEFM.
Relation to unconventional singularities. The results in Eqs. (3)(4) were obtained by assuming that the conventional slip velocity singularity vðXÞ $ 1= ffiffiffi ffi X p can be treated as the leading order solution of the frictional rupture problem (behind the propagating rupture front). That is, it was implicitly assumed that in some sense the singularity order of rate-and-state frictional rupture only slightly deviates from the conventional À 1 2 singularity. Yet, it remains unclear at this stage how Eqs. (3)-(4) are related to the singularity order, and in particular how these are related to the smallness of the deviation from the conventional singularity and to ½τðXÞ À τ res =τ res ( 1. The key to answering these questions is Δξ(c r ) in Eqs. (3)-(4) and its physical meaning. While it is common to assume that the conventional square root singularity of LEFM remains approximately valid for frictional rupture, and while this assumption is a posteriori supported by some observations (see refs. 1,23 ), for the effective linear viscous friction in Eq. (2) nothing should be assumed, the singularity order can be explicitly derived in light of the linearity of the problem. That is, we have where ξ(c r ) does not necessarily and a priori equal À 1 2 , i.e., it may correspond to an unconventional singularity emerging from the intrinsic rate dependence of the frictional stress 45 .
Using Eq. (2) and specializing here for anti-plane conditions, one can show that ξ(c r ) satisfies cot½π ξðc r Þ ¼ À2 η c r =½μ α s ðc r Þ (see "Methods"). This relation shows that the singularity order is not a constant, but rather a dynamic quantity that varies with the rupture velocity c r . Moreover, assuming that ξ(c r ) indeed deviates from À 1 2 only slightly, we obtain where surprisingly Δξ(c r ) is the same one given in Eq. (4). Consequently, Eq. (3) is indeed valid when Δξ(c r ) ≪ 1, i.e., when the deviation from the conventional LEFM singularity is small. Equations (5)- (6), together with Eqs. (3)-(4), constitute the major results of this work. Equation (6) shows that frictional rupture is in fact characterized by an unconventional singularity, yet that the deviation from the conventional À 1 2 singularity is small when Δξ is small. According to Eq. (4), the latter is small when the rate dependence of friction is weak, i.e., when the properly nondimensionalized dτ ss ðv res Þ=dv is small. This is generically the case for rate-and-state frictional interfaces. In such cases, the excess breakdown energy ΔE BD in Eq. (3) is proportional to the very same small quantity Δξ of Eq. (4), but ΔE BD is not necessarily small because the smallness of Δξ may be compensated by an accumulated spatially extended contribution. Therefore, we identify Δξ as a hidden small parameter in rate-and-state frictional failure dynamics. The origin of this small parameter is the rate dependence of the frictional stress, which in turn implies that the strict scale separation assumed in LEFM is only approximately valid in frictional rupture dynamics (manifested in the slow decay of τ(X) toward τ res , while satisfying ½τðXÞ À τ res =τ res ( 1). Moreover, some physical quantities (e.g., ΔE BD ) may be more strongly affected than others (e.g., ξ) by the lack of strict scale separation. Finally, we note that for other interfacial constitutive relations Δξ may not be small and additional new physics may emerge. Such situations will not be extensively discussed here but will be mentioned below in relation to seismological observations. But first, we set out to quantitatively test the predictions of the theory against detailed cutting-edge computer simulations 23 .
Testing the theory. In order to test the theoretical predictions in Eqs. (3)(4)(5)(6), we consider the recent computer simulations of ref. 23 , where generic rate-and state-dependent frictional dynamics have been studied. In these 2D anti-plane simulations, frictional rupture fronts spontaneously emerge, allowing accurate calculations of all of the physically relevant quantities discussed above. In particular, it has been shown that the near rupture edge fields (e.g., those shown in Fig. 2) follow the conventional LEFM À 1 2 singularity to a very good approximation 23 . That is, while Δξ in Eq. (6) has not been explicitly calculated, the results of ref. 23 clearly indicate that Δξ ≪ 1, which is precisely the validity condition of Eqs. (3)(4).
In Fig. 3a, we present E BD (X(t); x i ) for four different observation points x 1−4 along with the fault/interface. It is observed that the E BD (X(t); x i ) curves for all x i 's overlap on a small lengthscale and then branch out. The curves then keep on increasing and appear to saturate at x i -dependent values that are substantially larger than the value of E BD at the branching out point. This behavior is qualitatively different from the one of ordinary cracks, cf. Fig. 1. On the other hand, it appears to be consistent with the theoretical predictions of Eq. (3), not considering for the meantime the x i dependence and the saturation (to be discussed later).
A clear signature of the analytic prediction in Eq. (3) is the logarithmic dependence of ΔE BD on X in the intermediate range ℓ ≪ X ≪ L. To test this prediction, we need to identify G c , which is nothing but the value of E BD at the branching out point. Indeed, it is shown in ref. 23 that this value of G c is exactly the one that balances the elastic energy flux G into the edge region, as determined from the extracted stress intensity factor K. Consequently, the LEFM edge-localized energy balance G ≈ G c is satisfied to a very good approximation, which in turn determines the rupture velocity c r 23 . Moreover, G c allows to explicitly extract the localization length ℓ. Having at hand both G c and ℓ, we plot in Fig. 3b ΔE BD (corresponding to the data of Fig. 3a) against ln ðX='Þ, for the lowest and largest x i 's. It is observed that ΔE BD depends logarithmically on X in an intermediate range (for the two extreme values of x i ), as predicted analytically in Eq. (3), lending strong support to the theory.
The slope/prefactor of the logarithmic relation depends on the observation point x i , which in turn implies that Δξ in Eqs. (3-4) depends on x i . Equation (4) predicts, assuming that the effective linear viscous-friction coefficient η is independent of x i , that the observed x i dependence is attributed to the rupture propagation velocity c r , in particular to the combination c r /α s (c r ). Indeed, frictional rupture in the numerical simulation corresponding to Fig. 3 continuously accelerated 23 , i.e., _ c r ðtÞ > 0, where c r ðt x 1 Þ ¼ 0:94c s and c r ðt x 4 Þ ¼ 0:983c s . Using these values inside Eq. (4), the theory predicts the ratio of the slopes in Fig. 3b to be 1.94. This prediction is in reasonably good quantitative agreement with the observed ratio, which equals 1.76. Finally, the individual slopes satisfy Δξðc r Þ $ Oð10 À1 Þ, which is in agreement with a direct estimation of Δξ(c r ) according to Eq. (4), using the constitutive parameters of the numerical simulations 23 .
Taken together, these results provide direct support to the theoretical predictions. In particular, the results show that ΔE BD can be, and in fact is, quite significantly larger than Δξ(c r ) ≪ 1. This happens due to accumulated spatial contribution associated with the huge difference between X-that can reach the fault/ interface size-and the localization length ℓ (and despite the logarithmic dependence on their ratio). We thus conclude that for rate-and-state frictional interfaces, the non-edge-localized excess breakdown energy in Eq. (3) is a product of a typically small number, given by Eq. (4), and an accumulated spatially extended contribution that can compensate the smallness of Δξ(c r ). Consequently, the breakdown energy E BD can in general deviate significantly from the fracture energy G c .
The position dependence of the breakdown energy and its saturation level. How large can the deviation of E BD from G c be? What determines the magnitude of the deviation? In the example shown in Fig. 3a, the deviation can be as large as~100%, but more importantly it is observed that the saturation value of E BD (X; x i ) depends on the observation point x i . That is, in addition to the x i dependence discussed above in relation to nonsteady rupture propagation, the simulational results indicate an intrinsic relation between the observation point and the saturation value of E BD . It is clear that the ln ðX='Þ dependence of ΔE BD in Eq. (3), which was discussed and validated in the intermediate range ℓ ≪ X ≪ L in Fig. 3b, cannot persist indefinitely. This is simply the case because the logarithmic dependence is directly related to the singular part of v(X), which is no longer dominant at large X.
To understand the behavior of ΔE BD ðXÞ at large X, note that the during crack propagation the relation L(t) = x i + X(t) holds for L(t) ≥ x i , where both L(t) and X(t) increase, while x i is fixed. At short propagation times, measured relative to the time at which L(t) = x i , we have X ≪ x i . At intermediate propagation times, ΔE BD ðX; x i Þ varies logarithmically with X, as demonstrated in Fig. 3b, and finally, at long propagation times, we have X(t) → L(t), which implies X ≫ x i and for which the logarithmic law is not valid anymore. If τ(X) approaches τ res for large X in a way that leads to the saturation of ΔE BD ðXÞ, then we expect the logarithmic behavior of ΔE BD ðX; x i Þ to cross over to a plateau on a scale X~x i , i.e., when ln ðX='Þ roughly equals ln ðx i ='ðx i ÞÞ. This prediction is tested and supported in Fig. 3b, demonstrating that ΔE BD ðX; x i Þ indeed crosses over to a plateau on a scale X~x i . This is a surprising and somewhat non-intuitive result that shows that it is not the rupture size L per se that determines the magnitude of the deviation of E BD from G c but rather the observation point x i . Obviously, larger ruptures generally allow larger x i , so in general these can feature larger deviations of E BD (X; x i ) from G c .
This physical insight can be used to quantitatively predict the E BD (X; x i ) curves, shown in Fig. 3a, over the full range of X's and different x i 's. To that aim, we need to include higher order, non-singular contributions to v(X; c r , L). This is done by using Broberg's full-field solution for a self-similar crack propagating at a constant velocity c r , which takes the form vðX; c r ; LÞ ¼ c r ffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 8 G c L π μ α s ðc r Þ q = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi L 2 À ðL À XÞ 2 p 30 , and is valid for ℓ ≤ X ≤ L. Furthermore, as the full-field expression anyway requires numerical integration in Eq. (1), we can relax the assumption that τ ss ðvÞ À τ res can be linearized around v res . While this assumption appears plausible, the smallness of ðτ ss ðvÞ À (1). It is important to note that as Broberg's full-field solution is obtained in the framework of LEFM, the proposed procedure is still perturbative in nature; that is, it employs an LEFM solution as the leading order contribution in order to calculate the relevant dissipation integral in the main approximation 45 . ΔE BD ðX; c r ; L; 'Þ of Eq. (1) is a spatial integral over a snapshot of the rupture fields, which in itself is independent of the observation point x i . The observation point dependence is introduced in two steps, corresponding to different pieces of physics. First, as the integral in Eq. (1) extends up to X = L, the insight about the saturation at X~x i can be captured by setting L = x i . This saturation is totally unrelated to the additional x i dependence introduced by the nonsteadiness of rupture propagation. The latter, as already discussed earlier, is captured by setting c r = c r (x i ) and ℓ = ℓ(x i ). Consequently, we calculated ΔE BD ðX; c r ðx i Þ; L ¼ x i ; 'ðx i ÞÞ of Eq. (1) using the fully nonlinear N-shaped τ ss (v) (see "Methods") and Broberg's full-field solution v(X; c r (x i ), L = x i ), and compared it to the simulational results. The comparison, which is presented in Fig. 4, reveals reasonable quantitative agreement between the theory and the simulations, lending strong support to the former.
Possible implications for seismological observations. As explained above, the breakdown energy constitutes an important contribution to the total frictional rupture energy budget, which includes also the background frictional dissipation (heat) and the radiated energy. In the context of earthquake physics, these three contributions sum up to the potential energy release during an earthquake 10 . It would be interesting to discuss whether, and if so to what extent, our theory might have some implications for seismological observations. The latter typically aim at using source spectra to obtain coarse-grained average estimates of the following quantity 9,12,[46][47][48][49] Note that G f (δ) differs from the breakdown energy E BD ðδ; x i Þ ¼ R δ 0 ½τðδ 0 ; x i Þ À τ res dδ 0 in two respects. First, it makes no reference to a fault observation point x i . Second, the reference stress used in it is τ(δ), rather than the constant residual stress τ res .
Before discussing seismological observations, let us calculate G f (δ) in the framework of the theory developed in this work. As explained above, the dissipation in the spatial range 0 ≤ X ≤ ℓ near the rupture edge gives rise to a well-defined effective fracture energy G c , marked in Fig. 3a. The edge-localized dissipation G c is related to a strong strength reduction (cf. Fig. 2) over a characteristic slip displacement δ c , such that G f (δ c ) = G c (note that G c of Fig. 3a, which is based on E BD , slightly differs in its value from the one associated with G f due to the difference in the definition of these quantities). This strong frictional strength reduction is associated in the rate-and-state constitutive framework with the evolution of the internal state field ϕ. It has been shown 50,51 that while the rate-and-state constitutive framework does not make explicit reference to δ, the strength reduction from τ 0 -reached after the very initial increase in slip velocity near the rupture edge-to τ c at δ = δ c (where τ c is close to, but still larger than, τ res ) follows an effective linear slip-weakening law of the form τ(δ) ≃ τ 0 − (τ 0 − τ c )δ/δ c . Plugging the latter into Eq. (7), we obtain According to Eq. (8), G f (δ) follows a quadratic power law for δ ≤ δ c , i.e., for G f (δ) ≤ G c . For δ > δ c , where the frictional strength slowly reduces from τ c to the residual stress τ res (cf. Fig. 2, τ c is not marked), G f (δ) is associated with the dissipation in the extended spatial range X > ℓ. Our viscous-friction theory predicts that the excess dissipation in this range-on top of G c -is intimately related to the emergence of unconventional singularities in frictional rupture, which in turn mainly depends on the rate dependence of friction (and not on the internal state field ϕ). To obtain G f (δ) in this regime, we use the slip velocity in Eq. (5), the viscous-friction relation in Eq. (2) and the steady-state relation v = c r dδ/dX. These yield τðδÞ À τ res $ δ ξ 1þξ , which upon substitution inside Eq. (7) leads to For the conventional singularity, ξ ¼ À 1 2 , Eq. (9) predicts that G f is independent of δ for δ > δ c , as expected (in fact, the prefactor also vanishes in this case, implying G f = G c ). In cases in which unconventional singularities emerge, Eq. (9) predicts a power law that depends on the unconventional singularity order ξ.
Our theory thus predicts that G f (δ) follows a quadratic power law for δ ≤ δ c , cf. Eq. (8), which is associated with the strong frictional strength reduction taking place near the rupture edge. The quadratic law is a signature of an effective linear slipweakening characterizing this strength reduction process. At G f (δ c ) = G c , the theory predicts a crossover to another power law, valid for δ > δ c (cf. Eq. (9)), which is associated with spatially extended dissipation and is determined by the unconventional singularity order ξ. These predictions are being tested in Fig. 5 for the smallest x i data presented earlier in Figs. 3-4 (green curves). The numerical results quantitatively agree with the theoretical predictions, revealing a quadratic power law at small δ as predicted by Eq. (8), and a weaker power low (here with an exponent 0.172) for G f (δ) > G c , as predicted by Eq. (9). Using the latter, together with the transformation ξ ¼ À 1 2 ð1 À ΔξÞ (cf.  Fig. 3b. Note that the theoretical curves, by construction, extend up to X = x i (recall that ℓ(x i ) ≪ x i ) and that the simulational ones are limited by the fault/interface half-length W, the smaller x i the larger the maximal X.
Eq. (6)), one obtains Δξ = 0.094, which is in perfect quantitative agreement with Δξ extracted in Fig. 3a (see also Eqs. (3)(4)). Note that Δξ exhibits some dependence on the observation point x i (cf. Fig. 3a), which is not discussed here since seismological observations-to be considered next-completely lack the spatial resolution required to reveal this dependence.
As explained above, seismological observations aim at using source spectra to obtain coarse-grained average estimates of G f in Eq. (7), e.g., see refs. 12,[46][47][48][49] . Yet, such seismological observations completely lack the spatiotemporal resolution to probe the slip δ as a function of time at a given observation point x i on the fault. Instead, it is common to plot the seismological estimate of G f as a function of the total average slip δ in an earthquake, making no explicit reference to the spatiotemporal evolution of slip during rupture. Moreover, it is common to superimpose the seismological estimates of G f vs. δ for many earthquakes (including both crack-like and pulse-like events) occurring on different faults in a single plot, while it is not a priori clear that the data should at all collapse on a master curve. Finally, natural faults exhibit richer constitutive behaviors at high slip velocities (e.g., related to flash weakening and thermal pressurization 47 ) compared to the rateand-state framework used in this work, feature geometrical complexity and are 3D in nature. Yet, with these caveats in mind and following other authors 12,46-48 , we identify δ with δ and discuss the qualitative salient features of the theoretical predictions in Eqs. (8)(9) in relation to the available G f vs. δ seismological observations. Various authors compiled seismological observations from many earthquakes on different faults, spanning a broad range of total average slip δ, ranging from the micron-scale to the scale of tens of meters 12,46-48 . Several authors 46,47 reported G f ð δÞ $ δ 2 for relatively small δ, consistently with Eq. (8), i.e., with an effective linear slip-weakening behavior near the rupture edge. Others, e.g. 12,48 , suggested a weaker-than-quadratic small δ power law and interpreted it in terms of a sub-linear slipweakening behavior near the rupture edge. None of these, to the best of our knowledge, managed to single out G c -i.e., the part of G f that is balanced the edge-localized energy flux G and that in turn controls the rupture propagation velocity-from the data, as our theory allows.
Probably most relevant for our theoretical predictions are the data compiled in ref. 47 , where a quadratic power law G f ð δÞ $ δ 2 at small δ appears to cross over to a weaker power law G f ð δÞ $ δ 2=3 at large δ. This behavior appears to be in qualitative agreement with the theoretical predictions in Eqs. (8)(9), suggesting that different physics controls the two power-law regimes, and in particular that the latter is associated with an unconventional singularity and a dissipation contribution from a spatially extended region behind the rupture edge. Interestingly, the two power laws suggested in ref. 47 have been interpreted in terms of a thermal pressurization constitutive model, where the fluid pore pressure plays a central role. The quadratic power law has been interpreted to correspond to an effective linear slip-weakening behavior associated with undrained conditions and the 2 3 power law with drained conditions 47 . Most interestingly, the 2 3 power-law regime has been related to an unconventional singularity of order ξ ¼ À 1 4 , associated with the thermal pressurization model under drained conditions and corresponding to large slips accumulated far behind the rupture edge. In fact, substituting ξ ¼ À 1 4 in Eq. (9), one obtains a 2 3 power law, even though the linear viscous-friction approximation of Eq. (2) does not seem to be directly relevant to the analysis of ref. 47 .

Discussion
In this work, we developed a theory that elucidates the interrelation between unconventional singularities, scale separation, and energy balance in frictional rupture. We have shown that the intrinsic rate dependence of friction, dτ ss (v)/dv ≠ 0, generically leads to deviations from the conventional LEFM near-edge singularity. It is this rate dependence, which in turn implies that the frictional stress is self-selected, that leads to the emergence of singular fields different from those of LEFM. For the widespread rate-and-state-friction constitutive law, these deviations can be small, yet they are accompanied by non-edge-localized breakdown energy that significantly deviates from the edge-localized dissipation.
The developed theory sheds basic light on frictional rupture energy balance and the underlying lengthscales. The crux of the theory is the identification of a hidden small parameter Δξ that quantifies the deviation from the conventional LEFM singularity and that is intrinsically related to the rate dependence of friction. The theory quantitatively explains recent puzzling observations in cutting-edge numerical simulations and offers predictions that are amenable to laboratory testing using available techniques 1 . Finally, the theory offers tools and concepts that can be used to interpret seismological estimates of earthquake breakdown energies.
The concepts and ideas developed in this work are applicable to more complicated interfacial constitutive relations, incorporating even richer multiphysics of frictional systems. These can include healing, pore fluid effects, thermal pressurization, flash heating, off-fault damage and plasticity, and more. The developed theory remains valid as long as the frictional stress continues to evolve behind the rupture front over scales much larger than the localization (cohesive/process zone) length ℓ, implying that strict scale 3). G f (δ) features two power laws, marked by the two triangles, in agreement with the theoretical predictions in Eqs. (8)- (9). The crossover between the two power laws, marked by the vertical dashed line, occurs at (δ c , G c ), as predicted theoretically. G c and δ c are used to normalize G f and δ, respectively (the value used for the former, G c = 0.51J/m 2 , is slightly smaller than the one observed in Fig. 3a due to the difference in the definition of G f and E BD ). Finally, comparing the power-law exponent in the δ > δ c regime, 0.172, to the analytic prediction in Eq. (9) and using the transformation ξ ¼ À 1 2 ð1 À ΔξÞ (cf. Eq. (6)), one obtains Δξ = 0.094. The latter is in perfect quantitative agreement with Δξ extracted in Fig. 3a (see also Eqs. (3)(4)).
separation does not hold. In the most general case, the power-law exponent ξ(c r ) in Eq. (5) is not necessarily close to À 1 2 (i.e., Δξ is not necessarily small). This generalized theory is addressed elsewhere 52 , and applications to specific interfacial constitutive relations are expected to emerge in the future.

Methods
Determining the singularity order. This work is analytic in nature and all of the derivations are detailed in the text, except for the solution for the unconventional singularity order, which is provided below. The theoretical predictions are compared to numerical results that have been published in ref. 23 based on 2D antiplane spectral boundary integral method simulations [53][54][55] . These numerical simulations employed a rate-and-state-friction constitutive relation τ ¼ σsgnðvÞf ðjvj; ϕÞ, where σ is the normal stress and f ðjvj; ϕÞ ¼ ½1 þ blog ð1 þ ϕ=ϕ Ã Þ½f 0 = ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ ðv Ã =vÞ 2 q þ a log ð1 þ jvj=v Ã Þ. The internal state field ϕ satisfies _ ϕ ¼ 1 À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 1 þ ðv Ã =vÞ 2 q jvjϕ=D and the values of the parameters appear in Table I in ref. 22 . Under steady-state conditions, _ ϕ ¼ 0, the frictional strength τ ss (v) follows an N-shaped curve, as plotted in Fig. 2a of ref. 22 and in Fig. 1b of ref. 23 , and as supported by numerous experiments 41 . The numerical results of ref. 23 have been presented in this work in different forms, depending on the theoretical predictions being tested, as detailed in the text.
The unconventional singularity order ξ can be obtained by considering the interfacial boundary condition for 2D anti-plane steadily propagating rupture 56 τðXÞ ¼ μ α s ðc r Þ 2πc r where the left-hand side is the frictional strength and the right-hand side is the shear stress at the interface, as obtained from bulk elastodynamics 56 . Using τ ss (v) of Eq. (2) for τ(X) and invoking the asymptotic power-law ansatz in Eq. (5), Eq. (10) implies that the unconventional singularity order ξ satisfies cotðπ ξÞ ¼ À2 η c r =ðμ α s ðc r ÞÞ, as reported in the text. Finally, plugging into the last relation Eq. (6) and expanding to the leading order in Δξ, the latter is calculated and is shown to identify with Eq. (4), as stated in the text.
Data availability