These continuum models typically consist of a seismogenic region along the fault, in which properties exhibit a tendency to weaken with increasing slip rate and hence potentially susceptible to earthquake-nucleating instability, regions in which strengthen with increasing slip rate and where slip may be stably accommodated, and a mechanical representation of far-field tectonic loading. In all of the continuum models, a fundamental concern is how large must the rate-weakening region be to give rise to an earthquake or aseismic transient and what length scale must be resolved by numerical solutions [e.g., Ben-Zion and Rice, 1995]. An elasto-frictional lengthscale determined by early linear stability analysis [Ruina, 1983;Rice and Ruina, 1983] has provided an approximate measure of susceptibility of a rate-weakening region to instability [Rice, 1993;Liu and Rice, 2007;Rubin, 2008]. We take a closer look at this linear stability analysis and its connection with fault models.

Ruina [1983 ] introduced a general rate-and state-dependent formulation for sliding friction, accompanied by a linear stability analysis of a spring-block model, identifying a critical threshold stiffness below which steady sliding is unstable to perturbations. Rice and Ruina [1983] utilized a mapping between perturbations to a spring-block system and periodic perturbations along a frictional interface within a continuum, considering the specific cases of in-plane (mode-II) or outof-plane (mode-III) rupture of a fault plane dividing two elastic half-spaces or blocks of fixed and equal height. In either case, the fault is presumed to slide infinitely along the two directions of its plane and the perturbations are about uniform sliding at steady state; the mapping of Rice and Ruina [1983] translated the critical stiffness of Ruina [1983] to a critical spatial perturbation wavelength above which perturbations may grow unstably. While this analysis provided an early indication of the stability of model faults and the critical wavelength has been a canonical reference for a nucleation length, the basis of this continuum analysis is not representative of the wide variety of model conditions and neglects the non-linear nature of the friction law. Seismic cycle models typically have variations in frictional properties with down-dip or along-strike distance; the fault may be subject to an increasing stress due to remote loading; and the fault typically has a finite domain, with its boundary either satisfying a locked condition, a prescribed slip rate intended to model tectonic forcing, or a free-surface condition. While possibly providing an initial indication of model stability, periodic perturbations along an unbounded fault cannot capture this range of properties, loading, and boundary conditions.

To address this first deficiency, in section 3 of this work, we revisit the problem of linear stability and consider two elementary cases that exhibit the effects of fault finiteness and different loading configurations. First, we consider a finite fault patch sliding at steady state while surrounded by sliding that remains at a fixed rate. Second, we examine a finite fault whose boundaries are locked and the fault loaded by a uniform, constant rate of shear stress. As also discussed by Ciardo and Viesca [2025], we show that the first case can be considered by the same mapping to a spring-block system pursued by Rice and Ruina [1983]. The second case requires special treatment owing to the non-uniformity of slip rate and state, about which linear perturbations are considered. In both cases we identify a dimensionless parameter, above which the fault is linearly unstable. This dimensionless parameter can be interpreted as identifying a critical size of the rate-weakening fault or the rate-weakening patch on a fault. Notably we find that inherent non-uniformity of fault slip leads to a distinct scaling of the critical size with the friction rate-weakening parameter, b -a.

To address the second deficiency of neglecting frictional non-linearity, in section 4 we revisit and refine prior analyses in which solutions for self-sustaining non-linear instability are found [Rubin and Ampuero, 2005 ;Viesca, 2016a,b;Viesca, 2023]. In contrast with the linear stability analysis, these solutions describe quasi-static acceleration towards dynamic rupture taking into account the full non-linearity of rate-and-state friction. A remarkable feature of these solutions is their near-independence of the loading and boundary conditions, or universality. We highlight one exception to this independence and differentiate between so-called “free” and “pinned” solutions: the former qualification applies when faults are sufficiently large that this non-linear instability can develop without encountering any locked, steadily creeping, or terminal boundaries of a fault, and the latter applies to faults that are sufficiently small to constrain the lateral extent of instability development. In either case, we identify a second dimensionless parameter that can, for example, be used to identify a second critical size of a fault patch or fault.

In the subsequent sections, we use the results of linear and non-linear instability development to make a priori predictions of fault model behavior on the basis of the choice of fault frictional properties, rock elastic properties, fault dimensions and loading. In section 5, we highlight the relevance of our results to understand the transition from linear to non-linear stages of instability development. In section 6 we discuss in greater detail several aspect of the results. Foremost, is the implications for quantitatively and precisely predicting the behavior of seismic cycle models. Additionally, we will show extensions of the results of sections 3 and 4, considering the twodimensional elastic response to a single-mode of slip with one-dimensional variations along the fault, to three-dimensional elasticity, involving mixed-mode slip with variations occurring across two dimensions of the fault; and finally, we will discuss consequences of changing the so-called state-evolution law. The techniques are evidently generalizable to more complex cases, such as spatially varying frictional properties and effective normal stress.

We consider a planar fault lying in the x -y plane and sliding such that displacement varies within the x -z plane either as an in-plane or anti-plane state of deformation. This implies that the fault-resolve shear stress can be written as

where δ is the distribution of relative displacement across the fault plane, or simply, slip, τ b is the fault-resolved shear stress in the absence of slip, and L is a linear operator that determines the change in shear stress along the fault due to a distribution of slip. The operator has units of stress per unit distance. In the example of a fault dividing two elastic half-spaces [e.g., Rice, 1968],

where µ ′ = µ/[2(1 -ν)] or µ/2 for in-plane or anti-plane deformation.

For we may write the shear strength as τ s (x, t) = σf (x, t) where σ is the fault-normal stress and f is a rate-and state-dependent friction coefficient that may be written as

Each law is a function of the dimensionless parameter V θ/D c where D c is a characteristic distance.

The three laws share the steady-state point V θ/D c = 1 and a common linearization about this point.

We now look to linearize the pair of non-linear evolution equations for V and θ that result from the combination of (1), (3), and one of (4). We will be concerned with the evolution of departures Ṽ and θ from a regime of steady-state sliding under slip rate V o and state

V o and θ o may not necessarily be uniform in space but are such that

The three state evolution laws hence linearize to

and the friction evolution equation linearizes to

Also requiring that ∂τ /∂t = σ∂f /∂t, implies that, to within linearization,

and that V o and ∂τ b /∂t must satisfy

The latter requirement follows from considering that the unperturbed regime is steady such that both state θ o and slip rate V o do not evolve.

Robert C. Viesca In this section we analyze linear stability of faults under different geometries and loading conditions. We begin by considering faults that slide at a uniform, steady rate. In section 3.1, we show that linear stability has some common features, regardless of whether the fault is unbounded or if we focus on a finite-sized patch of the fault (e.g., a finite asperity loaded by surrounding creep), and we highlight the distinguishing features of these two cases. Finally, in section 3.2, we examine a distinct problem in which a finite fault is loaded at a constant rate of shear and for which the steady-state solution for slip rate and state is non-uniform. We show how this problem leads to scalings that are qualitatively different than cases of steady-state uniform sliding.

Revisiting the linear stability analysis of Rice and Ruina [1983], except here for a generic form of L, we consider a fault for which ∂τ b /∂t = 0 and we perturb about uniform values of slip rate and state

We pass to non-dimensional variables in the manner

such that the linearized governing equations are

We look for perturbations in the form

where s, V and Θ are to be determined. The governing equations reduce to

the second of which allows us to eliminate Θ from the first since Θ = -V/(1 + s), leaving us with

which has the form of an eigenvalue problem where we identify the terms in parentheses on the right hand side as the eigenvalue -k, such that

and identify the occurrence of a Hopf bifurcation, when s takes on a complex conjugate pair of values that cross the imaginary axis when the terms in parentheses are zero. This defines a critical value for k

which can be considered as a dimensionless critical stiffness for the system, below which the real part of ( 19) is positive and perturbations are unstable. The corresponding conjugate values of s are

As will be illustrated in the following two subsections, the eigenvalue problem ( 17) is pivotal and its solution enables the translation of the critical value k c to a critical spatial dimension along the fault. That translation is determined by the particular form of the elastic stress-transfer operator L, which reflects the fault geometry, mode of sliding, and presence of free surfaces. Ciardo and Viesca [2025] recognized this in the context of a study of spring-block slider stability and solved this eigenvalue problem analytically and numerically for a number of elementary configurations, including those presented here. (The reader should be aware of a typographical error of an omitted minus sign in their equation (1.8): there is no consequence of this omission except that L therein is minus that considered here.)

We use the above to rederive the results for the specific case considered by Ruina and Rice [1983] of an unbounded fault, for which the dimensionless expression of the operator L is

where µ ′ = µ/2 for anti-plane (mode-III) sliding and µ ′ = µ/[2(1 -ν)] for in-plane (mode-II) slip. We non-dimensionalize distance as

We identify the left hand side of (23) as the negative of the Hilbert transform

x-s ds and can write for brevity

The Hilbert transform and the spatial derivative have the commutative property (H(f )) ′ = H(f ′ ), provided f (x) decays sufficiently fast as x → ±∞. In addition H has the property of being its own inverse, to within a sign change: H[H(f )] = -f . Keeping these two points in mind, we can reduce (24) to

Linear and nonlinear stability of rate-and-state faults Page 7 of 30 20 Dec 2025 which has the two linearly independent solutions cos(kx) and sin(kx) and we find that k corresponds to a wavenumber, or a perturbation wavelength λ with k = 2π/λ. The critical value k c in this case hence implies a corresponding critical wavelength λ c , above which perturbations are unstable,

Reverting to dimensional form for the critical wavelength via λ c ⇒ λ c σb/(µ ′ D c ) and using the lengthscale

introduced by Rubin and Ampuero [2005], λ c can be expressed as

or more simply,

The results of the preceding subsection establishing a critical wavelength for an unbounded fault are known and were established by Rice and Ruina [1983]. We now look to examine a novel case, in which a finite region |x| < L may freely slide but is surrounded by a region |x| ≥ L that slides at a rate V pl . This configuration may be thought of as a creeping region that encloses an asperity of size 2L, for which are interested in the response to perturbations about steady-state sliding (11) within |x| < L. For this case the dimensionless form of L is

and we may subsequently non-dimensionalize distance as x/L ⇒ x such that the eigenvalue problem 17 becomes 1 2π

in which L remains dimensional and the ratio L/L b can be considered as a dimensionless fault half-length. We recognize (31) as the eigenvalue problem examined by Uenishi and Rice [2003], for which the combination kL/L b corresponds to the eigenvalue λ in that work. This eigenvalue has a discrete spectrum with a countably infinite number of values, the lowest three of which are (kL/L b ) 1 = 0.57888694…, (kL/L b ) 2 = 1.3773774…, and (kL/L b ) 3 = 2.1584005… (see Appendix A for numerical solution). To determine the critical value of L above which the fault is linearly unstable, L c , the salient eigenvalue is the smallest one. Since linear instability is incipient when k = k c , L c can be expressed as

The eigenfunction corresponding to the leading eigenvalue can be simply and accurately approximated as [Viesca and Dublanchet, 2019]

Linear and nonlinear stability of rate-and-state faults Page 8 of 30 20 Dec 2025

A non-uniform sliding rate is an inherent feature of faults featuring “strongly” locked ends and uniform shear stress loading, as would be the case for a fault of finite extent and examined here, or a fault patch which experiences non-uniform loading, such as loading by steady creep on one end and “strong” or “weak” locking at the other. The former, “weak” locking may be due to a fault boundary and the latter due to a temporary, post-seismic seizing of fault regions that slipped seismically and may eventually slide again.

Here we consider a fault undergoing in-plane or anti-plane deformation and locked at it ends x = ±L such that slip only occurs on |x| < L, while being uniformly loaded at a constant rate, which we denote as ∂τ b /∂t = τb . In this case, we will consider perturbations about a steady-state condition, which for this case, will not have a uniform slip rate, but one that must satisfy (10), or

Linear and nonlinear stability of rate-and-state faults Page 9 of 30 20 Dec 2025

The solution to which is

where we identify for this problem the characteristic slip rate V pl as

and the corresponding distribution of the state variable for steady-state sliding is

We non-dimensionalize ( 7) and ( 9) in the manner of 12, and rearrange to arrive to the pair of coupled, linear evolution equations for Ṽ and θ

We again look for perturbations of the form ( 14), non-dimensionalizing distance as x/L ⇒ x, introducing (30), such that (37) becomes

where we also require that V(±1) = 0 to meet the locked boundary condition. We recognize (38) as an eigenvalue problem with eigenvalue s, eigenfunctions V(x) and Θ, and parameters a/b and L/L b . Unlike (15), we see that eliminating Θ in (38) does not simplify the problem and we instead to solve the system (38) numerically, using a procedure detailed in Appendix B.

We find that stability is lost at a critical value of the fault length, L c /L b , the value of which value is a function of a/b. We find this critical length by fixing the value of a/b, increasing the parameter L/L b from 0, and solving the spectrum of eigenvalues for each value L/L b . We observe that at the critical value, stability is lost as a pair of complex conjugate eigenvalues s c cross the imaginary axis. For L < L c , the fault is linearly stable to perturbations. As L/L b is increased beyond L c /L b a sequence of Hopf bifurcations occur at discrete values of L/L b : additional unstable modes are presented for longer faults. In Figure 1a,b we show the numerical solution for the dependence of L c /L b on a/b, found by a bisection method, as well as the approximation

which has < 1% relative error for 0 < a/b < 0.998. While this has the appearance as an asymptotic expansion in the limit 1 -a/b → 1, we have not conducted an asymptotic analysis of eigenvalue problem at the Hopf-bifurcation point. However, the numerical solutions do suggest the leadingorder asymptotic term O((1 -a/b) -3/2 ) in the limit 1 -a/b → 0. The remaining terms in (39) are chosen to provide an accurate, low-parameter fit of the numerical results over the entire range of a/b examined. We emphasize here that the apparent asymptotic scaling of L c /L b changes when Dmitry I. Garagash 2 shifting from a uniform to non-uniform baseline rate of sliding: the critical fault dimension

) as a/b → 1 here, representing an abrupt departure from the uniform (1 -a/b) -1 scaling seen in section 3.1 for uniform sliding, ( 28) and (32).

In Figure 1c,d we show the value of imaginary part of the eigenvalue s at the Hopf bifurcation point, as well as its approximation

which has < 2% relative error for 0 < a/b < 0.1 and < 1% relative error for 0.1 < a/b < 0.999. The approximation (40) makes an explicit attempt at matching an apparent leading-order b/a asymptotic behavior of Im(s c ) as a/b → 0, an asymptotic scaling also seen in ( 21). However, the following terms in (40) are chosen to optimize fit over the entire range of a/b explored in numerical solutions and are not to be taken to be an asymptotic analysis.

There is a temptation to use the critical lengths L c calculated in sections 3.1 and 3.2 as a nucleation length, but this is a fine point: one should be clear as to what one is considering the nucleation of. Here we define the nucleation length as the initial size, or nucleus, of a dynamic, earthquake-generating rupture. By this definition, L c cannot be qualified as a nucleation length.

The linear stability analysis only indicates whether spontaneous acceleration beyond a specific state of steady-state sliding can initiate provided the fault dimensions are above a critical size. This initially unstable acceleration will not necessarily progress to seismic slip speeds. The stability analysis assumes slip rate and state are within a linear regime, in which perturbations from steadystate are of a few percent of steady-state sliding. However, seismic slip speeds of O(m/s) are ten orders of magnitude above background plate rates of O(cm/yr). Thus in order to consider the earliest phase of a dynamic regime, one must be able to examine the quasi-static acceleration of the fault well into a non-linear regime.

A series of works [Rubin and Ampuero, 2005;Ampuero and Rubin, 2008;Viesca, 2016a,b] has examined this non-linear phase of acceleration to determine nucleation lengths. Viesca [2016aViesca [ ,b, 2023] ] demonstrated that an analysis commonly used for non-linear, time-dependent problems can be applied to unstable acceleration of faults with rate-and-state friction with evolution laws ranging from aging to slip. Specifically, that work sought solutions in which slip rate diverges within a finite time t f as

in which the spatial distribution W(x) is solved for and found to be a compact distribution within |x| ≤ L f . For in-plane or anti-plane sliding between two elastic half-spaces on a fault with uniform friction properties in the rate-weakening range 0 < a/b < 1, Viesca [2016a,b] showed that the spatial distribution of diverging slip W and the distance L f must satisfy

where we have nondimensionalzied distance as x/L f ⇒ x. This problem is a free-boundary problem: the size of the region over which slip rate diverges is given by L f , an unknown to be solved for. Equation ( 42) is not sufficient to determine L f and an additional constraint is required. This constraint is the requirement that there be no stress-rate singularity just ahead of the tips of the compact region at x = ±L f . The physical origin of this requirement is the presumption that unstable acceleration of slip is occurring on a region contained within a larger fault whose finite strength cannot support such a singularity. This requirement amounted to imposing the condition that W ′ (±1) = 0. Viesca [2016a,b] solved (42) for W and the dimensionless size L f /L b , as they both depend on the parameter a/b. The solution for L f may be approximated by the following expression

where the first case is precise to as many significant figures shown and the latter case has a relative error less than 1%. Remarkably, these solutions retrieved two lengthscales proposed earlier by

The first lengthscale was found by Rubin and Ampuero [2005] and appears as a special case of the general analysis pursued by Viesca [2016a,b] and the second by the former authors approximating instability development via a fracture analogy. Viesca [2016a,b] provided a consistent analysis of nucleation over the entire range of a/b and showed that the second lengthscale is, in fact, the appropriate asymptotic limit for L f as a/b → 1, deriving the equivalent fracture problem L f must satisfy. L f can be considered as the expected nucleation length on an aging-law fault provided the nucleation patch is smaller than the fault itself or than the region of the fault not imposed to slip at a fixed rate.

The finiteness of faults may play a role in determining the nucleation length. The results of section 4.1 presumed that a singularity in stress rate was inadmissible and that the rate-weakening portion of the fault extended beyond the nucleation zone. However, for sufficiently small faults, nucleation and dynamic rupture may be contained within the fault if the surroundings are able to support a singularity at fault boundaries: such cases may occur for a rate-weakening fault that is locked at its edges, or as a finite rate-weakening patch surrounded by steady creep. In the former case, the fault terminates and is surrounded by rock with finite toughness capable of sustaining a singularity. In the latter, the rate-weakening asperity lies within a steadily creeping and the surrounding region may be strongly rate-strengthening (i.e., in the limit a/b → ∞) or the fault exterior to the asperity is simply restricted to slide at a constant rate by model design. In any case, we may search for so-called “pinned” non-linear solutions of the form (41) again solving (42) (with L f therein replaced with L p and x implicitly scaled by L p ) and removing the requirement imposed by Viesca [2016a,b] that W(±1) = 0. Here L p /L b is no longer a solution variable determined by a boundary condition but now can be considered as a free parameter. Restricting ourselves to solutions for which W (and hence V ) is positive within a region |x| ≤ 1, we find that the pinned solutions occupy a range L min p < L p < L f , where L min p is the lowest possible value for which a solution can be found numerically. Since this is the lower bound is the most relevant value in future discussions, we will drop the “min” superscript and simply dub this minimum “pinned” nucleation length L p , letting it be understood that there exist a family of pinned solutions starting from this minimum value and up to L f . L p may be approximated, to less than 0.1% relative error over the rate-weakening range 0 < a/b < 0.999, by

When the direct effect is neglected (a/b = 0), (42) reduces to the Uenishi-Rice eigenvalue problem in which the physically admissible solution is that of the lowest eigenmode: i.e., L p /L b = 0.57888694… in that case.

2L p represents the size of the smallest asperity within a creeping region or smallest fault that can sustain a non-linear instability towards seismic slip speeds, if the the fault has uniform rateweakening properties. This may be considered a nucleation length in that it marks the starting size for a dynamic rupture. However, under the conditions described above-locked edges, or ratestrengthening or creeping surroundings-the dynamic rupture would be expected to be contained within the nucleation area, unless another weakening mechanism is enabled by the high slip and stress rates of the dynamic rupture to allow for dynamic rupture to penetrate the surrounding barriers.

4.3 “Free” and “pinned” nucleation from aging to slip laws

The nucleation lengths for pinned (43) and unpinned (46) faults were found assuming that the fault followed the aging state evolution law (4a), but the analysis and results are easily extended to friction laws ranging between aging and slip laws. Using the intermediate law (4c), whose parameter ϵ allows a continuous transition between aging (ϵ = 1) and slip (ϵ → 1) laws, Viesca [2023] Viesca showed that non-linear blow-up solutions of the form (41) continue to exist for the intermediate law. This is seen by repeating the derivation of Viesca [2016a,b] leading to (42), except with introducing a factor of ϵ to the right hand side of (41) at the start. Doing so, Viesca [2023] found that the equation governing W is exactly that for the aging law, (42), except the lengthscale L f (or L p ) is reduced by a factor ϵ. This means that the expressions for nucleations lengths L f and L p also apply for state evolution laws ranging from aging to slip, provided the right hand side of each equation is multiplied by ϵ.

An immediately apparent consequence of this result is that nucleation lengths vanish as ϵ → 0. This is born out in numerical solutions of quasi-static instability on slip-law faults [Ampuero and Rubin, 2008;Rubin and Ampuero, 2009], in which the slip rate simultaneously accelerates and contracts towards a point in space in finite time. The slip-pulse nature was shown by Viesca [2023] to emerge as the consequence of the partial loss of stability of the blow-up solution (41) such that migration occurs, but the point-like nature of the blow-up is preserved. In fully dynamic models of the rupture process, quasi-static acceleration is ultimately limited by inertia, which prohibits continued pointlike acceleration and induces outward propagation of a dynamic rupture. How close the quasi-static acceleration approaches this point-like nucleation ultimately depends on the strength of the perturbation on the fault that initiates the sliding instability. This can be understood by considering a non-linear stability analysis of a spring-block slider model following the intermediate evolution law [Ciardo and Viesca, 2025], in which the value of the critical stiffness delimiting stable and unstable trajectories was found to be dependent on the perturbation size for all evolution laws and that dependency was most sensitive in the slip-law limit (ϵ → 0). Taking the critical stiffness to be inversely proportional to a fault lengthscale (e.g., perturbation wavelength λ or fault half-length L), this result implies a perturbation-dependent lengthscale for nucleation. In light of this, the nucleation lengths ϵL p or ϵL f must be considered as lower bounds on the pinned or unpinned non-linear nucleation lengths for a given ϵ.

Linear stability analysis is loading-specific. As shown in section 3, the critical length L c is strongly determined by loading condition. Foremost, the loading condition is reflected in the form of the underlying distribution of slip rate and state, about which perturbations are considered. For instance, the asymptotic scaling of L c /L b in the limit a/b → 1 changes with loading. Denoting ϵ = 1 -a/b, that scaling changes from O(ϵ -1 ) for a finite rate-weakening patch driven by boundary creep to O(ϵ -3/2 ) for uniform rate of shear stress loading on a finite fault. The only commonality of the stability results for the two loading configurations occurs when the direct effect is neglected (a/b = 0), in which case the critical length L c /L b for both configurations is given by the first Uenishi and Rice [2003] eigenvalue (0.57888694…). An additional signature of the initial and loading conditions is the characteristic slip velocity (V pl ) and timescale (D c /V pl ) for evolution of the slip rate following perturbation. For example V pl in the case of constant rate of shear stress, is given by ( 36), whereas for the case of loading by steady creep, V pl reflects the the rate of slip in the boundary regions.

In contrast, in a non-linear regime, finite-time divergence of slip rate in the form ( 41) is independent of loading conditions. The characteristic timescale and slip rate contain no information on the process that triggered instability. Rather, the relevant timescale is given by the time leading up to a finite time instability t f (t) = t o -t and the relevant slip rate scale is D c /t f (t), which diverges as t → t o . The “free” and “pinned” nucleation lengths L f and L p , as well as the slip rate distribution W also contain no information of initial conditions. This independence implies that late-stage acceleration of fault sliding exhibits a type of universality: there exists a pathway for slip rate to continue accelerate unstably regardless of what triggered the instability How does an instability on the fault transition from an unstable acceleration determined by loading conditions to one independent of it? In other words, what ensures that this universal pathway is taken? For example, consider the initial growth of perturbations from an initial steadystate: V = V o (x)+ Ṽ (x, t) with Ṽ (x, t) = V(x) exp(st). Let’s follow the fastest growing mode V 1 (x) with corresponding s 1 . We recall that this mode shape and exponential growth rate presumed small departures such that the governing equation could be reasonably linearized. We also recall that the form of V o (x) and V 1 (x) and the growth rate and s 1 are determined by the manner of loading, the fault size L/L b , and ratio a/b. For linear instability, these three factors are such that Re(s 1 ) > 0. However, as acceleration continues, non-linear effects become important and the distribution of slip rate and acceleration will depart from this modal exponential growth. Why might slip rate transition to diverge as V (x, t) = W(x)D c /t f ? That slip rate would favor approaching this form follows from a dynamical-systems analysis of its asymptotic stability [Viesca, 2016a,b]. This analysis indicated that under most cases, solutions of the form (V (x, t) = W(x)D c /t f ) are attractive and asymptotically stable, meaning that they are expected to be approached in the limit t → t o . This asymptotic stability is readily observable under the aging-law for a/b < 0.74. In limited cases, for which a/b is in the vicinity of unity, these solutions are not asymptotically stable, but may be considered to be Lyapunov stable and attractive nonetheless; this is most noticeable in numerical solutions for aging-law with a/b > 0.9, in which the presence of several unstable modes leads to alternation between slip rate approaching the solution (41) and large, chaotic deviations. This analysis reconciled the apparent dichotomy in behavior of earlier numerical solutions that were observed to exhibit asymptotic stability for a/b < 0.3781… and a lack of it for a/b near 1 [Rubin and Ampuero, 2005].

However, the following question remains to be answered: under what conditions will a linear instability progress beyond its linear phase and continue to accelerate in a non-linear regime until a dynamic rupture is nucleated? Posed alternatively, does a system that is linearly unstable necessarily continue to accelerate? This question might be most easily discussed considering fault length L as a variable, with all other parameters held fixed. In this case, we can simply compare L to the critical linear stability length L c and “free” or pinned nucleation lengths L f and L p .

As a starting point, we compare L c among the two cases of linear stability involving a finite fault dimension L: (i) a fault uniformly creeping at steady-state and loaded by a constant slip rate at the boundaries; (ii) a finite fault loaded with a constant rate of shear stress. Using here a superscript to denote the L c corresponding to each case, the relative values of L c follow a clear order,

c , for all 0 < a/b < 1 (Figure 4a).

We now compare the nucleation lengths, which have their own clear order, L p < L f , with L c for the first two of these three cases. In Figure 4b we show that L

c , this ordering would follow for L (i) c as well. We now return to the problem of considering a variable L with all other parameters fixed. For an L < (L

), the fault is linearly stable to perturbations, meaning any small perturbations to the initial steady state. It is then also necessarily the case that L < L p , the smallest possible nucleation length. Thus at this fault length L, there exists no solution for finite-time divergence of slip rate accounting for the full nonlinearity of the evolution equations. This implies that the fault is also nonlinearly stable: there is no finite (large) perturbation that will provoke the system to unstably accelerate. In other words, for L < (L (46). While L c is independent of the chose state evolution law (aging, slip or an intermediate law), L p is dependent on that choice; here L p is drawn for aging-law state evolution. its loading. There will be no spontaneous, transient acceleration, other than that directly imposed by an applied perturbation, and such an acceleration will ultimately decay back towards steady state. For fault lengths L in the window (L

c ) < L < L p , the fault can be considered to be linearly unstable but non-linearly stable. This means that small perturbations to steady state can grow with temporary excursions into the non-linear regime. However, because the fault is not sufficiently large to admit solutions that diverge in this non-linear regime, these excursions cannot transition into spontaneous acceleration and diverging slip rate. A likely outcome is that the fault will ultimately decelerate back towards steady-state, before another perturbation repeats this process.

How do the results of sections 3-5 inform models of the seismic cycle coupling elastic deformation of the host rock with rate-and state-dependent friction on the fault? Here we focus our discussion on the specific case of a rate-weakening fault patch embedded within an otherwise stably creeping fault, though the discussion would be similar for the other two cases for which we considered linear stability.

In Figure 5 32) and ( 46), respectively. The expression of L p is assuming aging-law state evolution; a prefactor of ϵ must be introduced if assuming the intermediate evolution law (4c). When L < L c , steady sliding of the rate-weakening patch at a uniform rate V p is both linearly and non-linearly stable: any perturbations to the fault will ultimately decay and the fault will return to steady creep.

For L in the range L < L c < L p we may expect the spontaneous emergence of aseismic transients: the fault is linearly unstable meaning small perturbations will grow exponentially quickly, but the fault is also non-linearly stable meaning that growth cannot continue. On the basis of results for a spring-block slider [Gu et al., 1983], we know that linear stability is lost as a supercritical Hopf bifurcation [e.g., Strogatz, 1994]. This implies that, as L is increased beyond the phase boundary L = L c , a stable limit cycle emerges in which state and slip rate may periodically vary about steady-state, and the amplitude of this stable limit cycle grows with increasing L. When the fault size L is finitely larger than L c , the linear description of this bifurcation loses relevance as the growth of slip rate following a small perturbation enters the non-linear regime. Exponential growth of perturbations does not continue indefinitely, and provided L remains below L p , there does not exist an attractive manner for slip to continue to accelerate. This means that any transient, finite acceleration induced by being linearly unstable is ultimately limited as the fault is non-linearly stable, and the fault will decelerate towards steady-state and the process repeats.

When L > L p , the fault can spontaneously accelerate towards dynamic rupture. In this case, the fault is linearly unstable, small perturbations spontaneously grow and because the fault is now also nonlinearly unstable, those perturbations may continue to grow in the non-linear regime. This is the fault is sufficiently large to admit attractive pathway for the slip rate to diverge in this nonlinear regime, in the form of the solution (41). For L in the range L p ≤ L ≤ L f , that acceleration towards inertial sliding and dynamic rupture occurs over the entire fault length L. This means that sliding reaches seismic slip velocities occurs simultaneously over the entire fault and there is no outward propagation of this dynamic rupture, since in this model, sliding beyond x = ±L occurs at a fixed rate V pl ; the rate-weakening patch slips will slip seismically, lock, and will subsequently be reloaded by the boundary creep to repeat the process. For L above L f (not shown on the phase diagram of Figure 5, but asymptotically approaching L p in the limit a/b → 1, as shown in Figure 3), the dynamic rupture can propagate before arresting at the boundaries x = ±L.

Both the linear stability and non-linear nucleation analyses here can be adapted to account for heterogeneity of frictional properties and effective normal stress. Seismic cycle models typically include variations in the coefficients a and b and also normal stress. Variations of the former typically takes the shape of a transition from rate-weakening behavior in the to rate-strengthening behavior down dip. Such variations were originally in early models of the seismic cycle [Tse and Rice, 1986;Rice, 1993;Lapusta et al., 2000], justified on the basis of hydrothermal dependence of frictional properties [Blanpied et al., 1991[Blanpied et al., , 1995]]. One may repeat the linear stability analyses of section 3 accounting for the heterogeneity of a and b, solving for the eigenfunctions V and Θ and eigenvalues s numerically. We defer this to subsequent work. Regarding unstable acceleration in he non-linear regime, Ray and Viesca [2017] considered the problem of finding solutions of the form (41) and the nucleation length L f considering heterogeneity of frictional properties a, b, and D c . They demonstrated that the equation governing the slip-rate distribution W and L f is equivalent to a classical fracture problem of a slip-weakening crack in equilibrium with an applied load. Subsequent work [Ray and Viesca, 2019], demonstrated the possibility of homogenization in which property variations occur over wavelengths disparate from elasto-frictional lengthscales (such as L b appearing in this work). While these works focused on the so-called “free-boundary” nucleation L f , the governing equations remain the same for pinned nucleation L p lengths under heterogeneity. Solving for L p can be done as pursued in 4.2: simply by removing the condition that W ′ vanish at its endpoints.

6.2 Linear stability and nucleation lengths in 3D (planar faults)

Up to now, in sections 3-5 we considered a single mode of slip, in-plane (mode-II) or anti-plane (mode-II), in which variations occurred along a single dimension and boundaries were represented by points. How do results change if a fault is undergoing mixed-mode slip, in which variations can occur in both dimensions of a two-dimensional fault plane and the boundaries are curves? Here we are interested in determining critical dimensions under such conditions. The simplest geometry to consider under mixed-mode conditions is a circular region of radius R. We will present a solution for the critical radius, R c , of a rate-weakening patch at which point its linear stability is lost. For brevity, we will do so for the specific loading case analogous to the single-mode problem considered in (3.1.2). In addition we will present solutions for the free and pinned nucleation radii R f and R p . The former corresponds to the size of a circular patch that can nucleate dynamic rupture within a larger, rate-weakening fault, which itself is not not necessarily circular. The latter corresponds to the smallest circular rate-weakening region that may nucleate dynamic rupture while confined at its edges, because the region exterior to the rate-weakening patch is either locked or undergoing an imposed slip rate.

For a circular region with an axisymmetric distribution of uni-directional slip within r < R, the traction on the fault plane is only parallel to that slip provided the Poisson ratio ν = 0. In this case, the operator L in (1) is given by [Salomon andDundurs, 1971, 1977;Bhattacharya and Viesca, 2019]

where µ ′ = µ/2, E and K are the complete elliptic integrals of the second and first kind, and k(u) = 2 √ u/(1 + u) is the elliptic modulus. While ν = 0 may appear to be a peculiar case to be concerned with, it can be considered as the leading-order approximation in a perturbation series solution to a problem in which ν ̸ = 0 and for which ν is the perturbation parameter. For shear cracks, this was first exhibited for the problem of a singular crack in equilibrium with a uniform loading and having a uniform energy release rate (matching a uniform fracture energy) along its circumference, analytically treated by Gao [1988], following Rice [1985]. There the leading-order (ν = 0) solution to this mixed-mode problem is a circular crack. The next-order correction to the boundary shape and slip distribution is of the order ν. This is also seen for non-singular ruptures. An example of problem in this vein is a fault with constant Coulomb friction, the rupture of which is driven by the migration of fluid along the fault away from a point source injecting at constant rate [Sáez et al. 2023;Viesca, 2025]. Closed-form expressions for the leading-order (ν = 0) circularcrack solution to this problem were provided as a perturbation expansion by Viesca [2025] and comparisons of these with numerical solutions to the problem in which ν > 0 [Sáez et al., 2023] show that the perturbations to the circular shape and peak slip are, to leading order, linear in ν. Remarkably while some terms require correction, the authors found, to within the accuracy of the For R < R c , the fault is linearly and non-linear stable, and steady creep is an unconditionally stable solution. For R c < R < R p the fault is linearly unstable, but non-linearly stable: the fault exhibits spontaneous emergence of aseismic transients, but acceleration is ultimately limited. For R > R p the fault is both linearly and non-linearly unstable and acceleration can spontaneously progress towards seismic slip speeds and the nucleation of dynamic rupture.

Robert C. Viesca 1, 3 Dmitry I. Garagash 2

Linear and nonlinear stability of rate-and-state faults Page 20 of 30 20 Dec 2025 numerical solutions, that the fracture area is preserved and independent of ν. Thus the solution for which ν = 0 is an informative and often adequate basis to consider the three-dimensional effects on rupture.

We now look to extend linear stability results for a single mode of slip to the mixed-mode problem, taking as an example the mixed-mode equivalent to the problem of section 3.1.2: a circular rate-weakening fault patch or asperity that has radius R and is surrounded by steady creep at the rate V pl (Figure 6a, left). The steady-state solution is uniform, steady sliding at the rate V pl within r < R. As part of the class of problems of section 3.1, the critical radius R c is determined by first solving the eigenvalue problem (17) with L therein as given by (47). Nondimensionalizing and solving the eigenvalue numerically, we find that the leading order eigenvalue is (kR/L b ) 1 = 1.003059165… . Linear stability is lost when the stiffness k = k c and thus the critical radius

The eigenfunction corresponding to the leading eigenvalue is accurately approximated, to within an arbitrary factor, by

For comparison, we also look to find the size of free and pinned nucleation patches that correspond to slip rate diverging in the manner of (41). For nucleation occurring on a uniformly rate-weakening patch of the fault away from fault boundaries (Figure 6a, bottom), the nucleation is unconstrained and its shape and size must be determined. The condition ν = 0 and the translation and rotational symmetry of the frictional properties implies that the slip-rate distribution of (41), W, will likewise be axisymmetric and translationally invariant within the rate-weakening region of the fault. Thus nucleation occurs over a circular region the radius of which is a freeboundary problem to be determined. The slip-rate distribution W and that radius R f must satisfy [Viesca, 2016a,b]

in which distances have been normalized as r/R f ⇒ r, L b = µD c /σb, and the boundary condition for non-singularity of the stress rate is W ′ (1) = 0. This last conditions determines the free-boundary problem for R f . Solving this problem numerically, we find that R f may be approximated as

) where the first expression in ( 51) is precise to as many significant figures as shown and the expression is accurate to within 1% relative error up to a/b = 0.99. The asymptotic behavior of R f as a/b → 1 can be derived by the considering equivalence of the problem (51) to that of a slip-weakening fracture in equilibrium with a uniform background stress [Viesca, 2016a,b]. Doing so yields that, as a

This nucleation length was also expressed by Chen and Lapusta [2009], citing a private communication with Rubin, presumably following the line of argument of Rubin and Ampuero [2005] based on identifying an apparent similarity between numerical solutions and a small-scale-yielding fracture.

Robert C. Viesca In contrast to the free-boundary problem above, the rate-weakening patch may be circular and small enough such that nucleation, if possible, occurs within a finite, pinned radius (Figure 6)a, right). We are interested to find smallest radius, R p , that can support the spontaneous nucleation of dynamic rupture, by admitting solutions of the form (41) satisfying (50). We solve (50) numerically, without the boundary condition W ′ (1) = 0 and instead imposing a value for R f /L b as a problem parameter. The minimum value of that parameter that admits a numerical solution is found by bi-section and is dubbed R p /L b . The dependence of R p on a/b is well approximated by the expression

which is to within 0.1% relative error for 0 < a/b < 0.999

In Figure 6 we draw comparisons between the critical radius R c from linear stability and the nucleation radii R f and R p . In Figure 6b we highlight that R p < R f for 0 < a/b < 1. In Figures 6c,d, we show the phase diagram of the behavior of a circular rate-weakening patch of radius R driven by creep of the surrounding fault. The boundaries of the phase diagram are given by R = R c and R = R p . Since R c < R p , there are three regions in the phase space, as in the equivalent problem for a single mode of slip (Figure 5). At small sizes, the circular asperity is unconditionally stable when its radius R < R c , such that it can only be expected to creep steadily at the rate V pl . At moderate sizes, R c < R < R p , the asperity is linearly unstable such that steady, uniform creep of the asperity can give way to spontaneous acceleration; however, as there is no non-linear pathway to indefinite acceleration towards dynamic rupture, such initial acceleration is ultimately limited to an aseismic transient and the fault will be expected to subsequent decelerate back towards steady state, before presumably repeating the process. For R > R p , the asperity may proceed to seismic slip speeds and dynamic rupture.

In section 2 we presented three state evolution laws proposed by Ruina [1983]: aging, slip, and an intermediate law bridging the first two. The results from linear stability analysis apply to all three forms, as their linearizations about steady state are the same. However, as discussed in section 4.3, the results concerning non-linear phase of nucleation require careful consideration of the form of state evolution. The expressions for the nucleation lengths L f , L p , R f , and R p presented in sections 4.1, 4.2, and 6.2 assumed an aging-law state evolution. Remarkably, these expressions may also be adapted to yield nucleation lengths under state-evolution laws ranging from aging to slip. As mentioned in section 4.3, such nucleation lengths also exist and are a factor ϵ smaller than their aging-law counterparts [Viesca, 2023], where ϵ is the parameter of the intermediate state evolution law. For the aging law, ϵ = 1 and the slip law is retrieved in the limit ϵ → 0. In other words, for a given value of ϵ we may rewrite the free and pinned nucleation lengths as ϵL f , ϵL p , ϵR f , and ϵR p for the single mode or mixed-mode problem configurations.

To highlight the consequence of this reduction in the nucleation length, we consider how the phase diagram of fault behavior changes when considering a state evolution other than aging. Specifically, we consider a finite rate-weakening patch of length 2L loaded by steady creep at a rate V pl beyond the patch. The critical fault size L c at which linear stability is lost is given by 32 and shown in Figure 7a. Additionally in Figure 7a, we show the minimum nucleation length L p for an aging-law fault reduced by a factor ϵ to give the correct corresponding length ϵL p for this 20 Dec 2025 intermediate evolution law. For this particular plot, we consider a state evolution law “halfway” between aging and slip, with ϵ = 0.5. The consequence is the addition of a fourth possible phase of fault behavior, in addition to the previously known phases of steady creep, aseismic transients, and dynamic rupture. This fourth phase occurs when the fault length is such that L p < L < L c . In this case, the fault is linearly stable, but non-linearly unstable. This implies that a steady uniform creep of the rate-weakening patch is stable to small perturbations, but unstable to sufficiently large perturbations: the fault will not spontaneously nucleate a dynamic rupture under the given loading conditions, can be driven to accelerate towards dynamic rupture in the manner of eq:bu if the fault experiences a large-amplitude perturbation. In Figure 7b, we highlight that decreasing the value of ϵ towards the slip law limit enlarges the region in phase space occupied by this fourth phase. In the limit ϵL p → 0 in the slip-law limit, implying that any L < L c is potentially nonlinearly unstable. However, stability analyses of spring-block model with slip and intermediate-law state evolution indicate that the size of the perturbation required for instability increases with the equivalent stiffness of the system, which is expected to scale with the inverse of the fault dimension [Gu et al., 1983;Ciardo and Viesca, 2025]: faults much smaller than L c will need a substantially larger amplitude perturbation to trigger dynamic rupture.

We examine the linear stability of faults under conditions beyond the original analysis of an unbounded fault undergoing uniform sliding. Our analysis takes into account the finiteness of a fault and the manner in which the fault is loaded. We find that loss of linear stability occurs as a Hopf bifurcation when the dimensionless size of the fault L/L b reaches a critical value, L c /L b , which is a function of the ratio of friction parameters a/b. This critical value diverges in the limit a/b → 1 and we find its asymptotic scaling depends on the manner of loading.

In parallel, we also characterize the loss of non-linear stability in terms of non-linear solutions for diverging slip rate introduced previously in the literature [Rubin and Ampuero, 2005;Viesca, 2016a,b;Viesca, 2023] We present an additional solution in which slip rate diverges but is laterally confined, or “pinned”, by fault boundaries. These pinned solutions represent the minimum nucleation length of a fault. We provide accurate, approximate expressions for the dependence of pinned L p on the ratio a/b. We also provide such expressions for larger, so-called “free” nucleation length L f as it depends on a/b, based on previous numerical solutions for L f over the entire range of rate-weakening a/b [Viesca, 2016a,b]. L f is the expected nucleation length when the fault size is larger than 2L f . Unlike the results for linear stability, L p and L f are independent of the manner of loading, with the exception of the issue of spatial confinement, and can be considered universal.

We show that L c and L p may be used to construct a phase diagram to determine the behavior of fault models a priori on the basis of the dimensionless parameters L/L c and L/L p . If L/L c < 1 the system is linearly stable to perturbations. If L/L p < 1 the system is also nonlinearly stable, meaning that no perturbations may lead to a subsequently self-sustaining instability. If model parameters are such that both conditions are met, steady-state condition can be considered to be unconditionally stable. However, if L/L c > 1 and L/L p < 1 the system is linearly unstable but nonlinearly stable, meaning aseismic transients may spontaneously emerge, but cannot continue to accelerate without bound. If both L/L c > 1 and L/L p > 1, the fault is both linearly and non-linearly unstable and small perturbations can lead to a slip rate that diverges until dynamic rupture is nucleated.

For conceptual simplicity, discussions in sections 5 and 6 were framed with fault length L as a model parameter to be adjusted, with all other parameters fixed, and comparing L to critical linear (L c ) and non-linear (L p ) lengths. However, one might consider fault length L to be fixed and adjust one or more of the other problem parameters, such as µ, ν, σ, D c , a, or b. The more general consideration to be given is that of the value of the dimensionless parameters L/L c and L/L p , and specifically, whether they are greater than or less than 1. To calculate these dimensionless parameters, expressions for L c and L p in terms of a/b and L b have been provided, bearing in mind that the elasto-frictional lengthscale L b is itself a function of model parameters, (27). We provide an extension of the results to 3D, identifying a critical radius R c for loss of linear stability of a circular asperity surrounded by steady creep, and the pinned/free nucleation lengths R p and R f . We find the impact the form of state evolution law has on the phase diagram. To do so, we use results from Viesca [2023], who considered an intermediate state evolution law of Ruina [1983], which lies somewhere on a continuum of behavior between aging to slip laws. For the aging law, L c < L p and we have the three possible scenarios discussed above. For deviations from the aging law, L p is reduced by a known factor leading to instances in the phase diagram in which L p < L c . This leads to a fourth possibility that L/L p > 1 and L/L c < 1, meaning the fault is linearly stable but nonlinearly unstable. The implication of this that the fault is susceptible to instability under finitely larger perturbations. The window in phase space for which this can occur grows and is the largest under the slip law. The linear stability results are independent of particular form of state evolution law, among those considered here.

We use the approximation of the integral transform, (A.6), as well as that of V(x), A.8, and discretely evaluate Θ(x) and V (x) at points x i to arrive to a system of 2n -2 equations for 2n -1 unknowns ϕ(s j ) and Θ(x i )

Department of Civil and Environmental Engineering, Tufts University

Department of Civil and Resource Engineering, Dalhousie University

Visitor, Earthquake Research Institute, University of Tokyo