# DYNAMICS OF ELASTIC EXCITABLE MEDIA

## 1 Introduction

The Burridge–Knopoff model \[Burridge & Knopoff, 1967\] mimics the interaction of two plates in a geological fault as a chain of blocks elastically coupled together and to one of the plates, and subject to a friction force by the surface of the other plate, such that they perform stick–slip motions — Fig. 1a. This simple system reproduces some statistical features of real earthquakes \[Carlson & Langer, 1989\] such as the Gutenberg–Richter power-law distribution \[Gutenberg & Richter, 1956\], considered an example self-organized criticality, \[Carlson et al., 1994\] which obviously involves a large number of degrees of freedom. However, recent laboratory experiments \[Rubio & Galeano, 1994\] that attempt to reproduce these dynamics in a real stick–slip dynamical system consisting of an elastic gel sliding around a metallic cylinder, have shown that low-dimensional phenomena are more robust in reality, and, although proven to be unstable in the Burridge–Knopoff model, do show up in the laboratory.

In a different realm, excitable media are usually studied using the model of van der Pol, FitzHugh, and Nagumo \[van der Pol & van der Mark, 1928, FitzHugh, 1960, FitzHugh, 1961, Nagumo et al., 1962\]. This model normally includes only diffusive coupling. Originally from physiology and chemistry, excitable media have also captured the attention of physicists and mathematicians working in the area of nonlinear science because of the apparent universality of many features of their complex spatiotemporal properties \[Meron, 1992\].

We have shown \[Cartwright et al., 1997\] that a Burridge–Knopoff model with a lubricated creep–slip friction force law showing viscous properties at both the low and high velocity limits (Fig. 1b) is a type of van der Pol–FitzHugh–Nagumo excitable medium in which the local interaction is elastic rather than diffusive. We have investigated the dynamics of the model and have shown that its behaviour is dominated by low-dimensional structures, including global oscillations and propagating fronts. Here we investigate further the dynamics of elastic excitable media, and focus on the behaviour of the fronts.

## 2 The Model

Our elastic excitable medium model may be written \[Cartwright et al., 1997\]

$$\frac{^2\chi }{t^2}=c^2\frac{^2\chi }{x^2}(\chi \nu t)\gamma \varphi \left(\frac{\chi }{t}\right),$$
(1)   
where, in the language of frictional sliding, $`\chi (x,t)`$ represents the time-dependent local longitudinal deformation of the surface of the upper plate in the static reference frame of the lower plate, $`\varphi (\chi /t)=(\chi /t)^3/3\chi /t`$ is the friction function, as the dashed line in Fig. 1b, $`\gamma `$ measures the magnitude of the friction, $`c`$ is the longitudinal speed of sound, and $`\nu `$ represents the pulling velocity or slip rate. Compare this with the discrete Burridge–Knopoff model from whence Eq. (1) may be derived in the continuum limit

$$m\frac{d^2x_i}{dt^2}=k_c(x_{i+1}2x_i+x_{i1})k_p(x_iVt)F_f\left(\frac{dx_i}{dt}\right),$$
(2)   
where $`x_i`$ is the departure of block $`i`$ from its equilibrium position. It has been noted \[Carlson et al., 1994\] that in some cases discrete Burridge–Knopoff models fail to attain a well-defined continuum limit. This is not a consequence of any numerical instability in computer simulations, but of the hyperbolic nature of the equation and of the shape of the friction force commonly used. In such cases, Eq. (1) should be considered as a symbolic representation of the well-defined discrete dynamics of Eq. (2). From Eq. (1) we can obtain an expression for the local velocity $`\psi =\chi /t`$ of the interface that gives us the model written as a couple of differential equations of first order in time

$`{\displaystyle \frac{\psi }{t}}`$ $`=`$ $`\gamma (\eta \varphi (\psi )),`$ (3)   
$`{\displaystyle \frac{\eta }{t}}`$ $`=`$ $`{\displaystyle \frac{1}{\gamma }}\left(\psi \nu c^2{\displaystyle \frac{^2\psi }{x^2}}\right).`$ (4)   
## 3 Front Dynamics

The dynamical behaviour of the elastic excitable medium model has been reported in detail in a previous paper \[Cartwright et al., 1997\]. It is notable that much of the chaotic behaviour shown by the discrete versions of the Burridge–Knopoff models built on the basis of a monotonic velocity weakening friction law \[Rice, 1993, Schmittbuhl et al., 1993, Xu & Knopoff, 1994\] becomes more organized in our case. Global oscillations and propagating fronts dominate a large proportion of the relevant parameter space. Moreover, the global oscillations show interesting instability mechanisms leading to the appearance of the propagating fronts. Among these instabilities, the most interesting is perhaps the occurrence of a period-doubling bifurcation at a finite spatial wavelength. As a result of this bifurcation, the globally synchronized oscillatory medium breaks into a finite number of equidistant pacemaker zones from which pairs of counterpropagating fronts are emitted. Fronts emitted from neighbouring pacemakers annihilate upon collision and the annihilation point then becomes an emitting centre for the next generation of fronts; the whole process repeats after two iterations. These results are graphically summarized by Fig. 2, where the transient evolution of a nearly uniform initial state is shown for $`\nu `$ just above the period-doubling instability. In the first stages of the evolution the dynamics is dominated by synchronized global oscillations. The instability then grows, giving rise during a certain time interval to a transient period-doubled structure. Finally this structure decays into a set of propagating fronts (or propagating pulses, since a pair of neighbouring fronts can be considered a pulse). This bifurcation was discussed in more detail in a previous work \[Cartwright et al., 1997\].

Here we focus on the properties of the propagating front regime with special emphasis on the selection mechanisms for the front velocity and spatial configuration. We suppose a solution of the type $`\psi (x,t)=f(\stackrel{~}{z})`$, where $`\stackrel{~}{z}=x/v+t`$, and $`v`$ is the front velocity. This together with the further rescaling $`z=\stackrel{~}{z}/\sqrt{1c^2/v^2}`$ leads to

$$\frac{d^2f}{dz^2}+\mu (f^21)\frac{df}{dz}+f=\nu ,$$
(5)   
which is the van der Pol equation with the nonlinearity rescaled by $`\mu =\gamma /\sqrt{1c^2/v^2}`$. The propagating fronts are then periodic solutions of the van der Pol equation. The parameter $`\mu `$ is undefined until the value of the front velocity $`v`$ is chosen. However, we know that the period of the solution is a function $`T=T(\mu )`$ of $`\mu `$: in the limit of large $`\mu `$, $`T`$ behaves as $`T=k\mu +O(\mu ^1)`$, where $`k=3+(\nu ^21)\mathrm{ln}[(4\nu ^2)/(1\nu ^2)]`$ \[Feingold et al., 1988\]. Since this period should be commensurate with the system size $`S`$, we have the condition $`nT(\mu (v))=S/(v\sqrt{1c^2/v^2})`$, where $`n`$ is an integer, to select the allowed front velocities, which in the large $`\mu `$ limit gives us the quantizing condition $`v=S/(nk\gamma )`$. The integer $`n`$ can be interpreted as the total number of pulses propagating in the system. Because Eq. (5) has bounded solutions only if $`v^2>c^2`$, the propagating fronts are supersonic. Since our analytical estimations predict that $`v`$ becomes smaller than $`c`$ and approaches zero as the number of pulses in the system is increased, there should be a maximum number of pulses allowed in the system. However, numerical simulations show that it is possible to find solutions composed of an arbitrary number of propagating pulses with a proper choice of the initial conditions. A question then arises as to what is the long term behaviour of these solutions when the number of pulses exceeds the maximum allowed by the restrictions on the front velocity. There are three simple scenarios logically compatible with the analysis:

1. Some of the pulses are annihilated by the dynamics.
2. Since the configuration cannot propagate rigidly at a compatible uniform speed, nondecaying fluctuations of the velocity of individual pulses should be observed.
3. The system evolves into a uniformly propagating solution with speed approaching the singular limit $`v=c`$. In this limit the pulses become discontinuous, so that Eq. (1) becomes ill-defined. The original system of blocks and springs is thus no longer well represented by the partial differential equation Eq. (1), to which our analytic estimations pertain, but the discrete effects present in Eq. (2) prevail.

The first scenario occurs, but only during a transient phase: at long times states with arbitrary numbers of pulses are reached and maintained in time. Examples of how the second scenario could be realized will be briefly discussed below. We first focus on the third scenario, for which we find abundant numerical evidence. Figure 3 summarizes this. In upper frame of the first column we see a propagating solution composed of only one pulse travelling around the system to the right. Using the spatial coordinate $`x`$ as a parameter we have drawn below the phase-space diagram of such a solution. This diagram shows geometrically that the spatial dependence of the solution is accurately described as the periodic solution of a van der Pol–FitzHugh–Nagumo equation. The second column shows a two-pulse solution propagating to the left at a speed approximately one half that of the single-pulse solution, and thus closer to the sound velocity $`c`$, as expected from the analytical estimations. The phase-space projection below shows general agreement with the prediction of Eq. (5), although a small discrepancy is already visible. Finally, the third column in Fig. 3 shows a breakdown of the continuum model. A train of six pulses travelling at constant speed has been generated, which is forbidden within the framework of the continuum analysis. They travel to the left at speed essentially $`c`$. Notice the sharp gradients at both the leading and trailing edges of the pulses. If we refine the numerical spatial discretization, this just increases the gradients, so that we may expect discontinuities in the continuum limit. The discrete model Eq. (2) introduced in the computer should better represent the physical system of springs and blocks than the partial differential equation Eq. (1). The phase-space diagram clearly shows a strong departure from the solutions of the van der Pol–FitzHugh–Nagumo equations. A well-defined phase-space structure appears, which should be understood as a property of the purely discrete system, different in this case from the van der Pol–FitzHugh–Nagumo phase-space structure. Discreteness effects have been studied in several models with velocity weakening friction \[Rice, 1993, Schmittbuhl et al., 1993, Xu & Knopoff, 1994, Galeano et al., 1998\].

When $`\nu `$ has been taken very far from the slipping threshold the behaviour locally is of autonomous relaxation oscillations, and a slightly different kind of propagating front dynamics occurs. By setting $`\nu =0`$ we recover the symmetric form of the van der Pol equation, supplemented in our case with the elastic spatial term. Now a Floquet stability analysis of the homogeneous oscillatory state shows the development of a long wavelength instability as $`\nu `$ decreases. The evolution of this instability is displayed in Fig. 4. One can roughly describe this behaviour as the formation of spatial regions where the medium undergoes almost synchronous relaxation oscillations. Neighbouring regions, however, oscillate in antiphase and a phase-change front defining the border between them travels around the system. On the left we can see some of these fronts moving at different velocities during a transient stage, while on the right is displayed an asymptotic state where the fronts have reached an equilibrium configuration in which the phase distribution travels rigidly. This behaviour is reminiscent of the phase dynamics found in the complex Ginzburg–Landau equation \[Montagne et al., 1996, Montagne et al., 1997\].

A closer look at Fig. 4 indicates that the spatial regions that we describe as synchronized, are in fact regions where the excitation or relaxation propagates extremely fast. This is visualized in the picture as a very slight tilt of the bands representing the oscillation. Notice that the projection of these bands shrinks considerably near the places where the phase jumps. By taking spacelike slices of these pictures we obtain the instantaneous configuration of the excitation field. Such configurations look like travelling pulses varying wildly in size and velocity. This behaviour can be considered an extreme manifestation of scenario 2 above.

## 4 Discussion

A lubricated friction law in the Burridge–Knopoff model can be justified both by theory \[Persson, 1995, Persson, 1997\], and experiments \[Heslot et al., 1994, Brechet & Estrin, 1994, Kilgore et al., 1993, Demirel & Granick, 1996, Budakian et al., 1998\]. Moreover, studies of peeling adhesive tape \[Hong & Yue, 1995\], of Saffman–Taylor fracture in viscous fingering \[Kurtze & Hong, 1993\], and of the Portevin–Le Châtelier effect \[Kubin & Estrin, 1985, Lebyodkin et al., 1995\] lead to the same form of friction law as we use here. Certainly our model well represents the qualitative characteristics of the laboratory stick–slip dynamics experiments \[Rubio & Galeano, 1994\] referred to above, and might have relevance to the present debate on shear stress and friction in real geological faults \[Sleep & Blanpied, 1992, Melosh, 1996, Cohen, 1996\]. In the context of excitable systems elastic coupling has been left aside, because in the chemical and biological systems studied up to now the coupling is diffusive. However, an elastic excitable medium can be realized as an active transmission line or optical waveguide \[Cartwright et al., 1997\].

## Acknowledgements

The authors acknowledge the financial support of the Spanish Dirección General de Investigación Científica y Técnica, contracts PB94-1167 and PB94-1172.