|Home | About | Journals | Submit | Contact Us | Français|
We examine the development of a frictional instability, with diverging sliding rate, at the interface of elastic bodies in contact. Evolution of friction is determined by a slip rate and state dependence. Following Viesca (2016 Phys. Rev. E 93, 060202(R). (doi:10.1103/PhysRevE.93.060202)), we show through an appropriate change of variable, the existence of blow-up solutions that are fixed points of a dynamical system. The solutions show self-similarity of the simple variety: separable dependence of time and space. For an interface with uniform frictional properties, there is a single-problem parameter. We examine the linear stability of these fixed points, as this problem parameter is varied. Specifically, we consider two archetypical elastic settings of the slip surface, in which interactions between points on the surface are either local or non-local. We show that, independent of the nature of elastic interactions, the fixed-points lose stability in the same matter as the parameter is increased towards a limit value: an apparently infinite sequence of Hopf bifurcations. However, for any value of the parameter, the nonlinear development of the instability is attraction, if not asymptotic convergence, towards these fixed points, owing to the existence of stable eigenmodes. For comparison, we perform numerical solutions of the original evolution equations and find precise agreement with the results of the analysis.
The asymptotic behaviour of spatially continuous time-evolution equations may approach self-similarity at finite or long time (e.g. ). Instances in which the asymptotic behaviour is that of a finite-time instability have a long, pervasive and somewhat recent history, with some of the earliest solutions found for the nonlinear heat and Schrödinger equations [2–6]. Subsequent and renewed attention was found in other fields, a prominent example being the self-similar pinchoff of a fluid drop [7–12]. This was followed and paralleled by solutions for self-similar, finite-time instabilities in a variety of other systems: e.g. thin-film rupture [13–17], solidification , solid or liquid surface diffusion [19,20] and problems of aggregation and gravitational collapse [21–23]. Further review on the issues and examples of such problems is found in Eggers & Fontelos [24,25].
A similarity variable transformation reduces the problem of finding such self-similar solutions to one of the findings a fixed point of the transformed evolution equation [2,3]. That such a self-similar solution is ever realized corresponds to the attraction to this fixed point. A stability analysis may be performed to determine whether fixed points are attractive and asymptotically stable, and requires special attention to space and time translational symmetry modes; furthermore, there may exist a family of self-similar fixed points, and a stability analysis often indicates that only one member of the solution family is asymptotically stable [12,15,17,19,22,26–28].
Here, we examine in greater detail, following Viesca , the problem of a sliding instability developing at the interface of similar elastic bodies in contact where frictional strength is governed by a slip rate- and state-dependent friction . We look for and find self-similar solutions with diverging slip rate and determine whether these solutions are asymptotically approached within the finite-time limit of instability. To do so, we find the appropriate change of variables, which involves both scaling considerations to determine similarity variables as well as a change in the nominal state variable. We show that solutions with diverging slip rate are the fixed points of the transformed evolution equations. We then find fixed-point solutions for two elastic configurations: e.g. slip between elastic half spaces or near a boundary, which embody non-local and local interactions between points on the slip surface. The stability of these fixed points is determined by a linear stability analysis. For homogeneous frictional properties, there is a single parameter: the ratio of direct and evolution effect parameters 0<a/b<1. We find that a subset of fixed points are asymptotically stable and lie within a range 0<a/b<(a/b)crit and that above (a/b)crit stability is lost as a cascade of Hopf bifurcations. We compare the results of analysis with those from numerical solutions to the original evolution equations and find agreement.
Laboratory experiments on rock friction at the slow slip rates by Dieterich [31,32] and Ruina [30,33] led to their development that the friction at a position along an interface may be determined by the rate of slip V and a state variable θ encapsulating information on the history of sliding. On the basis of experimental and theoretical considerations, Ruina [30,33] proposed that the dependence may take a form equivalent to
where fi is the coefficient of friction at a reference sliding velocity Vi, a is the direct effect coefficient relating changes in velocity to instantaneous changes in friction, and b is the coefficient of the non-instantaneous evolution effect resulting dependent on an evolving variable state. Dc is a characteristic slip distance over which the state variable evolves. The coordinates of the slip surface are denoted by x: e.g. in its place we may use either the coordinate x or the pair x,y depending on the problem considered. Differentiating the above in time,
The evolution of the state variable was expressed in two currently common forms by Ruina , one of which is referred to as the ageing (or Dieterich–Ruina) law
as evolution is allowed to occur in stationary (V =0) contact. The other common expression, known as the slip (or Ruina–Dieterich) law,
requires a non-zero slip rate for state evolution. Both forms yield the same steady-state dependence of friction on slip velocity
where a<b implies velocity-weakening behaviour of steady state.
The frictional shear strength of a slip surface is τs=σf, where σ is the local surface normal stress (positive in compression). The rate of strength evolution is τs/t=fσ/t+σf/t. Here, we consider cases where σ/t=0 and treat the case where σ is also spatially uniform. The rate of strength is then simply
During incipient or active slip, τs is equal to the magnitude of the shear stress on the slip surface τ. For a bulk material that responds in a linearly elastic manner to slip, τ along the slip surface may be expressed as
where τo is the magnitude of the shear stress in the absence of slip and τel is the elastic shear stress change due to a non-uniform slip distribution. We may write τel as
where is a functional of δ that depends on position, and, for quasi-static elasticity, only on time insofar as δ depends on it. More appropriately for flat surfaces accommodating slip, should be considered as operating not on δ directly, but rather its spatial gradient, because a uniform distribution of slip is a simple rigid translation that induces no stress change; however, we will use the form (2.8) to simplify notation. In many cases of interest, may be written as a spatial convolution over the region of slip that involves spatial derivatives of δ and kernels that depend on the (possibly heterogeneous) elastic properties of the surrounding medium. By the time-independence and assumed linearity of elasticity, has the property that
Of interest is when also has the property that its time derivative is , like for the case of unidirectional slip on a planar surface occurring over a fixed or unbounded region. In such cases, the rate of τel is then
Rice & Ruina  promptly established that the above system is unstable to perturbations when the interface is rate-weakening at steady state (i.e. a<b): a linear stability analysis of steady-state sliding at a uniform rate showed that perturbations above a critical wavelength grow unboundedly with time. That critical wavelength is the same for both evolution laws, given their identical linearizations about steady state. Here, we are interested in the progression of the instability beyond these linear stages, and, more specifically, whether or not there is an asymptotic development of the instability that is independent of the initial conditions and external forcing that provoke the instability.
A quasi-static slip instability is indicated by a divergence in the slip velocity, implying that τel/t will come to dominate any bounded external forcing rate τo/t over the region of divergence such that
We consider state evolution following the ageing law (2.3). To examine behaviour during a slip instability, we may combine the constitutive equations (2.2) and (2.3) with the condition (3.1) to arrive at a system of autonomous evolution equations of V and θ; however, we find that replacing θ with a dimensionless quantity
results in an autonomous system that will be more convenient for analysis:
where Φ can be interpreted as an indicator of the regime of the state of slip, with steady-state sliding corresponding to Φ=0 and Φ=1 corresponding to sliding far from a steady state; Φ>1 is not achievable without backslip.
Given the form of (3.3a) and that Φ is bounded when positive, we find a scaling in which V/t~V 2. Thus, we may expect that a slip instability may occur with V diverging at a finite-time tin and diverging with the inverse of the time from instability, tf(t)=tin−t, motivating a search for solutions of the form
where we will be interested in the behaviour of W in the limit of instability (), and to aid in examining its asymptotic behaviour, we define a new independent variable s as an alternative to time that diverges as . We may simply define s as
For the purposes here, the choice of the particular characteristic time Dc/Vi is arbitrary and inconsequential, given that we may as well have defined s to within an arbitrary constant implicitly by the relation ds/dt=1/tf(t).
Substituting (3.4) into (3.3) and making the change of independent variable from t to s, for which ()/t=(1/tf)()/s, we arrive to the new system
with the notation .
The fixed points of (3.6), and satisfy
To satisfy (3.7b), the fixed points may provide the root of one or both terms in brackets at each point of x within the region of instability. Because (by constraint on range of permissible values of Φ), the root of the latter term can only exist in regions on x where ; however, the root of the former term can exist on regions of x where provided there as well. In other words, we may look for fixed-point solutions with
which requires that, after a rearrangement of (3.7a), satisfy
Looking for such fixed-point solutions over a compact region S (over which the diverging velocity is continuous with for x within S and for x outside and along the boundaries of S), we remarkably find that the resulting problem of determining these solutions has full equivalence to a simple slip-weakening fracture problem. The equivalent problem is that of solving for the slip distribution of an equilibrium slip-weakening crack under a uniform background shear stress τb. Specifically, in the equivalent problem, the strength decays linearly with slip from a peak value τp to residual level τr over an amount of slip δc. The corresponding equation of that problem is
In the former rate-and-state context, there is the implicit requirement that the only singularity in stress rate arrives by the finite-time divergence of V (i.e. there is no singularity in the stress-rate distribution owing to the distribution ). This requirement is provided that the compact region S is not externally constrained by the absence of a continuing slip surface (or one otherwise locked). In the latter slip-weakening context, the equivalent implicit requirement is a non-singular stress distribution on the slip surface, in keeping with the underlying concept of a finite slip-weakening strength. Both requirements impose the same conditions on the distributions of δ(x) and . The problems (3.9) and (3.10) are free boundary problems, in which the boundary to the region S is part of the problem solution; these additional requirements provide the necessary additional conditions for its solution.
We first consider a relatively simple elastic configuration, that of an elastic half-space in which a sliding surface lies parallel to and a distance h away from the free surface. We consider the slip on the surface to follow either in-plane (modeII) or anti-plane (modeIII) conditions. Both conditions presume slip varies only on a single direction (denoted here by x), and the direction of sliding is, respectively, either along or perpendicular to the direction of variation. When the distances over which variations slip on the surface take place are much larger than h, the elastic response of the half-space to slip on the surface tends to that of slip being accommodated entirely by the displacement of material above the slip surface, the deformation of which is uniform in depth. This amounts to expressing the elastic stressing rate as 
The fixed-point solutions are determined by substituting the above operator into (3.9). Immediately following the substitution, we find an apparent characteristic, elasto-frictional length scale
The free boundary problem then reduces to solving for the distribution and the domain over which it is defined. Owing to the problem spatial symmetry, we may define a solution variable L that represents the half-width of the domain of .
Given the simple form of the operator, we arrive to closed-form expressions for the fixed-point solutions (see appendix A for solution details). This involves imposing the non-singular stress-rate condition , and additionally imposing the continuity of and across c (which amount to, respectively, imposing the continuity of rates of slip and stress). The solutions are then, for 0.5≤a/b<1,
where the half-length L and the distance c are given by
The distance c is defined as the point at which . For 0<a/b<0.5, on |x|≤L and the point c does not exist. In this range, the solution is simply that which applies for a/b=0.5. We illustrate these solutions in figure 1.
From the above, we find that the behaviour at ends of the range of 0<a/b<1 are given by
and as the limit to which the end-zone size R(L−c) tends is (π/2)Lbh.
As an aside, we note that the above family of solutions above (parametrized by a/b) do not provide a unique solution for any value of a/b in that range. The above solutions, while varying in form, all have the characteristic of monotonically increasing from at |x|=L to a global maximum at x=0. However, it may be readily shown that a countably infinite number of solution families exist, with each family having an integer number more peaks than the family of solutions detailed above. Specifically, taking a spatial derivative of the problem (3.9) with operator (4.1) results in an eigenvalue problem whose eigenfunctions represent potential solutions for and eigenvalues represent potential solutions for (L/Lbh)2. Of these potential solutions, only a subset will satisfy the boundary conditions on . However, as will be borne out of numerical solutions to the evolution equations (and may also be borne out by a linear stability analysis of these fixed points), the families of multiply peaked solutions are not observable. These multiply peaked fixed points are unstable and unrealized in numerical solutions to the evolution equations, in deference to their singly peaked counterpart. This feature is not unique to the specific problem considered here: the existence of a countable set of similarity solutions frequently appears in other nonlinear problems with finite-time instability, and likewise, it is often the case that a single one of these fixed-point solutions is stable and realizable [12,15,19,22,28].
When the slip surface is far from any material boundaries, the surface may be treated as an interface between two elastic half-spaces for which, under a single mode of slip we may write 
The free boundary problem then again reduces to solving for the distribution and the domain over which it is defined. We again let L denote the half-length of the domain of , which, given that slip rate is diverging on |x|≤L implies that the boundaries of the actively slipping region reduce to L+=L−L.
Defining a dimensionless coordinate and , we may write the problem as
and the problem is to determine L/Lb and for a given a/b. This problem was solved by Viesca , and the solutions are presented therein. The endmember behaviours of the one-parameter family of solutions are
where the latter approximation for was found via a truncated Chebyshev polynomial expansion of and is accurate to within 3.5%. As discussed in , we retrieve the nucleation size scalings of Rubin & Ampuero  as the above endmember scalings of L in the analysis pursued here. Furthermore, similar to the brief aside in the previous section, we may likewise find a countably infinite number of solution families (each parametrized by a/b) that correspond to multiply peaked distributions , but we defer an extended discussion of this point to subsequent work.
Here, we look to determine whether the fixed point solutions found above are asymptotically stable as . We examine stability by looking at the behaviour of small perturbations from the fixed point:
Substituting (5.1) into system (3.6) and linearizing, we arrive to
Looking for solutions to (5.2) of the form and , we arrive to a linear eigenvalue problem to determine the eigenvalues λ and corresponding eigenfunctions ω(x), ϕ(x). The sign of the real part of λ indicates whether that mode grows or decays and the existence of an imaginary part implies that this growth or decay occurs as an oscillation in time t with a diverging frequency chirp as .
Among the resulting eigenvalues, two correspond to perturbations of the time and location of instability [15,17,19,27,28]. A perturbation of the time of instability (by the small value β) leads to a perturbed time from instability tf,ϵ=tf+β, which, recalling that , we may also write as tf,ϵ=tf[1+β/(Dc/Vi)es]. A fixed point solution with this perturbed time to failure can be seen as diverging away from a fixed point solution with an unperturbed time to failure
Specifically, the divergence occurs as a perturbation to the fixed point with an eigenmode, whose shape is itself given by and with the corresponding eigenvalue λ=1. We may likewise consider perturbations to the location of instability by the shift xϵ=x+α
where here no divergence from an unperturbed solution corresponds to a perturbation with eigenvalue λ=0, and when x only consists of one dimension, the corresponding eigenfunction is simply the derivative of with respect to that dimension.
For the class of fixed points where and everywhere in S, (5.2b) simplifies to
Because , the perturbations decay asymptotically, and the asymptotic stability of the fixed point requires that perturbations do not grow as well. Determining this is made somewhat simpler, given the asymptotic decay of (independent of ), which implies that if perturbations of the form satisfying a reduced version of (5.2a)
result in eigenvalues λ with a negative real part, then the fixed point can be said to be asymptotically stable.
We begin with examination of the stability of the fixed points for the limiting case of shallow subsurface slip. Initially, we examine the stability of the fixed points for 0<a/b≤1/2, for which and within |x/Lbh|≤π. For convenience, we use the notation , where we define the dimensionless coordinate and write . Stability of these fixed points hinges on the solution of the eigenvalue problem set by (5.6), which in this case is the singular eigenvalue problem:
with the boundary conditions . We immediately find by inspection the two eigenvalues and eigenfunctions corresponding to perturbations to the instability location and to the time to instability, respectively
In addition to these two discrete eigenmodes, we find the continuous spectrum of eigenvalues λc≤−1/8, with each value of λc having two corresponding eigenmodes: one odd and one even, given by
where and 2F1 is the hypergeometric function. We note that and : i.e. (5.10) generate all eigenfunctions for the discontinuous spectrum . In the absence of any other remaining modes and with all of the eigenvalues being negative (excluding those associated with spatial or temporal translational symmetries), we can conclude that these fixed-point solutions are asymptotically stable.
To determine the fixed-point stability for a/b>1/2 requires solution of the full eigenvalue problem (5.2). We pursue that numerically here, discretizing the closed-form solutions of and (§4a) and the second-order derivative in and looking for the discrete eigenvectors and and corresponding eigenvalues λm of the resulting linear eigenvalue problem. Specifically, we discretize , , and along the set of points for i=0,…,n, and spatial derivatives (i.e. as occur in ) at these points via the spectral Chebyshev differentiation matrix (e.g. ). We then arrive to the system of equations
where the set of points (j=1,…,n) coincide with and the components of Mij are determined following the discretization of (5.2). Solving the linear matrix eigenvalue problem numerically, we arrive to the set of n eigenvalues λm and n eigenvectors .
As an important aside, we note that we may have pursued this discretization to discretize the eigenvalue problem (5.7) to assess the stability of fixed points in the range a/b≤1/2. In this case, discretization would break up the continuous spectrum into a band of discrete eigenvalues with no apparent convergence of individual eigenvalues with increasing resolution (i.e. an increase in grid resolution n leads to a continued ‘filling-in’ of the semi-infinite range of the continuous spectrum). While this would adequately resolve the true discrete modes and , there would be poor convergence towards the even and odd modes at the end of continuous spectrum, and , where the true eigenfunctions exhibit sharp decay to 0 at the interval endpoints. Furthermore, we could not hope for any commonly used discretization (e.g. spectral Chebyshev, local finite difference, etc.) to adequately resolve the even and odd modes for λ<−1/8 owing to oscillatory behaviour at the endpoints of the type (found, for example, by considering the power series expansion method of Frobenius for these singular endpoints). However, we would find that the discretization of the continuous spectrum does not lead to eigenvalues exceeding the upper bound of λc such that the conclusion of fixed point stability for a/b≤1/2 would be robust, if somewhat naively reached.
With this in mind, we continue the pursuit for numerical solutions to the full eigenvalue problem, which retains a characteristically singular behaviour at the endpoints. We anticipate that there will continue to exist a continuous spectrum (or spectra) that will be poorly treated by the discretization; however, we will find that the most relevant eigenvalues to the linear stability analysis will be discrete ones whose corresponding eigenfunctions can be adequately resolved. A prior linear stability analysis for self-similar fixed points of a problem of fluid pinch-off  also observed highly oscillatory eigenmodes and lack of convergence for discrete eigenvalues with increasing discretization, but dismissed this behaviour as artefacts of discretization. Our results here suggest that at least some of these modes may not be artefacts and may be inherently part of the eigenvalue problem, possibly indicative of the existence of a continuous spectrum of eigenvalues with highly oscillatory eigenmodes (although reaching this conclusion requires a closer inspection of that eigenvalue problem).
We find that the fixed points lose stability for a/b>(a/b)c=0=1/2. Stability is first lost as the bifurcation of an eigenvalue from λ=0 with movement along the positive real axis with increasing a/b. The corresponding mode is odd with respect to position x, implying that within a finite range above a/b>1/2, the fixed point may remain asymptotically stable provided symmetry is imposed. However, upon further increase of a/b, we find that a complex–conjugate pair of eigenvalues moves along the complex plane and crosses the imaginary axis as a Hopf-bifurcation at a/b≈0.94. The corresponding eigenmode are even functions and . The trajectories of these eigenvalues appear to have their origin as a bifurcation from the edge of the continuous spectrum at λ=−1/8 near a/b=1/2. Remarkably, we find that this marks the first of a cascade of Hopf bifurcations whose number increase with further increases to a/b, with each successive bifurcation alternating between even and odd modes. The numerical solutions indicate a countably infinite number of distinct modes undergoing Hopf bifurcations as .
What are the implications of a cascade of Hopf bifurcations? Following the first Hopf bifurcation, a limit cycle exists about the fixed point. For a supercritical Hopf bifurcation, that limit cycle is asymptotically stable and its amplitude is expected to grow as a/b is increased beyond its value at bifurcation (e.g. ). Here, the limit cycle would be manifested as the oscillation of the spatial modes and with a period in s determined by the imaginary part of the eigenvalue. A second Hopf bifurcation suggests that the existence of an additional oscillatory cycle. A linear superposition of the two cycles maintains the periodicity of a limit cycle if the frequency of this second mode is rationally related (commensurate) with the first, and quasi-periodic behaviour if it is not. In this latter case, the limit cycle is replaced with a limit torus. In the linear stability analysis conducted here, we find that the imaginary parts of the numerically obtained, Hopf-bifurcating eigenvalues appear to be incommensurate. Here, we may naively expect that additional Hopf bifurcations simply add dimensions to the torus and that behaviour is quasi-periodic with known underlying frequencies; however, it may be expected that such tori may not be stable and that the dynamical system may give rise to chaotic, aperiodic behaviour [39,40]. Therefore, it is of interest to determine the values of a/b at which the Hopf bifurcations occur, and particularly, the value of a/b at which the third bifurcation occurs.
We show the trajectories of the eigenvalues of the first six modes to undergo Hopf bifurcation in figure 2a. The order in which the modes undergo Hopf bifurcations follows the outward sequence of the trajectories. The trajectory of each mode has its origin along the real line, with the bifurcation of each trajectory onto the complex plane occurring at distinct and increasing values of a/b, as seen in figure 3a,b: e.g. the trajectory of the second Hopf-bifurcating mode begins with a/b≈0.67, that of the third begins with a/b≈0.72, etc. However, the Hopf bifurcations occur for a/b close to 1. Specifically, the additional second and third even-mode bifurcations occur at a/b≈0.98 and a/b≈0.994, and the first and second odd-mode bifurcations occur at a/b≈0.96 and a/b≈0.99. In figure 4, we show the eigenfunctions for the first even eigenmode to undergo a Hopf bifurcation, ωH(x) and ϕH(x), at fixed values of a/b.
Investigations of the stability of self-similar gravitational collapse, pinchoff of one inviscid fluid within another, and buoyant thermal plumes have shown a loss of stability by a single Hopf bifurcation as the sole parameter of each problem (a coupling constant, a density contrast and the Reynolds number) is decreased or increased in the latter two [27,41,42]. Leppinen & Lister  observe a single Hopf bifurcation in their linear stability analysis, while Peng & Lister  additionally observe the existence of a limit cycle in solutions to the original evolution equations within a parameter range following bifurcation that is in agreement with the unstable mode of their linear stability analysis. Hirschmann & Eardley  also allude to indications for the appearance of other unstable modes, possibly through one or more subsequent Hopf bifurcations. However, as far as the author is aware, the problem examined here and in  is the first demonstration of a dynamical system whose asymptotic self-similar behaviour becomes unstable through successive Hopf bifurcations.
As briefly discussed in Viesca , the loss of fixed-point stability as may have been anticipated considering that, in this limit, the fixed-points yield solutions of diverging slip under steady-state conditions and that uniform steady-state slip is known to be unstable to perturbations greater than a critical wavelength . As , over |x|≤L, except within a boundary layer near the endpoints |x|=L, where vanishes. The boundary layer width scales with Lbh: the size of the region near |x|=L over which is RL−c=πLbh in this limit. Our measure of distance from steady-state sliding, Φ, for the fixed-points is given by . Because when , over |x|≤L, excluding the boundary layer, as . Repeating the stability analysis of Rice & Ruina  using the local elastic operator (4.1), we find that the critical wavelength for stability of steady-state sliding is
Given the asymptotic behaviour of L as in the preceding section, in this limit. This evaluates to be for a/b≈0.94, the point at which the first Hopf bifurcation occurs.
We now consider the stability of fixed points for the elastic configuration of in-plane or anti-plane slip between elastic half-spaces, as discussed in Viesca . The key variation introduced from the previous configuration is the non-locality of the elasticity, and this makes the fixed-point solutions more precarious: i.e. loss of stability occurs at lower values of a/b.
The loss of stability occurs in the same manner observed in the preceding section. Namely, stability is first lost as a/b increases beyond (a/b)c=0 as an eigenvalue traverses into, and subsequently travels along, the positive real axis (and again, the corresponding eigenmode is odd). Further increases in a/b again result in a cascade of Hopf bifurcations. The eigenvalue trajectories of the first six Hopf-bifurcating modes are shown in figure 2b. Here, the order of bifurcations is slightly altered: the innermost trajectory is now the second mode to undergo a Hopf bifurcation and is preceded by the second innermost trajectory. Consequently, the first Hopf-bifurcating mode is now odd and the second and third are even. Subsequent bifurcations return to the alternation between odd and even modes. The values of a/b at which these bifurcations occur are significantly further from the rate-neutral limit of a/b=1 than the corresponding values of a/b in the preceding section: the first Hopf bifurcation (odd-mode) occurs at a/b≈0.72, the second and third Hopf bifurcations (even modes) at a/b≈0.74 and a/b≈0.79, respectively. The third even-mode bifurcation occurs at a/b≈0.87. The second and third odd-mode bifurcations occur at a/b≈0.84 and a/b≈0.91. An interesting departure in behaviour of the trajectories from that observed in the previous section is the return of the first Hopf-bifurcating pair of eigenvalues to the real axis. The complex-conjugate pair are observed to meet at the real axis when a/b≈0.86. Subsequently, the two eigenvalues follow opposite directions along the real axis. In figure 5, we show the eigenfunctions for the first even eigenmode to undergo a Hopf bifurcation, ωH(x) and ϕH(x), at various values of a/b. Here, numerical solution of the full eigenvalue problem (5.2) is pursued in a similar manner as in the preceding section: i.e. the problem is reduced to one of the form (5.11) by discretization. Here, the discretization is a piecewise-constant approximation of ω and ϕ over evenly spaced intervals between x=±L.
Likewise, the onset of stability loss reflects L exceeding λcr, the critical wavelength for unstable perturbations to steady-state slip. The elastic configuration considered here was one considered by Rice & Ruina , for which
Comparing this to the asymptotic behaviour for L, we find that L/λcr=1/[π2(1−a/b)], which is for a/b≈0.74, the value at which the first Hopf bifurcation occurs.
We find numerical solutions to the quasi-static evolution equations for slip rate and state with initial conditions and external forcing that spontaneously induce a slip instability (methods detailed in the electronic supplementary material). We focus on the configuration of single-mode slip between elastic half-spaces, but also consider single-mode slip on a planar surface below and parallel to a free surface. We find that the evolution of slip rate and state during a spontaneously developing instability is remarkably well described by considering the dynamical system’s fixed points and their asymptotic stability.
For simplicity, we take as initial conditions steady-state slip under uniform slip rate Vi. Instability is induced by imposing a compact rate of external loading τo/t, held constant in time beginning at t=0+ with the spatial distribution
for |x|≤Lτ and 0 otherwise. Lτ is the half-length of the loaded region, and may be chosen, for example, to have the length Lb/(1−a/b) or Lb/(1−a/b)2. However, we emphasize that the linear stability analysis yields the expectation that the system, if provoked to develop an instability, will be attracted to the fixed points (though not necessarily with asymptotic convergence), regardless of the specifics of the initial conditions or external forcing.
Given the symmetry of the problem posed to be solved numerically, fixed points are expected to be asymptotically approached for a/b below 0.74. We first examine slip rate and state evolution for a/b in this range. On the basis of our linear stability analysis, fixed point stability is lost to symmetric perturbations as a series of Hopf bifurcations, the first three of which occur at the approximate a/b values of 0.74, 0.79 and 0.87. We subsequently examine the change in the character of slip instability development following each bifurcation.
The solutions require numerical evaluation of the Hilbert transform and we use two methods to do so. For initial comparisons with fixed point solutions in the range 0<a/b≤0.7 (figure 6), we use a commonly used spectral method relying on the convolution property of the Fourier transform and the fast Fourier transform to numerically evaluate the Hilbert transform on a regularly spaced grid (e.g. ). However, this method has the deficiency of an implicit periodic replication of a finite domain, whereas the fixed-point solutions correspond to the development of a single slip instability on the real line. For more precise comparisons in the vicinity of Hopf bifurcations (a/b≥0.7) (figures 7 and and8),8), we use a lesser known but similarly advantageous and implementable spectral method that evaluates derivatives and the Hilbert transform on the entire real line [44,45]. Further details are available in the electronic supplementary material.
Our solutions to the evolution equations confirm our expectations that the fixed points are attractive and asymptotically stable over the range : i.e. when an instability does develop, it will ultimately develop in the manner prescribed by the fixed points. As examples, we present the evolution of a normalized slip rate, stress rate, and distance from steady state for three numerical solutions in figure 6 (grey scale lines), with values a/b=0.5,0.6,0.7. We find that fixed point is approached asymptotically. We present the solutions as a series of snapshots in time where the fixed-point distributions are overlain, red-dashed lines. Pointwise convergence to a fixed point can be seen, for example, in fig. 4a in . Among these solutions, we find that convergence is slower (i.e. requires the development to a larger peak slip rate) as the value of a/b increases. This slow convergence corresponds to an eigenmode approaching the first Hopf bifurcation (i.e. the real part of its eigenvalue is decreasingly negative, resulting in a slower, oscillatory decay towards the fixed point as ).
Realistically, inertia will ultimately limit slip rate, and a natural question is how closely will the system approach the fixed point by the time a dynamic rupture is finally nucleated. As cyan-dashed lines in figure 6, we overlay the numerical solution at the instant when , a rough estimate of a slip rate when dynamic stress changes become comparable to those developed quasi-statically (given a characteristic shear modulus, shear wave speed and normal stress, and taking inter seismic plate rates for Vi). We find that the characteristic nucleation length, slip-rate distributions and distance from steady state of the fixed-point solutions are quite close to those observed in the numerical solutions at the predicted onset point for dynamic effects. The conclusions that the fixed points are asymptotically approached for a/b in this range and that they may be representative of behaviour in moments just prior to dynamic rupture is consistent with a retrospective evaluation of numerical solutions performed by others under a variety of initial conditions and forcing.
Under the imposed symmetry, the stability of the fixed points is lost following the first even mode Hopf bifurcation occurring at a/b≈0.74 and we now move to study the behaviour approaching and beyond this critical value where a succession of such bifurcations occur. For the numerical solutions of the evolution equation in this part, we make use of the spectral method that does not resort to a periodic replication of a finite domain to make precise comparisons of behaviour with expectations set by the fixed point and our analysis of its stability.
Figure 7 shows the evolution of spatial slip rate distributions along the slip surface (normalized by the initial value Vi) at snapshots in time for a/b within the interval 0.7≤a/b<1. The instants are chosen as the peak slip rate reaches incrementally larger orders of magnitude. For a/b=0.7, figure 7a, we again confirm the expectation that the fixed-point is asymptotically stable. For a/b=0.75, figure 7b, which follows the first Hopf bifurcation, we observe periodic oscillations in the distribution of the diverging slip rate. The periodic behaviour is soon explored in greater detail, but can be seen in oscillations apparent in the fringes of the velocity distributions for |x| near L. We may similarly identify nearly quasi-periodic behaviour for a/b=0.8, figure 7c, which follows the second Hopf bifurcation. For a/b=0.85 and 0.9, figure 7d,e, the behaviour appears to be aperiodic, with the acceleration of slip rate varying between occurring over an extended region of size 2L to occurring over a more localized region. I
The deviations from the fixed point solutions are controlled by the unstable eigenmodes. For example, the periodic behaviour observed for a/b=0.75 is described by the eigenmode which underwent the first Hopf bifurcation and whose components for perturbations from and are denoted ωH(x) and ϕH(x), respectively. In figure 8a, we show the deviations of Φ from the fixed point distribution through one period of the limit cycle, which agrees with the oscillations of the eigenmode ϕH through one period, as shown in figure 8b. Likewise, we note that the most prominent variations in slip rate seen in figure 7c–e are also consistent with the unstable mode shapes ωH(x) of figure 5 for the corresponding values of a/b: i.e. the expanding crack-like behaviour observed by Rubin & Ampuero  and Kaneko & Ampuero , and similarly shown in in77d,e, is consistent with the emergence of the corresponding unstable eigenmode in figure 7c,d.
We now numerically solve for the evolution of slip rate and state considering the local elasticity for slip below a free surface in the limit that hLb, in which elasticity reduces to a local operator. We use the same Fourier spectral approach used for the non-local operator adapted for the local operator, although we may as well have used a finite-difference discretization. While there is an implicit replication of the domain, this has no influence in the asymptotic development of the instability owing to the local nature of the elastic interactions (unlike the elastic operator considered in the preceding subsection). Our initial conditions are again of uniform initial slip rate V =Vi and steady state Φ=0 and the instability is provoked by a constant external stress rate distribution imposed in the form (6.1) with Lτ=Lbh/(1−a/b).
Here, we make here only a cursory comparison with the behaviour observed in the preceding behaviour to highlight a key difference: the rapid convergence of the numerical solutions to the fixed-point. In figure 9, we show the results of the numerical solutions on the spatial development of state, slip and stress rates for values of a/b in which the fixed points are stable under symmetric conditions. We show results up to a maximum value of the slip rate eight orders of magnitude larger than its initial value. Owing to the relative stability of the fixed points, they are approached relatively more quickly at the same values of a/b when compared with the results of the preceding section.
We examined the development of a sliding instability at the interface of linearly elastic bodies in contact, whose frictional strength is determined by a nonlinear slip rate- and state-dependent frictional law. We find that such a system has the following generic features:
While we find that the above results are independent of the elastic configuration, the specifics of some features (e.g. critical values of a/b, as well as the form of the distributions of diverging slip rate) do depend on the nature of the elastic bodies in contact and the geometry of sliding (the mode of slip). We find the blow-up solutions for two archetypical configurations. We examine frictional contact between elastically identical bodies: specifically, two half-spaces in contact or a layer sliding above an elastic half-space. These embody, respectively, non-local and local elastic interactions. For each case, we find solutions for single mode slip (in-plane or anti-plane).
We perform numerical solutions to the original pair of evolution equations, in which an instability is provoked; we compare these solutions to the expectations set from the above analysis of nonlinear stability development as fixed-point attraction. We focus on single-mode slip between elastic half-spaces. We find that the behaviour of the numerical solutions to the equations conforms precisely to the expectations set by the fixed points and their stability. Below a critical value of a/b, the solutions asymptotically converge to the blow-up solution; we observe the expected limit cycle for a/b increased just beyond the value at which the first Hopf bifurcation occurs; and we find aperiodic behaviour for solutions following the second and third Hopf bifurcations.
We examine the asymptotic limit when variations in slip occur over distances much greater than the depth h of the slip surface below the free surface (expected when h/Lb1). In this limit, material overlying the crack undergoes approximately uniform deformation with depth and the material below is effectively rigid. Equilibrium of the overlying layer then requires that the in-plane horizontal normal stress (mode II) or out-of-plane shear stress (mode III) σh(x) satisfies
where τo is the shear stress acting on the plane of the rupture in the absence of slip and ϵ is the extensional (mode II) or out-of-plane shear strain in the layer (mode III). Given that displacements of the material below the slip surface are negligible, gradients in slip δ(x) are wholly accommodated by strain in the overlying layer
Combining the above, we find that the slip along the crack must satisfy
where τ(x,t) is determined by the strength of the slipping region and the left-hand side is . This result may also be deduced by examining the long-wavelength (shallow-depth) limit of an integral transform . The reduction of the problem to such a simple state of deformation has been leveraged in models of landslide and snow-slab avalanche initiation [47–49], and in a similar form with inertial terms as a model elastic system for earthquake rupture [50–52].
Using the above operator for in (3.9),
we find that the distribution has the form
where the first expression was found using symmetry and the condition and the latter found by imposing and at x=L (implying non-singular stress rate beyond the patch of instability). The lengths of L and c remain to be solved for. Imposing on the latter expression and rearranging, we find that
where we abbreviate the length scale .
The author has no competing interests.
The author acknowledges support from the National Science Foundation (grant no. EAR-1344993).