## Abstract

The rupture fronts that mediate the onset of frictional sliding may propagate at speeds below the Rayleigh wave speed or may surpass the shear wave speed and approach the longitudinal wave speed. While the conditions for the transition from sub-Rayleigh to supershear propagation have been studied extensively, little is known about what dictates supershear rupture speeds and how the interplay between the stresses that drive propagation and interface properties that resist motion affects them. By combining laboratory experiments and numerical simulations that reflect natural earthquakes, we find that supershear rupture propagation speeds can be predicted and described by a fracture mechanics–based equation of motion. This equation of motion quantitatively predicts rupture speeds, with the velocity selection dictated by the interface properties and stress. Our results reveal a critical rupture length, analogous to Griffith’s length for sub-Rayleigh cracks, below which supershear propagation is impossible. Above this critical length, supershear ruptures can exist, once excited, even for extremely low preexisting stress levels. These results significantly improve our fundamental understanding of what governs the speed of supershear earthquakes, with direct and important implications for interpreting their unique supershear seismic radiation patterns.

## INTRODUCTION

The onset of frictional sliding occurs through the dynamic propagation of rupture fronts (*1*). Local slip occurs only once these fronts rupture the microcontacts that define the rough contacting surfaces that separate sliding bodies. Propagation speeds of idealized rupture fronts that are driven by singular stresses at their tips are limited by the Rayleigh wave velocity (*2*, *3*). However, when regularization of these singularities is taken into account, theoretical (*4*) and numerical (*5*) work has predicted the existence of supershear cracks, a class of cracks that propagate beyond the shear wave speed and may approach the longitudinal wave speed. While supershear ruptures have since been observed in laboratory friction experiments (*6*–*10*) and inferred in natural earthquakes (*11*–*14*), it is not completely understood what determines their speed. Understanding what determines the speed of supershear earthquakes has important implications, as their propagation speeds strongly influence the structure of their radiated waves. The supershear radiation structure, which is fundamentally different from sub-Rayleigh earthquake radiation (*15*), has important consequences for resulting seismic hazards (*16*).

While the conditions for the transition of ruptures from sub-Rayleigh to supershear regime have been studied extensively (*5*, *17*–*22*), the question of what determines the evolution of supershear rupture speed is still entirely open. In the sub-Rayleigh regime, using the framework of linear elastic fracture mechanics (LEFM), the equation of motion for accelerating ruptures has been constructed via the superposition of solutions for constant-velocity (singular) shear cracks (*23*). This solution has recently been validated experimentally for frictional ruptures (*2*, *3*). In the supershear regime, however, this approach fails (*24*, *25*), as the elastic fields in the crack tip vicinity are coupled to the crack velocity history.

Here, we construct an approximate equation of motion for supershear cracks by using a particular solution given by Broberg (*26*). This equation of motion builds on previous work on energy flux for supershear cracks (*26*–*28*) but provides a direct link between crack length and rupture speed. It therefore enables us to predict supershear crack speeds for various nonuniform systems. We will show that this equation of motion provides a good description of both experiments and numerical observations of frictional supershear ruptures that have been proposed to model natural earthquakes. Our results demonstrate the broad applicability and the scope of these approximate solutions. These results demonstrate how predictions based on fracture mechanics can be successfully extended to the supershear regime.

## RESULTS

### Sustained supershear propagation

Our experimental setup consists of two poly(methyl methacrylate) (PMMA) plates of same thickness that are pressed together by an applied normal load (Fig. 1A). We use both dry and boundary-lubricated frictional interfaces (*3*). Shear forces *F*_{S} are applied, and slip events are either nucleated spontaneously (*7*) or induced (*3*) by a slight out-of-plane shear perturbation at *x* ≈ 0. A high-speed camera at 580,000 frames/s records the dynamic changes of the real contact area *A*(*x*, *t*). Figure 1B shows an example of a frictional rupture event that nucleates at *x* ≈ 0, accelerates, and transitions at *x* ≈ 50 mm to supershear speeds (*C*_{f} > *C*_{S}), which it maintains until reaching the leading edge.

In addition, the stresses σ_{ij} are calculated from strains measured every 1 μs at multiple locations along (and slightly above) the frictional interface (*2*). Figure 1C shows the variation of the shear stress ( is the residual frictional resistance) and the dynamic stress drop τ_{0}, as the rupture passes by the strain measurement location. At the interface, we expect the shear stress to increase to the peak shear strength of the frictional interface (see definition of τ_{p} in Fig. 1C). Although this highly localized increase of stress is not directly measured by our (slightly) off-fault measurements, τ_{p} as well as the fracture energy Г (the dissipated energy per unit area), were obtained independently (*2*, *3*) by analyzing the sub-Rayleigh crack propagation regime for this system. Briefly, Г was extracted from the 1/r LEFM singular form of the stresses, and τ_{p} was obtained from measurements of the cohesive zone size (details are provided in Materials and Methods). While inhomogeneities along the frictional interface are reflected in spatial variations of both and Г and τ_{p}, we find that Г and τ_{p} are independent of the rupture velocity. We therefore use these inferred values in analyzing the behavior of the system in the supershear regime.

The rupture front velocity *C*_{f}(*l*), inferred from the measurements of *A*(*x*, *t*) for the event presented in Fig. 1B, is shown in Fig. 1D (blue curve), along with an additional rupture event (red curve) whose high propagation speed approaches the longitudinal wave speed *C*_{L}. Here, *l* denotes the rupture length (as ruptures are nucleated at *x* ≈ 0), and . What dictates *C*_{f}(*l*) and how it relates to τ_{0} and the frictional properties of the interface is the central question investigated in this article.

We now briefly describe how LEFM can be used to provide quantitative predictions for the speeds of supershear frictional rupture fronts. Within the wake of frictional ruptures, the bodies are always in partial contact [σ_{yy}(*x*) ≠ 0], and the frictional resistance, , opposes sliding. The problem of frictional rupture fronts, propagating within an interface separating identical materials, can be mapped to the stress-free conditions that define the mode II crack problem (*29*) by using the linearity of the governing equations and defining the shear stress variations, (see definition of τ, τ_{0}, and τ_{p} in Fig. 1C). The stress fields near the tip of an idealized supershear crack with no cohesive zone are singular (*26*, *30*, *31*), τ ~ *K*/(*x* − *l*)^{g}, where *K* is the stress intensity factor and where, in contrast to sub-Rayleigh ruptures, the singular exponent *g* depends on *C*_{f} [with *g*(*C*_{f}/*C*_{L}) ≤ 1/2]. *K* has been calculated explicitly only for a handful of particular examples, such as a semi-infinite crack subject to a pair of concentrated shear forces (*27*) or for a symmetric bilateral crack subject to uniform remote shear stress (*26*). In both cases, cracks were assumed to expand at constant velocity. A consequence of the singular description of supershear cracks is the vanishing energy flux into the crack tip for *C*_{f} ≠2*C*_{S}. It was shown (*31*), however, that when a cohesive zone is introduced, a finite region where these singularities are regularized to a finite τ_{p} value, the requirement for a positive energy flux to the crack tip is fulfilled for any *C*_{f} > *C*_{S}. By combining the cohesive zone model with the constant-velocity singular crack solution, the energy flow to the crack tip per unit area (the energy release rate) *G* can be calculated (*26*–*28*). Here, we follow Broberg’s bilateral crack solution (*26*), where *G* was expressed in terms of *C*_{f}, τ_{0}, *l*, and the cohesive zone size *x*_{c}. Here, motivated by recent experimental results (*2*), which showed that the interface shear strength τ_{p} is roughly independent of *C*_{f} (in contrast to the *C*_{f} dependence of *x*_{c}), we rederive *G* in terms of and *C*_{f}/*C*_{L}, *l*, τ_{0}, and τ_{p}. Balancing *G* and Г yields a crack propagation criterion(1)where μ is the shear modulus and (*C*_{f}/*C*_{L}) and are known functions. The shape of the stress distribution within the cohesive zone is carried within [see Materials and Methods for the derivation and definition of (*C*_{f}/*C*_{L}) and ].

Studies of supershear crack arrest (*24*) showed that steady-state singular fields are not established instantaneously. Therefore, a solution of an accelerating crack cannot be constructed by the superposition of constant-velocity crack solutions, as had been performed previously for sub-Rayleigh propagation (*23*). In this sense, while Eq. 1 was strictly derived to describe cracks propagating at a constant velocity, we hypothesize that its applicability can be extended for “slowly” accelerating cracks, where velocity history can be neglected. As Eq. 1 implicitly links *C*_{f} to *l*, it will provide an equation of motion for supershear cracks.

The salient property of the supershear crack growth criterion (Eq. 1) is its dependence on the cohesive zone characteristics, in striking contrast to sub-Rayleigh brittle cracks (*9*, *30*). In addition, note the emergence of a characteristic length scale , which is proportional to the static cohesive zone size. Finally, we note that the construction of Eq. 1 describes velocity response to spatial changes in interface properties [Г(*x*) and τ_{p}(*x*)], as for sub-Rayleigh cracks. The rupture response to nonuniform spatial profiles, τ_{0}(*x*), is not captured by Eq. 1. We will consider this question in a later section of this article.

To investigate the applicability of the LEFM equation of motion, we supplement the experiments with two-dimensional (2D) numerical simulations of dynamic bilateral shear cracks using a spectral boundary integral method (*32*). The modeled system consists of two linear elastic semi-infinite half spaces with an interface governed by a linear slip weakening friction law (see Materials and Methods) and τ_{0} applied along the interface. Poisson ratio *v* = 0.35 and plane strain are assumed (*C*_{S}/*C*_{L} = 0.48). We use a weak patch (τ_{0}/τ_{p} = 0.8) around the nucleation center to promote a direct transition (*20*) to supershear propagation. To test the predictions of Eq. 1 for various τ_{0}/τ_{p} values, we fix the value of τ_{0}/τ_{p} to 0.2 to 0.7 after a finite propagation distance (see Fig. 2A, top) by increasing τ_{p} and Г (*l*_{0} remains constant). This leads to an instantaneous response of the crack.

Figure 2 demonstrates that Eq. 1 provides an excellent description of both the numerical simulations and the experiments with no adjustable parameters. Lower τ_{0} are associated with slower ruptures, as observed in early numerical work (*33*) and in sub-Rayleigh propagation (*3*). In our simulations, we observe transient discrepancies for the τ_{0}/τ_{p} = 0.2 case due to the strong localized increase of the interfacial strength (*l*/*l*_{0} = 140 to 170). We expect these discrepancies when strong crack accelerations are present, as was previously discussed. We note that other types of transition to the supershear regime, for example, supershear seed crack or transition through a secondary crack, have been applied and result in equivalent observations (see Materials and Methods). In our experiments, rupture nucleation, acceleration, and the supershear transition are sensitive to the high stress gradients in the block edge vicinity and are beyond the scope of the current study. Instead, we concentrate on the sustained supershear propagation within a central region where these edge effects are negligible. Note the spatial variation of the τ_{p} obtained profile and the accompanying response of the equation of motion (Fig. 2B).

Figure 2C demonstrates the generality of the experimental results, where we consider both dry and boundary-lubricated interfaces, each with significantly different values of Г and τ_{p} (*3*, *34*). Here, we compare the predictions of Eq. 1 at a fixed spatial position for multiple rupture events driven at different levels of τ_{0}/τ_{p}. Slight discrepancies for the boundary lubrication case might result from nonuniformity of τ_{0} along the interface, which is not taken into account here.

In systems with uniform τ_{0}, τ_{p}, and Г, the sub-Rayleigh–to–supershear transition [the Burridge-Andrews mechanism (*4*–*6*, *9*, *18*)] occurs only along weak interfaces above a critical value of τ_{0}/τ_{p}, and no supershear crack can exist for lower values of τ_{0}/τ_{p}. These critical values depend on *C*_{S}/*C*_{L}; τ_{0}/τ_{p} = 0.31 and 0.37 for *C*_{S}/*C*_{L} = 0.48 (simulations) and 0.57 (experiments), respectively. Both our simulations and experiments reveal strong hysteretic behavior; once a supershear crack has been excited [for example, the transition is facilitated through spatial nonuniformity (*17*, *19*, *35*)], supershear cracks may propagate in regions of strong interfaces, well below the critical τ_{0}/τ_{p} values given by the homogeneous Burridge-Andrews transition mechanism. For example, the green and red simulated examples (Fig. 2A) correspond to τ_{0}/τ_{p} = 0.3 and 0.2, respectively. Similarly, the particularly slow laboratory supershear rupture, *C*_{f} ≈ (blue example in Fig. 2B), corresponds to τ_{0}/τ_{p} < 0.37, suggesting that, in this case, the transition has not been induced by the Burridge-Andrews mechanism but rather by favorable heterogeneities (*19*) at the block edge.

### Minimal length for supershear propagation

The LEFM-based equation of motion describes two distinct shapes. For high prestress levels, for example, τ_{0}/τ_{p} = 0.8 (blue example in Fig. 2A), *C*_{f} monotonically increases with increasing *l*. As τ_{0}/τ_{p} is decreased, however, the solution bifurcates into two branches, for example, τ_{0}/τ_{p} = 0.2 (red example in Fig. 2A), one of which appears to be unstable (dashed line in Fig. 2A). The accelerating branch (the physical solution) only appears at a finite crack length *l*_{m}, predicting that no supershear cracks can exist for *l* < *l*_{m}. Thus, *l*_{m} represents a critical length for supershear propagation, similar to Griffith’s length for the sub-Rayleigh regime (*31*). We have confirmed this prediction for several values of τ_{0}/τ_{p} by varying the length of the weak patch in our simulation (Fig. 3A). If the weak patch exceeds *l*_{m}, then the crack obeys the supershear equation of motion in the following stronger region. However, in cases where the weak patch ends at *l* < *l*_{m}, the crack transitions instantaneously to the sub-Rayleigh regime. From this point, the crack remains at the sub-Rayleigh regime even for *l* > *l*_{m}, where the supershear solutions exist. This demonstrates the coexistence of the sub-Rayleigh and supershear propagation solutions and the hysteretic transition between the two.

Figure 3B describes the dependence of *l*_{m} on τ_{0}. No critical length exists for sufficiently high stress levels. For τ_{0}/τ_{p} ≲ 0.45, a finite value of *l*_{m} emerges and *l*_{m} increases rapidly with decreasing τ_{0}. We find empirically that *C*_{f}(*l*_{m}) ≈ ; hence, by using Eq. 1, *l*_{m} can be approximated by(2)Figure 3B shows a comparison of Eq. 2 with the exact solution of Eq. 1. This approximation reveals that when it exists, the scale of *l*_{m} is determined by the Griffith length for crack nucleation *l*_{c} ~ μГ/ (*l*_{m} ≈ 3.7*l*_{c} for *C*_{S}/*C*_{L} = 0.48). We emphasize that the minimal length presented here defines whether the supershear solution does or does not exist, which is fundamentally different from the question of what determines the sub-Rayleigh–to–supershear transition and whether it actually occurs (for example, Andrew’s transition length, which is, in a uniform setup, always larger than *l*_{m}).

### Supershear propagation along interfaces with nonuniform loading

Both changes in frictional strength (barriers) and nonuniform prestress (asperities) affect fault ruptures. These heterogeneities influence not only rupture propagation (*17*) but also produce different ground motions (*36*). In Fig. 4A, we show simulations that model nonuniform loading to illustrate the differences in the effect of barriers versus asperities for supershear propagation; changes in τ_{0}/τ_{p} are now induced by varying τ_{0}, in contrast to the simulations shown in Fig. 2A, where τ_{0}/τ_{p} profiles were varied by spatial changes in τ_{p} and Г.

Note that in the context of LEFM, spatial variations in τ_{0} are fundamentally different from those in τ_{p} (see Materials and Methods). Equation 1 assumed uniform prestress and cannot be strictly applied to nonuniform loading configurations. We therefore need to generalize Eq. 1 for spatially varying τ_{0}. To accomplish this, the information about the τ_{0}(*x*) profile should be incorporated within *K*. We assume *K* = κ(*C*_{f}/*C*_{L})*K*_{s}(τ_{0}, *l*, *g*) and hypothesize that *K*_{s} is a weighted integral functional of τ_{0}(*x*) that can be calculated in a manner analogous to the static stress intensity factor used for sub-Rayleigh cracks (see Materials and Methods). This yields(3)which coincides with Eq. 1 for uniform loading. Figure 4A shows predicted rupture speeds based on Eq. 3. Equation 3 generally compares well to simulations with minor discrepancies that appear in regions immediately following nonuniform areas. Note that even though the τ_{0}/τ_{p} profiles appear to be very similar to the nonuniform strength configurations in Fig. 2A, the dynamics are significantly different, reflecting the influence of the spatial variation of τ_{0}. For instance, while Eq. 2 predicts that at τ_{0}/τ_{p} ≈ 0.15 no supershear cracks can exist for *l*/*l*_{0} < 200 (see Fig. 3), Eq. 3 predicts sustained supershear propagation at τ_{0}/τ_{p} = 0.1 (see violet example in Fig. 4).

Only for extremely low background stress levels do we observe the transition to sub-Rayleigh propagation (see yellow and brown curves in Fig. 4A). At the crack length *l*_{a} of this transition, similar to the minimal supershear crack length *l*_{m} (Fig. 3), the energy flux to the crack tip is insufficient to drive supershear cracks. Figure 4B shows a systematic comparison between predicted supershear to sub-Rayleigh transitions and measured results from simulations for 50 different prestress configurations. To evaluate the criticality of prestress configurations, we computed the additional background stress Δτ_{0} required to prevent the transition to sub-Rayleigh (shown in color in Fig. 4B) and sustain supershear propagation. Simulations with lower prediction accuracy have lower Δτ_{0} and hence exhibit an extreme sensitivity to values of the prestress driving the transition. The prediction works generally well, although low values of Δτ_{0}/τ_{p} result in underpredicted *l*_{a} for *l*_{a}/*l*_{0} > 180. Thus, small changes in the background prestress, for instance due to wave radiation from the crack, which have been neglected in the current theoretical description, should have considerable effects on the rupture dynamics and could explain the prediction discrepancies.

## DISCUSSION

Our measurements and simulations have shown that the energy balance at the rupture tip (Eqs. 1 and 3) provides quantitative predictions for the evolution of *C*_{f}(*l*). These results are general so long as several necessary assumptions are satisfied. For instance, a region in the rupture tip’s vicinity should exist where the stresses are singular τ ~ 1/(*x* – *l*)^{g}. While typically this is assumed to be true when *l* ≫ *x*_{c}, this condition is hard to meet for τ_{0} → τ_{p}, where *g* → 0. Furthermore, in all rupture events considered here, shear waves that originate at the rupture nucleation (*9*, *18*) trail behind the supershear rupture tip. However, one might consider more complex scenarios. Rough faults, for example, may lead to rapid acceleration and deceleration of sub-Rayleigh ruptures (*35*) that would result in stress wave radiation. If the sub-Rayleigh–to–supershear transition were to occur at later times, then the supershear rupture will eventually catch up with these waves, and our assumption of time-independent loading will be violated. Finally, in our experiments, there is no significant variation of after the passage of the rupture tip (*37*), and frictional ruptures can be mapped to simple cracks with traction-free faces. Slip pulses (*15*, *38*), which may result because of strong velocity weakening of friction (*39*) and rapid healing of the interface, are fundamentally different from crack-like slip fronts and would require a different predictive theoretical model.

Natural faults are complex entities that include significant heterogeneity of fracture energy, friction laws, and stresses. Furthermore, complicated fault geometry, together with significant challenges in measuring the speeds of supershear earthquakes (*11*, *12*, *14*) and determining fault properties, often hinder comparisons with modeling. In some cases, qualitative and sometimes quantitative features of natural earthquakes can be successfully compared with idealized numerical models (*13*) and laboratory experiments (*40*). This fracture mechanics–based picture presented in this article is one of these, as it provides new and significant fundamental understanding of supershear rupture dynamics. These ideas may be potentially generalized to further account for the complexities associated with natural earthquakes.

Our results not only highlight how key parameters such as τ_{0}, τ_{p}, and Г control supershear earthquakes but also provide a tool for quantitative predictions of supershear rupture speeds. These results have important implications, as the propagation speeds of supershear earthquakes dictate the amplitudes of the stress fields along the radiated Mach fronts (*15*). These amplitudes, while disappearing at *C*_{f} = , greatly increase as *C*_{f} → *C*_{L}. These subtle changes of the supershear earthquake speeds, therefore, have significant impact on off-fault damage and provide a key parameter for predicting the serious seismic hazards that these earthquakes impose.

Finally, while supershear earthquakes are typically associated with highly stressed faults, our results support suggestions (*19*) that supershear earthquakes may exist even at extremely low stress levels, if favorable heterogeneities for supershear transition are present. The proposed equation of motion provides exact conditions for which these low stress earthquakes may or may not sustain supershear propagation and can explain supershear–to–sub-Rayleigh transitions such as those in (*19*). This work, furthermore, provides fundamental insights for understanding how the crucial interplay between fault roughness and stress levels (*35*) may govern supershear rupture.

## MATERIALS AND METHODS

### Experimental system

Our experiments were conducted using PMMA plates (ρ ≈ 1170 kg/m^{3}) of dimensions 200 mm × 100 mm × 5.5 mm (top block) and 240 mm × 100 mm × 5.5 mm (bottom block) in the *x*, *y*, and *z* directions, respectively (see Fig. 1A). The blocks were pressed together by applying ≈5 MPa of nominal pressure. The contacting surfaces were cleaned with distilled water and isopropyl alcohol and then dried for about 2 hours (termed here dry interfaces). We also conducted experiments in the boundary-lubrication regime (*3*, *34*), where contacting surfaces were coated by a thin layer of lubricant (silicon oil with kinematic viscosity *v* ≈ 100 mm^{2}/s). Material shear, *C*_{S}, and longitudinal, *C*_{L}, wave speeds were obtained by measuring the time of flight of ultrasonic pulses, yielding *C*_{S} = 1345 ± 10 m/s and *C*_{L} = 2700 ± 10 m/s, corresponding to the plane strain hypothesis (ε_{zz} = 0). Using these measured values, we calculated that for plane stress conditions (σ_{zz} = 0), *C*_{L} ≈ 2333 m/s (*C*_{S}/*C*_{L} ≈ 0.58) and *C*_{R} ≈ 1237 m/s. The corresponding dynamic shear modulus is μ = 2.1 GPa (*9*). The three components of the 2D strain tensor were continuously and simultaneously measured by Vishay 015RJ rosette strain gages every 1 μs at 16 to 19 spatial locations, ~3.5 mm above the frictional interface (*2*).

In previous work (*2*), it was found that the strains in the vicinity of rupture tips, propagating below *C*_{R}, are well described by the square-root singular LEFM solutions, originally derived to describe brittle shear cracks. This comparison of the solutions with the measured strains provides a measure of the fracture energy Г, which was found to be roughly independent of rate. Furthermore, it was suggested that measurements of the contact area *A* variations provide a direct measurement of the cohesive zone size, *x*_{c} (*2*). The theoretical interface strength τ_{p} is defined as the difference between the actual interface strength and the residual stress (see Fig. 1C). τ_{p} is obtained from measured values of Г and *x*_{c} using the simplest cohesive zone models (*29*, *41*). For the current set of experiments, we measured slight spatial inhomogeneities of the frictional interface. In the dry case, for 50 mm < *x* < 100 mm, Г ≈ 1.5 J/m^{2} and τ_{p} ≈ 1.3 MPa, while for 100 mm < *x* < 150 mm, Г ≈ 2.5 J/m^{2} and τ_{p} ≈ 1.6 MPa. This effective step increase of Г was used to compare our measurements to the predictions of Eq. 1 in Fig. 2B. In the boundary-lubrication regime, for 100 mm < *x* < 150 mm, Г ≈ 10 J/m^{2} and τ_{p} ≈ 2.66 MPa. Details can be found in the Supplemental Material of (*3*). In general, the measured values of Г and τ_{p} obtained for the sub-Rayleigh regime were used here for supershear ruptures.

### Numerical simulations

Our numerical results were generated by solving elastodynamic equations with the spectral boundary integral method (*32*, *42*). The propagation of a dynamic crack between two half spaces was modeled with a cohesive-type approach to describe the tractions along the weak interface. An explicit time integration was applied. The spectral formulation of the tractions and displacements at the interface results in a periodic setup. In our simulations, we used a replication length of 1.2 m with a discretization of 8192 nodes. The half-space material is linear elastic. To compare with the experiments, we applied a dynamic value of the elastic modulus *E* = 5.65 GPa, Poisson’s ratios *v* = 0.35, and density ρ = 1170 kg/m^{3} and used a plane strain assumption (*C*_{S}/*C*_{L} ≈ 0.48). The interface tractions were governed by a linear slip-weakening cohesive law, τ(δ) = τ_{p}(1 – δ/*d*_{c}) for 0 < δ < *d*_{c}, which imposes a strength that decreases linearly with slip δ from a peak value τ_{p} to zero over a characteristic slip weakening distance *d*_{c}. Reference values applied in the uniform setup were τ_{p} = 1.0 MPa and *d*_{c} = 2 μm, which led to a fracture energy of Г = 1.0 J/m^{2}. Rupture nucleation was triggered via slowly propagating a seed crack of imposed velocity 0.1*C*_{R} [following (*9*)] through the nucleation zone. In the nucleation zone, the value of τ_{p} was gradually reduced to zero (over a length of ≈6 mm). Once the seed crack reached a critical distance *l*_{c} (the Griffith length), rupture acceleration initiated and the ruptures propagated dynamically.

While *v* = 0.35 occurs in serpentinized mantle material (*43*) as, for instance, along the Denali fault (*44*) that hosted a supershear earthquake (*12*, *13*, *45*), *v* = 0.25 is a common value for granite. Additional numerical results were obtained with *v* = 0.25 and plane strain boundary conditions (*C*_{S}/*C*_{L} = 0.577). Comparisons with the theoretical predictions are shown in fig. S1. The results are qualitatively identical to Fig. 2. The main quantitative difference is that the ruptures accelerate over shorter distances to considerably higher speeds, that is, *C*_{f}/*C*_{L} > 0.9. Therefore, the differences in rupture speed for very different prestress levels are already relatively small shortly after transition.

Here, we used spatially nonuniform τ_{0}/τ_{p} profiles. Direct supershear transition was triggered by using high τ_{0}/τ_{p} levels for very short crack lengths. Supershear propagation was then studied as cracks entered into regions with lower values of τ_{0}/τ_{p} (for example, brown curve in Fig. 5). We have also tested the quality of the theoretical prediction for supershear ruptures after transition through the well-known Burridge-Andrews mechanisms (orange curve in Fig. 5). The rupture first propagates in the sub-Rayleigh regime until a radiated shear wave ahead of the crack tip nucleates a secondary crack that then propagates at supershear speeds. A third simulated rupture was initiated by a seed crack that propagated at supershear speeds already during the nucleation procedure. All three tested supershear transition mechanisms lead to supershear ruptures with propagation speeds that are, after brief transient differences, quantitatively well described by our theoretical model. The observed discrepancies within this transient period, we believe, are related to the history dependence, discussed by Huang and Gao (*24*), which was neglected in the derivation of Eq. 1.

### Theoretical framework

In what follows, we briefly describe how LEFM can be used to provide quantitative predictions for supershear rupture propagation following Broberg (*26*). The self-similar solution (*26*, *31*) for a symmetrical bilateral singular crack expanding at a constant speed under uniform remote shear stress τ_{0} provides(4)where *l* is the crack half length and κ(*C*_{f}/*C*_{L}) a known function.

A nonrealistic consequence of the singular description is the vanishing energy flux into the crack tip for *C*_{f} ≠ . It was shown (*31*), however, that when a cohesive zone is introduced, a finite region where these singularities are regularized, the requirement for a positive energy flux to the crack tip is fulfilled for any *C*_{f} > *C*_{S}. The energy flux is given by(5)Here, μ is the shear modulus, τ_{p} (the peak shear strength) is defined in Fig. 1, *x*_{c} is the cohesive zone size, and (*C*_{f}/*C*_{L}) is a known function containing information about the shape of the stress distribution within the cohesive zone. Furthermore, as small-scale yielding is assumed (that is, *l* ≫ *x*_{c} and τ_{p} ≫ τ_{0}), the shear stress of these cohesive zone models recovers the singular form at large distances from the crack tip (*31*)(6)where κ_{c}(*g*) is another known function. The complete form of Eqs. 4, 5, and 6 are given, respectively, by Eqs. 52, 66, and 64 in (*26*), with slight changes in the notations: τ_{p}, *x*_{c}, τ_{0}, and *g* denote τ_{D}, *d*, , and γ.

*G* can be expressed in terms of τ_{0} and *l*, by comparing the two independent forms for τ (Eqs. 4 and 6). This intermediate asymptotic matching was also used by Huang and Gao (*24*) and Antipov *et al*. (*28*). In (*26*), *G* was determined as a function of *C*_{f}/*C*_{L}, *l*, τ_{0}, and *x*_{c} [Eq. 68 in (*26*)] and given by(7)where *B*(*C*_{f}/*C*_{L}) and Г_{D}(*g*) are known functions where the shape of the stress distribution within the cohesive zone is contained within Г_{D}(*g*). Here, motivated by recent experimental results (*2*), which showed that τ_{p} is roughly independent of *C*_{f} (in contrast to the *C*_{f} dependence of *x*_{c}), we rederived *G* in terms of *C*_{f}/*C*_{L}, *l*, τ_{0}, and τ_{p}(8)(*C*_{f}/*C*_{L}) and _{D}(*g*) are given by(9)(10)where *I*_{0}, *N*, *f*_{1}, and ω_{D} are given explicitly in (*26*) and can be calculated once *C*_{f}/*C*_{L} and *C*_{S}/*C*_{L} are specified. Note that the cohesive zone properties are only contained within _{D}(*g*). *D*(ω) defines the shape of the stress distribution within the cohesive zone; for ξ = *x* – *x*_{tip} < 0, the shear stress gradually decreases according to a prescribed shear stress profile, τ(ξ) = τ_{p}*D*(ξ/*x*_{c}). We chose to describe the experimental results. Our numerical simulations use a more simply implemented linear slip-weakening cohesive law τ(δ) (see “Numerical simulations” section). For simplicity, we modeled simulated ruptures by a linear spatial stress distribution within the cohesive zone through *D*(ξ/*x*_{c}) = 1 + ξ/*x*_{c} for –1 < ξ/*x*_{c} < 0. Although some deviations between the models exist [for example, linear spatial profiles of *D*(ξ) result in nonlinear slip laws τ(δ), see Fig. 16 in (*15*)], we explicitly verified that these deviations have negligible effects on the resulting equation of motion.

Figure 6 shows *g*(*C*_{f}/*C*_{L}) and the product of (*C*_{f}/*C*_{L}) and _{D}(*g*) [_{D}(*g*) was calculated for a linear cohesive zone profile], so that the predictions of the supershear equation of motion (Eq. 1 in the main text) can be easily reproduced. At *C*_{f} = , the singularity exponent *g* = 1/2 and Г_{D}(1/2) = _{D}(1/2) = 1. Therefore, Eqs. 7 and 8 coincide, and *G*() is independent of the cohesive zone characteristics. Using these observations and explicitly verifying that *C*_{f}(*l*_{m}) ≈ , the approximated form for a critical length for supershear propagation, *l*_{m}, given by Eq. 2 in the main text, is easily obtained.

A crucial difference between Eqs. 4 and 6 exists. While Eq. 6 adapts rather quickly to spatial changes of interface properties, Eq. 4 (and therefore Eq. 1) does not account for nonuniform τ_{0} due to the underlying assumption of the solution. In what follows, we consider supershear propagation along interfaces with nonuniform spatial profiles of τ_{0}(*x*).

We assume for supershear cracks the condition of small-scale yielding, τ ~ *K*/(*x* – *l*)^{g}, where *K* is the stress intensity factor. The value of *K* explicitly incorporates the information about the τ_{0}(*x*) profile. The fundamental solution for semi-infinite supershear cracks propagating at a constant speed and subjected to a pair of suddenly applied concentrated shear forces is given by Eq. 30 in (*27*)(11)For a general τ_{0}(*x*) profile, *K* = (*C*_{f}/*C*_{L})*K*_{s} is obtained by integrating the fundamental solution following Freund’s superposition method (*30*), which results in (12)Note the analogy with semi-infinite sub-Rayleigh cracks, where *K*_{s}(*l*) ~ .

As far as we know, an explicit solution analogous to Eq. 12 (which was derived for semi-infinite cracks) for propagating supershear bilateral cracks in our simulations has not been formulated. Therefore, we suggest that the sub-Rayleigh static stress intensity factor *K*_{s}(*l*) ~ can be extended to supershear propagation by(13)The normalization factor _{s}(*g*) = ensures that the generalized form of the shear stresses given by(14)recovers Eq. 4 for spatially uniform τ_{0}(*x*) profiles.

It is worth emphasizing again that Eqs. 11 and 12 are derived for cracks propagating at a constant speed. Previous work (*24*) has shown that a rigorous solution for accelerating cracks cannot be constructed by the superposition of these constant-velocity solutions, as had been performed previously for sub-Rayleigh propagation (*23*). In this sense, we hypothesized that these equations can be applied to problems with slowly varying crack speed. As the kernel in Eqs. 12 and 14 is singular, we assumed that the main contribution to the integral comes from the region *s* ≈ *l*. Therefore, assuming that *C*_{f} varies slowly, we calculated the integral with a value of *g*(*C*_{f}) that corresponds to the local crack speed at the crack position *l*. Finally, *C*_{f}(*l*) is determined numerically by Eq. 3: For any crack length *l*, *C*_{f} is evaluated so that Eq. 3 is satisfied. Although this is in no way a rigorous derivation and some of the approximations here are uncontrolled, our conjecture is in good agreement with the simulations and provides a useful tool for describing supershear cracks subjected to nonuniform loading.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/4/7/eaat5622/DC1

Fig. S1. Comparison of theoretical predictions of supershear crack velocities with numerical simulations for various shear strength levels for *v* = 0.25 (plain strain).

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank M. Adda-Bedia for fruitful discussions and his help.

**Funding:**We also acknowledge the support of both the US-Israel Binational Science Foundation (grant no. 2016950) and Israel Science Foundation (grant no. 1523/15).

**Author contributions:**D.S.K. performed the numerical work. I.S. performed the experimental and theoretical work. All authors contributed to the analysis and writing of the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The computer code used for simulations, as well as the additional data related to this paper, may be requested from the authors.

- Copyright © 2018 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).