Search tips
Search criteria 


Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
Proc Math Phys Eng Sci. 2016 August; 472(2192): 20160254.
PMCID: PMC5014108

Self-similar slip instability on interfaces with rate- and state-dependent friction


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.

Keywords: sliding friction, fracture, self-similar nonlinear instability, rate- and state-dependent friction

1. Introduction

The asymptotic behaviour of spatially continuous time-evolution equations may approach self-similarity at finite or long time (e.g. [1]). 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 [26]. Subsequent and renewed attention was found in other fields, a prominent example being the self-similar pinchoff of a fluid drop [712]. This was followed and paralleled by solutions for self-similar, finite-time instabilities in a variety of other systems: e.g. thin-film rupture [1317], solidification [18], solid or liquid surface diffusion [19,20] and problems of aggregation and gravitational collapse [2123]. 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,2628].

Here, we examine in greater detail, following Viesca [29], 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 [30]. 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.

2. Governing equations

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 [33], one of which is referred to as the ageing (or Dieterich–Ruina) law

θ(x,t)t=1V(x,t)θ(x,t)Dc(ageing law)

as evolution is allowed to occur in stationary (V =0) contact. The other common expression, known as the slip (or Ruina–Dieterich) law,

θ(x,t)t=V(x,t)θ(x,t)Dcln[V(x,t)θ(x,t)Dc](slip 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 [partial differential]τs/[partial differential]t=f[partial differential]σ/[partial differential]t+σ[partial differential]f/[partial differential]t. Here, we consider cases where [partial differential]σ/[partial differential]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 L 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, L 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, L 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, L has the property that


Of interest is when L[ x;g( x,t)] also has the property that its time derivative is L[ x;g( x,t)/t], 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 [34] 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.

3. Slip instability as a fixed point of a dynamical system

A quasi-static slip instability is indicated by a divergence in the slip velocity, implying that [partial differential]τel/[partial differential]t will come to dominate any bounded external forcing rate [partial differential]τo/[partial differential]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 [partial differential]V/[partial differential]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)=tint, motivating a search for solutions of the form


where we will be interested in the behaviour of W in the limit of instability (tf0), and to aid in examining its asymptotic behaviour, we define a new independent variable s as an alternative to time that diverges as ttin. 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 [partial differential]([center dot])/[partial differential]t=(1/tf)[partial differential]([center dot])/[partial differential]s, we arrive to the new system




with the notation g~( x,s)g[ x,t(s)].

The fixed points of (3.6), W( x,s)=W( x) and Φ~( x,s)=P( x) 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 P1 (by constraint on range of permissible values of Φ), the root of the latter term can only exist in regions on x where W( x)1; however, the root of the former term can exist on regions of x where W( x)1 provided P( x)=1 there as well. In other words, we may look for fixed-point solutions with


which requires that, after a rearrangement of (3.7a), W satisfy


Looking for such fixed-point solutions over a compact region S (over which the diverging velocity is continuous with W>0 for x within S and W=0 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 W). 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 W( x). 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.

4. Fixed-point solutions for two archetypical elastic settings

(a) Slip below a free surface

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 (mode II) or anti-plane (mode III) 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 [29]

τel(x,t)t=L[x;V(x,t)]E¯h2V(x,t)x2E¯={2μ(1ν)mode IIμmode III.

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 W(x) 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 W.

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 W(±L)=0, and additionally imposing the continuity of W and W 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,

W(x)={1+1a/b2c2x2Lbh2for |x|<cab[1cos(L|x|Lbh)]for c<|x|L,

where the half-length L and the distance c are given by




The distance c is defined as the point at which W=1. For 0<a/b<0.5, W<1 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.

Figure 1.
Fixed-point solutions which yield a diverging slip rate (a) with distribution W and (b) occurring with a fixed measure of distance from steady-state sliding, P. The solutions shown are for the elastic configuration of a thin layer sliding over an elastically ...

From the above, we find that the behaviour at ends of the range of 0<a/b<1 are given by


and as a/b1 the limit to which the end-zone size R[equivalent](Lc) 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 W=0 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 W and eigenvalues represent potential solutions for (L/Lbh)2. Of these potential solutions, only a subset will satisfy the boundary conditions on W. 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].

(b) Slip between elastic half-spaces

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 [35]

τel(x,t)t=L[x;V(x,t)]μ¯2πLL+V(s,t)/ssxdsμ¯={μ(1ν)mode IIμmode III,

where L+ and L are the boundaries of the slipping region. Introducing the operator (4.4) to the problem of (3.9), we find the elasto-frictional length scale of [36]


The free boundary problem then again reduces to solving for the distribution W(x) and the domain over which it is defined. We again let L denote the half-length of the domain of W, which, given that slip rate is diverging on |x|≤L implies that the boundaries of the actively slipping region reduce to L+=L[equivalent]L.

Defining a dimensionless coordinate x^=x/L and W^(s^)=W(s^L), we may write the problem as


and the problem is to determine L/Lb and W^(x^) for a given a/b. This problem was solved by Viesca [29], and the solutions are presented therein. The endmember behaviours of the one-parameter family of solutions are


where the latter approximation for W was found via a truncated Chebyshev polynomial expansion of W/1x2 and is accurate to within 3.5%. As discussed in [29], we retrieve the nucleation size scalings of Rubin & Ampuero [36] 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 W, but we defer an extended discussion of this point to subsequent work.

5. Fixed-point stability

Here, we look to determine whether the fixed point solutions found above are asymptotically stable as s. 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 Wϵ( x,s)=ω( x)eλs and Pϵ( x,s)=ϕ( x)eλs, 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 tf0.

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 tf=(Dc/Vi)es, 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 W with an eigenmode, whose shape is itself given by W 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 W with respect to that dimension.

For the class of fixed points where P=1 and W1 everywhere in S, (5.2b) simplifies to


Because W<1, the perturbations Pϵ decay asymptotically, and the asymptotic stability of the fixed point requires that perturbations Wϵ do not grow as well. Determining this is made somewhat simpler, given the asymptotic decay of Pϵ (independent of Wϵ), which implies that if perturbations of the form Wϵ=ω( x)eλs 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.

(a) Slip below a free surface

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 P=1 and W1 within |x/Lbh|≤π. For convenience, we use the notation W~(x~)=(a/b)[1+cos(x~)], where we define the dimensionless coordinate x~=x/Lbh and write f~(s~)=f(s~Lbh). 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 ω~(±π)=0. 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 κ(λc)=(1+1+8λc)/4 and 2F1 is the hypergeometric function. We note that ω~c,o(x~;0)=ω~0(x~)/2 and ω~c,e(x~;1)=ω~1(x~): i.e. (5.10) generate all eigenfunctions for the discontinuous spectrum λ(,1/8],0,1. 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 W and P (§4a) and the second-order derivative in L and looking for the discrete eigenvectors ω~m and ϕ~m and corresponding eigenvalues λm of the resulting linear eigenvalue problem. Specifically, we discretize W, P, ω~m and ϕ~m along the set of points x~i=cos(iπ/n)(L/Lbh) for i=0,…,n, and spatial derivatives (i.e. as occur in L) at these points via the spectral Chebyshev differentiation matrix (e.g. [37]). We then arrive to the system of equations


where the set of points x~j (j=1,…,n) coincide with x~i 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 {ω~mϕ~m}T.

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 ω~0 and ω~1, there would be poor convergence towards the even and odd modes at the end of continuous spectrum, ω~c,e(x~;1/8) and ω~c,o(x~;1/8), 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 cos[ln(L/Lbh|x~|)] (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 [27] 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 ω~H(x~) and ϕ~H(x~). 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 a/b1.

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. [38]). Here, the limit cycle would be manifested as the oscillation of the spatial modes ω~H(x~) and ϕ~H(x~) 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.

Figure 2.
Trajectories in the complex plane of the eigenvalues λ of the first six eigenmodes to undergo Hopf bifurcations as the parameter a/b is increased from (a/b)c=0 towards a value of 1. These trajectories are found by numerically solving the fixed-point ...
Figure 3.
Evolution with a/b of the real and imaginary parts of the eigenvalues shown in figure 2: (a,b) correspond to the eigenvalues shown in figure 2a, the case in which elastic interactions are local, and (c,d) correspond to the eigenvalues in figure 2b, the ...
Figure 4.
(Solid lines) real and imaginary parts of the first even eigenmode mode to undergo a Hopf bifurcation for the near-surface elastic configuration of figure 2a. This mode consists of ωH(x), ϕH(x), which are the mode shapes for the perturbations ...

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 [27] observe a single Hopf bifurcation in their linear stability analysis, while Peng & Lister [42] 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 [41] 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 [29] is the first demonstration of a dynamical system whose asymptotic self-similar behaviour becomes unstable through successive Hopf bifurcations.

As briefly discussed in Viesca [29], the loss of fixed-point stability as a/b1 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 [34]. As a/b1, W1 over |x|≤L, except within a boundary layer near the endpoints |x|=L, where W vanishes. The boundary layer width scales with Lbh: the size of the region near |x|=L over which 0<W<1 is R[equivalent]Lc=πLbh in this limit. Our measure of distance from steady-state sliding, Φ, for the fixed-points is given by P. Because P=1/W when W>1, P0 over |x|≤L, excluding the boundary layer, as a/b1. Repeating the stability analysis of Rice & Ruina [34] 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 a/b1 in the preceding section, L/λcr1/(2π1a/b) in this limit. This evaluates to be O(1) for a/b≈0.94, the point at which the first Hopf bifurcation occurs.

(b) Slip between elastic half-spaces

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 [29]. 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 xL.

Figure 5.
Real and imaginary parts of the first even eigenmode mode to undergo a Hopf bifurcation (consisting of ωH(x), ϕH(x)) for the contacting half-spaces elastic configuration of figure 2b. This mode is that which follows the innermost trajectory ...

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 [34], for which


Comparing this to the asymptotic behaviour for L, we find that Lcr=1/[π2(1−a/b)], which is O(1) for a/b≈0.74, the value at which the first Hopf bifurcation occurs.

6. Comparisons with numerical solutions to evolution equations

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 [partial differential]τo/[partial differential]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.

(a) Slip between elastic half-spaces

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. [43]). 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.

Figure 6.
Plot of normalized slip rate (a(i),b(i),c(i)), stress rate (a(ii),b(ii),c(ii)) and dimensionless distance from steady-state (a(iii),b(iii),c(iii)) for three values of a/b that are below the critical value at which fixed-point stability to symmetric perturbations ...
Figure 7.
Normalized slip rate distributions at increments in time from numerical solution of the evolution equations for slip rate and state, with a stress-rate perturbation of the form (6.1) with Lτ=Lb/[π(1−a/b)2] for each a/b. Time increments ...
Figure 8.
(a) Plots of the difference between the distribution of Φ given by the fixed point solution P(x) (for a/b=0.75) and that resulting from the numerical solution of the problem corresponding to figure 7b and fig. fig.4(b)4(b) in [29]. The ...

Our solutions to the evolution equations confirm our expectations that the fixed points are attractive and asymptotically stable over the range 0<a/b0.74: 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 [29]. 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 s).

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 max[V(x,t)]/Vi=108, 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 W and P are denoted ωH(x) and ϕH(x), respectively. In figure 8a, we show the deviations of Φ from the fixed point distribution P(x) 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 7ce 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 [36] and Kaneko & Ampuero [46], and similarly shown in in77d,e, is consistent with the emergence of the corresponding unstable eigenmode in figure 7c,d.

(b) Slip below a free surface

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 h[double less-than sign]Lb, 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.

Figure 9.
(grey scale) A version of figure 6 now for the elastic configuration for a slip surface near (h[double less-than sign]Lb) and parallel to a free surface. Plot of normalized slip rate (a(i),b(i),c(i)), stress rate (a(ii),b(ii),c(ii)) and dimensionless distance from ...

7. Conclusion and summary of results

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:

  •  there exist solutions to the governing equations in which slip rate diverges in a finite time with a fixed spatial distribution. We consider a transformation of the initial pair of nonlinear evolution equations for slip rate and state, on the basis of considerations of scaling and the state of sliding. These blow-up solutions are fixed points of the transformed dynamical system. The development of the nonlinear instability of this elastic system is in common with that of other nonlinear, time-dependent problems with asymptotic self-similar solutions, including several instabilities in fluid dynamics;
  •  we show that finding the blow-up solutions (the fixed points) reduces to solving an equivalent cohesive-zone fracture mechanics problem, that of a piecewise linearly slip-weakening fracture. This problem similarity implies that stress and slip distributions during instability may appear to be consistent with linear slip-weakening, despite the actual dependence of frictional strength being on slip rate and its history. There is a single problem parameter 0<a/b<1;
  •  as fixed points of an infinite-dimensional dynamical system, we find that these blow-up solutions are attractive by linear stability analysis; however, they become less so as problem parameter a/b is increased. The fixed points lose stability as discrete modes become unstable;
  •  the linear stability analysis of the fixed points reduces to a singular eigenvalue problem with both discrete and continuous spectra. The fixed points are asymptotically stable for a/b below a critical value. As the parameter a/b increases, discrete eigenvalues are observed to branch off of the continuous part of the spectrum as complex conjugate pairs. A loss of stability follows as each pair undergoes a Hopf bifurcation in an apparently infinite succession as a/b1 from below; and
  •  the results of the linear stability analysis implies that the instability development is universal provided the value of a/b for an interface is below a critical value (when the fixed points are asymptotically stable). Above a critical value of a/b, possibly as early as that value of a/b at which the third Hopf bifurcation occurs, the development will be chaotic. Despite the appearance of several unstable modes, the fixed point will still be attractive (if not asymptotically so) in that there will exist a family of stable modes by way of which the system may be drawn towards the fixed point, before being expelled away by way of excitation of the unstable modes. In this case, the fixed point solutions still provide useful information (e.g. characteristic lengths for slip rate variations), though the precise evolution of slip rate is unpredictable and highly sensitive to initial conditions.

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.

Supplementary Material

Supplementary material for: Self-similar slip instability on interfaces with rate- and state-dependent friction:

Appendix A. Fixed point solutions for shallow slip below free surface

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/Lb[double less-than sign]1). 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

hσh(x,t)x=τ(x,t)τo(x,t),withσh(x,t)=E¯ϵ(x,t)E¯={2μ(1ν)mode IIμmode III,
A 1

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

A 2

Combining the above, we find that the slip along the crack must satisfy

A 3

where τ(x,t) is determined by the strength of the slipping region and the left-hand side is τel(x,t)=L[x;δ(x,t)]. This result may also be deduced by examining the long-wavelength (shallow-depth) limit of an integral transform [29]. The reduction of the problem to such a simple state of deformation has been leveraged in models of landslide and snow-slab avalanche initiation [4749], and in a similar form with inertial terms as a model elastic system for earthquake rupture [5052].

Using the above operator for L in (3.9),

A 4

we find that the distribution W has the form

W(x)={σbE¯hDca/b12(c2x2)+1for |x|<cab(1cos[(L|x|)σbE¯hDc])for c<|x|L,
A 5

where the first expression was found using symmetry and the condition W(c)=1 and the latter found by imposing W=0 and dW/dx=0 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 W(c)=1 on the latter expression and rearranging, we find that

A 6

which restricts the set of solutions 0≤c/L<1 to 0.5≤a/b<1. Requiring that the stress rate be continuous across x=c implies that W(c)=W(c+), which along with (A.5) and (A.6) can be used to find that

A 7

where we abbreviate the length scale E¯hDc/(σb)Lbh.

Competing interests

The author has no competing interests.


The author acknowledges support from the National Science Foundation (grant no. EAR-1344993).


1. Barenblatt GI. 1996. Scaling, self-similarity, and intermediate asymptotics. Cambridge, UK: Cambridge University Press.
2. Giga Y, Kohn RV 1985. Asymptotically self-similar blow-up semilinear heat equations. Commun. Pur. Appl. Math. 38, 297–319. (doi:10.1002/cpa.3160380304)
3. Giga Y, Kohn RV 1987. Characterizing blowup using similarity variables. Indiana Univ. Math. J. 36, 1–40. (doi:10.1512/iumj.1987.36.36001)
4. Zakharov VE, Kuznetsov EA, Miusher SL 1985. Semiclassical regime of a three-dimensional wave collapse. Pisma Zh. Eksp. Teor. Fiz. 41, 125–127.
5. McLaughlin DW, Papanicolaou GC, Sulem C, Sulem PL 1986. Focusing singularity of the cubic Schrödinger equation. Phys. Rev. A 34, 1200–1210. (doi:10.1103/PhysRevA.34.1200) [PubMed]
6. Rypdal K, Rasmussen JJ 1986. Blow-up in nonlinear Schroedinger equations-II similarity structure of the blow-up singularity. Phys. Scr. 33, 498–504. (doi:10.1088/0031-8949/33/6/002)
7. Eggers J. 1993. Universal pinching of 3D axisymmetric free-surface flow. Phys. Rev. Lett. 71, 3458–3460. (doi:10.1103/PhysRevLett.71.3458) [PubMed]
8. Papageorgiou DT. 1995. On the breakup of viscous-liquid threads. Phys. Fluids 7, 1529–1544. (doi:10.1063/1.868540)
9. Almgren R, Bertozzi A, Brenner MP 1996. Stable and unstable singularities in the unforced Hele–Shaw cell. Phys. Fluids 8, 1356 (doi:10.1063/1.868915)
10. Day RF, Hinch EJ, Lister JR 1998. Self-similar capillary pinchoff of an inviscid fluid. Phys. Rev. Lett. 80, 704–707. (doi:10.1103/PhysRevLett.80.704)
11. Lister JR, Stone HA 1998. Capillary breakup of a viscous thread surrounded by another viscous fluid. Phys. Fluids 10, 2758–2764. (doi:10.1063/1.869799)
12. Zhang WW, Lister JR 1999. Similarity solutions for capillary pinch-off in fluids of differing viscosity. Phys. Rev. Lett. 83, 1151–1154. (doi:10.1103/PhysRevLett.83.1151)
13. Bertozzi AL, Brenner MP, Dupont TF, Kadanoff LP 1994. Singularities and similarities in interface flows. In Trends and perspectives in applied mathematics (eds L Sirovich, F John, JE Marsden), pp. 155–208. Applied Mathematical Sciences, vol. 100. New York, NY: Springer.
14. Zhang WW, Lister JR 1999. Similarity solutions for van der Waals rupture of a thin film on a solid substrate. Phys. Fluids 11, 2454–2462. (doi:10.1063/1.870110)
15. Witelski TP, Bernoff AJ 1999. Stability of self-similar solutions for van der Waals driven thin film rupture. Phys. Fluids 11, 2443–2445. (doi:10.1063/1.870138)
16. Witelski TP, Bernoff AJ 2000. Dynamics of three-dimensional thin film rupture. Physica D 147, 155–176. (doi:10.1016/S0167-2789(00)00165-2)
17. Bernoff AJ, Witelski TP 2002. Linear stability of source-type similarity solutions of the thin film equation. Appl. Math. Lett. 15, 599–606. (doi:10.1016/S0893-9659(02)80012-X)
18. Bernoff AJ, Bertozzi AL 1995. Singularities in a modified Kuramoto-Sivashinsky equation describing interface motion for phase-transition. Physica D 85, 375–404. (doi:10.1016/0167-2789(95)00054-8)
19. Bernoff AJ, Bertozzi AL, Witelski TP 1998. Axisymmetric surface diffusion: dynamics and stability of self-similar pinchoff. J. Stat. Phys. 93, 725–776. (doi:10.1023/
20. Wong H, Miksis MJ, Voorhees PW, Davis SH 1998. Universal pinch off of rods by capillarity-driven surface diffusion. Scr. Mater. 39, 55–60. (doi:10.1016/S1359-6462(98)00127-4)
21. Brenner MP, Witelski TP 1998. On spherically symmetric gravitational collapse. J. Stat. Phys. 93, 863–899. (doi:10.1023/B:JOSS.0000033167.19114.b8)
22. Brenner MP, Constantin P, Kadanoff LP, Schenkel A, Venkataramani SC 1999. Diffusion, attraction and collapse. Nonlinearity 12, 1071–1098. (doi:10.1088/0951-7715/12/4/320)
23. Huang Y, Bertozzi AL 2010. Self-similar blowup solutions to an aggregation equation in Rn. SIAM J. Appl. Math. 70, 2582–2603. (doi:10.1137/090774495)
24. Eggers J, Fontelos MA 2009. The role of self-similarity in singularities of partial differential equations. Nonlinearity 22, R1–R44. (doi:10.1088/0951-7715/22/1/R01)
25. Eggers J, Fontelos MA 2015. Singularities: formation, structure, and propagation. Cambridge Texts in Applied Mathematics Cambridge, UK: Cambridge University Press.
26. Brenner MP, Lister JR, Stone HA 1996. Pinching threads, singularities and the number 0.0304…. Phys. Fluids 8, 2827 (doi:10.1063/1.869086)
27. Leppinen D, Lister JR 2003. Capillary pinch-off in inviscid fluids. Phys. Fluids 15, 568–612. (doi:10.1063/1.1537237)
28. Bernoff AJ, Witelski TP 2010. Stability and dynamics of self-similarity in evolution equations. J. Eng. Math. 66, 11–31. (doi:10.1007/s10665-009-9309-8)
29. Viesca RC. 2016. Stable and unstable development of an interfacial sliding instability. Phys. Rev. E 93, 060202(R) (doi:10.1103/PhysRevE.93.060202) [PubMed]
30. Ruina A. 1983. Slip instability and state variable friction laws. J. Geophys. Res. 88, 10 359–10 370. (doi:10.1029/JB088iB12p10359)
31. Dieterich JH. 1978. Time-dependent friction and the mechanics of stick-slip. Pure Appl. Geophys. 116, 790–806. (doi:10.1007/BF00876539)
32. Dieterich JH. 1979. Modeling of rock friction 1. Experimental results and constitutive equations. J. Geophys. Res. 84, 2161–2168. (doi:10.1029/JB084iB05p02161)
33. Ruina AL. 1980. Friction laws and instabilities: a quasistatic analysis of some dry frictional behavior. PhD thesis, Brown University.
34. Rice JR, Ruina AL 1983. Stability of steady frictional slipping. J. Appl. Mech. 50, 343–349. (doi:10.1115/1.3167042)
35. Rice JR. 1968. Mathematical analysis in the mechanics of fracture. In Fracture: an advanced treatise (ed. H Liebowitz), pp. 192–311. New York, NY: Academic Press, Inc.
36. Rubin AM, Ampuero J-P 2005. Earthquake nucleation on (aging) rate and state faults. J. Geophys. Res. 110, B11312 (doi:10.1029/2005JB003686)
37. Berrut J-P, Trefethen LN 2004. Barycentric lagrange interpolation. SIAM Rev. 46, 501–517. (doi:10.1137/S0036144502417715)
38. Strogatz S. 1994. Nonlinear dynamics and chaos. New York, NY: Westview Press.
39. Ruelle D, Takens F 1971. On the nature of turbulence. Commun. Math. Phys. 20, 167–192. (doi:10.1007/BF01646553)
40. Newhouse S, Ruelle D, Takens F 1978. Occurrence of strange axiom A attractors near quasiperiodic flow on Tm, m[greater, double equals] 3. Commun. Math. Phys. 64, 35–40. (doi:10.1007/BF01940759)
41. Hirschmann EW, Eardley DM 1997. Criticality and bifurcation in the gravitational collapse of a self-coupled scalar field. Phys. Rev. D 56, 4696–4705. (doi:10.1103/PhysRevD.56.4696)
42. Peng GG, Lister JR 2014. The initial transient and approach to self-similarity of a very viscous buoyant thermal. J. Fluid Mech. 744, 352–375. (doi:10.1017/jfm.2014.75)
43. King FW. 2009. Hilbert transforms. Encyclopedia of Mathematics and its Applications 124 Cambridge, UK: Cambridge University Press.
44. James RL, Weideman JAC 1992. Psuedospectral methods for the Benjamin-Ono equation. In Advances in computer methods for partial differential equations VII (eds R Vichnevetsky, D Knight, G Richter), pp. 371–377. International Association for Mathematics and Computers in Simulation.
45. Weideman JAC. 1995. Computing the Hilbert transform on the real line. Math. Comput. 64, 745–745. (doi:10.1090/S0025-5718-1995-1277773-8)
46. Kaneko Y, Ampuero J-P 2011. A mechanism for preseismic steady rupture fronts observed in laboratory experiments. Geophys. Res. Lett. 38, L21307 (doi:10.1029/2011GL049953)
47. Palmer AC, Rice JR 1973. Growth of slip surfaces in progressive failure of over-consolidated clay. Proc. R. Soc. Lond. A 332, 527–548. (doi:10.1098/rspa.1973.0040)
48. Puzrin AM, Germanovich LN 2005. The growth of shear bands in the catastrophic failure of soils. Proc. R. Soc. A 461, 1199–1228. (doi:10.1098/rspa.2004.1378)
49. McClung DM. 2009. Dry snow slab quasi-brittle fracture initiation and verification from field tests. J. Geophys. Res. 114, F01022 (doi:10.1029/2007JF000913)
50. Carlson JM, Langer JS 1989. Mechanical model of an earthquake fault. Phys. Rev. A 40, 6470–6484. (doi:10.1103/PhysRevA.40.6470) [PubMed]
51. Shaw BE. 1995. Frictional weakening and slip complexity in earthquake faults. J. Geophys. Res. 100, 18 239–18 251. (doi:10.1029/95JB01306)
52. Bouchbinder E, Brener EA, Barel I, Urbakh M 2011. Slow cracklike dynamics at the onset of frictional sliding. Phys. Rev. Lett. 107, 235 501–235 505. (doi:10.1103/PhysRevLett.107.235501) [PubMed]

Articles from Proceedings. Mathematical, Physical, and Engineering Sciences are provided here courtesy of The Royal Society