J E M Morales

G James

A Tonnelier

Solitary waves in the excitable Burridge-Knopoff model

Keywords: Solitary waves, Burridge-Knopoff model, Excitable media, Lattice dynamical systems Ondes solitaires, modèle de Burridge-Knopoff, milieux excitables, systèmes dynamiques

The Burridge-Knopoff model is a lattice differential equation describing a chain of blocks connected by springs and pulled over a surface. This model was originally introduced to investigate nonlinear effects arising in the dynamics of earthquake faults. One of the main ingredients of the model is a nonlinear velocity-dependent friction force between the blocks and the fixed surface. For some classes of non-monotonic friction forces, the system displays a large response to perturbations above a threshold, which is characteristic of excitable dynamics. Using extensive numerical simulations, we show that this response corresponds to the propagation of a solitary wave for a broad range of friction laws (smooth or nonsmooth) and parameter values. These solitary waves develop shock-like profiles at large coupling (a phenomenon connected with the existence of weak solutions in a formal continuum limit) and propagation failure occurs at low coupling. We introduce a simplified piecewise linear friction law (reminiscent of the McKean nonlinearity for excitable cells) which allows us to obtain analytical expression of solitary waves and study some of their qualitative properties, such as wavespeed and propagation failure. We propose a possible physical realization of this system as a chain of impulsively forced mechanical oscillators. In certain parameter regimes, non-monotonic friction forces can also give rise to bistability between the ground state and limit-cycle oscillations and allow for the propagation of fronts connecting these two stable states.

Introduction

A significant body of work has been devoted to elucidating nonlinear mechanisms of earthquakes [START_REF] Kagan | Observational evidence for earthquakes as a nonlinear dynamic process[END_REF][START_REF] Scholz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Kaproth | Slow earthquakes, preseismic velocity changes, and the origin of slow frictional stick-slip[END_REF]. Almost fifty years ago, Burridge and Knopoff [START_REF] Burridge | Model and theoretical seismicity[END_REF] introduced a nonlinear lattice model to investigate the generation of earthquakes along faults, or more generally the occurrence of dynamical instabilities at frictional interfaces. The Burridge-Knopoff (BK) model formally describes a chain of blocks connected by springs and pulled over a surface, each block being attached to a spring pulled at constant velocity and subject to a friction force with the surface. When considering two plates under friction, the blocks can either correspond to a discretization of a plate or account for an existing microstructure [START_REF] Huang | Biomimetic design of elastomer surface pattern for friction control under wet conditions[END_REF]. It has been shown that this simple slider-block model is able to reproduce some statistical features of earthquakes generated by fault dynamics [START_REF] Carlson | Properties of earthquakes generated by fault dynamics[END_REF].

A key feature of the BK model lies in the friction force between the blocks and the fixed surface, which depends nonlinearly on sliding velocity. Experiments performed with a broad range of materials have revealed that the steady-state kinetic friction coefficient is non-monotone versus sliding velocity (see [START_REF] Bar-Sinai | On the velocity-strengthening behaviour of dry friction[END_REF] for a review). Friction is velocity-strengthening (i.e. frictional resistance increases with sliding velocity) for high enough velocities and becomes velocity-weakening in a regime of low velocities, a behaviour intimately linked with the occurrence of stick-slip instabilities and earthquake phenomena. Such velocity-weakening or velocity-strengthening regimes can be described by different types of friction laws. Rate-and-state friction laws incorporating state variables for the frictional interface are frequently used, see e.g. [START_REF] Putelat | Regimes of frictional sliding of a spring-block system[END_REF] for a review and a detailed bifurcation analysis in the case of a single block, and [START_REF] Wang | One-dimensional dynamic simulations of slip complexity of earthquake faults[END_REF][START_REF] Ohmura | Rate-and state-dependent friction law and statistical properties of earthquakes[END_REF][START_REF] Erickson | Periodicity, chaos and localization in a Burridgeâ-Knopoff model of an earthquake with rate-and-state friction[END_REF] for dynamical simulations of the BK model. However, the mathematical analysis of such models can be delicate due to the additional degrees of freedom introduced to describe the state of the interface. Another class of widely used friction laws is given by generalized Coulomb laws with velocity-dependent kinetic friction coefficient [START_REF] Wang | One-dimensional dynamical modelling of earthquakes: a review[END_REF]. In that case, the set-valued character of the friction laws requires an adapted numerical treatment [START_REF] Acary | Numerical methods for nonsmooth dynamical systems[END_REF][START_REF] Bastien | Theoretical study of a chain sliding on a fixed support[END_REF][START_REF] Xiong | Implicit Euler simulation of one-dimensional Burridge-Knopoff model of earthquakes with set-valued friction laws[END_REF] and leads to analytical complications in dynamical studies [START_REF] Lebellego | Phénomènes ondulatoires dans un modèle discret de faille sismique[END_REF]. Alternatively, single-valued laws can be used to approximate set-valued friction laws [START_REF] Carlson | Properties of earthquakes generated by fault dynamics[END_REF] or account for a velocity-strengthening regime at small enough velocities [START_REF] Burridge | Model and theoretical seismicity[END_REF][START_REF] Vorvolakos | The effects of molecular weight and temperature on the kinetic friction of silicone rubbers[END_REF][START_REF] Gong | Friction and lubrication of hydrogels -its richness and complexity[END_REF].

The dynamics of the BK model has been extensively studied in the case of steady-state velocity-weakening friction. Depending on the choice of parameters and system size, this regime can lead to chaotic dynamics or to the propagation of nonlinear wavetrains [START_REF] Carlson | Mechanical model of an earthquake fault[END_REF][START_REF] Schmittbuhl | Propagative macrodislocation modes in an earthquake fault model[END_REF][START_REF] Erickson | Periodicity, chaos and localization in a Burridgeâ-Knopoff model of an earthquake with rate-and-state friction[END_REF]. In particular, periodic travelling waves close to solitary waves (with highly localized slipping events propagating at constant velocity) have been reported in numerical and analytical studies [START_REF] Schmittbuhl | Propagative macrodislocation modes in an earthquake fault model[END_REF][START_REF] Español | Propagative slipping modes in a spring-block model[END_REF][START_REF] Schmittbuhl | A dissipation-based analysis of an earthquake fault model[END_REF][START_REF] Erickson | Periodicity, chaos and localization in a Burridgeâ-Knopoff model of an earthquake with rate-and-state friction[END_REF][START_REF] Lebellego | Phénomènes ondulatoires dans un modèle discret de faille sismique[END_REF].

In this paper, we consider a different situation corresponding to 'spinodal' friction laws [START_REF] Burridge | Model and theoretical seismicity[END_REF][START_REF] Estrin | On a model of frictional sliding[END_REF][START_REF] Putelat | Regimes of frictional sliding of a spring-block system[END_REF], where steady-state friction is velocity-strengthening both for small and large enough velocities and an intermediate velocity-weakening region exists. This situation has been reported in a number of friction experiments performed with rocks, rubber and hydrogels [START_REF] Shimamoto | Transition between frictional slip and ductile flow for halite shear zones at room temperature[END_REF][START_REF] Wu-Bavouzet | Stick-Slip: wet versus dry[END_REF][START_REF] Gong | Friction and lubrication of hydrogels -its richness and complexity[END_REF]. In that case, the dynamics of a single block is described by a Van der Pol type equation, a situation frequently encountered in the modelling of excitable media [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF][START_REF] Zeeman | Differential equations for the heartbeat and nerve impulse[END_REF][START_REF] Feingold | Phase locking, period doubling, and chaotic phenomena in externally driven excitable systems[END_REF][START_REF] Izhikevich | Dynamical systems in neuroscience: The geometry of excitability and bursting[END_REF]. Numerical studies of the BK model with spinodal friction laws (either generalized Coulomb laws or regularizations thereof) have revealed different types of wave patterns, ranging from synchronous oscillations, periodic travelling waves and phase fronts [START_REF] Cartwright | Burridge-Knopoff models as elastic excitable media[END_REF][START_REF] Cartwright | Dynamics of elastic excitable media[END_REF] to chimera-like states (see Fig. 8 of reference [START_REF] Xu | Periodicity and chaos in a one-dimensional dynamical model of earthquakes[END_REF] for an early observation of this phenomenon). These different studies were focused on the oscillatory regime where the pulling velocity lies within the velocity-weakening region.

With regard to spatially localized travelling waves, the existence of fronts has been established in a continuum limit of the BK model [START_REF] Muratov | Traveling wave solutions in the Burridge-Knopoff model[END_REF] (see also [START_REF] Bouchbinder | Slow cracklike dynamics at the onset of frictional sliding[END_REF][START_REF] Bar-Sinai | Instabilities at frictional interfaces: Creep patches, nucleation, and rupture fronts[END_REF][START_REF] Bar-Sinai | Velocity-strengthening friction significantly affects interfacial dynamics, strength and dissipation[END_REF] for numerical studies of rupture fronts in other continuum models based on rate-and-state laws). However, the existence of localized waves was not established so far for the spatially discrete BK model with spinodal friction laws. The analysis of finite amplitude travelling waves in the discrete system is more delicate because it leads to a nonlinear advance-delay differential equation for the wave profiles. In order to tackle this problem, an interesting inverse approach was described in reference [START_REF] Comte | Discrete Burridge-Knopoff model, with exact solitonic or compactlike traveling wave solution[END_REF], where particular friction laws were computed in order to fit prescribed explicit wave profiles. An analytical moving kink solution (with block displacements given by an odd function of the moving frame coordinate z) was proposed but turns out to be erroneous (the velocity dependent friction force induced by the kink is even in z, and cannot be balanced by the odd inertial and stress terms in the dynamical equation [START_REF] Cahn | Traveling wave solutions for systems of ODEs on a two-dimensional spatial lattice[END_REF] of [START_REF] Comte | Discrete Burridge-Knopoff model, with exact solitonic or compactlike traveling wave solution[END_REF]).

In the present paper, we analyse in detail the existence of localized waves in the discrete BK model with spinodal friction laws, using both extensive numerical simulations and analytical techniques. We restrict our attention to solitary waves and will treat the case of fronts in a forthcoming publication. We consider the excitable regime where the pulling velocity lies within a velocity-strengthening domain of the spinodal friction laws (above velocity-weakening). In that case, each single block admits a stable state of continuous slip but displays a large response to perturbations above some threshold. When blocks become coupled, our numerical simulations reveal that this response generates a solitary wave for a broad range of friction laws (smooth or nonsmooth) and parameter values. In the case of regularized Coulomb laws, the solitary wave consists of a moving localized region where sticking occurs, in contrast to the propagation of localized slipping events previously reported for velocity-weakening friction [START_REF] Schmittbuhl | Propagative macrodislocation modes in an earthquake fault model[END_REF].

We analyse in detail the influence of parameters on the existence and qualitative properties of solitary waves. For low enough coupling between blocks, we observe a phenomenon of propagation failure corresponding to the rapid extinction of initially propagating pulses. We find that solitary waves develop shock-like profiles in the opposite limit of large coupling, a phenomenon connected with the existence of weak solutions in a formal continuum limit. When the pulling velocity is decreased towards the boundary of the velocity-strengthening domain, the dynamics of the blocks becomes underdamped and we observe solitary waves with oscillatory tails, while propagation failure takes place above some critical pulling velocity. For certain friction laws (near the transition to a velocity-weakening law), we also observe bistability between continuous slip and limit-cycle oscillations and the existence of propagating fronts connecting these two stable states.

In order to obtain analytical expressions for solitary waves and explain some of their qualitative properties, we introduce a simplified piecewise linear friction law with two velocitystrengthening regions separated by a (negative) jump discontinuity. The discontinuity in the friction force can be interpreted as a rough approximation for the existence of a small intermediate velocity-weakening region. Alternatively, we propose a possible physical realization of this piecewise-linear system as a chain of impulsively forced mechanical oscillators. Such types of nonlinearities have been extensively used for the mathematical study of travelling waves in different types of PDEs [START_REF] Mckean | Nagumo's equation[END_REF][START_REF] Rinzel | Traveling wave solutions of a nerve conduction equation[END_REF][START_REF] Zemskov | Stability of travelling fronts in a piecewise-linear reaction-diffusion system[END_REF] and spatially discrete systems [START_REF] Flytzanis | High velocity dislocation motion and interatomic force law[END_REF][START_REF] Cahn | Traveling wave solutions for systems of ODEs on a two-dimensional spatial lattice[END_REF][START_REF] Tonnelier | McKean caricature of the FitzHugh-Nagumo model: traveling pulses in a discrete diffusive medium[END_REF][START_REF] Elmer | Spatially discrete FitzHugh-Nagumo equations[END_REF][START_REF] Truskinovsky | Kinetics of martensitic phase transitions: Lattice model[END_REF][START_REF] Schwetlick | Existence of dynamic phase transitions in a one-dimensional lattice model with piecewise quadratic interaction potential[END_REF][START_REF] Atkinson | Motion of a Frenkel-Kontorova dislocation in a onedimensional crystal[END_REF]. Following this approach allows us to obtain explicit solitary waves in the form of oscillatory integrals and numerically compute their profiles and wave speed. Moreover, this method provides analytical expansions of the solitary waves at small coupling and for pulling velocities near the jump discontinuity [START_REF] Tonnelier | McKean caricature of the FitzHugh-Nagumo model: traveling pulses in a discrete diffusive medium[END_REF]. In this limit, we explain the occurrence of propagation failure below some critical coupling by the existence of a saddle-node bifurcation of solitary waves.

It is worthwhile to stress that we obtain fully localized solitary waves, i.e. blocks lie in the state of stable slip at infinity on both sides of the chain (this constitutes another important difference with the previous works [START_REF] Schmittbuhl | Propagative macrodislocation modes in an earthquake fault model[END_REF][START_REF] Lebellego | Phénomènes ondulatoires dans un modèle discret de faille sismique[END_REF]). Similar solutions have been previously obtained for other types of excitable lattice dynamical systems where the coupling is diffusive rather than elastic. These systems correspond to spatially discrete FitzHugh-Nagumo equations with either smooth [START_REF] Hupkes | Travelling pulses for the discrete FitzHugh-Nagumo system[END_REF] or piecewise-linear [START_REF] Tonnelier | McKean caricature of the FitzHugh-Nagumo model: traveling pulses in a discrete diffusive medium[END_REF] bistable nonlinearities. While our analysis of solitary waves at small coupling closely follows the lines of [START_REF] Tonnelier | McKean caricature of the FitzHugh-Nagumo model: traveling pulses in a discrete diffusive medium[END_REF], the properties of solitary waves are quite different at large coupling with the occurrence of shocks in the BK model.

The outline of this paper is as follows. In section 2, we introduce the excitable BK model with different types of spinodal friction laws, summarize several applications of the model (frictional interfaces, mechanical oscillators, nonlinear transmission lines) and illustrate the excitation of solitary waves. In section 3, we use numerical simulations to study the existence and shape of solitary waves depending on parameters. We also consider a continuum limit of the model and illustrate the bistable dynamics near the transition to a velocity weakening law. Section 4 provides analytical results for the idealized piecewise linear friction law. Section 5 summarizes the main findings, describes potential applications of the present work and points out interesting open problems.

2 The Burridge-Knopoff model and its applications

Dynamical equations and solitary waves

The BK model describes the interaction of two sheared surfaces, one being pulled at a constant speed V and the second is discretized as a chain of N identical blocks with mass m resting over a surface with frictional interaction described by a nonlinear friction function F (v). The blocks are connected to each other through Hookean springs of stiffness k c , and to the upper plate by springs of stiffness k p (see Fig. 1). The motion of an isolated block is given by

m d 2 x dt 2 + k p x + f 0 F V + dx dt v 1 = 0, (1)
where x(t) is the displacement of the block relative to the point of attachment on the moving plate, f 0 and v 1 are two scaling parameters of the friction law.

Using the change of variables y 1) can be reformulated as the dimensionless dynamical system on the plane:

= x/ρ, τ = ρ/v 1 , t = t/τ , γ = mv 2 1 k p /f 2 0 , V = V /v 1 and ρ = f 0 /k p , Eq. (
dy dt = u, γ du dt = -F (V + u) -y, (2)
where the rescaled time t and pulling velocity V have been rewritten as t and V . The qualitative dynamics of (2) highly depends on the shape of the nonlinear function F (v). For a spinodal friction law, similar to the one plotted in Fig. 2A, the system shows excitability: a strong enough perturbation of the steady state (u s , y s ) = (0, -F (V )) produces a large trajectory in the phase plane whereas small perturbations monotonically decay to the resting state. In addition, there exists a refractory period where the block is unresponsive to further perturbations. The geometrical illustration of this threshold effect is shown in Fig. 2B. This excitability property contrasts with the one due to static friction where a critical amount of stress has to be accumulated in order to reach the unstable slipping mode. In our case, for a sufficiently large slip rate V , the slip solution is stable and a temporary incursion in the velocity strengthening creep region is performed following a strong enough perturbation. When the blocks are interconnected, the resulting network is a discrete excitable medium and the perturbation of one block may propagate to its neighbours producing complex spatiotemporal dynamics, including periodic wave trains, propagating fronts, global oscillations, chaotic dynamics or complex transitions [START_REF] Cartwright | Burridge-Knopoff models as elastic excitable media[END_REF][START_REF] Cartwright | Dynamics of elastic excitable media[END_REF]. In the original one-dimensional BK model, the blocks are connected together by springs and the resulting force on the n-th block is given by k c ∆ d x n where ∆ d x n = x n+1 -2x n + x n-1 is the discrete Laplacian and k c is the linear elastic coupling Inria parameter. Using the dimensionless coupling parameter k = k c /k p , the governing equations for the n-th block read

F (v) v 1
dy n dt = u n , γ du n dt = k∆ d y n -F (V + u n ) -y n , (3)
where n ∈ Z.

We report here results on the existence of solitary wave solutions where each block returns to the stable state of continuous slip after wave passage (see Fig. 2C,D). By a solitary wave solution we mean a solution of (3) such that

u n (t) = ϕ(n -ct), y n (t) = ψ(n -ct), (4)
with (ϕ(±∞), ψ(±∞)) = (0, -F (V )). The function (ϕ, ψ) defines the waveform, and c is the wave speed that has to be determined. Substitution of ( 4) into (3) yields the following advancedelay differential equation:

-cψ (ξ) = ϕ(ξ), -cγϕ (ξ) = k[ψ(ξ + 1) -2ψ(ξ) + ψ(ξ -1)] -F (V + ϕ(ξ)) -ψ(ξ), (5)
where ξ = nct ∈ R represents the travelling wave coordinate.

The spatiotemporal plot of a solitary wave for the cubic friction law is shown in Fig. 2C and the corresponding waveforms are depicted in Fig. 2D. In the forthcoming sections, we show that the BK model exhibits solitary wave solutions for a large class of nonlinearities F of spinodal type. We also review different physical contexts where such nonlinearities arise.

Friction dynamics

Several non-monotonic friction laws have been proposed in the context of earthquake fault dynamics. Many studies consider generalized Coulomb friction laws with a multivalued part at the origin (see Fig. 3A,B). This friction law involves a sticking condition when a block achieves a zero velocity with respect to the lower plate. Single-valued regularized Coulomb laws are also often considered for numerical purpose, or to account for the existence of a velocity-strengthening region for low sliding velocities (this can represent e.g. lubrication effects). In addition, friction laws that switch from a velocity-weakening to a velocity-strengthening regime for large sliding velocities are also relevant for applications, see [START_REF] Heslot | Creep, stick-slip, and dry friction dynamics: Experiments and a heuristic model[END_REF] and references therein. Combining these different effects results in spinodal friction laws [START_REF] Cartwright | Burridge-Knopoff models as elastic excitable media[END_REF][START_REF] Cartwright | Dynamics of elastic excitable media[END_REF] of the type depicted in Fig. 2A.

In this paper, we investigate the dynamics of the BK model using the prototypical smooth spinodal friction law given by the cubic function F c (v) = 3.2v 3 -7.2v 2 + 4.8v, depicted in Fig. 2A. We shall denote by v = a the local minimum of the spinodal friction force, the above cubic friction law corresponding to a = 1. In order to analyse the effects of a larger class of spinodal friction laws (in particular with different degrees of smoothness or stiffness), we also introduce the following class of piecewise linear functions

F P L (v; a, b, µ) =        v b , 0 ≤ v < b α v-b b-a + 1, b < v < a µ(v -a) + 1 -α, a < v (6)
RR n°8996 that have been normalized so that F P L (b) = 1. For v < 0 we set F P L (v) = -F P L (-v). The parameter b is the threshold where the transition from velocity strengthening to velocity weakening occurs. If µ ≥ 0, velocity strengthening is recovered above the sliding velocity a and F (a) = 1-α (α ∈ (0, 1)) corresponds to a local minimum of the friction force. We shall consider more specifically four piecewise linear functions defined as follows and depicted in Fig. 3C-F. A crude approximation of spinodal friction laws is given by F 0 (v) = F P L (v; a, a, 1/a) (Fig. 3C). In that case, the damping ratio in the two velocity-strengthening regimes are identical and the velocity-weakening region is reduced to a jump discontinuity. This friction function can be rewritten as F 0 (v) = v/a -αH(va) where H is the Heaviside step function. The friction law defined by F 1 = F P L (v; a, a/2, 2/a) is obtained from F 0 with the addition of an intermediate velocity-weakening region with finite thickness (Fig. 3D). To regularize a multivalued generalized Coulomb law, we also consider F 2 (v) = F P L (v; a, b, 2/a) with b = 0.01 (Fig. 3E). For this parameter value, the friction law is close to the multivalued case depicted in Fig. 3A. It would be also interesting to investigate the limit case of Fig. 3B (considered e.g. in [START_REF] Muratov | Traveling wave solutions in the Burridge-Knopoff model[END_REF]), but the numerical simulation of the BK model is more delicate in that case and will be investigated in future work. Finally, to investigate the transition from a velocity-strengthening law to a velocity-weakening law at large sliding velocities we introduce F 3 (v) = F P L (v; a, a/2, µ) (Fig. 3F) where the slope µ may be negative.

F (v) v 1 1 -α a A F (v) v 1 1 -α B F 0 (v) v 1 1 -α a C F 1 (v) v 1 1 -α b a D F 2 (v) v 1 1 -α b a E F 3 (v) v 1
Note that the present study examines different types of friction laws but is limited to purely linear nearest-neighbours coupling between blocks. Other types of couplings would be relevant for applications, such as long-range elastic couplings or contact interactions (unilateral springs). We have checked that solitary waves can be generated as well in the case of Hertzian interactions or for a linear contact law (see Fig. 4A,B), but a detailed investigation of such systems lies beyond the scope of the present study.

Oscillator chain with impulsive excitations

In this section we describe a possible physical implementation of the BK model with piecewiselinear friction function F 0 (Fig. 3-C). This system corresponds to a chain of coupled linear oscillators, where each element is subject to an impulsive force (applied e.g. using piezo actuators) when reaching a critical deflection z = ζ. The dynamical equations read

m zn + ν żn + k p z n = k c ∆ d z n + λ t k ∈ R, z n (t k ) = ζ sign( żn (t - k )) δ t k , (7)
where z n (t) is the deflection of the nth oscillator in the chain. In the right side of (7), δ t k denotes the Dirac distribution at t = t k , sign the usual sign function (odd and equal to unity on (0, +∞)), and we use the notation żn (t - k ) (resp. żn (t + k )) for the left (resp. right) limit of żn at t = t k . For the series of impulses to be well defined as a distribution, it is assumed that each component z n (t) crosses the critical value z = ζ for a countable set of times t k (generally depending on n) without finite accumulation point (i.e. Zeno behaviour does not occur). The external (state-dependent) impulsion has a fixed magnitude λ > 0 and is always applied in the direction of motion. Each oscillator is damped (ν > 0 denotes the associated damping constant) and the other physical parameters in [START_REF] Biktashev | Solitary waves in excitable systems with crossdiffusion[END_REF] are the same as in section 2.1.

The components z n correspond to continuous piecewise-differentiable functions, displaying jumps of derivatives at t = t k . One obtains from [START_REF] Biktashev | Solitary waves in excitable systems with crossdiffusion[END_REF] [ żn ]

t k = λ m sign( żn (t - k )), (8)
with

[f ] t0 = f (t + 0 ) -f (t - 0 )
denoting the jump discontinuity of a function f at t 0 . In order to map [START_REF] Biktashev | Solitary waves in excitable systems with crossdiffusion[END_REF] to the piecewise-linear BK model, we need the following technical lemma. Lemma 1. Let u denote a continuous piecewise-differentiable function such that u(0) = 0, u(0 -) = 0. The following properties are equivalent :

i) sign([ u] 0 ) = sign( u(0 -)), ii) sign([ u] 0 ) = [H(u)] 0 .
Proof. Let us first assume i) and show that ii) holds true. Since u(0

+ ) = u(0 -) + [ u] 0 , i) implies sign( u(0 + )) = sign( u(0 -))
. This equality implies [H(u)] 0 = sign( u(0 ± )) (since u(0 ± ) = 0), which leads to ii) using i). Now let us show that ii) implies i). We prove this statement for u(0 -) > 0, the proof for u(0 -) < 0 being similar. The above assumptions imply sign(

[ u] 0 ) = [H(u)] 0 ≥ 0, hence u(0 + ) ≥ u(0 -) > 0, leading to [H(u)] 0 = 1.
From ii) we then get sign([ u] 0 ) = 1 and thus i) holds true.

In the sequel we restrict our attention to solutions of (7) satisfying the following transversality condition for all n and all crossing times t k :

żn (t - k ) = 0. (9)
Applying lemma 1 to

u(t) = z n (t + t k ) -ζ, property i) is satisfied thanks to (8) (with λ, m > 0) and thus sign( u(0 -)) = [H(u)] 0 , meaning that sign( żn (t - k )) = [H(z n -ζ)] t k . (10)
Consequently, equation ( 7) and the transversality condition (9) lead to

m zn + ν żn + k p z n = k c ∆ d z n + λ d dt H(z n -ζ). (11)
Similarly, applying lemma 1-ii) to any solution of ( 11) satisfying [START_REF] Burridge | Model and theoretical seismicity[END_REF], one establishes the equivalence between equations ( 7) and [START_REF] Carlson | Mechanical model of an earthquake fault[END_REF] under the above transversality condition.

Assuming further z n (0) = ζ for all n and integrating (11) on [0, t], one obtains the equivalent formulation

m żn + ν (z n -ζ - λ ν H(z n -ζ)) = (k c ∆ d -k p ) t 0 z n (s) ds + p n (12)
with

p n = m żn (0) + ν(z n (0) -ζ) -λ H(z n (0) -ζ).
Let us either consider an infinite lattice with periodic or free end boundary conditions, or bounded sequences (z n , żn ) on an infinite lattice. In each case the linear map k c ∆ dk p is invertible (k c , k p > 0) and one can introduce the new variable

x n (t) = t 0 z n (s) ds + (k c ∆ d -k p ) -1 p n - f 0 k p ,
where f 0 denotes an arbitrary constant. Substitution in [START_REF] Carlson | Properties of earthquakes generated by fault dynamics[END_REF] leads to the piecewise-linear BK model

m ẍn + F ( ẋn ) + k p x n = k c ∆ d x n , (13)
where

F (u) = ν (u -ζ) -λ H(u -ζ) + f 0 . ( 14
)
This family of piecewise-linear functions is parametrized by the slope ν, position ζ and magnitude λ of the jump discontinuity and local maximum f 0 .

Inria

Nonlinear transmission lines

The mechanical BK system (3) has an electrical equivalent as an active transmission line made of electronic L-C circuits coupled to a diode. The corresponding two-node circuit has been already introduced and studied [START_REF] Venturi | Marginal stability and chaos in coupled faults modelled by nonlinear circuits[END_REF].

C L C L C L L c L c V n-1 V n V n+1 V d
Application of Kirchoff's current law at the node n (see Fig. 5) gives

C d 2 Y n dt 2 + 1 L + 1 L c Y n - 1 L c (Y n+1 + Y n-1 ) + F V d + dY n dt = 0, (15)
where

Y n = (V n -V d ), V n is
+ L c ), C = γ, k = 1/L c , y n (t) = Y n (t)
we recover the BK system (3). The tunnel diode has a spinodal current-voltage characteristic as the one depicted in Fig. 2A (in the I-V plane) where a negative resistance occurs at intermediate voltage values. A sufficiently strong voltage perturbation may initiate the propagation of a solitary wave along the transmission line. Potential applications of nonlinear electrical wave propagation in electronic devices include broadband circuits, microwaves, instrumentation and the electronic generation of optical signals [START_REF] Remoissenet | Waves called solitons -Concepts and experiments[END_REF].

3 Qualitative properties of solitary waves

Regimes of solitary wave generation

Unless stated otherwise, simulations are done with a network size and for a time duration such that a stationary solution is reached. The network is initialized at its resting state (u s , y s ) n and free boundary conditions are taken. For the numerical integration, the adaptive scheme Lsoda is employed with a minimal error tolerance set at 1.5e-8. A time step ∆t = 0.001 is used to follow the numerical solution. All the numerical results are illustrated with the fixed parameter γ = 0.15 that sets the simulations in the overdamped regime. The underdamped case will be briefly discussed. The existence of solitary waves is explored in the excitable regime, i.e. for a pulling velocity such that a < V < V max , where V max > a is defined through F (V max ) = 1. The two limiting regimes, V → a + and V → V - max , will be referred to as the small pulling velocity regime and the maximal pulling velocity regime, respectively. To explore intermediate velocity regimes, we also consider a medium pulling velocity V med > a defined through F (V med ) = 1-α/2. For ensuring a stable localized wave propagation, typical network size is of order one hundred for small pulling velocities, whereas medium velocities require around ten thousands blocks. Such large network simulations have to be performed because the increase of V also requires increasing k in order to excite solitary waves (see below). This parameter regime results in higher wave speed, convergence time towards steady state, and pulse width (simulations for these network sizes can be efficiently addressed using Python programming language).

Propagating waves are initiated using a shock-like initial condition: a strong perturbation is applied on the resting state of a subset of blocks in the network, possibly reduced to a single block perturbation near the left edge of the network. The shock excitation is applied to the position or to the velocity. The magnitude of the excitation required to generate a propagation depends on V . In the regime of small pulling velocity, a velocity perturbation of magnitude in the interval [0.125, 2.5] is sufficient to excite a solitary wave for all friction forces. An excitation of the position is more efficient to produce trains of solitary waves and a perturbation magnitude in the interval [0.01, 0.5] is used (this higher sensitivity to perturbations of positions is clearly illustrated for a single block by Fig. 2-B). For a medium pulling velocity (V ∼ V med ), a stronger perturbation on the velocity is needed, up to 10, whereas perturbations in the position remain similar. In the maximal velocity regime, excitations of the network fail to produce travelling waves. Localized solutions can propagate stably along the network in the excitable regime for a sufficiently large coupling value. The wave speed and width increase with the coupling value when the other parameters are fixed, while the increase of γ decreases the wave speed. For a given pulling velocity, there exists a critical coupling value k * below which no propagation occurs. This phenomenon has been coined propagation failure and is well documented for excitable lattices of diffusive type [START_REF] Keener | Propagation and its failure in coupled systems of discrete excitable cells[END_REF]. As V increases, the critical value k * increases and a vertical asymptote is observed near V = V max . This is shown in Fig. 6B where the regime of existence and stability of solitary waves is indicated in the Vk plane for the cubic friction law. Qualitatively similar results are obtained for the other spinodal friction functions (results not shown). Table 1 presents computations of k * values for different friction laws and pulling velocities.

Inria

Friction function Small pulling velocity Medium pulling velocity

V = 1.025 V med Fc(v) k * ≈ 1.4 k * ≈ 360.8 F0(v) k * ≈ 0.43 k * ≈ 1915.3 F1(v) k * ≈ 4.5 k * ≈ 99 F2(v) k * ≈ 48 k * ≈ 602
Table 1 -Critical coupling k * for different friction functions and different pulling velocities.

Parameter values for a, b and α are as in Fig. 3.

The increase of k * versus V can be intuitively explained by the geometry of excitability of a single block: for small pulling velocities, the initial state of the network is near the right knee of the spinodal friction function and a small perturbation is sufficient to reach the threshold of excitability. As V increases, the steady state gets away from the threshold separatrix and the perturbation has to be stronger in order to initiate a large response of the block. Above V max the excitability property of the block disappears. The waveforms obtained for different friction laws, pulling velocities V and coupling strengths k share common features. The displacement profile presents a hump whereas the velocity profile shows a rapid downstroke followed by an upstroke reminiscent to the action potential of spiking neurons. This is illustrated in Fig. 7A-D using the different friction laws. Note that R ϕ dξ = c (ψ(+∞)ψ(-∞)) = 0. Moreover, the discontinuous law defined by F 0 produces a rather symmetric velocity profile where the positive and negative peaks are quasi-equidistant from the resting state. This property is more pronounced in the small pulling velocity regime and tends to break down as V increases. This symmetry is broken for the other friction laws (F c , F 1 , F 2 and F 3 ) where the positive peak of the velocity profile is less marked. When the velocity weakening range increases, the amplitude of the negative peak of the velocity profile becomes larger while the wave speed decreases. The tail of the velocity profile presents either a monotonic or an oscillatory return to the resting state. This feature is determined by the modes around the resting state. Linearisation of (3) around the steady state (u n , y n ) = (0, -F (V )) reads

γ ÿn = k∆y n -y n -F (V ) ẏn . ( 16
)
In order to analyse the relaxation towards the rest state, we look for normal modes y n (t) = Ae iqn+σ(q)t + c.c. where A = 0 and σ(q) are two complex numbers, q ∈ [0, π] and c.c. stands for complex conjugate. This leads to

σ(q) = -F (V ) ± ∆(q) 1 2γ , (17)
where ∆(q) = F (V ) 2 -4γ 1 + 4k sin 2 (q/2) is a decreasing function of q on [0, π]. The asymptotic decay rate is controlled by the slowest decaying mode (obtained for q = 0), which also describes the relaxation of the uncoupled system. When V approaches the local minimum of the cubic friction force F c , the value of F c (V ) tends towards 0 and all modes become underdamped (σ being imaginary). The resulting feature is an oscillatory tail of the wave as shown in Fig. 7D where the envelope of the oscillations becomes larger as V decreases. Due to the constant derivative of the piecewise linear friction forces F 0 , F 1 , F 2 for sliding velocities above a, the modes do not depends on V . In that case the occurrence of oscillatory behaviour can be analysed as above and only depends on γ, k and the slope µ in (6).

In the following two subsections we explore two limiting regimes of particular interest: the continuum limit and a bistability regime.

Continuum limit

We are interested in the study of solitary waves when the discrete BK model approaches a continuum limit, i.e. for k 1. Near the continuum limit, system (3) takes the form:

γ ∂ 2 y ∂t 2 + F V + ∂y ∂t + y = ∂ 2 y ∂x 2 + h 2 12
∂ 4 y ∂x 4 + O(h 4 ) (18)
where we have set 18) is the so-called Boussinesq approximation. We look for solitary wave solutions of ( 18) and set y(x, t) = Y (s) with s = x/ct (for system [START_REF] Bar-Sinai | Instabilities at frictional interfaces: Creep patches, nucleation, and rupture fronts[END_REF], this corresponds to solitary waves with large velocity c = c√ k). The above Ansatz leads to the ODE

y n (t) = y(x, t), x = n/ √ k and h = 1/ √ k 1. Equation (
ΛY -F (V -Y ) -Y = - h 2 12c 4 Y (4) + O(h 4 ) (19)
where Λ = 1/c 2γ.

The fourth derivative at the right side of ( 19) may be useful to account for the appearance of high frequency oscillations (see below), but this problem will not be analysed in the present Inria study. In this work we restrict to the second order model obtained by neglecting O(h 2 ) terms, which reads

ΛY (s) -F (V -Y (s)) -Y = 0. (20)
This equation admits the unique equilibrium (Y, Y ) = (-F (V ), 0). In what follows we study the case of the cubic friction law F = F c . The equilibrium point is a saddle point for Λ > 0, an unstable point (source) for Λ < 0 and ∆ ≥ 0, an unstable focus for Λ < 0 and ∆ < 0, where ∆ = F (V ) 2 + 4Λ. The existence of a smooth solitary wave corresponds to the existence of an homoclinic orbit for equation [START_REF] Estrin | On a model of frictional sliding[END_REF]. In the unstable configuration (Λ < 0) it is clear that such an orbit cannot exist. In the saddle case (Λ > 0), an additional equilibrium inside the homoclinic orbit has to be present because the phase space is two-dimensional. Since there is only one single equilibrium point, smooth solitary waves cannot exist in the continuum model [START_REF] Estrin | On a model of frictional sliding[END_REF]. However for large k values, numerical simulations of the discrete BK model reveal the existence of solitary waves. Typical profiles of these waves are shown in Fig. 8. The velocity component displays two shocks (resulting in two slope discontinuities for the displacement component), and we observe long-lived (i.e. metastable) oscillations in the vicinity of the rear shock. An interesting open question concerns the existence of exact travelling wave solutions of this type (i.e. with an oscillatory shock) in the BK model (3) or in the higher order continuum model obtained by neglecting O(h 4 ) terms in [START_REF] Español | Propagative slipping modes in a spring-block model[END_REF]. These problems will be addressed in future works using numerical continuation techniques. After a very long transient, the fast oscillations near the rear shock disappear, which leads to the shock-wave profiles displayed in the right column of Fig. 8. Given the long time scale of simulations, it is not clear whether this slow drift towards a non-oscillatory shock is a dynamical property of the BK model or is due to numerical integration errors.

The profile of the asymptotic non-oscillatory pulse can be analytically captured by considering equation ( 20) with Λ = 0, which corresponds to fixing c = 1/γ. In that case, the solitary wave of (3) propagates at the sound velocity c = k/γ. This case leads to the differential-algebraic equation

Y = -F (V -Y ). (21)
Solutions of ( 21) are plotted in Fig. 8 (blue lines). They correspond to weak solutions y(x, t) = Y (x/ct) of the nonlinear PDE obtained by setting h = 0 in [START_REF] Ermentrout | Mathematical foundations of neuroscience[END_REF]. In particular, the solution of [START_REF] Feingold | Phase locking, period doubling, and chaotic phenomena in externally driven excitable systems[END_REF] shown the right column of Fig. 8 perfectly matches the stationary profiles shown in Fig. 8C,E. For the parameter values provided in Fig. 8, the wave speed of the numerically observed solitary wave is c ≈ 2587 and the analytical approximation c = k/γ = 2582 is quite accurate (relative error < 2.10 -3 ). In the case of oscillatory shocks, solutions of (21) provide reasonably accurate approximations of block displacements (see Fig. 8D), but they only approximate velocity oscillations on average near the rear shock (Fig. 8B).

Bistability regime

The excitable nature of system (3) can give rise to an oscillatory behaviour, similar to the selfsustained oscillations produced by excitable membranes in biology [START_REF] Ermentrout | Mathematical foundations of neuroscience[END_REF]. The state of continuous slip loses its stability when the friction law operates a transition from a velocity-strengthening regime to a velocity-weakening regime. For the friction law F 3 defined in section 2.2, this case occurs when V > a and the slope parameter µ switches from positive to negative. The dynamics of a single block is illustrated in the top panels of Fig. 9. In the velocity-weakening regime µ < 0, there exists a large-amplitude stable limit cycle, which persists in the range µ > 0 provided µ is not too large (see the bifurcation diagram in Fig. 9A). For the state of continuous slip, a subcritical Hopf bifurcation occurs at µ = 0, leading to the existence of an unstable A Spatiotemporal plot of un(t). When the pulse reaches n > 15000, a quasi-stationary state is reached. B Plot of un(t) (black line) for n = 17500 and D, plot of yn(t) (black line) where the inset is an enlargement that shows the small fluctuations of yn(t) in the highly oscillatory region of un(t). In C and E, we plot the asymptotic solitary wave profiles for n > 95000 (time has been shifted for convenience). In F and G, the blue lines are different solutions of ( 21 periodic orbit when µ is positive and not too large. If µ further increases, the two periodic orbits disappear through a saddle-node bifurcation. Consequently, there exists a narrow interval of positive values of µ where the model exhibits bistability between the state of continuous slip and a limit cycle (a typical phase portrait is shown in Fig. 9B). The two attractors are separated by an unstable periodic orbit, and a change of position (or velocity) can switch the system from rest to oscillation and vice versa.

In this parameter regime, it is an interesting problem to examine the possible spatiotemporal patterns which can be generated by coupling the bistable units. The stability of steady sliding can be analysed through the linear system [START_REF] Elmer | Spatially discrete FitzHugh-Nagumo equations[END_REF], with eigenvalues given by [START_REF] Erickson | Periodicity, chaos and localization in a Burridgeâ-Knopoff model of an earthquake with rate-and-state friction[END_REF] for F (V ) = µ. For µ > 0, the equilibrium state (u s , y s ) is locally asymptotically stable, while for µ < 0 it becomes unstable. For µ > 0, different types of dynamics can be observed depending on the initial perturbation, as exemplified by Fig. 9C-D. In this example, velocity perturbations of one block with large magnitudes produce a front wave with oscillations in the back (periodic travelling waves are clearly visible), while much smaller perturbations produce complex oscillatory patterns. The dynamics of the BK model in the bistable regime is likely to be very rich and needs to be explored in more detail in future works.

Construction of solitary waves for the discontinuous piecewiselinear friction force

In this section we study analytically and numerically the solutions of the following travelling wave equation derived from (5):

c 2 γ d 2 ϕ dξ 2 = k∆ d ϕ + c d dξ F (V + ϕ) -ϕ, (22)
with ϕ(ξ) → 0 as ξ → ±∞. We consider the friction force given by the piecewise linear function:

F 0 (v) = v a -αH(v -a), v > 0,
where the velocity weakening regime is assumed to be instantaneous. A strict velocity threshold occurs at the discontinuity ϕ = -β = a -V < 0 unlike the smooth velocity threshold induced by the cubic friction law. Numerical simulations suggest that this piecewise linear model captures the qualitative features of the localized solitary waves observed for a broader class of friction forces.

The strict velocity threshold makes it possible to classify travelling waves based on the threshold hitting. Let us define the decreasing (finite or infinite) sequence

(ξ k ) k through ϕ(ξ k ) = -β.
Translation invariance of the solution allows us to fix ξ 0 = 0. Different wave types can be obtained according to their signature (σ j ) j defined as (see Fig. 10)

σ j = [H(ϕ + β)] ξj .
We distinguish between threshold crossing from above, threshold crossing from below, and threshold reaching without crossing that read σ j = -1, σ j = 1, and σ j = 0, respectively. The solitary wave solutions are classified according to the changes of sign of ϕ+β that is given by the sequence (σ j ) j . It is worth nothing that travelling sliding solutions are excluded, i.e. the solution is not allowed to remain on the discontinuity line of the friction force. Fixing (σ j ) j , we have:

d dξ F 0 (V + ϕ) = d dξ F 0 (a + ϕ + β) + i [ F 0 (a + ϕ + β) ] ξi δ ξi , = {F 0 (a + ϕ + β)} {ϕ } + [F 0 ] a i σ i δ ξi , = ϕ a -α i σ i δ ξi ,

Inria

where δ ξi is the Dirac distribution at ξ i . Travelling waves in the class (σ j ) j satisfy the linear non-homogeneous differential equation:

c 2 γ ϕ - c a ϕ + ϕ = k ∆ d ϕ -c α i σ i δ ξi . (23)
A solution of ( 23) is a travelling wave solution of ( 22) if the crossing conditions at (ξ j ) j are satisfied together with the sign conditions given by the signature (σ j ) j . This leads to the following admissibility conditions:

   ϕ(ξ 2j ) = ϕ(ξ 2j+1 ) = -β ϕ(ξ) < -β if ξ ∈]ξ 2j+1 , ξ 2j [ ϕ(ξ) > -β otherwise, (24)
where we assume for clarity that there is no threshold reaching without crossing (no i such that σ i = 0). By Fourier transform, we construct a localized solution ϕ(ξ) of system [START_REF] Flytzanis | High velocity dislocation motion and interatomic force law[END_REF]. A travelling wave solution can be expressed as

0 ξ ξ -1 ξ 1 β σ 1 = -1 σ 0 = 0 σ -1 = 1 ϕ(ξ) + β
ϕ(ξ) = p R e i2πλξ g c,ξp (λ)dλ (25)
where

g c,ξp (λ) = αcσ p e -2πiλξp 4π 2 c 2 λ 2 γ -4k sin 2 (πλ) + 2πicλ a -1 -1
.

The analytical expression (25) only provides trial solutions. Indeed, the admissibility conditions [START_REF] Gong | Friction and lubrication of hydrogels -its richness and complexity[END_REF] have to be fulfilled in order to define a solution of [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF]. The simplest solitary waves cross the threshold only two times, at ξ 0 = 0 and at ξ 1 < 0, which gives two nonlinear equations

ϕ(0) = -β, ϕ(ξ 1 ) = -β. (26)
The crossing points of the two level curves [START_REF] Heslot | Creep, stick-slip, and dry friction dynamics: Experiments and a heuristic model[END_REF] in the plane (c, ξ 1 ) determine the speed and the width of trial solutions. We numerically compute the two level curves using a grid with mesh size 0.001. A Gauss-Konrod quadrature formula in a truncated interval [-10 6 , 10 6 ] is used to evaluate [START_REF] Mckean | Nagumo's equation[END_REF]. Approximate solutions for [START_REF] Heslot | Creep, stick-slip, and dry friction dynamics: Experiments and a heuristic model[END_REF] are found with a trust-region dogleg algorithm. As shown in Fig. 11A different crossing points exist.

The inequality constraints of the admissibility conditions have to be checked in order to exclude spurious solutions. Plot of the waveform (25) allows to determine if the sign conditions are fulfilled (see Fig. 11B,C). It can be seen that one of the two crossing points is not associated to a travelling wave solution (see the enlargement in Fig. 11C).

In Fig. 12, we plot the (c, k) curve of the travelling wave solutions. The speed curve is Ushaped with a fold at (c * , k * ). For c > c * stable solitary wave solutions are obtained. For c < c * solitary wave solutions are conjectured to be unstable (they are not observed in the simulations of the network) and become spurious solutions as k increases. A border collision bifurcation is reached below which no solution to (22) can be found with exactly two threshold crossings. We suspect that a travelling sliding solution appears or an additional threshold crossing occurs (both situations are not captured by the present analysis).

Discussion

In this paper, we have reported the existence of solitary waves and oscillatory patterns in the excitable Burridge-Knopoff model. We have studied the influence of the friction law and parameters on the robustness and qualitative properties of solitary waves. Wave profiles have been obtained analytically for a piecewise-linear caricature of the friction law. In that case, we have shown that two branches of solitary waves coexist: one corresponding to stable waves whereas the other consists of presumably unstable waves or spurious solutions.

Despite the very different nature of the coupling, the waves observed in the excitable BK model share many similarities with the taxi waves recently reported in excitable systems with cross-diffusion [START_REF] Biktashev | Solitary waves in excitable systems with crossdiffusion[END_REF][START_REF] Tsyganov | Waves in systems with cross-diffusion as a new class of nonlinear waves[END_REF][START_REF] Tsyganov | Classification of wave regimes in excitable systems with linear cross-diffusion[END_REF]. In particular, reflection from boundaries and penetration of solitary waves through each other are observed in both cases, and annihilation of two solitary waves can be also observed in the discrete BK model (results not shown).

A key problem left open in this paper concerns the theoretical justification of the trial solutions ( 25)- [START_REF] Heslot | Creep, stick-slip, and dry friction dynamics: Experiments and a heuristic model[END_REF], which requires to check the inequality constraints in [START_REF] Gong | Friction and lubrication of hydrogels -its richness and complexity[END_REF]. This problem will be addressed in a next paper in the small coupling limit.

Inria

Fig. 1 -

 1 Fig. 1 -The mechanical interpretation of the Burridge-Knopoff model, where V is the pulling velocity or slip rate of the loader plate, m is the mass of a block, k c and k p are the stiffness of the Hookean springs.

Fig. 2 -

 2 Fig.2 -AThe spinodal friction force defined by the cubic polynomial F c (v) = 3.2v 3 -7.2v 2 +4.8v for v > 0. In B, the responses of an isolated block to perturbations of different magnitudes are shown in the phase plane. The nullclines are represented (dashed lines). For a strong enough perturbation of the resting state, a large trajectory is elicited (solid black), whereas a small perturbation does not produce a significant response (grey line). C When the blocks are connected, an initial perturbation of the network elicits the propagation of a stable solitary wave with constant shape and speed. D Travelling wave profiles of the velocity (filled line) and the displacement (dashed line) in the moving frame coordinate. Parameters are: k = 3, V = 1.05, γ = 0.15. In B, the two perturbations are of magnitude 0.11 and 0.05. In C, an initial perturbation (u 0 (0), y 0 (0)) = (u s -1, y s ) is imposed on the first block (left edge) of the chain.

FFig. 3 -

 3 Fig. 3 -Typical piecewise-linear friction laws for the Burridge-Knopoff model. A-E display different spinodal friction laws (velocity-weakening for an intermediate range of sliding velocities) whereas F is velocity-weakening for large enough velocities. A, B correspond to multivalued laws (generalized Coulomb laws). Velocity-weakening is assumed instantaneous in cases B, C (the velocity-weakening domains are reduced to a jump discontinuity in the friction force). The single-valued laws C-F are used in the present work, with a = 1 and α = 0.2.

Fig. 4 - 2 + 2 +

 422 Fig. 4 -Solitary waves profiles ϕ(nct) = u n (t) (solid) and ψ(nct) = y n (t) (dashed) obtained with (A) Hertzian and (B) linear contact interactions in the Burridge-Knopoff model (cubic friction law). The Hertzian coupling is obtained by replacing the discrete linear Laplacian by k[(y n-1y n ) 3/2 + -(y ny n+1 ) 3/2 + ], and the linear contact law corresponds to k[(y n-1y n ) + -(y ny n+1 ) + ], where (x) + = max(0, x). Parameters are: A, k = 20, V = 1.025 and B, k = 150, V = 1.05. The wave speeds are A, c ≈ 6.36 and B, c ≈ 32.47, respectively. The wave is generated from the chain at rest except for the last block, where the initial condition is A (u s -2, y s ) and B (u s -1, y s ).

Fig. 5 -

 5 Fig. 5 -Electrical analogue of the excitable Burridge-Knopoff model as a nonlinear transmission line.

  the voltage at node n, V d is the driving voltage, C is the capacitance, L and L c are inductances, and the function F defines the current-voltage relationship of the nonlinear component. For L = L c /(1

Fig. 6 -

 6 Fig.6-A Wave profiles of the velocity ϕ(nct) (solid) and the displacement ψ(nct) (dashed) for the cubic friction law F c (v), V = 1.183 and k = 400. The wave speed of the localized travelling wave is c ≈ 50.05. B Regimes of existence and stability of solitary waves in the Vk plane (grey shading) for the cubic friction law. The solid line indicates the critical coupling values k * below which propagation failure occurs. Regimes are computed from numerical simulations in which the wave is initiated by imprinting a velocity u p to the first block in the chain (magnitude up to u p = -10). Note that here V max = 1.25.

Fig. 7 -

 7 Fig. 7 -Plot of the waveforms ϕ(nct) (solid) and ψ(nct) (dashed) for, A F 0 (v), V = 1.025, k = 2, B F 1 (v), V = 1.025, k = 5, C F 2 (v), V = 1.025, k = 70 and, D F c (v), V = 1, k = 2. The local minimum of each spinodal friction force is located at a = 1. The computed wave speed values are c = 4.39, c = 5.36, c = 21.23, c = 4.04 for A-D, respectively. Simulations are done for a strong perturbation (u p = -10) of the velocity of the first block.

Fig. 8 -

 8 Fig. 8 -Simulation of the BK model (3) in the high coupling limit (k = 10 6 ), for V = 1.025 and friction law Fc(v).

1 DFig. 9 -

 19 Fig.9-Simulations of the BK model with friction law F 3 (v) in a bistable regime, with µ 0.2, V = 1.025 and k = 5. A Bifurcation diagram of a single block that shows the amplitudes of the stable (solid line) and unstable (dashed line) periodic orbits. The amplitude is defined by max t>0 (u(t)u s , y(t)y s ) 2 , where (u s , y s ) is the equilibrium state. The latter is stable (solid straight line) or unstable (dotted straight line) depending on µ. B Phase plane of a single block. The nullclines (dashed lines), the stable limit cycle (bold solid circle) and unstable limit cycle (dashed line) are represented. C,D Wave patterns generated in the chain of blocks. A perturbation (u 0 (0), y 0 (0)) = (u s + u p , y s ) with u p = -10 is used in C, whereas u p = -0.25 is used in D. For the simulation shown in C, the asymptotic states of two blocks (n = 0, n = 380) are also plotted in figure B, where the trajectory of the first block converges towards the limit cycle.

Fig. 10 -

 10 Fig. 10 -Signature (σ j ) j for the different threshold crossings.

Fig. 11 -

 11 Fig. 11 -In A plot of the two level curves ϕ(ξ 1 ) = -β (black) and ϕ(0) = -β (grey) in the (c, ξ 1 ) parameter plane. Crossing points define trial solutions for solitary waves of the Burridge-Knopoff model. The crossing point (B) ≈ (4.39, -6.21) is associated with the waveform shown in B (black line). We also plot the solution obtained by the numerical simulation of the network (grey dashed). A perfect match is realized with the analytical solution. The crossing point (C) ≈ (1.65, -0.8) defines a spurious solution as shown in C where the enlargement (inset) reveals the violation of the inequality constraints. Parameters are k = 2, V = 1.025.

Inria

Fig. 12 -

 12 Fig.12-A and B (zoom), curves of the coupling strength k with respect to the wave speed c obtained by solving the non linear system[START_REF] Heslot | Creep, stick-slip, and dry friction dynamics: Experiments and a heuristic model[END_REF]. k * is the critical coupling where propagation failure occurs and c * is the related wave speed. The full line corresponds to stable solitary wave solutions. The dotted line represents solitary waves that are conjectured to be unstable (not observed in dynamical simulations). The grey line denotes spurious solutions. The red crosses are obtained from network simulations. The fold of the speed curve is located at c * ≈ 1.34135, k * ≈ 0.42585. As k increases, the solitary waves with c < c * disappear through a border collision bifurcation at (c s , k s ) ≈ (1.2598, 0.4360). Parameter V = 1.025.

2 2 . 1 2 . 2 2 . 3 2 . 4 11 3 15 4

 212223241115 The Burridge-Knopoff model and its applications 5 Dynamical equations and solitary waves . . . . . . . . . . . . . . . . . . . . . . . 5 Friction dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 Oscillator chain with impulsive excitations . . . . . . . . . . . . . . . . . . . . . . 9 Nonlinear transmission lines . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Qualitative properties of solitary waves 11 3.1 Regimes of solitary wave generation . . . . . . . . . . . . . . . . . . . . . . . . . 11 3.2 Continuum limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 3.3 Bistability regime . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Construction of solitary waves for the discontinuous piecewise-linear friction force Inovallée 655 avenue de l'Europe Montbonnot 38334 Saint Ismier Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

RR n°8996


D J Voneshen

K Refson

E Borissenko

M Krisch

A Bosak

A Piovano

E Cemal

M Enderle

M J Gutmann

M Hoesch

M Roger

L Gannon

A T Boothroyd

S Uthayakumar

D G Porter

J P Goff
email: jon.goff@rhul.ac.uk

Suppression of thermal conductivity by rattling modes in thermoelectric sodium cobaltate

The need for both high electrical conductivity and low thermal conductivity creates a design conflict for thermoelectric systems, leading to the consideration of materials with complicated crystal structures 1 . Rattling of ions in cages results in low thermal conductivity 2-5 , but understanding the mechanism through studies of the phonon dispersion using momentumresolved spectroscopy is made difficult by the complexity of the unit cells 6 . We have performed inelastic X-ray and neutron scattering experiments that are in remarkable agreement with our first-principles density-functional calculations of the phonon dispersion for thermoelectric Na 0.8 CoO 2 , which has a large-period superstructure 7 . We have directly observed an Einstein-like rattling mode at low energy, involving large anharmonic displacements of the sodium ions inside multi-vacancy clusters. These rattling modes suppress the thermal conductivity by a factor of six compared with vacancy-free NaCoO 2 . Our results will guide the design of the next generation of materials for applications in solid-state refrigerators and power recovery.

New thermoelectric materials often have large unit cells and, therefore, calculations of the thermal conductivity, κ, using acoustic phonon branches within the usual Debye framework are inadequate. Instead, it is necessary to take account of the multitude of optical modes and, therefore, detailed knowledge of the phonon dispersion is required. The ability to measure phonon dispersions using inelastic X-ray scattering (IXS) and inelastic neutron scattering (INS) provides unprecedented opportunities to understand thermal transport in these systems. These experiments on samples with a complicated unit cell require a model calculation to assign spectral features. For complex unit cells ab initio calculations are also demanding, and still require experimental validation 6 .

The efficiency of operation of a thermoelectric material is usually measured by its figure-of-merit, ZT = TS 2 σ /κ, a dimensionless number proportional to the square of the Seebeck coefficient, S, and the ratio (σ /κ) of the electrical to thermal conductivities. The phonon-glass electron-crystal model 2 has been proposed, where cages form regular periodic crystal lattices in which electrons or holes can move freely, but the rattling of loosely bound atoms in the cage is believed to disrupt the vibrational modes of the cage so that heat conduction is comparable to glasses 3-5 . However, other theoretical studies have called into question the validity

of the phonon-glass electron-crystal and rattler concepts [START_REF] Feldman | Lattice dynamics and reduced thermal conductivity of filled skutterudites[END_REF][START_REF] Dong | Theoretical study of the lattice thermal conductivity in Ge framework semiconductors[END_REF] and, moreover, INS measurements demonstrate that phonon excitations are sharp for these materials and are certainly not glass like [START_REF] Koza | Breakdown of phonon glass paradigm in La-and Ce-filled Fe 4 Sb 12 skutterudites[END_REF][START_REF] Christensen | Avoided crossing of rattler modes in thermoelectric materials[END_REF] . The presence of a low-energy rattling mode can lead to an avoided crossing with acoustic phonons. Although such an avoided crossing is suggestive of the presence of a rattler, it is not proof, as is confirmed by the fact that it is observed on some occasions, but not others [START_REF] Christensen | Avoided crossing of rattler modes in thermoelectric materials[END_REF] .

Particularly good thermoelectric properties have been observed for sodium cobaltate [START_REF] Terasaki | Large thermoelectric power in NaCo 2 O 4 single crystals[END_REF][START_REF] Wang | Spin entropy as the likely source of enhanced thermopower in Na x Co 2 O 4[END_REF][START_REF] Foo | Charge ordering, commensurability, and metallicity in the phase diagram of the layered Na x CoO 2[END_REF][START_REF] Lee | Large enhancement of the thermopower in Na x CoO 2 at high Na doping[END_REF] . Theoretical studies have considered the electronic contributions to ZT using density functional theory (DFT) and Boltzmann transport theory [START_REF] Xiang | Suppression of thermopower of Na x CoO 2 by an external magnetic field: Boltzmann transport combined with spin-polarized density functional theory[END_REF] , explaining the suppression of S in large magnetic fields [START_REF] Wang | Spin entropy as the likely source of enhanced thermopower in Na x Co 2 O 4[END_REF] . The phonon dispersion for stoichiometric NaCoO 2 has been calculated using DFT (refs 17,18) and the results are broadly consistent with Raman [START_REF] Iliev | Raman phonons and ageing-related surface disorder in Na x CoO 2[END_REF][START_REF] Lemmens | Effect of Na content and hydration on the excitation spectrum of the cobaltite Na x CoO 2 • yH 2 O[END_REF][START_REF] Shi | Raman spectroscopy study of Na x CoO 2 and superconducting Na x CoO 2 • yH 2 O[END_REF] , infrared [START_REF] Wang | Infrared probe of the electronic structure and charge dynamics of Na 0.7 CoO 2[END_REF] and IXS (ref. 23) experiments. However, thermoelectric Na x CoO 2 contains multi-vacancy clusters that order to form longranged superstructures [START_REF] Roger | Patterning of sodium ions and the control of electrons in sodium cobaltate[END_REF] (Fig. 1a) and those calculations omit the low-energy rattling modes caused by the multi-vacancy clusters.

Here we have included the multi-vacancy clusters in the calculation of the phonon dispersion and compared it with measurements of momentum-resolved IXS and INS spectra from single-crystal Na 0.8 CoO 2 with a well-defined superstructure. We are able to unambiguously identify an Einstein-like rattler mode at low energy and, as a result, to quantitatively account for the suppression of κ for this class of materials.

IXS measurements were performed on a single crystal of Na 0.8 CoO 2 of dimension a few hundred micrometres using the ID28 spectrometer at the European Synchrotron Radiation Facility. The sample was cooled to a temperature T = 200 K, where it was found to be predominantly in the square lattice tri-vacancy phase [START_REF] Roger | Patterning of sodium ions and the control of electrons in sodium cobaltate[END_REF] . We also performed measurements at elevated temperature, at T = 320 K in the disordered stripe phase [START_REF] Morris | Crystal-to-stripe reordering of sodium ions in Na x CoO 2 (x ≥ 0.75)[END_REF] , and at T = 400 K in the disordered phase. High-resolution INS measurements were performed on a larger crystal using the triple-axis spectrometer IN8 at the Institut Laue-Langevin, and phonon lifetimes were determined as a function of temperature in the square tri-vacancy phase.

Figure 1 shows the X-ray intensity as a function of energy transfer with a fixed wavevector transfer, Q, at T = 200 K. The solid lines show the results of first-principles lattice dynamics calculations using the CASTEP code [START_REF] Clark | First principles methods using CASTEP[END_REF] , for the large-period square superstructure determined by X-ray diffraction and for a lattice of NaCoO 2 with no vacancies. The latter calculation used the primitive cell of space group P6 3 /mmc with CoO 2 layers of edge-sharing CoO 6 octahedra and two Wyckoff sites, 2b and 2d, available for Na, and we find consistency with previous calculations for this structure [START_REF] Li | First-principles lattice dynamics of NaCoO 2[END_REF][START_REF] Jha | Phonon properties of the intrinsic insulating phase of the cobalt oxide superconductor NaCoO 2[END_REF] . It is clear that the calculated spectrum of the square phase differs markedly from that of the stoichiometric compound at energies below E ∼ 20 meV; in particular there is an additional narrow mode at E ∼ 13 meV in the superstructure that is completely absent in stoichiometric NaCoO 2 .

Scans of variable energy transfer were performed at a number of settings of Q along the high-symmetry directions → K → M from (111) and for → M from (100), and the fitted peak positions are shown in Fig. 2 as red hexagons. The calculated scattering intensity is shown as colour contour diagrams for the square superstructure and the ideal compound NaCoO 2 . The agreement between the ab initio calculations for the square superstructure and the experimental data is remarkable. The mode at E ∼ 13 meV is captured by the model, and a typical scan is shown in Fig. 1b. In contrast, the ideal vacancy-free compound has no modes below E ∼ 20 meV away from the zone centre.

It is possible to trace the mode at E ∼ 13 meV back to the zone centre, and the relative amplitudes of the ionic displacements are shown in Fig. 1a. The blue spheres in Fig. 1a are the Na 2d sites, and these are fully occupied in the stoichiometric compound. The red spheres are on the inequivalent Na 2b sites and, because they are associated with clusters of three vacancies, they are expected to be more loosely bound. The atomic displacement pattern is

NATURE MATERIALS | VOL 12 | NOVEMBER 2013 | www.nature.com/naturematerials
characterized by a predominant in-plane amplitude of vibration of the Na 2b site with only very minor or negligible contribution from other sites. Thus, we find that the Na 2b sites inside tri-vacancy clusters sandwiched between CoO 2 layers are rattler sites, and the mode at E ∼ 13 meV is a rattling mode.

The Na 2b ions also have large-amplitude vibrations in other modes (Fig. 3a). The influence of the rattling mode on the rest of the dispersion can be investigated theoretically by replacing the Na 2b ions of atomic mass 23 with fictitious heavier ions of atomic mass 40. The calculated density of states for these two cases is shown in Fig. 3b. The largest differences between the two curves are found for energies below E ∼ 40 meV and, therefore, the rattling of the Na 2b ions in the tri-vacancy clusters is seen to affect mainly these low-energy modes. The density of states is also shown for NaCoO 2 , and it is clear that the transfer of density from sharp peaks above E ∼ 20 meV to a broad range at low energy will affect κ.

The atomic displacement parameters for Na are larger for Na 0.8 CoO 2 than for NaCoO 2 , confirming the larger amplitude of vibration in the superstructure (Supplementary Fig. S1). Frozen phonon calculations were performed for the rattler mode, and the potential energy is plotted as a function of displacement in Fig. 3c. The crystallographic asymmetry of the Na 2b site gives rise to a cubic term in energy and consequently a significant degree of anharmonicity. This contrasts with symmetric Ce guest rattler modes in skutterudites where the leading term is quartic and the anharmonicity smaller [START_REF] Feldman | Lattice dynamics and reduced thermal conductivity of filled skutterudites[END_REF] . The energy of this mode is E ∼ 13 meV and, therefore, anharmonicity is expected to become important above T ∼ 150 K.

Measurements of κ using a Quantum Design physical property measurement system on a crystal of the same composition are presented in Fig. 4a. The magnitude of κ, in the region of a few watts per Kelvin-metre, is similar to the values reported for crystals of comparable composition [START_REF] Foo | Charge ordering, commensurability, and metallicity in the phase diagram of the layered Na x CoO 2[END_REF][START_REF] Lee | Large enhancement of the thermopower in Na x CoO 2 at high Na doping[END_REF] . The temperature dependence indicates a contribution from phonon scattering.

To measure phonon linewidths at higher resolution, we have measured a region of the dispersion dominated by a single optical phonon at Q ∼ (1.01,1.01,1) and E ∼ 8 meV using a triple-axis neutron spectrometer. The energy line shape is resolution-limited below T ∼ 80 K, and above this temperature we fitted the data using a Lorentzian convolved with the instrumental resolution function. Figure 4b shows a typical energy scan and the phonon lifetimes, τ , related to the energy linewidths, ph , by τ ∼ h/ ph , as a function of temperature. These data are consistent with our measurements averaged over several modes at T ∼ 200 K using IXS. They are narrow in comparison with previous Raman [START_REF] Shi | Raman spectroscopy study of Na x CoO 2 and superconducting Na x CoO 2 • yH 2 O[END_REF] and IXS (ref. 23) measurements, presumably owing to the high crystalline quality of our sample. There is no evidence for the marked energy broadening expected for a phonon glass.

The lattice contribution to κ was estimated using the expression

κ L = j,Q v 2 jx (Q) τ j (Q,T ) c j (Q,T )
where v jx is the group velocity along the thermal gradient, c j is the heat capacity of the jth phonon mode, and τ j is the phonon lifetime. The values of v jx and c j were determined using the calculated dispersions assuming Bose statistics, and the average value of τ at each temperature was obtained from Fig. 4b. The estimated value of κ L calculated using our model for the square superstructure agrees with the data within experimental uncertainty. The electronic contribution to κ estimated using the Wiedemann-Franz law from the electrical conductivity [START_REF] Lee | Large enhancement of the thermopower in Na x CoO 2 at high Na doping[END_REF] is an order of magnitude lower.

For comparison, we have estimated κ L at T ∼ 200 K for the stoichiometric compound, assuming (artificially) that τ is the same (blue circle in Fig. 4a). The fact that this value is roughly six Significant differences occur below E ∼ 40 meV suggesting that these low-energy modes are strongly affected by the rattling ions. The calculation for NaCoO 2 shows that density has been transferred from sharp peaks above E ∼ 20 meV towards lower-energy modes that contribute strongly to the thermal conductivity. c, Frozen phonon calculation for the rattling mode (black). The displacement is scaled to normal mode amplitude, which corresponds to a maximum in-plane displacement of the Na 2b site of 0.3 Å. The asymmetry of the potential energy as a function of displacement of the rattling ions clearly demonstrates the anharmonicity of this mode.

times higher clearly indicates that rattling modes are important in lowering κ, because they contribute significantly to the specific heat but have small propagation velocities. The marked increase in κ for the stoichiometric compound is borne out experimentally [START_REF] Lee | Large enhancement of the thermopower in Na x CoO 2 at high Na doping[END_REF] , where there is an additional factor of two increase over our calculated value for NaCoO 2 . This could be attributed to weaker anharmonicity due to the absence of rattling modes, reduced electron-phonon scattering because it is an insulator, or fewer defects in a vacancy-free sample. However, it is clear that the changes to the dispersion from the rattling modes give the dominant contribution to the suppression of κ.

Multi-vacancy clusters exist over a wide range of concentrations 26 and, therefore, it is highly likely that they are also responsible for rattling modes throughout the range of compositions relevant to thermoelectric applications [START_REF] Lee | Large enhancement of the thermopower in Na x CoO 2 at high Na doping[END_REF] . A non-dispersive mode was also observed at a similar energy to our rattling mode in a previous IXS study of Na 0.7 CoO 2 focused on the softening of Co-related optical branches near the zone boundary and implications for the unconventional superconducting state [START_REF] Rueff | Phonon softening in Na x CoO 2 • yH 2 O: Implications for the Fermi surface topology and the superconducting state[END_REF] (Supplementary Fig. S2). The squashing mode for CoO 6 octahedra from E ∼ 30-50 meV along → M seen in Fig. 2 is ideal for focusing on electronic properties (Supplementary Fig. S3). The low-energy region of the spectrum containing rattling modes is not very sensitive to the hole concentration in the Co layers.

We have performed IXS measurements at T ∼ 320 K in the disordered stripe phase [START_REF] Morris | Crystal-to-stripe reordering of sodium ions in Na x CoO 2 (x ≥ 0.75)[END_REF] and at T ∼ 400 K, where long-range ordering of the superstructure has disappeared. The absence of long-range order prevents ab initio calculations, but we are able to draw some qualitative conclusions. There is renormalization of the phonon dispersion to slightly lower energies and broadening of the energy linewidth by about 50%. It is not clear whether this is due to the increase in disorder, or simply the lattice expansion and additional phonon scattering. The rattling mode persists at high temperature despite the disorder due to the continued presence of multi-vacancy clusters (Supplementary Fig. S4). Thus, the rattling modes are still important at the elevated temperatures relevant to power recovery applications.

Combining experimental studies of phonon dispersions with DFT provides new insights into the factors that control thermal conductivity. In Na x CoO 2 the phonon scattering is mainly intrinsic, and a further challenge would be to calculate the contributions from anharmonicity and electron-phonon scattering to predict phonon lifetimes. These spectroscopic studies directly link first-principles calculations to the design of better thermoelectric materials.

Methods

Experiment. The IXS crystal of Na 0.8 CoO 2 was grown using the floating-zone technique at Royal Holloway. A single crystal of size 400×400×20 µm 3 was cleaved from the boule and studied using X-rays at the Mo K-edge at room temperature using an Agilent Xcalibur X-ray diffractometer, and large volumes of reciprocal space were surveyed with a CCD (charge-coupled device) detector. Roughly 80% of the sample was found to be in the so-called disordered stripe phase [START_REF] Wang | Infrared probe of the electronic structure and charge dynamics of Na 0.7 CoO 2[END_REF] . The sample was attached to a glass capillary with silicone grease, and a thin coating of the grease covered the sample to prevent reaction with the atmosphere.

The IXS studies were performed using the ID28 spectrometer at the European Synchrotron Radiation Facility in Grenoble. The sample was placed in a stream of nitrogen gas from a 700 series Oxford Cryostream temperature controller, and at each temperature, the superstructure was determined using a CCD detector. IXS measurements were performed at T = 200, 320 and 400 K, where the sample was found to be predominantly in the square [START_REF] Nolas | Effect of partial void filling on the lattice thermal conductivity of skutterudites[END_REF] , disordered stripe [START_REF] Wang | Infrared probe of the electronic structure and charge dynamics of Na 0.7 CoO 2[END_REF] and disordered phases, respectively. The sample was aligned accurately using a point detector in two settings with (h0l) and (hhl) crystallographic planes horizontal. The lattice parameters were a = 2.844(2) and c = 10.803 (7) for the P6 3 /mmc space group at T = 200 K, and the mosaic spread was about 0.5 • .

A resolution in energy transfer (full-width at half-maximum) of 3.0 meV was achieved using the Si(999) reflection for both the monochromator and analyser crystals in backscattering geometry with a Bragg angle of 89.98 • . The monochromator beam is focused by a toroidal mirror into a spot size of 250 × 90 µm 2 (horizontal × vertical) full-width at half-maximum at the sample position. Energy transfers in the range -20 to 90 meV were measured at fixed Q by varying the incident energy through the temperature variation of the monochromator crystal, and keeping the final energy fixed. The use of nine spherical analyser crystals allowed simultaneous measurements at multiple values of Q. One of the analysers was chosen to give Q along high-symmetry directions.

A larger fragment from the same boule of size 7 × 2 × 0.5 mm 3 , with the short direction parallel to the hexagonal c-direction, was measured in a Quantum Design physical property measurement system. Gold contacts were bonded to the sample using silver-loaded epoxy adhesive, and the in-plane κ was measured over the temperature range T ∼ 2-300 K.

A larger crystal for INS of size 35 × 10 × 3 mm 3 and of mosaic spread 3 • was screened using the neutron Laue diffractometer SXD at ISIS, and it was found to be a single grain and to adopt the square phase below T ∼ 280 K. It was sealed inside an aluminium can with a helium atmosphere, and mounted in a cryostat with the (hhl) plane horizontal on the triple-axis spectrometer IN8 at the Institut Laue-Langevin. A fixed final energy of 14.7 meV was obtained using a graphite analyser and higher-order contamination was suppressed with a graphite filter after the sample. Energy scans were performed by varying the energy selected by a PG(002) analyser/monochromator, and source-monochromator, monochromator-sample, sample-analyser, analyser-detector collimation 30', 20', 30', 30' gave an instrumental energy resolution of full-width at half-maximum ∼1.45 meV for the dispersion at E ∼ 8 meV and Q ∼ (1.01,1.01,1).

Theory. The phonon calculations used DFT in the plane-wave pseudopotential approach as implemented in the CASTEP code [START_REF] Clark | First principles methods using CASTEP[END_REF] and the spin-polarized Perdew-Burke-Ernzerhof generalized gradient approximation to exchange and correlation [START_REF] Perdew | Generalized gradient approximation made simple[END_REF] . Vanderbilt-type ultrasoft pseudopotentials [START_REF] Vanderbilt | Soft self-consistent pseudopotentials in a generalized eigenvalue formalism[END_REF] CASTEP with small core radii to give very high accuracy and transferability at a plane-wave cutoff of 500 eV. Monkhorst-Pack meshes were used to sample the electronic Brillouin zone at a resolution of 0.025 Å -1 or finer with a Gaussian smearing of 0.1 eV for the metallic square phase Na 0.8 CoO 2 calculations. This resulted in a maximum error in forces of |δF | < 0.0024 eV Å -1 . Lattice dynamics was computed using the ab initio finite displacement/supercell method [START_REF] Frank | Ab initio force-constant method for phonon dispersions in alkali metals[END_REF] . A single 10.276 × 10.276 × 10.8 Å 3 cell of the square phase structure containing 114 atoms was sufficient to compute S (Q, ω) across the full Brillouin zone of the primitive cell in Fig. 2. Reference calculations excluding the rattler mode were performed in the primitive cell containing 2 formula units of NaCoO 2 , using a 5×5×2 supercell for interpolation across the Brillouin zone.

Figure 1 |

 1 Figure 1 | Einstein-like rattling modes in sodium cobaltate. a, The supercell of Na 0.80 CoO 2 showing the rattling mode at E ∼ 13 meV comprising mainly displacements of the (red) Na 2b ions inside tri-vacancy clusters. Na 2d ions are blue, Co is green and O is gold. The length of the arrow indicates the relative amplitude for each ion, not an absolute displacement. b, IXS at Q = (1.25,1.25,1) measured from Na 0.80 CoO 2 at T = 200 K, focusing on the rattling mode in a. The red line shows the CASTEP calculation for the Na 0.80 CoO 2 superstructure comprising a square array of tri-vacancy clusters, and the blue line shows the corresponding calculation for NaCoO 2 . c, IXS at Q = (1.17,0,0) from Na 0.80 CoO 2 at T = 200 K showing typical data for the full spectrum. Error bars are determined by Poisson statistics.

Figure 2 |

 2 Figure 2 | Phonon dispersion of sodium cobaltate. a,b, CASTEP calculations for the Na 0.8 CoO 2 superstructure comprising a square array of tri-vacancy clusters (a) and for NaCoO 2 (b), shown as a colour contour diagram on a hot scale (blue, low intensity; white, high intensity). The experimental data obtained for Na 0.8 CoO 2 at T = 200 K are shown as red hexagons. The superstructure is essential to capture the scattering at lower energies below E ∼ 20 meV.

Figure 3 |

 3 Figure3| Influence of rattlers on the rest of the dispersion. a, A typical zone-centre optical phonon, showing how the rattling ions couple strongly to and alter modes. b, The calculated density of states (g(E)) for Na rattlers of mass (m r ) 23 (red) and for test atoms of mass 40 AMU (black). Significant differences occur below E ∼ 40 meV suggesting that these low-energy modes are strongly affected by the rattling ions. The calculation for NaCoO 2 shows that density has been transferred from sharp peaks above E ∼ 20 meV towards lower-energy modes that contribute strongly to the thermal conductivity. c, Frozen phonon calculation for the rattling mode (black). The displacement is scaled to normal mode amplitude, which corresponds to a maximum in-plane displacement of the Na 2b site of 0.3 Å. The asymmetry of the potential energy as a function of displacement of the rattling ions clearly demonstrates the anharmonicity of this mode.

Figure 4 |

 4 Figure 4 | Thermal properties of sodium cobaltate. a, The measured thermal conductivity of Na 0.8 CoO 2 (black) is compared with the calculated lattice contribution to the thermal conductivity for the Na 0.8 CoO 2 superstructure comprising a square array of tri-vacancy clusters using the phonon lifetimes derived using neutrons (red) and X-rays (yellow). The calculated values for NaCoO 2 (blue) are a factor of six higher. b, The phonon lifetimes determined from the widths of strong single modes. The lifetimes obtained using INS (red) are consistent with the measurement using IXS (yellow). The solid black line shows the fitted inverse temperature dependence. A typical INS scan at fixed Q = (1.01,1.01,1) at T ∼ 275 K is shown in the inset. The red line shows the best fit of a Lorentzian convolved with the resolution Gaussian, and the dashed black line shows the best fit of the resolution Gaussian. Error bars in a and b are determined by the fits to phonon lifetimes, and in the inset by Poisson statistics.

  were generated within NATURE MATERIALS | VOL 12 | NOVEMBER 2013 | www.nature.com/naturematerials LETTERS NATURE MATERIALS DOI: 10.1038/NMAT3739

Received 15

 15 November 2012; accepted 24 July 2013; published online 25 August 2013

© 2013 Macmillan Publishers Limited. All rights reserved

NATURE MATERIALS | VOL 12 | NOVEMBER 2013 | www.nature.com/naturematerials © 2013 Macmillan Publishers Limited. All rights reserved

Acknowledgements

We thank the EPSRC for support through grants EP/J011150/1 and EP/J012912/1, and for use of the UK national supercomputing facility HECToR under grant EP/F036809/1. Other calculations used STFC's e-Science facility. We wish to acknowledge the assistance of R. Hanson in developing the crystal structure and normal mode visualizations, which were created using the Jmol software (http://www.jmol.org/).

Author contributions

Initial planning was by J.P.G., K.R., M.R., A.T.B., A.B. and M.H. The IXS experiment was performed by D.J.V., J.P.G., E.B., L.G. and M.K., and INS measurements were performed by D.J.V., J.P.G., A.P. and M.E. Single crystals were grown by S.U. The thermal conductivity was measured by D.G.P. and the structural characterization was by D.J.V., D.G.P., E.C. and M.J.G. Theoretical modelling was by K.R. and D.J.V. The manuscript was drafted by J.P.G, D.J.V. and K.R. and all authors participated in the writing and review of the final draft.

Additional information

Supplementary information is available in the online version of the paper. Reprints and permissions information is available online at www.nature.com/reprints. Correspondence and requests for materials should be addressed to J.P.G.

Competing financial interests

The authors declare no competing financial interests.