Search tips
Search criteria 


Logo of procaThe Royal Society PublishingProceedings AAboutBrowse by SubjectAlertsFree Trial
Proc Math Phys Eng Sci. 2017 July; 473(2203): 20160606.
Published online 2017 July 5. doi:  10.1098/rspa.2016.0606
PMCID: PMC5549563

A phase-plane analysis of localized frictional waves


Sliding frictional interfaces at a range of length scales are observed to generate travelling waves; these are considered relevant, for example, to both earthquake ground surface movements and the performance of mechanical brakes and dampers. We propose an explanation of the origins of these waves through the study of an idealized mechanical model: a thin elastic plate subject to uniform shear stress held in frictional contact with a rigid flat surface. We construct a nonlinear wave equation for the deformation of the plate, and couple it to a spinodal rate-and-state friction law which leads to a mathematically well-posed problem that is capable of capturing many effects not accessible in a Coulomb friction model. Our model sustains a rich variety of solutions, including periodic stick–slip wave trains, isolated slip and stick pulses, and detachment and attachment fronts. Analytical and numerical bifurcation analysis is used to show how these states are organized in a two-parameter state diagram. We discuss briefly the possible physical interpretation of each of these states, and remark also that our spinodal friction law, though more complicated than other classical rate-and-state laws, is required in order to capture the full richness of wave types.

Keywords: friction, rate-and-state, self-healing slip pulse, detachment front, global bifurcation

1. Introduction

Inhomogeneous frictional sliding between solid bodies is ubiquitous and characterized by multiple spatio-temporal scales, compelling and practical examples being earthquake mechanics [1] or brake squeal [2]. In a sense, friction may be seen as an additional fundamental force, an irreversible process by which kinetic energy is transferred from macroscales to the microscale and dissipated in the form of heat. In mechanical systems, it has been estimated that up to 1.4% of GDP could be saved by a better management of interfacial wear and frictional energy loses [3]. Not only a nuisance, friction is also beneficial in many situations, such as in mechanical dampers, brakes and clutches, and there are many works aimed at controlling friction from robotics to turbine blades [46].

Within engineering dynamics, it is common to use approximate friction laws such as those of Coulomb or Stribeck [7], which allow simple distinction between stick and slip, or more elaborate versions such as the LuGre model [8]. The validity domain of such simple models is at best questionable [9], and theoretical studies tend to be limited to either point contact or to where surfaces in contact are assumed to behave homogeneously, although see [1015] for notable exceptions. Nevertheless, such approximations have proved useful in situations where friction is regarded merely as a loss mechanism on an otherwise macroscopic motion. Using such an approach, quite sophisticated mathematical theories of non-smooth mechanics have developed [16,17], which can predict complex temporal behaviours at the macroscale (e.g. [18,19]).

Tribology, on the other hand, being the detailed study of the wear, friction, adhesion and lubrication processes in the vicinity of interacting surfaces in relative motion, is a much more subtle topic that requires understanding of thermo-chemo-mechanical processes that occur at the microscale, or even below. In general, temperature dependence, micro-fluidics, molecular structure and the precise arrangement of surface asperities can all play a role in understanding the mechanisms that occur inside the contact region (e.g. [2024]).

At a much larger length scale, theories of earthquake formation often rely on models of frictional contact between tectonic plates, see [1,2527] for reviews. One of the things these models try to understand is how travelling slip pulses can be generated after years of inactivity or creep [14,2831]. Moreover, the intricate nature and diversity of earthquake types which have been recorded over the past decade, including aseismic events, episodic tremors, slow and fast earthquakes [32,33], which together suggest a continuum of sliding modes, remain a challenge in terms of mathematical modelling and physical understanding. Finally, the question of how locked zones along creeping faults are created, the failure of which leads to classical earthquakes, is also very intriguing [34].

For example, Avouac and collaborators [35] made observations in the 2003 Chengkung earthquake that (i) ‘seismic slip occurred on a fault patch that had remained partially locked in the interseismic period, (ii) the seismic rupture propagated partially into a zone of shallow aseismic interseismic creep but failed to reach the surface, and (iii) the aseismic afterslip occurred around the area that ruptured seismically’. A possible interpretation of these observations could be in terms of a new frictional slip mode that we describe in the present paper as a stick pulse. A stick pulse corresponds to a narrow sticking zone travelling against a slipping background and which might be identifiable as the seed of an earthquake.

Another major issue in the wider field of friction modelling is that experimental characterization of dynamic interfacial friction forces is fraught with difficulty. Apart from the branch of nanotribology, most results are restricted to macroscale measurements, which are typically characterized by scatter of data and lack of repeatability. One approach used to gain a quantitative match with theory, is to use high-resolution imaging or computation of the asperity map of the surface, to understand its precise roughness characteristics [3638]. The predictive power of this approach, limited by the computing power available, is somehow questionable because asperities continually change with wear, temperature and loading conditions. A notable exception is the work of Fineberg and co-workers (e.g. [3941]), who use ingeneous optical techniques to measure the detailed patterns of slip and stick between two transparent plates being slowly sheared against each other via a constant tangential force. The results suggest the existence of both fast and slow fronts and slip pulses. Another interesting study is that of Behrendt et al. [42] which, although not experimental, uses high-resolution numerical simulations to study a two-dimensional elastic block in contact with a rigid surface using a local Stribeck-like friction law. Different small-scale stick, slip and lift-off waves are found to occur as the macroscale slipping velocity varies between 0.5 and 40 mm s−1.

The ethos of the present paper is to attempt to bridge the gap between the macroscopic non-smooth dynamical systems approach of friction modelling and the detailed tribological point of view. We are also motivated by phenomena that might be universal across time and force scales, from earthquakes interseismic periods of decades to squeal oscillations in the kilohertz range. To that end we study a canonical, dimensionless, problem at the next order of complexity from point or homogeneous regional contacting behaviour, namely the existence of travelling fronts and pulses of stick-like and slip-like behaviour between two homogeneous frictional surfaces.

In contrast to most studies, the main purpose of our article then is to consider a general theory allowing at once for the generation of propagating slip- and stick-pulse, together with propagating slip and stick fronts, in models for regional interfacial contact. We emphasize that this study shows that these localized slip patterns can exist over a continuum range of wavespeeds; the range of wavespeeds can be explicitly computed as a function of the driving stress. We propose that this theory sheds new light on the processes responsible for the large variability of earthquake duration and frequency spectrum [32,33].

Our analysis builds on previous work [43,44] where the full three-dimensional elastodynamic problem is reduced to that of the idealized situation of a long thin elastic plate sliding against a rigid substrate. In [43,44], detachment front solutions were computed and shown to result from the non-monotonic crossover of the steady-state friction characteristic from velocity-weakening to velocity-strengthening regimes (e.g. [4549]). Such failure modes were interpreted as slow slip events and argued to explain incipient sliding friction along spatially extended contact as experimentally monitored [39,40]. Going much further, using the classical travelling-wave reduction of partial differential equations (PDEs) posed on an unbounded domain, we shall give precise descriptions for these localized states in terms of different types of connecting orbits [50].

We finally note that a variety of solitary waves in the form of slip pulses or fronts were also numerically found within the one-dimensional continuum limit of the Burridge–Knoppoff model [51], either with velocity-dependent non-smooth friction laws [5256] or the Ruina rate-and-state law [31].

This paper is organized as follows. Section 2 introduces and justifies our choice of friction law, namely a non-monotonic rate-and-state law that we analysed in a number of different mechanical settings in previous work [48,57,58]. Such models are well-motivated experimentally, yet are analytically tractable and avoid the non-smoothness associated with strictly Coulomb-like models. Section 3 introduces the geometry and mechanics of the problem to be studied. It is shown how the equations of motion reduce through long-wave approximation to a one-dimensional nonlinear wave equation, and through non-dimensionlization we identify its key parameters. We perform linear stability analysis and show numerically that the system supports travelling waves. Such waves are described by a two-dimensional dynamical system having a slow–fast multiple time-scale structure that is analysed in the phase plane, through a combination of analytical and numerical bifurcation analysis. Section 4 presents a complete exploration of the existence regions of various kinds of response, including wavetrains, stick and slip pulses and fronts between regions of homogeneous stick or slip. Our results are summarized in terms of two-parameter diagrams involving dimensionless parameters that represent the speed of travelling waves, and the applied shear stress. Finally, §5 comments on the physical applicability of the work, shows how our extended rate-and-state model has the minimal complexity necessary to capture the phenomena at hand, and suggests avenues for future work.

2. Rate-and-state friction model

At the heart of this paper lies the use of a recently proposed model for rate-and-state friction that captured many experimental observed effects that are not captured by non-smooth Coulomb-like models. We give here only a brief motivation, and refer the reader to recent work by the first two authors, e.g. [58,59] and references therein, for further details and discussion.

Basic Coulomb-like models classify material contacts via a single constant, the static coefficient of friction μs which is the maximum ratio of tangential interfacial shear stress τ to normal stress σ that can sustain stick. Alternatively, the bodies are said to slip, with the ratio equal to some other kinetic coefficient of friction μk which may, in general, be a (non-monotonic) function of relative velocity v in general (in a so-called Stribeck’s Law [7]). However, engineers have long recognized the limitations of such an approximation. In the 1950s, Rabinowicz, Bristow and co-workers, see e.g. [6062] argued from experimental evidence that most (μ,v) curves have a finite positive slope for low v, with the stick state better being described as a slow creep and what is called μs being some quantity near the maximum of the curve. Rabinowicz [62,63] also highlighted dependence on past sliding history and introduced the idea of a critical slip distance before any effect of the previous slip rate history has faded away. Equivalently, these memory effects determine a persistence length over which friction preserves its static value before dropping to its kinetic value for impulsively accelerating frictional interfaces. In the absence of such memory effects, the empirical assumption of a time-dependent static friction coefficient and a velocity-dependent kinetic friction coefficient has been shown to fail to reproduce the basic physics of stick–slip oscillations [62].

In a seminal 1966 paper, Brace & Byerlee [64] later on proposed that, rather than by brittle fracture, earthquakes could arise from recurrent stick–slip instabilities along pre-existing fault planes [65]. This shift of paradigm has caused the study of rock friction to flourish, propelled in particular by the work of Dieterich [66] and Ruina [67,68], who introduced the so-called rate-and-state formulation for a frictional interface. In such a formulation, memory effects are encoded in an internal variable ϕ(t) that characterizes the state of the contact region and quantifies its resistance to slip, originally thought of as representing average lifetime of a population of interfacial contacts [66,69] even though different microphysical origins could actually be conjectured [59]. In the rate-and-state framework of friction, the time evolution of the interfacial state ϕ(t) is determined by an empirical evolution equation whose precise mathematical expression remains an open question. We note, however, that the use of a first-order kinetics for the interfacial state evolution is mathematically justified whenever ϕ is interpretable as an nth statistical moment of the friction force [70]. On pure mathematical grounds, it has also been argued that piecewise-smooth descriptions of friction such as Coulomb’s need to be smoothed out by introducing an evolution equation for some hidden variable in order to account for any departure from steady sliding and hysteretic effects [18]. Originally, ad hoc choices of state evolution laws, such as (5.1), were motivated by fitting the frictional stress relaxation observed in response to sudden changes in driving velocity within ‘slide-hold-slide’ experiments [66,68]. The physical and experimental foundations of rate-and-state models can found in references [69,7173].

Rate-and-state models of friction are thus phenomenological in essence and assume that the interfacial shear stress τ is determined by nonlinear equations of the form


where the interfacial slip rate v=d(δu)/dt corresponds to the time derivative of the slip jump δu:=[[u]] across the frictional interface, t* is the characteristic time scale over which the interfacial state relaxes to equilibrium, and σ represents the normal stress applied to the interface. Different models correspond to different choices for the functions F and G, models with several internal variables, say ϕi, being possible, at least conceptually [68].

Most realizations of (2.1) assume that F is simply proportional to the normal stress, as in Coulomb friction, so that F=μ(v,ϕ)σ with


where V* and ϕ*:=L/V*=t* are reference values of the slip rate and the interfacial state so that the value of the kinetic friction coefficient μ takes the value μ=μ* in steady sliding at v=V*. The dimensionless material parameters a and b, which are of the order of 10−2 for most materials, can be fitted to experimental measurements (e.g. [71,74]) and represent the amplitude of the frictional response to sudden velocity variations [75]. A characteristic memory length L, which is connected to Rabinowitz’s persistence length [58], is usually introduced in place of the time scale t* in connection with the definition of (2.1)2, as we show below. Note meanwhile that the phenomenological analytical form (2.2) is well justified both experimentally and theoretically from physical first principles, see [59] and references therein.

In the absence of any ϕ dynamics, e.g. in steady sliding for which v:=V , the interfacial state must reach an equilibrium ϕss so that G(V,ϕss,σ)=0, which in combination with (2.2), yields the steady-state friction curve μss(V):=μ[V,ϕss(V)].

In particular, ‘regularized generalizations’ of the ageing law, defined by (2.2) with (5.1), can be built with expressions of the form


where c is an additional material constant describing the residual strength of the interface at high slip velocities, which requires fitting. tϕ(v) and ϕss(v) are (dimensionless) functions of the sliding velocity which describe the interfacial time scale and steady interfacial state. As in our previous work [58], we propose a ‘spinodal’ version of the friction law, taking


where R[dbl greater-than sign]1 is the ratio of the much slower time scale t** to t*; t** is the characteristic time scale for relaxation of the interfacial state variable in stationary contact, allowing for interfacial state saturation [76]. With (2.4), the steady-state friction curve μss(V) is non-monotonic and has a local maximum and a minimum located at


We remark that models (2.3) and (2.4) are regularized both mathematically in the sense that the logarithmic singularity of (2.2) at v0 is avoided, and physically in the sense that unbounded weakening at moderate sliding speed (i.e. greater than or equal to 1 mm s−1) is prevented. In addition, if the sinh1 expression for μ can be justified from microphysical theories of friction featuring thermally activated Eyring rate processes (see [59] and references therein), the sinh expression for G is somehow arbitrary, however, physically justified in ensuring a strong interfacial healing that prevents unphysical and unbounded acceleration of the interface in quasi-stationary contact (see [58] for more details).

Note that, compared to simple Coulomb-like models, there are many rate-and-state models and they often contain many ad hoc parameters. See, for example, the recent studies [74,77,78] that discuss how results from a pin-on-disc experiment allow the determination of parameters and discrimination between different rate-and-state models. Indeed, in [74] it was shown that the extra complexity associated with some rate-and-state models was necessary in order to replicate the qualitative nature of the frequency response curves.

In §5a, we return to the question of choice of friction law and we show how simpler monotonic or unregularized models do not capture the rich diversity of solution types that have been observed experimentally, and which are produced with the above friction law defined by (2.3) and (2.4). Alternative non-monotonic rate-and-state models can be found in [43,44,48].

3. A model problem for a thin sliding slab

(a) Formulation

Consider an infinite elastic plate of thickness h, density ρ, Young’s modulus E, shear modulus S and Poisson’s ratio ν, that is subject to a spatially uniform constant normal stress σ¯ and is driven from the top by a constant shear stress τ¯:=μ¯σ¯. The bottom of the plate is assumed to slide undeformed but with friction on a flat and horizontal rigid foundation (figure 1). Provided the wavelength of the elastic longitudinal wave propagating in the plate is large compared with its thickness, Lamb [79] showed that the distribution of the longitudinal stress and displacement components σxx and u is uniform across the plate’s cross-section. Following [80], the plate equation of motion can then be derived from a consideration of force balance in a cross-section of infinitesimal width δx and of unit length in the transverse direction (figure 1), yielding


Furthermore, assuming uniformity of the vertical component of the stress across the plate, i.e. σyyσ¯, Hooke’s Law implies that the vertical strain component v,y=(λu,x+σ¯)/(λ+2S). Hence the longitudinal stress σxx=(λ+2S)u,xv,y reads


The first and second Lamé coefficients are denoted λ and S (i.e. the shear modulus). Combining (3.2) with (3.1), the dimensional equation of motion of the plate becomes


where the longitudinal wave speed cl is defined by cl2=E¯/ρ where E¯:=E/(1ν2).

Figure 1.
Definition sketch of the elastic plate and forces balance on a cross-section of infinitesimal width δx. (Online version in colour.)

There are, in fact, alternative ways of deriving the slab equation (3.3). For example, we can use a ‘shallow-layer’ approximation to the full three-dimensional elasto-dynamic equations in which the slab displacement corresponds to the vertical averaging u(x,t):=0hu~(x,y,t)dy [43]. Alternatively, an asymptotic thin plate theory has been developed which shows that (3.3) arises at leading order. The details will appear elsewhere.

The model (3.3) is completed by using a friction law of the form (2.1) to express τ in terms of the sliding velocity u,t and the state variable ϕ. By using the scales h/cl, h, V*h/cl and E¯(V/cl) for non-dimensionalizing time, space, displacement and stress, respectively, the rate-and-state nonlinear dynamics of the plate is then governed by the dimensionless system


where we use the regularized spinodal form so that F and G are given by (2.3) and (2.4). The two dimensionless parameters appearing in (3.4) are


The parameter ζ compares the wave impedance to the friction impedance, which estimates the order of magnitude (for μ≈1) of the kinetic friction force that resists sliding at speed V*. Alternatively, ζ can be thought of as the reciprocal of dimensionless normal pressure measured in units of the stress scale E¯V/cl. The parameter r compares the characteristic time scale of information propagation carried by the slab elastic waves, over a distance of order h, to the characteristic rejuvenation time scale of the interfacial state (of order 1 for most materials). In practice, we expect both parameters to be small: ζ,r[double less-than sign]1.

(b) Uniform sliding states and their stability

Spatially uniform and steady sliding states (u,t(x,t),ϕ(x,t))=(v0,ϕ0) of the plate correspond to solutions of the nonlinear system


and are fully determined by the value of the dimensionless driving stress μ¯, independently of ζ and r. Using (2.3) we can rewrite (3.6) as


The number of solutions to (3.7) depends on the shape of the steady-state friction curve μss(v). Unlike the case of a monotonic friction law, which allows for a unique uniform steady-state sliding, the case of the spinodal law (2.3) and (2.4) allows for three possible uniform steady sliding solutions whenever μMμ¯μm. We interpret these three steady states as representing ‘stick’, ‘creep’ and ‘slip’, in order of increasing steady-state velocity. Typical orders of magnitude for these three states are 10 nm s−1, 1 μm s−1 and 10 mm s−1, referring to the steady-state velocity curve sketched in figure 10.

Figure 10.
Numerical continuation of wavetrain solutions in γ for fixed μ¯ (a,b) and μ¯ for fixed V (c,d), with the maximum and minimum amplitude super-imposed on the μss versus v curve. In a,c (respectively, b,d) ...

It is straightforward to show that only the stick and slip states, which lie on the velocity-strengthening branches of the friction curve, are stable to spatially homogeneous perturbations. In general, to examine the possibility of wave-like instabilities we consider the behaviour of infinitesimal perturbations to a uniform state, writing


where kR+ is a spatial wavenumber and p=s−iω gives the corresponding growth rate. The dispersion relation is as usual given by the determinant of the linearization of (3.4):


where, similar to before, μ,v denotes the partial derivative of the function μ with respect to u,t evaluated at (v0,ϕ0). The dispersion relation takes the form


where the slope of the friction curve μss′:=dμss(v0)/dv0 appears by the chain rule.

In the spatially homogeneous case k=0, we obtain the quadratic equation ζp2+(rζG,ϕ+μ,v)p+rG,ϕμss′=0, where r,ζ[double less-than sign]1, which can be rewritten more clearly as


so that the leading root can be estimated to be p/r=−G,ϕμss′/μ,v+O().

Neutral stability curves can be obtained by investigating the two cases p=0 and p complex but with a zero real part. In the first case, we observe that p=0 occurs only when k=0, so that neutrally stable modes having k≠0 are always oscillatory. Therefore, we look for modes for which p=−iω, ωR+. The real and imaginary parts of the cubic (3.8) then lead to


Combining these two expressions for ω2 yields the wavenumber kc and frequency ωc for neutrally stable modes:


Applying the Routh–Hurwitz criterion, we conclude that steady-state sliding is stable (i.e. R(p)<0 for all k) to all perturbations with wavenumbers k for which k2>kc2. In summary, uniform states on the velocity-strengthening branches of μss(v) are stable to all wavenumbers, whereas uniform states on the velocity-weakening branch are unstable to sufficiently long-wavelength perturbations for which k<kc. As Ruina [67] remarks, ‘the growth of instabilities is, then, insensitive to small wavelength (stiff) perturbations and very sensitive to long wavelength (soft) perturbations’ (see also [81,82]).

Figures Figures22 and and33 illustrate these results using the material parameter values listed in table 1 and the following dimensionless parameter values:


which correspond to material that has a thickness of the order of a few millimetre that is subjected to a normal pressure of the order of 10 Pa.

Figure 2.
Dispersion curves for the spinodal law (2.3): (a) growth rates of perturbations with wavenumber k; (b) temporal frequencies of perturbations with wavenumber k. Parameters: v0=10, r=10−5, ζ=0.1.
Figure 3.
(a) Stability region of perturbation wavenumbers k as a function of slip rate v0 as defined by (3.11). (b) Neutral angular frequency ωc as defined in (3.11). (c) The corresponding neutral wavelength λc:=2π/kc. (d) The neutral phase ...
Table 1.
Material parameter values used in the friction laws (2.3) and (2.4). These parameter values fit the experimental data obtained for the frictional properties of Bristol paper board [46] and paper sheet [83].

The linear stability analysis above is useful also to estimate the validity of the long wavelength approximation of the wave equation (3.3). It is clear from figure 3 that the critical wavenumber for the monotonic Dieterich Law (cf. §5a), kc2=(r/ζ)(ba)(1+rζv2/a), provides a good upper bound for the instability domain. Hence the long-wave approximation, λc[dbl greater-than sign]1, remains valid whenever r/ζ[double less-than sign]4π2/(ba), i.e. a confining pressure σ¯[4π2/(ba)](L/h)E¯1GPa, which is the case for both engineering systems and geophysical contexts. Equivalently, the thickness of the layer must satisfy h[4π2/(ba)](E¯/σ¯)L106/σ¯. For typical values of the confining pressures σ¯ ranging from a few GPa up to of the order of MPa, the latter expression suggests that the thin slab approximation remains valid for thicknesses ranging from a few millimetre up to of the order of metre.

(c) Multiple time-scale analysis of travelling waves

In addition to uniform steady states, the differential equations (3.4) are expected to support travelling waves. In this section, we investigate the existence of travelling waves and their generation in Hopf bifurcations from the uniform steady-state solutions. Without loss of generality, we assume that waves travel to the left, so that we introduce a travelling wave co-ordinate z=r(t+x/V) with travelling wave velocity V >0. System (3.4) then reduces to a planar system of nonlinear ordinary differential equations for the unknown functions v(z)=ru,t and ϕ(z):


where a dot now denotes d/dz. We find that the parameters V , ζ and r combine naturally into a single dimensionless parameter γ:


which is a monotonically decreasing function of the travelling wave speed V . Note that the regime γ≥0 corresponds the interval 0≤V ≤1 of subsonic wave speeds (note the longitudinal wave speed cl is scaled to unity in the non-dimensionalization leading to (3.4), which means that V is measured in units of cl).

System (3.13) has a natural slow–fast asymptotic structure because r,ζ[double less-than sign]1 and hence γ[double less-than sign]1, provided V is considered to be fixed and non-zero. Alternatively, we can consider the limit γ[double less-than sign]1 to be equivalent to the limit V1, that is the limit of fast travelling waves. This naturally defines z as a slow ‘time scale’ so that v evolves quickly until μ(v,ϕ) is close to μ¯ while ϕ evolves slowly. In the limit γ0 the slow variable ϕ is governed by the slow subsystem or reduced problem


For small γ, the full system (3.13) evolves slowly provided it stays within a small neighbourhood of the critical manifold defined by


which couples the evolution of v to that of the state variable ϕ: inverting (3.16) defines the relation v=v(ϕ;μ¯).

By contrast, if we rescale the equations on the fast time scale z^=z/γ we obtain the fast system


The natural approximation of this fast system is then to consider that the evolution of the fast variable v, when ϕ is constant, is given by the fast subsystem or layer problem


The one-dimensional dynamics along the critical manifold of the slow subsystem (3.15), governed by


is of importance as it provides a first insight into the solution types of the full system (3.13). For a spinodal friction model such as ours, g(ϕ;μ¯) is of course a non-monotonic function of ϕ whenever μMμ¯μm and three equilibria describing homogeneous sliding exist. As shown in figure 4a, the changes of sign of ϕ˙ indicate that the creep equilibrium point is unstable, whereas the stick and slip equilibria are stable. In terms of travelling wave solutions, this means that travelling fronts connecting the creep equilibrium to either the stick or slip equilibria can exist over open intervals in both the driving stress μMμ¯μm, and the wave speed V , whose set of values cannot be determined at this level of analysis. The connection between the creep and slip (respectively, stick) states represents a sharp jump between the corresponding interfacial slip rates, which can be interpreted as a generic detachment front (respectively, attachment front), see figure 4c. From this computation, a rough estimate of the front width is δx~δzV/r[dbl greater-than sign]1, whereas δz=O(1) and r[double less-than sign]1, which agrees with long wavelength approximation of this slab model. Physically, the dynamics of these fronts, which can travel fast with V =O(1), is determined by the interfacial state evolution law in such a way that the interfacial slip rate is slaved to the interfacial state so that the friction stress remains uniform and equal to that of the driving.

Figure 4.
The one-dimensional qualitative dynamics of the slow subsystems (3.15) and (3.22)2 for the fast fronts (ac) and slow fronts (bd) (moving to the left). The solid and dashed lines depict the fixed point connection type corresponding to ...

An alternative distinguished limit exists for the case of very slow travelling waves, where Vrζ, so that instead of being small, γ[dbl greater-than sign]1. We define ϵ:=γ−1[double less-than sign]1, and then define a new slow time variable z~ with z=z~/ϵ, so that the fast and slow systems, respectively, are now as follows:


where now v represents the slow variable, ϕ the fast variable and the critical manifold being defined by


The fast and slow systems (3.20) are analogous to the slow and fast systems (3.13) and (3.17), and the corresponding layer and reduced problems then read


Similarly, the slow dynamics along the critical manifold (3.21) governed by (3.22)2 allows for the existence of generic slow detachment and attachment fronts, but now, connecting either the stick or slip states to the creep state (figure 4bd). These slow fronts travel with V=o(rζ) and hence are characterized by an interfacial state which has time to relax to its steady-state value, becoming in turn enslaved to the evolution of the interfacial slip rate controlled by inertia. In contrast to the fast fronts, these slow fronts are consequently associated with a pulse shape friction stress spatial distribution whose width is estimated to be δxδz~ζ/r1.

We finally note that these different scalings of the problem indicate that dynamics is likely to be very ‘stiff’, with trajectories squeezed together in regions of phase space, and occasionally very rapid changes in the shape of typical orbits as parameters vary. Such phenomena could perhaps be fruitfully analytically tackled using Fenichel’s geometric singular perturbation theory [84]; we leave this to be the subject of future work. In the following sections, we present a detailed phase-plane and bifurcation analysis which shows that the full system (3.13) can exhibit a rich variety of solution types between the two asymptotic regimes of travelling fronts described here.

4. Detailed phase-plane analysis

We now study the bifurcation structure of travelling wave solutions to the system (3.13) using a combination of analytical and numerical methods. All numerical computations were carried out using the continuation software AUTO [85,86]. Except where otherwise stated, we use the spinodal law (2.3)–(2.4) with the material parameters as given in (3.12). We investigate the solution structure using the shear stress μ¯ and the travelling wave speed V as bifurcation parameters.

(a) Overall bifurcation structure

Figure 5 illustrates different kinds of spatially localized solution to the travelling wave system (3.13). The left-hand panels of each part of figure 5 show profiles of three different kinds of localized wave: (i) slip pulses, (ii) stick pulses, and (iii) detachment fronts. The interpretation and properties of these spatially localized solutions are described later. The right-hand panels in each part show how these waves appear as trajectories in the phase plane. In this and subsequent figures, the critical manifold (3.16) is indicated by a dashed line which coincides with the v-nullcline, whereas the ϕ-nullcline is indicated by a solid line. Equilibria (corresponding to uniform sliding states) are indicated by the solid black and white circles.

Figure 5.
Examples of three different kinds of localized solution trajectory. (a) A slip pulse on a background of homogeneous stick (red curve containing a long section near zero), together with single periods of periodic slip waves of large amplitude and long ...

Figure 6 depicts the two-parameter bifurcation diagram of the system, computed numerically using AUTO, showing curves of codimension-one bifurcations as lines in the (γ,μ¯) and (V,μ¯) planes: these are equivalent representations due to relation (3.14). Figure 7a is a topologically equivalent version of the bifurcation diagrams shown in figure 6 that links the numerical bifurcation results to the theoretical analysis summarized in figure 7b. These two figures therefore are the core of our results and illustrate how the numerical and analytic investigations together enable us to provide a complete summary of the bifurcations associated with the existence of uniform states and travelling waves (figures (figures88 and and99).

Figure 6.
Fundamental phase diagram in the form of a two-parameter bifurcation diagram of the travelling wave system (3.13) under the spinodal rate-and-state friction model (2.4). As (a) γ versus applied shear stress μ¯ and (b) wave speed ...
Figure 7.
(a) Qualitative version of the bifurcation diagram in figure 6, unfolding the various bifurcation curves that appear to overlie each other due to the slow–fast nature of the true diagram. Note that the curves separating regions iii and iv, and ...
Figure 8.
Main examples of qualitative phase portraits and waveforms corresponding to the saddle–saddle connections between the equilibrium points of (3.13). We denote their corresponding uniform slab slip rates by (vs,vc,vf), solutions of μss( ...
Figure 9.
Qualitative phase portraits and sketch waveforms corresponding to the generic connections between the equilibrium points of (3.13). Acronyms are explained in figure 6.

Uniform states exist for all values of γ>0 and shear stress μ¯. For sufficiently high shear stress μ¯>μmax=μss(vmax)0.4 or sufficiently low shear stress μ¯<μmin=μss(vmin)0.3, there is a unique uniform state. Saddle-node bifurcations occur on the lines μ¯=μmax and μ¯=μmin associated with the local extrema of the friction characteristics.

For μmin<μ¯<μmax, there are therefore three distinct uniform states due to the spinodal nature of the friction model. The middle state on the velocity-weakening branch of the steady-state friction curve, i.e. the ‘creep’ sliding state, can become unstable along a Hopf bifurcation curve, which is precisely the neutral stability line that was computed in §3b (figure 3) as presented in §4b. The bifurcation produces a periodic travelling wave train that exists for values of γ below, or equivalently wave speeds V above, this curve (marked H in figures figures66 and and77a). With increasing V or decreasing γ, the branch of periodic orbits is destroyed in one of two kinds of homoclinic bifurcation (blue lines in figures figures66 and and7),7), leading to either localized slip or stick pulses that in phase space are homoclinic orbits to the slip or to the stick equilibrium points, respectively. See also the phase portraits in figure 5 and bifurcation diagrams in figure 10 for examples of such homoclinic orbits.

The Hopf bifurcation curve (H) terminates in Bogdanov–Takens (BT) bifurcation points where it meets one of the two saddle node lines, the location of which is indicated by a hollow circle in the figures. A standard unfolding of such a point (e.g. [87]) shows that the curves of homoclinic orbits must also originate from each of these points. In the interior of the spinodal region, the two homoclinic orbits meet at a codimension-two global bifurcation point (P), represented by a solid circle, where a complete heteroclinic cycle exists between the stick and the slip equilibria. An unfolding of that point is presented in §4e; although this kind of analysis is, in general, very well known, we could not find a direct reference in the literature to this specific case. In the following subsections, we comment in turn on each of the three main kinds of motion: wave trains of periodic travelling waves, pulses and fronts.

(b) Periodic orbits: travelling waves

Spatially periodic wave trains arise from Hopf bifurcation points located on the velocity weakening branch of the steady-state friction curve. Study of the roots of the Jacobian matrix of (3.13) shows that the bifurcation point occurs at a critical value μ¯h of the shear stress μ¯ that provides the external driving, where μ¯h solves the equation


which defines the locus (H) of the Hopf bifurcation points shown in figures figures66 and and7).7). Note that (4.1) naturally follows from the neutral stability curve (3.11)1 derived in §3b after substituting the definition V :=ωc/kc of the phase velocity of a linear wave. A good approximation to the critical stress is given by μ¯hμ+(ab)ln(vh), in which the critical slab slip rate vh=a/γ follows from (4.1) using the classical Dieterich–Ruina Laws (5.1).

To go beyond, numerical continuation of the branch of periodic orbits indicates that there is typically one of two kinds of behaviour observed on varying μ¯ for different wave speeds V . These behaviours are illustrated in figure 10.

For sufficiently low wave speeds V , the behaviour of the bifurcating branch is qualitatively always as illustrated in figure 10a. The Hopf bifurcation gives birth to a branch of small-amplitude limit cycles that leads to a canard explosion (e.g. [88]) over an exponentially small range of the parameter μ¯μ¯h. The limit cycles are known as canard solutions of (3.13) and their amplitude and period increase significantly over this small range of parameter values as the limit cycle follows close to a segment of the slow manifold.

A typical canard cycle can be divided into three consecutive phases, as a consequence of the slow–fast asymptotic structure of (3.13). The fast phase defines the front of the velocity spike and corresponds to the sharp acceleration of the slab material points, whose dynamics is governed by (3.18) to leading order. Such an acceleration is accompanied by a decrease in the value of state variable (dynamic rejuvenation) since ϕ[dbl greater-than sign]ϕss(v) until a maximum slip rate is reached. Subsequently, the slip rate decreases throughout an intermediate phase along the ϕ-nullcline during which the interfacial state is forced to remain close to its equilibrium value ϕss(v). As a result, the deceleration is associated with an increase of ϕ leading to the slow phase of the cycle along the critical manifold (the v-nullcline). The dominant contribution to the period of the limit cycle is from this final phase, ϕ, because ϕ[double less-than sign]ϕss(v). Eventually, the branch of relaxation oscillations terminates in a homoclinic bifurcation where the limit cycle collides with the low-velocity stick equilibrium point. For higher V , as in figure 10b the branch of limit cycles terminates in a homoclinic bifurcation in which it collides with the high-velocity slip equilibrium point, before any relaxation oscillation can occur.

(c) Homoclinic orbits: travelling localized pulses

Homoclinic orbits, such as those that exist on the curves of homoclinic bifurcations shown in the (γ,μ¯) plane in figure 6, correspond to localized travelling pulses in the original (x,t) coordinates [89].

The orbit homoclinic to the equilibrium point with the lower value of V corresponds to a slip pulse: the interface creeps slowly, at the rate given by that of the ‘stick’ equilibrium, except in a short spatial region where the interface sees a velocity spike. The homoclinic orbit for moderate to large values of V corresponds to a stick pulse. Here, the whole interface slides at the rate associated with the high velocity saddle point while a localized patch of creeping material points is travelling along the interface.

(d) Heteroclinic orbits: fronts

A different kind of localized travelling solution is a front solution that describes a transition in space between low- and high-speed regions of slip. In the phase plane for our ODE system, such a solution corresponds to a heteroclinic orbit, asymptotic to two different equilibrium points as z tends to ± [89]. These particular solutions correspond to heteroclinic orbits which connect two different saddle points. For the model (2.3) and (2.4), there are two front solutions of particular interest.

The first of these correspond to a heteroclinic orbit formed by the intersection of the unstable manifold of the low velocity saddle point and the stable manifold of the high velocity saddle point. We interpret this ‘upper’ heteroclinic orbit as a detachment front since it describes a division of the slab into a part behind the front which is almost at rest while far ahead of the front the slab is moving rapidly (figure 5c).

The second (lower) heteroclinic orbit is the reverse: and intersection between the unstable manifold of the high velocity saddle point and the stable manifold of the low velocity saddle point. The rear of the slab is now sliding rapidly while ahead of the front the slab is moving rapidly. We interpret such a front as an attachment front.

We remark that numerically it is difficult to continue paths of some of the homoclinic and heteroclinic orbits, due to numerical stiffness caused by the slow–fast nature of the system. Nevertheless, we have confirmed the location of these curves as depicted in figure 6 through a series of one-parameter continuation runs of periodic orbits and noting points of divergence to large period.

(e) The codimension-two point P

It is interesting to observe in figure 6 how all four bifurcation branches hetd, heta, homsl and homst emerge from the organizing centre P, marked by a solid red dot, which represents a codimension-two heteroclinic cycle. We analyse the dynamics near P by constructing return maps describing trajectories in neighbourhoods of the two-saddle points involved in the heteroclinic cycle. This methodology and general notation is standard, going back to Poincaré and Shil’nikov; examples of this kind of calculation can be found in [90,91]. Interestingly, despite its simplicity we are unaware of a reference to this specific calculation in the literature.

Consider a planar system having two hyperbolic saddle points p1 and p2 with local coordinates aligned along eigenvector directions (x1,y1) and (x2,y2), respectively. Let the eigenvalues be −c1,e1 and e2,−c2, respectively, where c1,2,e1,2>0 denote the contracting and expanding eigenvalues at the saddle points. We further assume that at parameter values (α,β)=(0,0) there exists a heteroclinic cycle composed of two connecting orbits between the saddle points: one heteroclinic connecting orbit leaves p1 tangent to the y1>0 direction and approaches p2 tangent to the y2>0 direction. The other heteroclinic orbit leaves p2 tangent to the x2<0 direction and approaches p1 tangent to the x1>0 direction.

In terms of our travelling wave ODE system, the saddle point p1 corresponds to the low-velocity (stick) equilibrium, whereas p2 corresponds to the high velocity (slip) equilibrium point. The analysis proceeds by constructing leading-order approximations to the dynamics in two regimes: within small boxes |xj|,|yj|≤d near each saddle point, where the flow can be assumed to be linear (after a suitable coordinate transformation, due to hyperbolicity), and near each unstable manifold where we use the underlying differentiability of the vector field to make a Taylor series expansion. These two regimes lead to ‘local’ and ‘global’ maps between boundaries of these small boxes.

Accordingly, we define the sections (box sides) Σ1:={x1=d}, Σ2:={y1=d}, Σ3:={y2=d} and Σ4:={x2=−d}.

Consider first the local map Π1:Σ1Σ2 near p1, where the flow is well approximated by the linear flow x˙1=c1x1, y˙1=e1y1. Consider the trajectory starting from (x1,y1)=(d,y^0) and arriving at (x1,y1)=(x^1,d) at time T1>0. Integrating the ODEs we obtain x1(t)=d ec1t, y1(t)=y^0ee1t. At time t=T1, we can compute that

T1=1e1ln(dy^0)and hencex^1=d(dy^0)c1/e1.

Similarly, consider the local map near p2 given by Π3:Σ3Σ4 for which x˙2=e2x2, y˙2=c2y2. We take the initial condition to be (x2,y2)=(x^2,d), so that x2(t)=x^2ee2t, y1(t)=d ec2t. At time t=T2, we define the trajectory to hit the point (x2,y2)=(d,y^2) on Σ4, so that

T2=1e2ln(dx^2)and hencey^2=d(dx^2)c2/e2.

We now turn to the global maps Π2 and Π4 between neighbourhoods of the saddle points. For the heteroclinic connection from p1 to p2, we now consider the map Π2:Σ2Σ3 which takes a point (x^1,d)(x^2,d) (where x^2<0). For 0<x^11 this takes the form


where the constant A>0 because the flow is two-dimensional, and the global bifurcation parameter α parametrizes the unfolding of the heteroclinic connection from p1 to p2. Similarly, to examine trajectories near to the heteroclinic connection from p2 to p1 we form the map Π4:Σ4Σ1 which takes a point (d,y^2)(d,y^3) such that


where B>0 is a constant and β is the bifurcation parameter. We now compose the maps to obtain Π:=Π4Π3Π2Π1:Σ1Σ1 which takes the form



where B~ is a rescaled constant. After rescaling to set as many inessential positive constants as possible to unity this is a one-dimensional map of the form


and the new variable xn is a new variable. Alternatively, we could write the map in two parts, in the form


In the usual way, fixed point of this map indicate periodic, homoclinic and heteroclinic orbits for the ODEs. The conditions for the existence of various orbits (and their interpretation in the ODE dynamics) are therefore as follows:

  • (i) there exists a periodic orbit (a periodic wavetrain) if y=−α+xδ1>0 and x=−β+ayδ2>0;
  • (ii) there exists an orbit homoclinic to p1 (a slip pulse) if x=0 and y>0, i.e. 0=−β+a(−α)δ2 and α<0;
  • (iii) there exists an orbit homoclinic to p2 (a stick pulse) if y=0 and x>0, i.e. 0=−α+(−β)δ1 and β<0;
  • (iv) there exists a heteroclinic orbit from p1 to p2 (a detachment front) when α=0; and
  • (v) there exists a heteroclinic orbit from p2 to p1 (an attachment front) when β=0.

Figure 7b summarizes the curves and regions on which these orbit types exist; this figure then explains the qualitative nature of the ‘state diagram’ of slip and stick pulses and detachment and attachment fronts, near the codimension-two point P, as indicated in figure 6.

5. Discussion

We first provide a few comments on the necessity of our spinodal law, contrasting it with the use of monotonic or unregularized friction laws. We then make a few tentative physical remarks and suggest avenues for future work, before concluding by summarizing the results of the paper in a slighter wider context.

(a) Results for simpler friction laws

A natural question to ask is whether all the ingredients of the spinodal friction law we chose are necessary to obtain the rich bifurcation structure of frictional waves we have found. For example, within rate-and-state friction theory, there is no a priori obvious choice for the state evolution function G, see [68,71]. The two most studied laws are known as the Dieterich ageing law and the Ruina slip law defined, respectively, by (with dimension)


Regularized generalizations of the Dieterich–Ruina models such as the one introduced in §2 were proposed in [58] and shown to be consistent with many more physical phenomena.

We have repeated the analysis in this paper with simpler friction laws (results not shown). Note that for the Dieterich and Ruina laws, the shape of the friction laws μss(v) is monotonic so that only one equilibrium point of (3.13) is possible, corresponding to the creep uniform sliding state. Hence the planar slow–fast structure of (3.13) indicates that such a model allows only for wavetrains that have a form close to slip pulses. This solution type lies at the heart of many previous studies of the slip pattern formation along spatially homogeneous frictional interfaces (e.g. [31,82,92]).

Alternatively, if we introduce a spinodal version of the friction law without the hyperbolic regularization of the state evolution law, then we introduce both uniform stick and slip equilibria. However, both states are only able to support singular canard-type wavetrains leading to homoclinic orbits either of slip pulse or stick pulse type. Without the hyperbolic regularization of the state evolution law, the stiffness of the slow–fast dynamics is exacerbated and the wavetrain domain ii.WTm is reduced to an exponentially thin region along the Hopf bifurcation locus H. The saddle–saddle connections of homoclinic and heteroclinic orbits are also confined along the H line.

We conclude that our regularized spinodal law, used here, is the most parsimonious model able to capture the range of dynamical phenomena we have uncovered in a single two-parameter state diagram.

(b) Physical interpretation and future work

Before embarking on a discussion of physical interpretation, we should state at the outset that we have not yet considered the dynamics of the PDE system (3.4), nor established the stability of any of the wave structures we have constructed. Although the problem can be thought of as a generalized damped nonlinear wave equation, it has several non-trivial features, not least its slow–fast nature and that the state u is cyclic in the travelling-wave problem. We are unaware of any mathematical theory that can be applied directly to make generic statements about the dynamics and stability. We also need to consider boundary conditions carefully in order to describe a physical problem with either rigid or dead loading. Even straightforward time integration of the equations of motion requires careful treatment due to the numerical stiffness of the problem. Thus, any investigation into the dynamics and stability of the wave-train, pulse and front states is beyond the scope of the present study and will be left to future work.

Nevertheless, some preliminary remarks can be made. Firstly, the linear stability analysis of the uniform slip or stick states suggests that the Hopf bifurcations result in small amplitude waves that we expect to be stable on an infinite domain. Second, it would seem intuitive based on general theory (e.g. [93]), that one can indeed expect to find stable solutions to the nonlinear wave problem that are represented by connecting orbits between saddle points. However, such arguments are likely to be subtle owing to the multiple-time-scale nature of the nonlinearity in this model. An approach via geometric singular perturbation theory [84] may well prove feasible. We also point out that commonly used tools such as exponential dichotomies and Evans’ functions developed for reaction–diffusion systems (e.g. [9397]) should be most useful to establish linear stability results. Finally, in close relation to our problem, we note that such results for wave-train solutions in the one-dimensional sine-Gordon and Klein–Gordon equations have recently been published [98,99].

Physical quantities, principally characteristic wave speeds and length scales, and hence time scales, can be extracted from our main results as presented in the (V,μ¯) phase diagram in figures figures66 and and7.7. To indicate these as clearly as possible, we focus on the curves of heteroclinic, homoclinic and Hopf bifurcations from figure 6 corresponding to the various localized and periodic structures. These curves also give the stress–velocity relations V(μ¯) of practical importance for each type of travelling wave (figure 11a).

Figure 11.
Replotting of figure 6 in the (μ¯,V) plane to show the characteristic wave speeds and length scales associated with the wave train, slip and stick pulse solutions. (a) The (μ¯,V) plane as in figure 6. (b) Length scales ...

Although existing over a range of values for V , it is useful to attempt to distinguish qualitatively between slow and fast propagating waves. An estimate of the order of magnitude of the wave speed can be obtained from the location of the BT points lying at the intersection of the Hopf (H) and saddle-node (sn) bifurcation curves. The location of these BT points is determined by the extrema of the friction curve and given by (4.1) and (2.5), where we set V=rζ/(rζ+γ) and use the appropriate values of γ, i.e.


Slip pulses and attachment fronts can be slow localized waves, whereas all types of localized waves can be fast. This does not contradict recent experimental and theoretical work reporting on the slow propagation of detachment fronts (e.g. [39,40,43,44,100]) as the order of magnitude of the wave speed V can be substantially reduced while increasing the confining pressure σ¯, i.e. reducing ζ.

The characteristic length scales of periodic wave trains can be estimated from the neutral wavelength λc associated with the Hopf bifurcation. By contrast, we can compute numerical estimates of the length scale Δ associated with the two cases of slip and stick pulses by defining the integrals


where the first definition holds for slip pulses and the second for stick pulses, and we denote by δv the interfacial velocity jump, defined as δv:=max(v)min(v). Figure 11bc illustrates these length scale estimates based on the stress-velocity curves in figure 11a. We see that characteristic sizes of the slip and stick pulses deviate from the neutral wave train wavelength λc by a minimum of one order of magnitude. Note that the lines for Δst and Δsl end where the computations break down due to the numerical stiffness of the problem in AUTO [85,86]. Further, we note that the orders of magnitude of these length scales are compatible with the long-wavelength assumption underlying this work.

Finally, we comment on the friction stress distribution along the interface for the slip and stick pulses as shown in figure 11d. For the slip pulse case, the spatial variations of μ are localized within the pulse itself, which contrasts with crack-like models characterized by a stress singularity at the front of the pulse [29,30]. Our friction law truly generates a self-healing pulse as the friction stress in front of, and behind, the pulse zone is just the applied shear stress μ¯, the slip rate being equal to that of the corresponding ‘stick’ equilibrium. This is a consequence of the spinodal nature of the friction model we use. For the stick pulse, the friction stress distribution mimics a crack solution with a stress singularity at the front of the ‘stick’ zone, which is associated with a diffuse drop in the slip rate.

It is worth pointing out that the stress–velocity relation of the detachment fronts, i.e. hetd, is qualitatively equivalent to the V(μ¯) curve computed in [44] (figure 4). In their somewhat different non-monotonic friction model, its local minimum sets a lower bound, equivalent to our Vmin, for a rupture front propagation speed, which is consistent with experimental measurements in [40]. Our analysis also attributes a clear definition of Vmin as the front speed of the codimension-two non-central saddle-node heteroclinic bifurcation point (open squares). Similarly, the local maximum of our friction law, gives a maximum Vmax<1 for the front speed which is below the slab speed of sound.

Further physical consequences of these results for earthquakes mechanics and engineering systems remain to be explored and are left for future work. However, this work reveals that the quantitative features of the stress-velocity curves strongly depend on the mathematical details of the state evolution law, whose experimental identification and theoretical derivation from first principles are still open questions.

(c) Summary

A comprehensive understanding of the physical mechanisms that determine the diversity of frictional slip pattern formation along extended contacts between solids that has been reported over the years [40,42,82] is still lacking. This lack would appear to be at least in part because of incomplete physical and mathematical modelling of friction and how it couples with elastic wave radiation. A key phenomenon we have sought to explain is that such interfaces may possess, depending on the stress scales involved, either waves of slip on a background of uniform stick or creep, or waves of stick within a background of uniform slip. We have also explained that such waves can exist as isolated pulses, as periodic waves or as fronts between regions of homogeneous slip or stick.

Specifically, this work has shown how introducing a smooth interfacial friction model can explain the origin of slip pulses, stick pulses, travelling wavetrains, detachment fronts and reattachment fronts, all in the same mathematical formulation of regional contact and within the well-established theory of smooth dynamical systems theory. The key ingredient of the mathematical model is that the friction law should present non-monotonic velocity-dependent steady-state characteristics. Then, by assuming deformation in the form of a steady travelling wave, it is argued that slip and stick pulses can be understood in terms of homoclinic global bifurcations of travelling periodic slip patterns, as well as the travelling fronts as heteroclinic connections.

We emphasize that slip pulses are anchored at the equilibrium saddle point lying on the low velocity strengthening branch of the steady-state friction curve, while the existence of a high velocity strengthening branch in spinodal friction also allows the existence of ‘stick pulse’ which corresponds to a narrow travelling ‘stick’ zone. Along the bifurcated branch, travelling wave trains of slip pulses develop from a canard explosion, which can lead to relaxation oscillations. We note that the qualitative features of these patterns are independent of their propagation direction with respect to the slab motion. More broadly, we have shown how this plethora of behaviours is shown to be a consequence of the spinodal character of friction, by showing that simpler friction models are unable to reproduce this behaviour.

Finally, we stress that the localized pulse solutions exist only along specific lines in parameter space, giving stress–velocity relations and separating domains of generic travelling fronts and wavetrains of various types. This may contribute to a better understanding of why friction experiments are so notoriously difficult to perform in an experimentally reliable and repeatable way (e.g. [9]). For instance, our analysis may provide a possible origin for the scattering in experimental measurements of the detachment fronts’ speed [40], the experimental apparatus not only sampling proper heteroclinic detachment fronts but also the generic detachment fronts belonging to the adjacent regions of the heteroclinic bifurcation curve. In this respect, in any specific experiment, initial and boundary conditions also will play a key role in the selection of solutions; we will return to this rich avenue for investigation in future work.


The authors thank two anonymous reviewers for their thorough reading of our manuscript, comments and suggestions. The open access processing charge was supported by the University of Bristol.

Data accessibility

All data are contained within the published paper.

Authors' contributions

T.P. conceived and developed the model and performed the research. J.H.P.D. provided the analysis in §4e. All three authors contributed to the general overview and interpretation of the results, helped draft the manuscript and gave final approval for publication.

Competing interests

We declare we have no competing interests.


T.P. and A.R.C. acknowledge support from the UK EPSRC programme grant ‘Engineering Nonlinearity’ (ref. EP/K003836/1). J.H.P.D. acknowledges support for this work from the Royal Society through a University Research Fellowship.


1. Kanamori H, Brodsky EE 2004. The physics of earthquakes. Rep. Progress Phys. 67, 1429 (doi:10.1088/0034-4885/67/8/R03)
2. Butlin T, Woodhouse J 2011. A systematic experimental study of squeal initiation. J. Sound. Vib. 330, 5077–5095. (doi:10.1016/j.jsv.2011.05.018)
3. Jost HP. 1966. Lubrication (Tribology): Education and Research: A report on the present position and industry’s needs. London, UK: Department of Education and Science. HM Stationary Office.
4. Armstrong-Hélouvry B. 1991. Control of machines with friction. New York, Ny: Springer.
5. Armstrong-Hélouvry B, Dupont P, Canudas De Wit C 1994. A survey of models, analysis tools and compensation methods for the control of machines with friction. Automatica 30, 1083–1138. (doi:10.1016/0005-1098(94)90209-7)
6. Schwingshackl CW, Petrov EP, Ewins DJ 2012. Measured and estimated friction interface parameters in a nonlinear dynamic analysis. Mech. Syst. Signal. Process. 28, 574–584. (doi:10.1016/j.ymssp.2011.10.005)
7. Sheng G. 2008. Friction-induced vibrations and sound: principles and applications. Boca Raton, FL: CRC Press.
8. Canudas de Wit CH, Olsson H, Astrom KJ, Lischinsky P 1995. A new model for control of systems with friction. IEEE Trans. Autom. Control 40, 419–425. (doi:10.1109/9.376053)
9. Woodhouse J, Putelat T, McKay A 2015. Are there reliable constitutive laws for dynamic friction? Phil. Trans. R. Soc. A 373, 20140401 (doi:10.1098/rsta.2014.0401) [PubMed]
10. Adams GG. 1995. Self-excited oscillations of two elastic half-spaces sliding with a constant coefficient of friction. ASME. J. Appl. Mech. 62, 867–872. (doi:10.1115/1.2896013)
11. Adams GG. 1996. Self-excited oscillations in sliding with a constant friction coefficient–a simple model. ASME. J. Tribol. 118, 819–823. (doi:10.1115/1.2831614)
12. Adams GG. 1998. Steady sliding of two elastic half-spaces with friction reduction due to interface stick-slip. ASME. J. Appl. Mech. 65, 470–475. (doi:10.1115/1.2789077)
13. Adams GG. 2000. An intersonic slip pulse at a frictional interface between dissimilar materials. ASME. J. Appl. Mech. 68, 81–86. (doi:10.1115/1.1349119)
14. Adda-Bedia M, Ben Amar M 2003. Self-sustained slip pulses of finite size between dissimilar materials. J. Mech. Phys. Solids 51, 1849–1861. (doi:10.1016/S0022-5096(03)00068-1)
15. Scheibert J, Prevost A, Debrégeas G, Katzav E, Adda-Bedia M 2009. Stress field at a sliding frictional contact: experiments and calculations. J. Mech. Phys. Solids 57, 1921–1933. (doi:10.1016/j.jmps.2009.08.008)
16. Brogliato B. 1999. Nonsmooth mechanics, 2nd edn London: Springer.
17. Bernardo M, Budd C, Champneys AR, Kowalczyk P 2008. Piecewise-smooth dynamical systems: theory and applications. Applied mathematical sciences, vol. 163 Berlin, Germany: Springer.
18. Jeffrey MR. 2015. On the mathematical basis of solid friction. Nonlinear Dyn. 81, 1699–1716. (doi:10.1007/s11071-015-2100-7)
19. Champneys AR, Várkonyi PL 2016. The Painlevé paradox in contact mechanics. IMA J. Appl. Math. 81, 867–887. (doi:10.1093/imamat/hxw027)
20. Israelachvili JN. 1991. Intermolecular and surface forces. New York, NY: Academic Press.
21. Hutchings IM. 1992. Friction and wear of engineering materials. In Tribology: friction and wear of engineering materials. London, UK: Edward Arnold.
22. Bhushan B, Israelachvili JN, Landman U 1995. Nanotribology: friction, wear and lubrication at the atomic scale. Nature 374, 607–616. (doi:10.1038/374607a0)
23. Gao J, Luedtke WD, Gourdon D, Ruths M, Israelachvili J, Landman U 2004. Frictional forces and amontons’ law: from the molecular to the macroscopic scale. J. Phys. Chem. B 108, 3410–3425. (doi:10.1021/jp036362l)
24. Williams J. 2005. Engineering tribology. Cambridge, UK: Cambridge University Press.
25. Daub EG, Carlson JM 2010. Friction, fracture, and earthquakes. Annu. Rev. Condensed Matter Phys. 1, 397–418. (doi:10.1146/annurev-conmatphys-070909-104025)
26. Bizzarri A. 2014. The mechanics of seismic faulting: recent advances and open issues. Rivista del Nuovo Cimento 37, 181–271. (doi:10.1393/ncr/i2014-10099-0)
27. Kanamori H. 2014. The diversity of large earthquakes and its implications for hazard mitigation. Annu. Rev. Earth. Planet. Sci. 42, 7–26. (doi:10.1146/annurev-earth-060313-055034)
28. Brune JN. 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes. J. Geophys. Res. 75, 4997–5009. (doi:10.1029/JB075i026p04997)
29. Heaton TH. 1990. Evidence for and implications of self-healing pulses of slip in earthquake rupture. Phys. Earth Planet Int. 64, 1–20. (doi:10.1016/0031-9201(90)90002-F)
30. Perrin G, Rice JR, Zheng G 1995. Self-healing slip pulse on a frictional surface. J. Mech. Phys. Solids 43, 1461–1495. (doi:10.1016/0022-5096(95)00036-I)
31. Erickson B, Birnir B, Lavallée D 2011. Periodicity, chaos and localization in a Burridge-Knopoff model of an earthquake with rate-and-state friction. Geophys. J. Int. 187, 178–198. (doi:10.1111/j.1365-246X.2011.05123.x)
32. Peng Z, Gomberg J 2010. An integrated perspective of the continuum between earthquakes and slow-slip phenomena. Nat. Geosci. 3, 599–607. (doi:10.1038/ngeo940)
33. Ide S. 2014. Modeling fast and slow earthquakes at various scales. Proc. Jpn Acad. Ser. B Phys. Biol. Sci. 90, 259–277. (doi:10.2183/pjab.90.259) [PMC free article] [PubMed]
34. Kaneko Y, Avouac J-P, Lapusta N 2010. Towards inferring earthquake patterns from geodetic observations of interseismic coupling. Nat. Geosci. 3, 363–369. (doi:10.1038/ngeo843)
35. Thomas MY, Avouac J-P, Champenois J, Lee JC, Kuo LC 2014. Spatiotemporal evolution of seismic and aseismic slip on the Longitudinal Valley Fault, Taiwan. J. Geophys. Res. 119, 5114–5139. (doi:10.1002/2013JB010603)
36. Dieterich JH, Kilgore BD 1994. Direct observation of frictional contacts: new insights for state-dependent properties. Pure Appl. Geophys. 143, 283–302. (doi:10.1007/BF00874332)
37. Berthoud P, Baumberger T 1998. Shear stiffness of a solid-solid multicontact interface. Proc. R. Soc. Lond. A 454, 1615–1634. (doi:10.1098/rspa.1998.0223)
38. Yastrebov VA, Anciaux G, Molinari J-F 2015. From infinitesimal to full contact between rough surfaces: evolution of the contact area. Int. J. Solids Struct. 52, 83–102. (doi:10.1016/j.ijsolstr.2014.09.019)
39. Rubinstein SM, Cohen G, Fineberg J 2004. Detachment fronts and the onset of dynamic friction. Nature 430, 1005–1009. (doi:10.1038/nature02830) [PubMed]
40. Ben-David G, Cohen O, Fineberg J 2010. The dynamics of the onset of frictional slip. Science 330, 211–214. (doi:10.1126/science.1194777) [PubMed]
41. Shlomai H, Fineberg J 2016. The structure of slip-pulses and supershear ruptures driving slip in bimaterial friction. Nat. Commun. 7, 11787 (doi:10.1038/ncomms11787) [PMC free article] [PubMed]
42. Behrendt AC, Weiss J, Hoffmann NP 2011. A numerical study on stick-slip motion of a brake pad in steady sliding. J. Sound Vib. 330, 636–651. (doi:10.1016/j.jsv.2010.08.030)
43. Bouchbinder E, Brener EA, Barel I, Urbakh M 2011. Slow cracklike dynamics at the onset of frictional sliding. Phys. Rev. Lett. 107, 235501 (doi:10.1103/PhysRevLett.107.235501) [PubMed]
44. Bar-Sinai Y, Brener EA, Bouchbinder E 2012. Slow rupture of frictional interfaces. Geophys. Res. Lett. 39, L03308 (doi:10.1029/2011GL050554)
45. Kilgore BD, Blanpied ML, Dieterich JH 1993. Velocity dependent friction of granite over a wide range of conditions. Geophys. Res. Lett. 20, 903–906. (doi:10.1029/93GL00368)
46. Heslot F, Baumberger T, Perrin B, Caroli B, Caroli C 1994. Creep, stick-slip, and dry-friction dynamics: experiments and a heuristic model. Phys. Rev. E 49, 4973–4988. (doi:10.1103/PhysRevE.49.4973) [PubMed]
47. Estrin Y, Bréchet Y 1996. On a model of frictional sliding. Pure Appl. Geophys. 147, 745–762. (doi:10.1007/BF01089700)
48. Putelat T, Dawes JHP, Willis JR 2010. Regimes of frictional sliding of a spring-block system. J. Mech. Phys. Solids 58, 27–53. (doi:10.1016/j.jmps.2009.09.001)
49. Bar-Sinai Y, Spatschek R, Brener EA, Bouchbinder E 2014. On the velocity-strengthening behavior of dry friction. J. Geophys. Res. 119, 1738–1748. (doi:10.1002/2013JB010586)
50. Champneys AR, Sandstede B 2007. Numerical computation of coherent structures. In Numerical continuation methods for dynamical systems: path following and boundary value problems (eds Krauskopf B, Osinga HM, Galán-Vioque J), pp. 331–358. Dordrecht, The Netherlands: Springer.
51. Burridge R, Knopoff L 1967. Model and theoretical seismicity. Bull. Seismol. Soc. Am. 57, 341–371.
52. Carlson JM, Langer JS 1989. Mechanical model of an earthquake fault. Phys. Rev. A 40, 6470–6484. (doi:10.1103/PhysRevA.40.6470) [PubMed]
53. Schmittbuhl J, Vilotte J-P, Roux S 1993. Propagative macrodislocation modes in an earthquake fault model. Europhys. Lett. 21, 375 (doi:10.1209/0295-5075/21/3/020)
54. Cartwright JHE, Hernández-García E, Piro O 1997. Burridge-Knopoff models as elastic excitable media. Phys. Rev. Lett. 79, 527–530. (doi:10.1103/PhysRevLett.79.527)
55. Cartwright JHE, Hernández-García E, Piro O 1999. Dynamics of elastic excitable media. Int. J. Bifurc. Chaos 09, 2197–2202. (doi:10.1142/S0218127499001620)
56. Morales JEM, James G, Tonnelier A 2016. Solitary waves in the excitable Burridge-Knopoff model. Technical Report RR-8996, INRIA Grenoble – Rhône-Alpes.
57. Putelat T, Dawes JHP, Willis JR 2007. Sliding modes of two interacting frictional interfaces. J. Mech. Phys. Solids 55, 2073–2105. (doi:10.1016/j.jmps.2007.03.004)
58. Putelat T, Dawes JHP 2015. Steady and transient sliding under rate-and-state friction. J. Mech. Phys. Solids 78, 70–93. (doi:10.1016/j.jmps.2015.01.016)
59. Putelat T, Dawes JHP, Willis JR 2011. On the microphysical foundations of rate-and-state friction. J. Mech. Phys. Solids 59, 1062–1075. (doi:10.1016/j.jmps.2011.02.002)
60. Bristow JR. 1950. Frictional relaxation oscillations. Proc. Phys. Soc. B 63, 964 (doi:10.1088/0370-1301/63/11/114)
61. Burwell JT, Rabinowicz E 1953. The nature of the coefficient of friction. J. Appl. Phys. 24, 136–139. (doi:10.1063/1.1721227)
62. Rabinowicz E. 1957. The intrinsic variables affecting the stick-slip process. Proc. Phys. Soc. London 71, 668–675. (doi:10.1088/0370-1328/71/4/316)
63. Rabinowicz E. 1951. The nature of the static and kinetic coefficients of friction. J. Appl. Phys. 22, 1373–1379. (doi:10.1063/1.1699869)
64. Brace WF, Byerlee JD 1966. Stick slip as a mechanism for earthquakes. Science 153, 990–992. (doi:10.1126/science.153.3739.990) [PubMed]
65. Scholz CH. 1998. Earthquakes and friction laws. Nature 391, 37–41. (doi:10.1038/34097)
66. Dieterich JH. 1979. Modeling rock friction: 1. Experimental results and constitutive equations. J. Geophys. Res. 84, 2161–68. (doi:10.1029/JB084iB05p02161)
67. Ruina AL. 1980. Friction laws and instabilities: a quasistatic analysis of some dry frictional behavior. PhD thesis, Brown University.
68. Ruina AL. 1983. Slip instability and state variable friction laws. J. Geophys. Res. 88, 10 359–10 370. (doi:10.1029/JB088iB12p10359)
69. Baumberger T, Caroli C 2006. Solid friction from stick-slip down to pinning and aging. Adv. Phys. 55, 279–348. (doi:10.1080/00018730600732186)
70. Hulikal S, Bhattacharya K, Lapusta N 2014. Collective behavior of viscoelastic asperities as a model for static and kinetic friction. J. Mech. Phys. Solids 76, 144–161. (doi:10.1016/j.jmps.2014.10.008)
71. Marone C. 1998. Laboratory-derived friction laws and their application to seismic faulting. Annu. Rev. Earth Sci. 26, 643–696. (doi:10.1146/
72. Persson BNJ. 2000. Sliding friction: physical principles and applications, 2nd edn. Nanoscience and Technology. Berlin, Germany: Springer.
73. Scholz CH. 2002. The mechanics of earthquakes and faulting, 2nd edn Cambridge, UK: Cambridge University Press.
74. Cabboi A, Putelat T, Woodhouse J 2016. The frequency response of friction: enhanced rate-and-state models. J. Mech. Phys. Solids 92, 210–236. (doi:10.1016/j.jmps.2016.03.025)
75. Rice JR, Lapusta N, Ranjith K 2001. Rate and state dependent friction and the stability of sliding between elastically deformable solids. J. Mech. Phys. Solids 49, 1865–198. (doi:10.1016/S0022-5096(01)00042-4)
76. Richardson RSH, Nolle H 1976. Surface friction under time-dependent loads. Wear 37, 87–101. (doi:10.1016/0043-1648(76)90183-6)
77. Wang SK, Woodhouse J 2011. The frequency response of dynamic friction: A new view of sliding interfaces. J. Mech. Phys. Solids 59, 1020–1036. (doi:10.1016/j.jmps.2011.02.005)
78. Woodhouse J, Wang SK 2011. The frequency response of dynamic friction: model comparisons. J. Mech. Phys. Solids 59, 2294–2306. (doi:10.1016/j.jmps.2011.08.006)
79. Lamb H. 1917. On waves in an elastic plate. Proc. R. Soc. Lond. A 93, 114–128. (doi:10.1098/rspa.1917.0008)
80. Kolsky H. 1953. Stress waves in solids. Monographs on the physics and chemistry of materials Oxford, UK: Clarendon Press.
81. Rice JR., Ruina AL 1983. Stability of steady frictional slipping. J. Appl. Mech. 50, 343–349. (doi:10.1115/1.3167042)
82. Horowitz F, Ruina AL 1989. Slip patterns in a spatially homogeneous fault model. J. Geophys. Res. 94, 10 279–10 298. (doi:10.1029/JB094iB08p10279)
83. Jones AR. 1967. An experimental investigation of the in-plane elastic moduli of paper. PhD thesis, The Institute of Paper Chemistry.
84. Kuehn C. 2015. Multiple time scale dynamics. Applied mathematical sciences (eds SS Antman, L Greengard, PJ Holmes). Berlin, Germany: Springer International Publishing.
85. Doedel EJ, Keller HB, Kernevez JP 1991. Numerical analysis and control of bifurcation problems. Int. J. Bifurc. Chaos 1, 493–520. (doi:10.1142/S0218127491000397)
86. Champneys AR, Kuznetsov YA, Sandstede B 1996. A numerical toolbox for homoclinic bifurcation analysis. Int. J. Bifurc. Chaos 06, 867–887. (doi:10.1142/S0218127496000485)
87. Kuznetsov YA. 2004. Elements of applied bifurcation theory. Applied Mathematical Sciences, vol. 112, 3rd edn. New York, NY: Springer.
88. Eckhaus W. 1983. Relaxation oscillations including a standard chase on french ducks. In Asymptotic analysis II—surveys and new trends (ed. F Verhulst). Lecture Notes in Mathematics, vol. 985, pp. 449–494. Berlin, Germany: Springer.
89. Homburg AJ, Sandstede B 2010. Homoclinic and heteroclinic bifurcations in vector fields. In Handbook of dynamical systems III (eds H Broer, F Takens, B Hasselblatt), pp. 379–524. Amsterdam, The Netherlands: North-Holland.
90. Wiggins S. 1988. Global bifurcations and chaos: analytical methods. Applied Mathematical Sciences Series, no. 73. New York, NY: Springer.
91. Shilnikov LP, Shilnikov AL, Turaev DV, Chua LO 2001. Methods of qualitative theory in nonlinear dynamics. Singapore: World Scientific.
92. Rice JR, Ben-Zion Y 1996. Slip complexity in earthquake fault models. Proc. Natl Acad. Sci. USA 93, 3811–3818. (doi:10.1073/pnas.93.9.3811) [PubMed]
93. Sandstede B. 2010. Stability of travelling waves. In Handbook of Dynamical Systems II (ed. B Fiedler), pp. 983–1055. Amsterdam, The Netherlands: North-Holland.
94. Evans J. 1975. Nerve axon equations: IV the stable and the unstable impulse. Indiana Univ. Math. J. 24, 1169–1190. (doi:10.1512/iumj.1975.24.24096)
95. Jones CKRT. 1984. Stability of the travelling wave solution of the Fitzhugh-Nagumo system. Trans. Am. Math. Soc. 286, 431–469. (doi:10.1090/S0002-9947-1984-0760971-6)
96. Kapitula T. 2005. Stability analysis of pulses via the evans function: dissipative systems, pp. 407–428. Berlin, Germany: Springer.
97. Kapitula T, Promislow K 2013. Spectral and dynamical stability of nonlinear waves. Applied Mathematical Sciences. New York, NY: Springer.
98. Jones CKRT, Marangell R 2012. The spectrum of travelling wave solutions to the Sine-Gordon equation. Disc. Cont. Dyn. Syst. 5, 925–937. (doi:10.3934/dcdss.2012.5.925)
99. Jones CKRT, Marangell R, Miller PD, Plaza RG 2014. Spectral and modulational stability of periodic wavetrains for the nonlinear Klein-Gordon equation. J. Differ. Equ. 257, 4632-–4703. (doi:10.1016/j.jde.2014.09.004)
100. Hawthorne JC, Rubin AM 2013. Laterally propagating slow slip events in a rate and state friction model with a velocity-weakening to velocity-strengthening transition. J. Geophys. Res. Solid Earth 118, 3785–3808. (doi:10.1002/jgrb.50261)

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