## Abstract

We consider a FitzHugh–Nagumo system of equations where the traditional diffusion terms are replaced with linear cross-diffusion of components. This system describes solitary waves that have unusual form and are capable of quasi-soliton interaction. This is different from the classical FitzHugh–Nagumo system with self-diffusion, but similar to a predator–prey model with taxis of populations on each other's gradient which we considered earlier. We study these waves by numerical simulations and also present an analytical theory, based on the asymptotic behaviour which arises when the local dynamics of the inhibitor field are much slower than those of the activator field.

## 1. Introduction

Recently we have reported a new class of nonlinear waves, which have some properties similar to autowaves in excitable systems, and some to solitons (Tsyganov *et al*. 2003). These waves were observed in a system of partial differential equations describing two spatially distributed populations in a ‘predator–prey’ relationship with each other. The spatial evolution is governed by three processes, positive taxis of predators up the gradient of prey (pursuit) and negative taxis of prey down the gradient of predators (evasion), yielding nonlinear ‘cross-diffusion’ terms, and random motion of both species, the usual ‘self-diffusion’. We considered the problem in one spatial dimension, *x*, using the equations(1.1)where the local predator–prey interactions are as suggested by Truscott & Brindley (1994). Here the spatial interactions are described by diffusion with coefficient *D*, which for simplicity was considered constant, uniform and equal for both species, and taxis terms and , where *h*_{−} was the coefficient of negative taxis of *P* on the gradient of *Z*, *h*_{+} was the coefficient of positive taxis of *Z* on the gradient of *P* (Keller & Segel 1971; Murray 2003).

We demonstrated a new class of propagating waves in the system (1.1). The mechanism of propagation of these pursuit-evasion taxis waves essentially depends on the taxis of the species and is entirely different from waves in reaction–diffusion systems. These new waves are similar to reaction–diffusion waves in that their eventual stationary form and amplitude do not depend on initial conditions, and are different from them in that establishment of this stationary form is very slow, and that these new waves often penetrate through each other and reflect from impermeable boundaries. This latter property is similar to solitons in conservative systems (Tsyganov *et al*. 2003).

An explanation of this quasi-soliton behaviour has been suggested in terms of the taxis of species to each other's gradients (Biktashev *et al*. 2004), which is mathematically described by nonlinear terms. This explanation restricts this phenomenon to population dynamics. However, as noted in (Tsyganov *et al*. 2003), certain aspects of the unusual behaviour of those waves, such as spatially oscillatory fronts, can be understood in terms of linear terms, which are akin to cross-diffusion of species. Cross-diffusion terms are found in mathematical models of various processes, and have been discussed not only in population dynamics context (Shigesada *et al*. 1979; Okubo 1980; Kuznetsov *et al*. 1994; del Castillo Negrete *et al*. 2002; Murray 2003), but also occur in Ginzburg–Landau type equations near critical transition points in systems of various physical origins (Haken 1978; Cross & Hohenberg 1993), notably oscillatory chemical reactions (Kuramoto & Tsuzuki 1975; Kuramoto 1984) and asymptotic description of water waves (Johnson 1976, 1997), and specifically in description of such varied processes as chemical and biological pattern formation (Almirantis & Papageorgiou 1991), movement of tectonic plates (Cartwright *et al*. 1997) and dynamics of electrolytic solutions (Jorné 1975).

Thus the motivation of the present work is to find out if cross-diffusion interaction of nonlinear fields can be enough to produce quasi-soliton behaviour. We choose the nonlinearity in the form of a FitzHugh–Nagumo system (FitzHugh 1961; Nagumo *et al*. 1962), which is a prototypical excitable system model. However, instead of traditional description of spatial spread of the components by diffusion, we include only cross-diffusion terms:(1.2)We have chosen excitable local kinetics, to allow solitary waves which are methodologically easier to study than periodic waves, and to allow comparison with the vast knowledge about solitary waves in FitzHugh–Nagumo system with self-diffusion. We consider *D*_{u}≥0, *D*_{v}≥0, where *D*_{u}+*D*_{v}>0, so the eigenvalues of the diffusion matrix are , i.e. lie on the imaginary axis. This choice of signs mimicks the pursuit-evasion interaction of predators and prey, as in equation (1.1). These signs also correspond to the Burridge–Knopoff model; Cartwright *et al*. (1997) considered a system very similar to equation (1.2), with *D*_{v}=0, and slightly different local kinetics. This is also similar to Ginzburg–Landau models, in the case of the negligible real part of the matrix element of the diffusion matrix corresponding to the critical mode.

The purpose of this paper is to report new phenomena observed in this system, and to provide some analytical theory. We do not intend to discuss here particular physical or biological applications of these results. In §2 we describe selected numerical simulations. In §3 we present the analytical theory. The results are summarized and discussed in §4.

## 2. Numerical results

### (a) Methods

We did numerical simulations of equation (1.2) for *x*∈[0, *L*], where *L* varied in different numerical experiments. We used no-flux boundary conditions for both variables. When we wanted to study the behaviour of waves in an infinite space, we moved the calculation interval as the wave progressed, ensuring that it never approached closer than 60 to the boundaries, to avoid their influence. Thus, the labeling of the spatial coordinate axes on all figures only represents the spatial scale, not the position, except on figures depicting wave collisions. The value of parameter *a* was *a*=0.3 in all cases, whereas *ϵ*, *D*_{u} and *D*_{v} varied. We used first-order time stepping, explicit in the reaction terms and fully implicit in the cross-diffusion terms, with a second-order central difference approximation for the spatial derivatives. We used steps *h*_{x}=0.2 and *h*_{t}=0.005. The beginnings and ends of excitation waves were defined as the points where *u*(*x*, *t*)=*u*_{f}=0.5. Initial conditions were set as *u*(*x*, 0)=*Θ*(*δ*−*x*), *v*(*x*, 0)=0, to initiate a wave starting from the left end of the domain, and *u*(*x*, 0)=*Θ*(*δ*−*x*)+*Θ*(*x*−(*L*−*δ*)), *v*(*x*, 0)=0 to initiate two waves from both ends. Here *Θ* is the Heaviside function, and the wave seed length was chosen as *δ*=2.

### (b) Evolution of the waveform

Figure 1 illustrates a typical profile of a propagating wave in our excitable cross-diffusion model. We notice the following features:

there is a rather long transient, much longer than a typial duration of the wave itself, during which the the wave grows longer up until a certain stationary length;

the front of a pulse, both during the transient and in the established state, has an oscillatory character;

oscillations of the wave profile around the front and the back of a pulse are observed both in the activator and in the inhibitor components;

during the plateau of the pulse (see

*t*=120 and*t*=620 in figure 1), both activator and inhibitor do not change very much.

All these features are similar to pursuit-evasion waves in predator–prey model that we described in (Tsyganov *et al*. 2004), and very different from waves in FitzHugh–Nagumo systems with the usual diffusion, where

the transient is usually over during a time comparable to the duration of the wave, whereas the wave itself is much longer than in the taxis system with the same kinetics

fronts and backs are monotonic and

involve mostly changes in the activator but not the inhibitor, and

the changes of activator and inhibitor during the plateau phase are substantial.

Let us consider in more detail the evolution of the length of the wave profile *w*, defined as the distance, at any given moment of time, between the points with *u*=*u*_{f} and ∂*u*/∂*t*>0 (the wavefront) and *u*=*u*_{f} and ∂*u*/∂*t*<0 (the waveback). Figure 2 shows that with the initial conditions used, the wave width typically increased gradually until reaching the equilibrium saturation level. Both the rate of the increase of the width and its equilibrium value are affected by all parameters. However, dependence on *ϵ* is most interesting, as at *ϵ*=0 the saturation effect disappears altogether and the wave width seems to go on growing for ever, see figure 2*d*.

Figure 3 illustrates the evolution of the shape of a wave in more detail. This evolution becomes more evident at small values of parameter *ϵ*. The wave becomes longer and, correspondingly, the wavefront and waveback are better separated. During the long transient shown in figure 3, the wavefront remains practically unchanged, as corresponding curves in figure 3*c* overlap so much they are indistinguishable, and all the change, although still relatively small, is observed only in the shape of the waveback.

Considering the movement of fronts and backs of waves separately, we have recorded how their velocities changed with time. The results are illustrated in figure 4. Typically, after a short transient, the wavefront speed was established at a constant value. The waveback speed, on the contrary, started at a much lower value, and gradually increased, reaching the front speed in the long run. Thus the wave profile dynamics at small *ϵ* can be described as the waveback lagging behind the wavefront until it is at such distance that its speed equals the speed of the wavefront. The waveback speed changes more slowly for smaller *ϵ*. At *ϵ*, the speed of the waveback does not change at all, and so the separation between the front and the back goes on growing for ever.

Figure 5 illustrates the dependence of the established wave profile on *ϵ*. At small *ϵ*, the dependence of the established width *w* is well approximated by .

### (c) Reflection/penetration of waves

A remarkable property of pursuit-evasion waves in equation (1.1) was their ability to penetrate through each other or reflect from impermeable boundaries, a ‘quasi-soliton’ behaviour, in a broad range of parameters. We have also observed this property in excitable systems with cross-diffusion. Here we illustrate this phenomenon, in the context of the slow evolution of the wave shape discussed in the previous subsection. Namely, we have found that the ability of the waves to reflect/penetrate depends on their ‘age’, i.e. their width at the moment of the collision of waves with each other or with the boundary. Figure 6 shows the two kinds of interaction at different ages of the waves in systems with the same parameters. Different ages were achieved by varying the geometric size of the system. In figure 6*a*, the medium is small and the waves do not have time to grow wide before they collide. The resulting collision is quasi-soliton. The waves newly emerging after the collision are thin again and, since they do not have much time to grow wide before colliding with the boundaries, their state is then much the same as before the first collision with each other and they are reflected from the boundaries. This process repeats again and again and so we have self-sustained activity in the system. In the bigger medium in figure 6*b*, the waves grow wider before they collide. As a result, they annihilate.

The profiles of waves during selected moments of cases shown in figure 6 are shown in figures 7 and 8. Unlike pursuit-evasion waves in equation (1.1), this time the evolution of the components cannot be interpreted in terms of directed movements of predator and prey populations. We see, however, that the potential for reflection exists in both cases and is related to the tendency of the front and back of the wave to spatial oscillations. This tendency is due to the character of the cross-diffusion terms, which is the same as in Schrödinger equation, and is not found in reaction–diffusion systems with the usual (self) diffusion. The difference between figures 7 and 8 is that, in the second case, the backward wave is too weak and decays instead of developing into a propagating wave.

The transition from the reflecting to the non-reflecting regime occurs as the wave grows wider during propagation. The positions of that transition are indicated in figures 2 and 4. As we have seen in the previous section, the growth of the wave width is related to the dynamics of the waveback, both in terms of its spatial position and of its structure. Thus it appears likely that the result of a collision depends mostly on the properties of the waveback. The interaction between the wavefront and the waveback should be negligible if the wave is very wide. On the other hand, the structure of the back of a wide wave will remain the same as that of a thin wave, i.e. capable of reflection, if the dynamics during the plateau is slow enough, which is determined by *ϵ*. Thus, we conclude that if *ϵ*=0 then the wave should remain quasi-soliton even when it grows very wide. This conclusion is confirmed by simulations, see figure 9.

## 3. Analytical theory

### (a) Asymptotic setting

In this section, we present an asymptotic theory, capturing the main phenomenological features reported in §2, in the limit of small *ϵ*. We consider a generalization of (1.2)(3.1)where *f*(*u*) is an *N*-shaped-graph function, *D*_{u} and *D*_{v} are both assumed positive in this section, and *ϵ* is a small parameter. The asymptotic theory for the limit *ϵ*→0 can be developed similarly to the classical case of excitable waves with self-diffusion, a review of which can be found in (Tyson & Keener 1988). At finite times and distances, *t*, *x*=*O*(1), we observe the wavefronts and wavebacks, which are trigger waves between two asymptotic states, say from (*u*_{1}, *v*_{1}) to (*u*_{2}, *v*_{2}). These asymptotic states are equilibria in fast time, and therefore satisfy *f*(*u*_{j})=*v*_{j}, *j*=1,2. The behaviour at large times and distances, , can be formally described by passing to independent variables *X*=*ϵx*, *T*=*ϵt*. For , we have(3.2)Based on the phenomenology described in §2 and by analogy with (Tyson & Keener 1988), we consider typical solutions that consist of large, , pieces along the slow manifold *v*≈*f*(*u*), corresponding to wave plateaux and the intervals between the waves, separated by fast, *x*, *t*=*O*(1), transition zones, which represent wavefronts and wavebacks. Correspondingly, the fast transition zones are described by the limit *ϵ*→0 in equation (3.1), and the slow pieces are described by the limit *ϵ*→0 in equation (3.2). For brevitiy, in this section we will refer both to wavefronts and to wavebacks as fronts.

### (b) Fast movements and the piecewise caricature

The motion of the fronts is described by the *ϵ*→0 limit of equation (3.1), i.e.(3.3)As the fronts are far from each other, , and change their speed only at large time scales, , in the limit *ϵ*→0 we consider steadily propagating, solitary front solutions, , , where *ξ*=*x*−*ct*−*ξ*_{0} and *ξ*_{0} is an arbitrary constant. For definiteness, we consider a wave propagating to the right, *c*>0. Then the profiles of and satisfy(3.4)The arbirary constant *ξ*_{0} can be fixed by an additional condition(3.5)where *u*_{f} is an appropriate threshold constant corresponding to a selected point on the unstable middle branch of equation *u*=*f*(*v*). This is analogous, but not necessarily identical, to the constant *u*_{f} used in the numerics.

Note that *v* is a conserved quantity of equation (3.3) with flux *D*_{u}∂*u*/∂*x* and, correspondingly, equation (3.4) has a first integral . As , we have *v*_{1}=*v*_{2}=*v*_{*} and then *u*_{1} and *u*_{2} are two stable roots of function *f*(*u*)−*v*_{*}. Using the first integral to eliminate from equation (3.4), we get(3.6)The well-posedness of this problem and the type of beginning and end of the front (oscillatory/monotonic) can be understood from the asymptotics of the solutions of this equation at *ξ*→±∞, which can be obtained by linearizations of at . However, the front solutions themselves are to be obtained numerically. Further analytical progress can be obtained for a piecewise linear function *f*(*u*). We consider a piecewise function similar to one used in Rinzel & Keller (1973) caricature of the FitzHugh–Nagumo system(3.7)where *a*∈(0, 1) and *Θ*( ) is the Heaviside step function. Then(3.8)where *u*_{−}=−*v*_{*}/*α*, *u*_{+}=1−*v*_{*}/*α*, so *u*_{+}>*a*>*u*_{−} as long as *v*_{*}∈(0, *α*(1−*a*)). Then a wavefront corresponds to *u*_{1}=*u*_{−}, *u*_{2}=*u*_{+}, and a waveback to *u*_{1}=*u*_{+}, *u*_{2}=*u*_{−}. Let *u*_{2}=*u*_{1}+*Δ*, *a*=*u*_{1}+*Δθ*. Then a wavefront and a waveback, corresponding to *the same value of v*_{*}, differ from each other by the simultaneous transformations *Δ*→−*Δ* and *θ*→1−*θ*.

As the ‘unstable branch’ of function (3.7) is vertical we choose the threshold *u*_{f}=*a*, so we have . On either side of the threshold, the general solution of equations (3.6) and (3.7) iswhere subscripts 1,2 correspond to the antegrade, *ξ*>0, and retrograde, *ξ*<0, parts of the front, respectively, and *λ*_{j}, *j*=1,2,3 are the roots of the characteristic equation(3.9)The fronts observed in numerics are characterized by oscillatory onset and monotonic end, so we expect *λ*_{1}>0 and *λ*_{2,3}=−*μ*±*ik*, *μ*>0. Then *λ*_{1}=2*μ*, and *μ* and *k* should satisfy(3.10)(3.11)Asymptotics of imply that . The solution of equation (3.6) for a piecewise continuous *f* should be continuous and have two continuous derivatives. This provides three matching conditions at *ξ*=0, leading toSolving this system gives(3.12)Note that this sort of behaviour is only possible when *θ*<8/9. Substituting *k*^{2} from equation (3.12) into the characteristic equation (3.9) leads to(3.13)(3.14)Excluding *μ* from this system, we get an equation for *c*,(3.15)(note that the choice of *c*>0 guarantees *μ*>0 in equation (3.14)), where(3.16)and *Χ*(*β*) is the positive root of the cubic(3.17)By considering the graphs of the two sides of equation (3.17) one can easily deduce that such a root exists and is unique for any *β*∈(−∞,+∞), and that *Χ*:(−∞,+∞)→(0,+∞) is a monotonically increasing function. Its analytical form is(3.18)where(3.19)Note that for *β*>27/4, the function *Χ*(*β*) remains real despite *q* becoming complex (*casus irreducibilus*).

Equations (3.15), (3.16),(3.18) and (3.19) define dependence of *c*^{2}/*D*_{u} on the dimensionless threshold *θ* and the dimensionless ratio *p*=*D*_{v}*α*^{2}/*D*_{u}. This dependence is illustrated in figure 10*a*.

An important feature is that positive propagation speed *c* is achieved at all values of the threshold *θ* in the range (0, 8/9), i.e. both above and below *θ*=1/2. Due to the transformation *u*_{1}↔*u*_{2}, *θ*↔1−*θ* mentioned above, this result means that if *θ*∈(1/9, 8/9), we have the possibility of propagating either a wavefront or a waveback, both with oscillatory onsets, at the same value of *v*_{*}. Recall that this is different from typical reaction-diffusion fronts with standard self-diffusion, where at *θ*<1/2 we would have only a wavefront and at *θ*>1/2 only a waveback.

Now let us represent the results in terms of the first integral *v*_{*}. By considering *u*_{1}=*u*_{−}, *u*_{2}=*u*_{+} for wavefront and *u*_{1}=*u*_{+}, *u*_{2}=*u*_{−} for waveback, we have(3.20)(3.21)which are to be substituted for *θ* in equation (3.16). Figure 10*b* shows the resulting dependence of the velocities of the wavefront, *c*_{f}, and the waveback, *c*_{b}, on *v*_{*}, for selected values of parameters.

Although these results are obtained for a piecewise linear caricature, we note that the antegrade and retrograde asymptotics of a front are the same for a generic N-shaped nonlinearity *f*(*u*), and so we may expect a qualitatively similar behaviour for the cubical FitzHugh nonlinearity (1.2). Figure 10*c* illustrates this similarity.

### (c) Slow movement

Evolution in slow time is obtained by the limit *ϵ*→0 of equation (3.2), i.e.(3.22)Let denote the two stable solutions of the equation *f*(*U*)=*V*, where the upper solution applies within the plateau of the wave and applies before and after the wave. Then the slow system reduces to(3.23)This equation depends on *X* only as a parameter and its smoothness on *X* is only due to initial conditions. Such initial conditions are to be obtained by matching with front solutions, if any fronts passed through point *X* before time *T* under consideration, or otherwise from initial conditions of the whole problem, with account being taken of the initial fast transient if appropriate.

### (d) Asymptotic matching

Let us consider a wavefront first. If we characterize the position of a wavefront by its coordinate *x*_{f}(*t*), so that *u*(*x*_{f}(*t*), *t*)≡*u*_{f}, then a steadily propagating wavefront is described by *x*_{f}=ξ_{0}+*c*(*v*_{*})*t* for an arbitrary *ξ*_{0} or, equivalently(3.24)where *v*_{*} is the first integral, and function *c*_{f}(*v*_{*}) is determined from the solution of the steady propagating front problem (3.4).

The first integral *v*_{*}, being an arbitrary constant in the fast space–time scale (*x*, *t*), becomes a function of slow space-time coordinates (*X*, *T*)=(*ϵx*, *ϵt*) at the slow scale. Thus the wavefront motion equation (3.24) at the slow scale is(3.25)where .

To determine the slow function *v*_{*}(*X*_{f}, *T*), we need to match the fast-scale front solution with the slow scale solution.

For given front movement *x*_{f}(*t*), the fast solution is(3.26)where , is the solution of the steady front problem (3.4) and (3.5) for a given value of *v*_{*}, and *v*_{*} has the value corresponding to the location of this piece of front on the long scale, *v*_{*}=*v*_{*}(*X*_{f}, *T*).

By van Dyke's matching rule, the outer limit of the inner solution, i.e (3.26), should be equal to the inner limit of the outer solution, i.e. the solution of equation (3.23). For the *v*-component this gives(3.27)Thus, the matching condition simply means that the first integral of the front coincides with the limit of *V* of the slow solution on the right-hand side of the front. A similar matching condition is found on the other side of the front, *v*_{*}(*X*_{f}, *T*)=*V*(*X*_{f}−0, *T*). So the slow field *V*(*X*, *T*) is continuous over the front line *X*=*X*_{f}(*T*). The fast field *U*(*X*,*T*) is found to obey different conditions, ahead of the front and behind it, and so exhibits a jump between and at *X*=*X*_{f}(*T*).

Matching for wavebacks is considered similarly, and the result is the same except that the fast field *U* switches in the opposite direction.

To summarize, the matching conditions are(3.28)where *X*_{f}(*T*) is defined by equation (3.25), and *X*_{b}(*T*) by a similar equation(3.29)

### (e) Evolution of a pulse

Now we can describe the evolution of a single wave pulse as in the simulations above, at the slow time and space scales, away from the fronts. We start with considering an ageing pulse, as exemplified by figures 1 and 3. It is described by the system of equations (3.23), (3.25), (3.28), (3.29).

*Before the wavefront*: the system is in the resting state, *U*=*V*=0, as these are the initial conditions, and this state is stationary under equation (3.23).

*The wavefront*: propagates according to equation (3.25), where according to matching conditions (3.28), we have *v*_{*}≡0 by matching with the pre-front state. Thus the wavefront propagates with a constants speed, d*X*_{f}/d*T*=*c*_{f}(*v*_{*}(*X*_{f}, *T*))≡*c*_{f}(0), and(3.30)where and time *T* is measured from the moment when *X*_{f}=0.

**The plateau phase** is described by equation (3.23) where the initial conditions are determined by equation (3.28) via matching with the wavefront(3.31)and its solution is therefore(3.32)where(3.33)**The waveback** moves according to equation (3.29) in which *v*_{*}(*X*, *T*) is determined by equation (3.28) via matching with the plateau solution, which gives(3.34)Taking into account equation (3.30), this gives an equation for the width of the wave, *W*(*T*)=*X*_{f}(*T*)−*X*_{b}(*T*), in the form(3.35)readily solved in quadratures,(3.36)The possibility of an explicit answer here depends on the concrete form of the functions *V*_{+}(*T*) and *c*_{b}(*v*_{*}). The exact form of the function *c*_{b}(*v*_{*}) is so complicated that even though *V*_{+}(*T*) can be calculated exactly, the integral (3.36) does not look promising. Thus we do some further approximations valid for the concrete numerical examples discussed above.

First, we notice that the relevant range for the function *c*_{b}(*v*_{*}) is the interval , which is covered while the waveback catches up with the wavefront. From figure 10*b* we see that, in this interval, *c*_{b}(*v*_{*}) is very close to a linear function. We linearized this dependence around the value *v*_{0} at which *c*_{b}(*v*_{0})=*c*_{f}(*v*_{0}), i.e. the middle of the relevant interval. So we put *v*_{*}=*v*_{0}+*σ*, . Then from (3.16) we have . Additional simplification comes from the fact that for the numerical examples considered, *p*=*α*^{2}*D*_{v}/*D*_{u}≈0.0029, i.e. is very small, and consequently, relevant values of *β* are small. Thus the cubic (3.17) can be approximately solved as *Χ*≈1+*β*^{1/3}+*O*(*β*^{2/3}) and, substituting this all into equation (3.15) and linearizing in *σ*, we ultimately obtain(3.37)where and . In particular, . The quality of this linear approximation is shown in figure 10*b*.

The numerically obtained dependencies *c*_{f,b}(*v*_{*}) for the cubical model (1.2) also are well approximated by linear ones in the range *v*_{*}∈(0, 2*v*_{0}), see figure 10*c*.

With linear dependence *c*_{b}(*v*_{*}), quadrature (3.36) for the wave width is easily calculated. As *f*(*u*) is piecewise linear, we have(3.38)This is exact for the piecewise linear model, but a similar linear dependence can be obtained, albeit approximately, for the cubic model, too. Then integral (3.33) gives(3.39)where *k*=(*α*+1)/*α*. Substituting this and (3.37) into (3.36) gives(3.40)where *A*=*k*/*c*_{1} and . Therefore(3.41)where *W*_{∞}=−ln(1−2*kv*_{0})/*C*, *μ*=(1−2*kv*_{0})*C*/*A*. Figure 10*d* illustrates this analytical solution for selected values of the parameters. Comparison with figure 2*b* shows that, although we have replaced a cubic nonlinearity *f*(*u*) with a piecewise linear function, the analytical theory gives an excellent qualitative agreement with the direct numerical simulations of the original problem.

### (f) On the reflection conditions

As the numerics have suggested, the transition from reflective, quasi-soliton mode, to annihilating mode occurs as the wave ages and widens, and this is characterized by slight changes of the waveback but hardly any changes of the wavefront; see figure 3*c*,*d*. Results shown in figure 9 also suggest that reflection of a wave from the boundary is an event happening mainly to the waveback. Although the process of reflection is essentially non-stationary and thus very difficult for analytical treatment, even in the piecewise linear formulation, we try to present here some preliminary considerations towards analytical understanding of this process.

These considerations are based on an analysis of the order parameter, known as the first integral *v*_{*} on the fast scale and as the slow field *V*(*X*, *T*) on the slow scale. Equivalence of the two is guaranteed by the matching conditions (3.28); henceforth we do not distinguish between them and use notation *v*_{*} for both.

The key observation is that for a range of *v*_{*}, which for the piecewise linear variant of the model includes at least the interval *θ*∈(1/9,8/9), the fast system admits both wavefront and waveback solutions, i.e. trigger waves between the same asymptotic states but propagating in opposite directions. Clearly, which way a trigger wave will actually propagate, depends on initial conditions. Let us see what implications this may have for reflection of a waveback from a boundary. As reflection happens on the fast time scale, it means that the front propagating in the retrograde direction develops, as a result of the non-stationary perturbation introduced by the boundary, at approximately the same value of *v*_{*} as the waveback. As the numerics suggest (see figure 8), a non-reflective mode may actually involve generation of the reflected wave which, however, does not survive and decays. Note that the collision generates not only a backward propagating wavefront, but also a backward propagating waveback. Obviously, a *necessary* condition for the survival of the backward wave is that its front propagates faster than its back.

This necessary condition can be expressed in terms of the above analytical results. As we have seen both for the piecewise linear and cubic *f*(*u*), see figure 10*b*,*c*, the speed *c*_{f}(*v*_{*}) decreases and the speed *c*_{b}(*v*_{*}) increases as *v*_{*} increases, and the backward wave can survive only if its front goes faster than its back,(3.42)which implies that *v*_{b}=*v*_{*}(*X*_{b})<*v*_{0}, i.e. the wave is young and thin enough. Notice that in the equation (3.36), the denominator of the integral is a linear function of the waveback velocity *c*_{b}, and ranges from 0 (for *ξ*=*W*_{∞}) to 2*c*_{1}*v*_{0} (for *ξ*=0). Thus, the moment *T*_{*} when *v*_{b}=*v*_{0} and *c*_{b}=*c*_{0} occurs when the denominator of the integral (3.36) or, equivalently, (3.40) is half of its maximal value, so that the waveback has a speed halfway between *c*_{b}(0), its initial value at *W*_{0}=0, and , i.e. the wavefront speed, its ultimate value. This condition leads to(3.43)and(3.44)Figure 10*d* shows these moments for the selected sets of parameters. Comparison with figure 2*b* shows a reasonable qualitative correspondence with the moments when the reflection property of the propagating pulses is lost in direct numerical simulations. Again, a quantitative difference is expected here because of the different nonlinearity. Besides, the criterion (3.42) is only heuristic, obtained through a series of approximations, and is thought of as necessary but not sufficient. An immediate check of the heuristic condition (3.42) is possible by analysing the dynamics of the speed of the waveback, such as those shown in figure 4. Assume that the order parameter of the waveback starts from *v*_{*}=0 and the corresponding minimal value of *c*_{b,min} and, for *ϵ*>0, rises until *v*_{*}=2*v*_{0} and the corresponding value of *c*_{b,max} equals *c*_{f}. Assume further that the speed of the waveback depends on *v*_{*} linearly. Then criterion (3.42) corresponds to the waveback speed exactly equalling the arithmetic mean between *c*_{b,min} and *c*_{b,max}. Analysis of data shown in figure 4 and from numerous other simulations not presented here suggests that whereas this criterion works well for not very small *ϵ*, say *ϵ*=0.01, it gives a consistent overestimation of the transition time for smaller *ϵ*, say *ϵ*=0.001. That is, for very small *ϵ*, the loss of reflection property happens well before the velocity of the back reaches the middle speed (see figure 4*c*). This is an evidence that condition (3.42) is not the only one, as it only provides the possibility of the survival of the backward wave, assuming that the perturbation inducing that backward wave is present as a result of the collision. Some numerical experiments suggest that the wavefront and/or the slope of the plateau play some role in creating such a perturbation, and at very small *ϵ*, the wavefront may be so far ahead of the waveback and the plateau may be so smooth that the perturbation is too small, which is why the loss of the reflection property happens in that case earlier than is predicted by (3.42). This interpretation is consistent with the interpretation of the mechanism of reflection suggested in (Biktashev *et al*. 2004) which was in terms of the excess of the population of prey going backwards through the predators and being depleted by them, before breaking to the predator-free space and initiating the reflected wave. Proper analytical description of this scenario is going to be much more complicated than the simple formulas presented above, as this is an essentially non-stationary process.

Argentina et al. considered nonlinear dissipative wave systems in which colliding waves reflected in some conditions but annihilated in others. This included a continuous chain of damped pendula subjected to a constant torque (Argentina *et al*. 1997), equivalent to a two-component reaction–diffusion system with zero self-diffusion and one of the cross-diffusion coefficients non-zero, and the classical cubic FitzHugh–Nagumo system with self-diffusion of both components (Argentina *et al*. 2000). In both cases, the authors have been able to identify a stationary ‘nucleation’ solution with a single unstable eigenmode. The co-dimension one centre-stable manifold of that solution has been identified with the boundary separating the attraction basins of the annihilation and reflection. This is an intuitively appealing scheme which could ease the analysis of the reflection conditions. However, it does not necessarily work in our case. The piecewise linear caricature (3.1) and (3.7) admits analytical study of a non-trivial steady-state solution. Tedious but straightforward calculations show that such a solution only exists if *a*<*α*/(2(*α*+1)), which is, for example, far from true for the parameters used in figure 10.

## 4. Discussion

We have shown that the solitary waves in FitzHugh–Nagumo system with cross-diffusion terms have the following properties

The waves have peculiar profiles. The wavefronts have oscillatory onset. If the waves are wide enough to tell the waveback from the wavefront, then the wavebacks have oscillatory onset, too.

The waves have long transients, when the wavelength slowly changes approaching its stationary value. We call this ‘ageing’ of the waves.

In appropriate conditions, waves demonstrate quasi-soliton behaviour, i.e. penetrate through each other or reflect from boundaries. The conditions include the ages of the waves. Namely, younger and thinner waves reflect more easily than older and wider waves.

The described properties of cross-diffusion waves are similar to pursuit–evasion waves in the predator–prey system (1.1), and different from the waves in excitable systems with standard self-diffusion terms. Thus we conclude that these unusual waves are not restricted to population dynamics and can be expected to be found in a wide variety of physical systems with nonlinear kinetics and cross-diffusion interaction between fields.

The advantage of the present work is that unlike our previous papers cited above, where the results were based on numerical simulations, here we have been able to present an analytical theory, which allows us not only to confirm the universality of the key phenomena, but also to throw some light on their mechanisms.

In particular, notice that the dependence of the reflection property on the wave age here is opposite to that observed in (Tsyganov & Biktashev 2004); however, its dependence on the wave width is the same. This is because waves in (Tsyganov & Biktashev 2004), due to the initial conditions used there, were shrinking with time while here they have been growing. This is in perfect correspondence with the heuristic criterion (3.42), which requires that the speed of the waveback should not be larger than the speed of a wavefront at the same *v*_{*} would be, and therefore is not satisfied when the wave is rapidly shrinking so that its back catches up with its front.

Thus, the quasi-soliton property, both here and in the system considered in (Tsyganov & Biktashev 2004), appears to be closely related with the existence of wavefront and waveback solutions at the same value of the slow parameter, denoted here by *v*_{*}. Remember that this is very different from the FitzHugh–Nagumo system with classical self-diffusion coefficients, where only one trigger wave solution, either a waveback or wavefront, exists for every value of the slow variable. And this difference should therefore be attributed to the cross-diffusion terms here, or the mutual taxis terms in (Tsyganov & Biktashev 2004).

Another phenomenological manifestation of the cross-diffusion and mutual taxis terms is the oscillatory onset of the trigger waves, i.e. wavefronts and wavebacks. This oscillatory character is similar to the solitons in the nonlinear Schrödinger equation (NLS), and the diffusion terms in our system, at *D*_{u}>0 and *D*_{v}>0, are equivalent to the diffusion terms of the NLS, so that the diffusion matrix is not positive definite, as its eigenvalues lie on the imaginary axis. This is the second property, in addition to reflection, uniting our present system with NLS and distinguishing it from the classical FitzHugh–Nagumo.

On the other hand, certain properties of the cross-diffusion excitation waves are similar to those in self-diffusion excitable systems. Indeed, despite the relatively long transient, the cross-diffusion waves approach a stable stationary shape, which is independent on the details of initial conditions (e.g. in terms of our theory, on the initial width of the wave *w*_{0}). Also, although the trigger wave solution is not unique for every value of the slow variable, we only have two rather than a continuous spectrum of such solutions.

Thus, we have in (1.2) a system which combines some features of conservative systems, such as a diffusion matrix with purely imaginary eigenvalues as in the nonlinear Schrödinger equation, and some features of dissipative systems, such as the local nonlinearity terms as in the classical FitzHugh–Nagumo system. Correspondingly, the solutions have some dissipative properties, such as independence of the established wave shape and speed on initial conditions, and some conservative properties, such as reflection of waves.

## Acknowledgments

This study was supported by EPSRC grants GR/S08664/01 and GR/S75314/01 (UK) and by RFBR grant 03-01-00673 (Russia).

## Footnotes

- Received June 2, 2004.
- Accepted June 10, 2005.

- © 2005 The Royal Society