Stable bound states of pulses in an excitable medium

Physica D 135 (2000) 154–174

Stable bound states of pulses in an excitable medium Michal Or-Guil∗,a , Ioannis G. Kevrekidis b , Markus Bär a a

Max Planck Institute for the Physics of Complex Systems, Nöthnitzer Strasse 38, 01187 Dresden, Germany b Department of Chemical Engineering, Princeton University, Princeton, NJ 08544 USA Received 30 November 1998; received in revised form 7 May 1999; accepted 7 May 1999 Communicated by J.P. Keener

Abstract The interaction of one-dimensional pulses is studied in the excitable regime of a two variable reaction–diffusion model. The model is capable of exhibiting long range attraction of pulses and formation of stable bound pulse states. The important features of pulse interactions can be captured by a combination of various analytical and numerical methods. A kinematic ansatz treating pulses as particle-like interacting structures is described. Their interaction is determined using the dispersion relation for pulse trains, which gives the dependence of the speed c(d) of the wavetrain on its wavelength d. Anomalous dispersion for large d, i.e. a negative slope of c(d), corresponds to long range pulse attraction. Stable bound pairs are possible if the medium exhibits long range attraction and there is at least one maximum of the dispersion curve. We compare predictions of the kinematic theory with numerical simulations and stability analysis. If the slope of the dispersion curve changes sign, branches of non-equidistant pulse train solutions bifurcate and may lead to bound pulse states. The transition from normal long range dispersion, typical in excitable media, to the anomalous dispersion studied here can be understood through a multiscale perturbation theory for pulse interactions. We derive the relevant equations, which yield an analytic expression for non-monotonic dispersion curves with a finite number of extrema. ©2000 Elsevier Science B.V. All rights reserved. PACS: 82.20.Mj; 47.54.+r; 02.30.−f Keywords: Stable bound states; Pulses; Excitable medium

1. Introduction One-dimensional excitable media exhibit nonlinear waves such as wavetrains and solitary pulses. Examples include concentration waves in chemical reactions in solution and on surfaces, as well as electrical impulse propagation in neurons and cardiac tissue. These processes are often modelled by two-variable reaction–diffusion equations [1–3]. These models describe the coupled dynamics of fast autocatalytic species (activator) and a slow variable introducing a negative feedback on the fast species (inhibitor) [4,5]. A pulse is a localized structure moving with ∗ Corresponding author E-mail addresses: [email protected] (M. Or-Guil), [email protected] (I. G. Kevrekidis), [email protected] (M. Bär)

0167-2789/00/$ – see front matter ©2000 Elsevier Science B.V. All rights reserved. PII: S 0 1 6 7 - 2 7 8 9 ( 9 9 ) 0 0 1 3 6 - 0

M. Or-Guil et al. / Physica D 135 (2000) 154–174

155

constant velocity c. It can be created by a suitable localized perturbation (“excitation”) of a linearly stable, homogeneous rest state. Such a perturbation may evolve into a solitary pulse traveling with constant speed and shape. If the perturbation is applied repeatedly or for a longer time, a train of pulses may result. For multipulse solutions with large interpulse distances, the dynamics can be described in terms of weak interactions between the pulses. Perturbation theory states that these long range interactions lead to corrections to the pulse velocities that depend exponentially on the distance between neighboring pulses [6,7]. The long range interaction can be computed based on the spatial profile of the single pulse; in particular, its strength is determined by the spatial exponential decay of the pulse profile into the stable rest state of the medium. Interaction laws were likewise derived for stationary spots in two-dimensional domains, explaining “molecule” structures [8]. Pulses in excitable media usually consist of a steep front, caused by a fast excitation of the activator variable, and a smooth back, reflecting the slow recovery of the inhibitor variable. For such a strongly asymmetric shape, the velocity of a pulse depends only on the distance d to its predecessor. This can be exploited for the formulation of a reduced dynamics description of the pulse position. The so-called kinematic theory states that dynamics of the pulse positions are governed by a function c(d) [9,10]. This function corresponds to the dispersion curve, which describes the relation between wavelength and velocity of a periodic wavetrain. Hence, a numerically computed dispersion curve for these pulse trains gives insight into pulse interactions and its implications. Depending on the shape of the single pulse tail, three cases of long range interactions can be distinguished: attractive, repulsive and oscillatory interactions. They correspond to negative, positive and oscillating slopes of the dispersion curve in the limit of large distances. We will refer to a positive (negative) slope of the dispersion curve as normal (anomalous) dispersion. When attraction and repulsion between two adjacent pulses compensate each other, a bound pulse pair is formed. This requires the two pulses to have the same velocity. Such states are not possible when the dispersion curve is monotonic, but they may exist when it is oscillatory [6–12]. Here, we show that stable bound pairs of pulses can also be found when the curve c(d) exhibits only a single maximum. Part of the motivation for our investigations stems from recent experimental findings in chemical and biological reaction–diffusion systems. Anomalous dispersion for large interpulse distances has been observed in experiments on the aggregation of slime mold cells exposed to caffeine [13], in surface waves in Marangoni convection [14] and in a chemical reaction on a catalytic surface [15]. Furthermore, stable moving bound states were predicted for a three component RD system in a two-dimensional domain at the onset to propagation [16]. This paper investigates the interaction of pulses as well as the possibility of stable bound states and non-equidistant pulse trains and their reduced description in an excitable reaction–diffusion model. This is achieved by using a combination of perturbation and kinematic theory with bifurcation results and numerical stability analysis. In our opinion, all of these techniques are essential in the understanding of pulse interactions and may also be used in related problems. The investigation is performed with a simple activator–inhibitor model, which we analyze in detail. The methods presented here also apply to more complicated realistic models of chemical reactions and biological systems. In Section 2 we describe the model equations as well as the numerical methods used in this study. Section 3 outlines the ideas of reduced dynamics and their application to our case. In Section 4 the results of reduced dynamics are compared to numerical simulations and stability analysis. In Section 5 an analytical expression for the dispersion curve is motivated by a two-scale perturbation theory. The complete derivation of this expression is contained in Appendix A. Section 6 relates the above results to homoclinic bifurcations in ordinary differential equations. Summarizing the results of this paper, anomalous dispersion leads to destabilization of equidistant pulse arrangements and results either in formation of non-equidistant arrangements or in the eventual annihilation of one or more of the pulses initially present. The emerging non-equidistant pulse trains still propagate with constant overall shape and speed. For large system length L, these non-equidistant pulses may approach bound states. The simplest bound state is a bound pair that can emerge from an initial condition with two pulses.

156

M. Or-Guil et al. / Physica D 135 (2000) 154–174

Maxima and minima of the dispersion curve separate regions of normal and anomalous dispersion. If the system length L is taken as a bifurcation parameter and the considered solution branch corresponds to N equidistant pulses on a ring, it is found that N − 1 branches of non-equidistant pulse trains bifurcate near extrema in c of the branch of equidistant pulse trains, and approach various combinations of bound states of pulses (pairs, triples, etc.) as L → ∞. Interaction eigenvectors and corresponding eigenvalues are identified and give information about the type and strength of the pulse interactions. An analytical expression gives a simple explanation for the transition from a monotonic dispersion curve to a curve with a single maximum and the anomalous long range dispersion observed for the model studied here.

2. Model We investigate equations of the FitzHugh–Nagumo type with a fast activator variable u and a slow inhibitor v:   1 b+v ∂t v = f (u) − v, (1) + ∂x2 u, ∂t u = − u(u − 1) u −  a  0 ≤ u < 1/3  0, 2 f (u) = 1 − 6.75u(u − 1) , 1/3 ≤ u ≤ 1,  1, 1
    b + v0 1 u0 (u0 − 1) + u0 − (2u0 − 1) ,  a

B(z) =

u0 (u0 − 1), a

D(z) = ∂u f (u0 ).

The linear stability problem for traveling waves amounts thus to the determination of the spectrum of the Jacobian M in Eq. (3). In general, solutions u0 (z), v0 (z) and the eigenfunctions of the Jacobian M are not available in closed form, so the problem has to be approximated numerically. We compute wave-trains as stationary solutions in the traveling frame by using a pseudospectral discretization of Eqs. (2) with periodic boundary conditions and Newton–Raphson

M. Or-Guil et al. / Physica D 135 (2000) 154–174

157

iteration. The velocity c, which is not known a priori, is formally viewed in the computation as an additional variable. A so-called pinning condition singles out one of the infinitely many solutions existing due to translational symmetry. With that, the velocity is likewise calculated numerically. The eigenfunctions and the spectrum of eigenvalues of the Jacobian are calculated in Fourier space resulting from a 200 Fourier mode discretization of the stationary solution. The zero eigenvalue, which always exists due to the translational symmetry of the problem, can be used as a numerical accuracy check and has been obtained with a precision of 10−4 to 10−3 . The time evolution starting close to unstable solutions was computed using an explicit finite difference scheme to solve Eqs. (1). Because traveling wave solutions fulfill the conditions ∂t u = ∂t v = 0 in the comoving frame, they can also be found by studying the so-called traveling wave ordinary differential equations (ODEs) du = w, dz

dw 1 = −cw + u(u − 1)(u − uth ), dz 

dv = (v − f (u))/c, dz

(4)

with uth = (b + v)/a. A pulse corresponds to a homoclinic orbit in these equations, while a wavetrain emerges as a limit cycle therein. In the parameter range considered here, there are three relevant fixed points of Eqs. (4): A = (0, 0, 0), B = (b/a, 0, 0), and a focus. They correspond to the homogeneous states of Eqs. (1): Fixed point A to the stable rest state and B to the saddle. We have studied wavetrain solutions in Eqs. (1) by numerical continuation of limit cycles, detection of period doublings and automatic branch switching with the continuation software AUTO 94 [22] in the above ODEs (4). In the following, the parameters of the model are fixed at a = 0.84 and b = 0.07, while the time scales ratio  > 0 is varied. We shall base our analysis on pulses which decay into the stable rest state.

3. Reduced dynamics: perturbation theory and kinematic theory The principal idea of reduced dynamics is to consider the pulses as localized objects which only weakly interact with other pulses. Then, the dynamical changes in the system can be described through changes in the positions of the interacting pulses, while the shape of the pulses stays constant. Consider a train of pulses with large interpulse distances. For this case, it was shown by perturbation theory [7] that the scaling behavior of the spatial decay of a single pulse into the rest state determines also the scaling behavior of the pulse interactions with the pulse distance. The spatial decay behavior can be computed by linearizing the traveling wave ODEs (Eqs. (4) for our model) around the fixed point corresponding to the rest state. For right-going pulses, which we consider here without loss of generality, the eigenvalues of the linearization with positive real part give the rate of the exponential decay of the back of the pulse far from the pulse position, while the eigenvalues with negative real part describe the corresponding behavior for the front. The most relevant decay rates are given by the positive resp. negative eigenvalue with the smallest absolute value. Note that for a general reaction–diffusion system with P diffusing and Q non-diffusing components, there are 2P + Q traveling wave ODEs and thus 2P + Q eigenvalues of the linearization around the fixed point. Our example corresponds to P = Q = 1 and gives three eigenvalues. As they are all real in the case studied here (see Section 5), we will in the following restrict ourselves to this case, excluding thus oscillatory interactions throughout the paper. In this case, the eigenvalues λ of the traveling wave ODEs give directly the spatial decay rates, and notation gets simpler. Let λa > 0 and λb < 0 be the eigenvalues, giving the relevant scales of the back resp. the front of a single pulse. It was shown by perturbation theory [7] that, if λa ≈ |λb | holds, the pulse dynamics for approximately equidistant pulse trains can be described by a set of reduced equations for the pulse positions pi : p˙ i = A1 e−λa (pi−1 −pi ) + B1 e−|λb |(pi −pi+1 ) + c∞ ,

(5)

158

M. Or-Guil et al. / Physica D 135 (2000) 154–174

where i = 1, . . . , N and di := pi−1 − pi > 0 for i = 2, . . . , N and d1 := L − p1 + pN are the distances between neighboring pulses. In the specific formulation used here, we assume from now on that the pulses travel in a system with periodic boundary conditions and length L. The rightmost pulse is labeled by i = 1. The first (second) term in Eq. (5) describes the interaction of the ith pulse with neigboring pulses traveling in front (behind), while c∞ is the velocity of a solitary, unperturbed pulse in the medium. The interaction coefficients A1 and B1 are obtained from the shape of the single pulse solution. The magnitudes of the quantities λa and |λb | are a measure of the interaction scales, while the signs of the coefficients A1 and B1 decide the type of interaction (attractive or repulsive). A special case is given when the single localized structure is stationary and symmetric, yielding λa = −λb and A1 = B1 . On the other hand, traveling pulses in excitable systems are strongly asymmetric, having a steep front and a spatially slowly decaying back. In general, the smaller the time scales ratio, the more asymmetric the pulse shape. The interaction of such pulses, where λa  |λb |, can be described by [10] p˙ i = c∞ + A1 e−λa (pi−1 −pi ) .

(6)

In this case, the dynamics of the pulse position depends only on the distance to the preceding pulse. In the model studied here, the interaction scales for stable pulses are given by the three eigenvalues λi of the linearization of Eqs. (4) around the fixed point (0, 0, 0). In the range of interest, λ1,2 > 0 and λ3 < 0, indicating that the steep front is described by λ3 , while λ1 and λ2 characterize the smooth tail of the pulse. For small values of the time scales ratio , the eigenvalues scale as λ1 ∝  1/2 and λ2,3 ∝  −1/2 . Thus, |λ1 |  |λ2,3 | and we identify λa = λ1 (see Section 5). Since λ1  −λ3 for all parameter values of interest, the pulses under consideration are always strongly asymmetric and their dynamics can be described by Eq. (6). For larger  however, λ1 and λ2 may become comparable and then a refined treatment including competing interaction scales is needed (see Section 5). As the rhs of Eq. (6) gives an expression for the dependence of the pulse velocity on interpulse distances, it gives consequently also an approximation for the dispersion relation for equidistant pulse trains c(d) with large interpulse distances d. Thus, we can replace the exponential form of the interaction by the dispersion relation c(d) [9,10]. The resulting equations of motion are p˙ i = c(di ),

(7)

where again i = 1, . . . , N and di := pi−1 − pi for i = 2, . . . , N and d1 := L − p1 + pN are the distances between the pulses. Thus, in the framework of the kinematic theory, the dispersion relation c(d) determines pulse–pulse interactions and at the same time represents the branch of equidistant pulse solutions. In the limit of small , the dispersion relation c(d) for Eqs. (1) can be derived analytically [18]. For large distances d, this expression can be cast in the form of Eqs. (6) and reads c(d) ≈ c∞ −

2(a − 2b) −d/c∞ , e √ a 2

(8)

√ where c∞ = 1/2a(1 − 2b/a) is the solitary pulse velocity for  → 0. In this limit case, we get a monotonic, normal dispersion curve, where the dominating interaction scale is given by λa = 1/c∞ . The prefactor of the exponential term on the rhs of Eq. (8) can be identified with the quantity A1 in Eq. (6). Similar expressions have been derived for the Oregonator model of the Belousov–Zhabotinsky reaction and compare well to experimental measurements [23]. In general, c(d) is not known analytically and it has to be computed numerically. Consider now N pulses traveling in a ring of length L. Wave solutions traveling with constant velocity c0 are obtained from Eqs. (7) when the interpulse distances di are such that (i) all pulses move with equal speed, c(di ) = c0 P for all i, and (ii) the sum of the distances yields the ring length: L = N i=1 di . For finite ring length, equidistant (ED) pulse train solutions satisfy the condition d = L/N where N is an integer. Non-equidistant (NED) pulses can only exist, i.e., requirements (i) and (ii) can only be fulfilled, if different interpulse distances correspond to the same

M. Or-Guil et al. / Physica D 135 (2000) 154–174

159

velocities on the dispersion curve. This can only occur if the dispersion curve is non-monotonic; it should possess at least one extremum in the validity range of reduced dynamics. By transforming the reduced dynamics, Eqs. (7), into the appropriate comoving frame, i.e., by considering the equations p˙ i = c(di ) − c0 ,

(9)

traveling solutions are transformed into stationary ones. The stability of these solutions can then be determined by linearizing this equation. For perturbations δpi of the original pulse positions this results in: δ p˙ i = c0 (di ) · (δpi−1 − δpi ),

δ p˙ 1 = c0 (d1 ) · (L + δpN − δp1 ),

(10)

where 0 denotes differentiation with respect to the argument. The eigenvalues ω of this eigenvalue problem can be obtained by solving the characteristic equation of the Jacobian of Eqs. (10) N N Y Y (ω + c0 (di )) − c0 (di ) = 0. i=1

(11)

i=1

The eigenvalues for an ED train of N pulses with spatial period d are thus given by ωk = c0 (d)(ei2πk/N − 1);

k = 0, . . . N − 1.

(12)

Thus, the stability depends only on the slope of the dispersion relation c(d). If c0 (dk ) = 0 for a single pulse, the eigenvalues are obtained as ω

N Y

(ω + c0 (di )) = 0.

(13)

i=1,i6=k

Eq. (13) also gives the eigenvalues of a train of N pulses in the unbounded domain, where the dynamics of the first pulse p˙ 1 = c∞ is not influenced by interactions with other pulses. In the following we will apply kinematic theory to our specific model. Fig. 1(a) shows a dispersion curve c(d) for  = 0.1. It exhibits anomalous behavior for large distances d and it displays a maximum at dmax . Notice that there exists a finite distance db with the velocity c(db ) = c∞ of the free pulse. Notice furthermore that the destabilization of a pulse train at d ≈ 10 (see Fig. 1(a)) is due to an eigenmode which does not shift the pulse positions, but changes the pulse shape (result of numerical stability analysis, not shown). Thus, solutions for d < 10 can no longer be captured by reduced dynamics. We will therefore renounce to show this range in the remaining figures. To better illustrate the predictions made by reduced dynamics, we consider at first the case of two pulses in the ring. In this case, the above requirements for stationary solutions in the reduced dynamics formulation predict a branch of NED pulse pairs that emerges from the maximum of the ED branch c(2d) and traces the mean value L(c) = d1 (c) + d2 (c) between the two intersections of the ED branch with the line c = const. The resulting c(L)-curve for the NED pulse pairs is depicted in Fig. 1(b) along with the ED branch. The stability of a stationary pulse pair in the reduced dynamics formulation is determined by the two equations of the eigenvalue problem equations (10). One of the eigenvectors, eE1 = (1, 1)T , corresponds to a shift of both pulses in the same direction. It is a Goldstone mode, corresponding to the translational invariance of the problem, and has thus the eigenvalue ω1 = 0. The second eigenvector, eE2 = (−c0 (d1 )/c0 (d2 ), 1)T ,

(14)

160

M. Or-Guil et al. / Physica D 135 (2000) 154–174

Fig. 1. Dispersion curve c(d) for  = 0.1 (a); solution curves for a ring with two (b) resp. three (c) pulses. Branches of non-equidistant pulse trains bifurcate near the maxima of the equidistant branch c(L). They approach bound two (b) resp. three (c) pulse solutions for large ring lengths. Thin lines show predictions of kinematic theory, while thick lines are numerically computed stationary solutions of Eqs. (2) and/or limit cycles of Eqs. (4). Solid (dashed) lines denote stable (unstable) solutions. The circle in (c) depicts a saddle-node bifurcation which leads to a change of stability of the leftmost non-equidistant three-pulse branch 1.

corresponds to a relative shift of the pulse positions. It shall thus be called interaction eigenvector. The respective eigenvalue, ω2 = −(c0 (d1 ) + c0 (d2 )),

(15)

determines the stability of the pulse pair solutions. In the case of ED pulse pairs, the interaction eigenvector

M. Or-Guil et al. / Physica D 135 (2000) 154–174

161

eE2 = (−1, 1) shifts the two pulses in opposite directions. The eigenvalue is ω2 = −2c0 (d), representing a special case (N = 2) of Eq. (12). This implies that ED pulse pairs are stable (unstable) when the slope of the dispersion curve is positive (negative). As for the NED pulse pairs, it can be drawn from Fig. 1(b) that, in our case, the sum of the slopes corresponding to the NED pulse pairs is always positive. Thus, Eq. (15) implies that these solutions should be stable. For very long systems, one interpulse distance approaches db (which is the distance where c(db ) = c∞ holds, see Fig 1(b)), while the other distance becomes roughly L − db . Thus the NED solution approaches a stable bound pair of pulses with velocity c ≈ c∞ . As the slope of the dispersion curve approaches zero for very large distances d, the interaction eigenvalue approaches ω2 = −c0 (db ). The bound pulse pair is thus only stable if c0 (db ) > 0. Consider now the general case of N pulses on a ring. For a non-monotonic dispersion curve with a single maximum, there is a range of velocities c consistent with two different interpulse distances c = c(d 0 ) = c(d 00 ) with d 0 < d 00 . Here, N-pulse solutions with up to 2N permutations of pulse distances will be found for a given velocity c. It is straightforward to see that all combinations with certain numbers of short (d0 ) and long distances (d00 ) will have the same speed c and, from Eq. (11), that they will all have equal stability properties in the reduced description. Thus, we are left with only N + 1 distinguishable combinations. Since two of them (the “all short” and the “all long”) already belong to the ED-branch, N − 1 new separate NED branches are thus possible in the (c, L)-plane for non-monotonic dispersion. They emerge from the ED solution at the bifurcation point, i.e., the point of maximum speed, as L is increased (see comment at the end of Section 5). The emergence of N − 1 NED branches is consistent with the fact that N − 1 different eigenvalues cross the imaginary axis at the bifurcation point in the kinematic description of an ED solution of N pulses, as stated by Eq. (12), yielding also N −1 different destabilizing eigendirections. Kinematic theory predicts that (a) all NED branches emerge at the velocity maximum of the ED branch and (b) approach bound solutions with velocity c ≈ c∞ , if there exist a stable ED solution with wavelength db and c(db ) = c∞ . In order to regard NED branches more closely, we consider pulse triples as an example. For this case, two distinct new NED branches exist. They correspond to solutions with the sequences of interpulse distances (d 0 , d 0 , d 00 , branch 1) and (d 0 , d 00 , d 00 , branch 2). These two branches are shown in Fig. 1(c). Kinematic theory predicts that branch 2 is always unstable. In contrast, branch 1 in Fig. 1(c) bifurcates subcritically towards small L, turns around via a saddle-node bifurcation and approaches for large L a stable bound triple of pulses. The latter observations are also true for the general case of N pulses: any solution with large enough L, N − 1 small distances d 0 , and one large distance d 00 , is stable, while all other combinations are unstable. This can be shown by estimating the eigenvalues given by the characteristic polynomial, Eq. (11). To accomplish that, we consider the case of N pulses distributed such that there exist p small and q large distances, and p + q = N . By introducing c0 (d 0 ) =: a and c0 (d 00 ) =:  we can write Eq. (11) as: (ω + a)p (ω + )q − a p  q = 0.

(16)

If  = 0, the solution is stable as we obtain a p-fold root ω = a and a q-fold degenerate eigenvalue ω = 0. Supposing d 00  d 0 ,  can be considered to be small and a to be of order one. Thus we can assume that the eigenvalues ω are of order . With that, and neglecting terms of order higher than  q , Eq. (16) becomes: a p ((ω + )q −  q ) = 0.

(17)

The eigenvalues are thus given by ω,k = (ei2πk/q − 1),

k = 0, . . . q − 1.

(18)

Since  is negative, all these eigenvalues have a non-negative real part. Hence, in the case q > 1, i.e., when more than one of the distances is large, the solution is unstable against interaction modes.

162

M. Or-Guil et al. / Physica D 135 (2000) 154–174

Altogether, the above reduced dynamics description provides a simple method to compute the interaction behavior of pulses either from explicit perturbation theory or from the numerically obtained dispersion curve. We remind the reader of its basic assumption: only one relevant degree of freedom (translation) is assigned to each pulse. This description works if the zero eigenvalue corresponding to the translational invariance of the problem is well separated from the (negative) eigenvalues of the continuous spectrum as well as from other eigenvalues of the discrete spectrum that may appear near bifurcations of the pulse solution. Explicit computation of the complete spectrum by solving Eqs. (3) will be presented in the next section. It enables us to check the validity of the perturbation and kinematic approaches for our model.

4. Comparison with numerical results In order to verify the predictions of the reduced dynamics, solutions with constant shape and constant velocity c were computed as stationary solutions of Eqs. (2) in a periodic domain with length L (see Section 2). This periodic ring length L was used as a bifurcation parameter. For N = 2, the NED branches can alternatively be computed by tracing branches born at period doublings of the periodic solution of Eq. (4) corresponding to ED wavetrains. Numerically computed branches of ED solutions for N = 2, 3 have already been shown in Figs. 1(b) and (c). In addition, the NED branches that approach stable bound pulse states for large L are also depicted in Figs. 1(b) and (c) for N = 2, 3. The bifurcations appear slightly right of the maximum of the ED branches in moderate contrast to the prediction of kinematic theory. To determine the stability of the pulse trains and obtain a systematic comparison with reduced dynamics, eigenvalues and eigenvectors of the Jacobian matrix M, Eqs. (3), have been computed using the numerical method described in Section 2. The spectrum shows always a zero eigenvalue, due to translational invariance. Its eigenvector, the Goldstone mode, is given by the spatial derivative of the wave profile and corresponds to an infinitesimal spatial shift of the wave. The eigenvalue spectrum of a single pulse in an infinite system contains a continuous part, which can be identified with the spectrum of the stable rest state, as well as a discrete part related to eigenfunctions localized near the pulse position. As the single pulse is stable, all these eigenvalues have a non-positive real part. If we consider multiple pulses, additional discrete eigenvalues appear in the spectrum. If the interpulse distances are large enough, the new localized eigenfunctions are roughly linear combinations of the discrete eigenfunctions of the single pulse. Eigenvalues are then slightly shifted and may therefore gain positive real parts, destabilizing the pulse arrangement. Eigenmodes consisting of combinations of the Goldstone mode of the single pulses correspond essentially to relative shifts of the pulses. These modes are the most probable candidates for providing an instability in the multipulse solution. They correspond to the eigenvectors captured by the reduced dynamics. A typical numerically obtained eigenvalue spectrum of a periodic pulse train is shown for N = 16 in the complex plane in Fig. 2. For sufficiently large interpulse distances, a clear separation between the major part of the spectrum and the N interaction eigenvalues is seen (see the inset in Fig. 2). The eigenfunctions of these interaction eigenvalues correspond to the modes used in the reduced dynamics; we shall compare them with the reduced dynamics eigenvectors. The computed interaction eigenmodes for the ED and NED pulse solutions depicted in Fig. 3(a) are shown in Fig. 3(b). While the NED solution is stable here, since all eigenvalues are smaller or equal to zero up to numerical accuracy, the ED solution is not: one eigenvalue is positive. For the ED pulse pair, the interaction eigenvector is an antisymmetric combination of the Goldstone modes of the single pulses, i.e., the pulses are shifted with the same strength in opposite directions. This corresponds to the eigenvector eE2 = (−1, 1)T of the reduced dynamics (Eq. (14)). In the NED pulse case, the corresponding interaction eigenvalue becomes negative and the pulse pair is stabilized in its NED position. The leading pulse is barely influenced by its predecessor because of the considerably

M. Or-Guil et al. / Physica D 135 (2000) 154–174

163

Fig. 2. Leading part of the eigenvalue spectrum of a train of 16 equidistant pulses on a ring computed through numerical discretization. The parameters are a = 0.84, b = 0.07 and  = 0.1, and the ring length is 290 dimensionless units. Most of the approximation to the continuous spectrum lies “far” left of the imaginary axis. Near zero we find sixteen eigenvalues, shown as a blowup in the inset, lying almost on a circle, according to predictions of kinematic theory (see Eq. (13)).

Fig. 3. Unstable equidistant pulse pair solution (dashed line) and stable non-equidistant pulse pair solution (solid line) in a ring of length L = 40 with  = 0.1 (a) and corresponding interaction eigenvectors (b). Perturbation of the pulse solutions along the interaction eigenmodes result in a shift of the pulse positions in opposite directions in the equidistant case, while in the non-equidistant case only one pulse is shifted. The non-equidistant pulse pair corresponds closely to a bound pulse solution.

larger distance (if we think of the domain as periodic). This can also be seen from the shape of the interaction eigenmode (full line in Fig. 3(b)): its amplitude near the rightmost pulse is almost zero. Consequently, the interaction eigenvector indeed corresponds to the prediction of the reduced dynamics eE2 ≈ (0, 1)T (see Eq. (14)) for large ring lengths L. In this limit, the velocity of the NED pulses approaches c∞ (Fig. 1). Thus, the NED solution evidently approaches, for large system lengths, a bound two pulse solution. Consider again two pulses in the ring and let us now vary the ring length L as a bifurcation parameter. The change of L implies a change of the interpulse distances. The stability of the pulse branches, computed from Eq. (2) for

164

M. Or-Guil et al. / Physica D 135 (2000) 154–174

Fig. 4. The two most significant eigenvalues of the linearization around pulse pair solutions (circles) compared with their reduced dynamics prediction, Eq. (10), as a function of the system length L for equidistant (a) and non-equidistant (b) pulse solutions.

different ring lengths, is visualized in Fig. 1(b) by solid (stable) resp. dashed (unstable) lines. The destabilization at the maximum of the ED branch is due to the interaction of the two pulses, while the instability for smaller lengths corresponds to a destabilization of the single pulse itself. The two largest eigenvalues of the Jacobian M, Eq. (3), are shown in Fig. 4 for different system lengths and are compared to the predictions of reduced dynamics (Eq. (15)). The agreement of the interaction eigenvalue with the prediction of kinematic theory in the case of ED pulse pairs (Fig. 4(a)) is good to the right of the maximum of the dispersion curve, getting progressively worse for smaller pulse distances d. Consequently, agreement with reduced dynamic predictions is less accurate in the case of NED pulses (Fig. 4(b)), as these solutions comprise one interpulse distance to the left of the dispersion maximum (see Fig. 1). The destabilization of the ED pulse train by interaction is also evidenced by direct numerical integration of Eqs. (1) using a two-pulse initial condition (Fig. 5(a)). The pulses attract and form a stable bound pair. As an example of multipulse solution branches, the ED and the two three-pulse NED branches have been depicted in Fig. 1(c). Fig. 5(b) shows the formation of a stable NED pulse triple out of an unstable ED initial condition through shift of the pulse positions. Rinzel et al. [9] show, in a case with similar dispersion curve, that unstable ED pulse trains with many pulses tend to cluster to arbitrary “bunches” of several pulses (what we term NED pulse trains). There are cases where pulse distances are too small for predictions based on reduced dynamics to be accurate. If the value of the time scale ratio  is increased, the maximum of the dispersion curve c(d) is shifted towards smaller distances d and the velocity c∞ decreases. Fig. 6 (a) shows the dispersion curve for  = 0.101. There is no stable ED solution left of the dispersion maximum which has the velocity c∞ , and a possible NED pulse pair branch is expected to disappear at the velocity cmin . Nevertheless, a branch of NED pulses still exists at large ring lengths and apparently approaches a stable bound pair. Fig. 6(b) shows a different dispersion curve for  = 0.103. Here, there is also no ED wave solution which has the velocity c∞ , and the NED pulse pair might be expected to disappear almost immediately at cmin , but it persists for much smaller velocities. In this case the branch of NED pulse pairs transforms into a branch of double-humped pulses which possess a different background state (corresponding to the unstable fixpoint B). We should point out

M. Or-Guil et al. / Physica D 135 (2000) 154–174

165

Fig. 5. Space-time plots representing results of numerical integration of Eqs. (1) with unstable equidistant pulse solutions as initial conditions. Black denotes regions of high values of the variable u. (a) Pulse pair  = 0.1 and L = 31.6; (b) pulse triple,  = 0.1 and L = 47.4; (c) pulse pair,  = 0.103 and L = 40.

Fig. 6. (a) Dispersion curve resp. branch of equidistant pulse pairs c(L) = c(2d) for  = 0.101 (thick line) and branch of non-equidistant pulse pairs (thinner solid line). (b) Dispersion curve resp. branch of equidistant pulse pairs c(L) = c(2d) for  = 0.103 (thick line), branch of non-equidistant pulse pairs (thinner solid line) and branch of double-humped pulses with rest state (b/a, 0) (thinner dashed line). Stability (instability) is represented by solid (dashed) lines.

that no single-humped pulses to B exist in this parameter regime (see Section 4). For increasing system lengths, our NED branch approaches this solution (see Fig. 5(b)). Integration of the model (1) for this case with L = 30, where no NED pair exists, shows that the unstable ED pulse pair, as initial condition, evolves through movement of the pulses relative to each other, to an irregular behavior with backfiring and eventually to annihilation of all pulses (see Fig. 5(c)). Clearly, this phenomenology lies beyond the scope of kinematic theory.

166

M. Or-Guil et al. / Physica D 135 (2000) 154–174

5. Extended perturbation theory: the two-scale approach to pulse interactions In the previous sections it has been shown that the reduced dynamics approach of Eqs. (7) can capture the main features of pulse dynamics which are introduced by the presence of a maximum in the dispersion curve. However, this maximum cannot be explained by the assumption of a single spatial scale or equivalently a monotonic approximation of the dispersion curve as in Eq. (6). In contrast, we conjecture that the dispersion curve and thus the pulse dynamics obey a two-scale behavior of the type: c(d) = A1 e−λ1 d + A2 e−λ2 d + c∞ ,

(19)

where A1,2 are suitable constants and λ1,2 are positive spatial decay rates of the single pulse. In other words, we consider two spatial decay rates of the tail of the pulse instead of only one. In principle, the pulse tail is going to influence the dynamics of all subsequent pulses, but, as the tail decays exponentially, this influence is likewise going to decay exponentially with the distance to the regarded pulse. Since we want to consider two decay scales, we have to make sure that the influence of the faster decaying scale is not weaker than influences of the slower decaying scale on the next nearest neighbor. Otherwise, we would have to consider also these interactions and would no longer obtain the above described reduced dynamics. Indeed, if the condition 1/2 < λ2 /λ1 < 2 is fulfilled, a two-scale perturbation approach yields reduced dynamics and a two-scale dispersion curve as in Eq. (19). A detailed derivation of this condition, of reduced dynamics and of Eq. (19) is carried out in Appendix A. The spatial decay of the single pulse with speed c∞ into the rest state (u, v) = (0, 0) is governed by the modulus of the real part of the eigenvalues of the linearization of the ODE system (4) around the fixed point (u, w, v) = (0, 0, 0). They are thus determined by the eigenvalue problem      u u 0 1 0  b/a −c∞ (20) 0 w  = λw  0 0 1/c∞ v v which yields the eigenvalues 1 , λ1 = c∞

λ2,3

c∞ ± =− 2

s

2 b c∞ + . 4 a

(21)

Note that the eigenvalues are always real for a, b,  > 0, and correspond to the decay rates already mentioned in Section 3. Fig. 7 shows these eigenvalues as a function of the parameter , computed from the respective single pulse velocities. The two positive eigenvalues λ1 and λ2 , which give the decay rate of the tail of the pulse, cross at  = L = 0.069. Notice that, despite their being equal at L , they stay real due to the special form of the matrix in Eq. (20). Throughout the parameter range shown in Fig. 7, the condition 1/2 < |λi |/|λj | < 2 applies for the combination i = 1, j = 2, but neither for i = 1, j = 3 nor for i = 2, j = 3. As |λ3 | is always at least twice as large as λ1 and λ2 , the corresponding influence can be neglected to leading order. The subsequent perturbation theory yields Eq. (19). The decay rates λ1 and λ2 in Eq. (19) are given by Eqs. (21), while the amplitudes A1 and A2 have to be computed from the shape of the solitary pulse (see Appendix A). Fig. 8 shows analytical approximations to the dispersion curves calculated with Eq. (42) derived in Appendix A, in comparison with the curve obtained with the methods described in Section 4 from the full partial differential equations (PDE) (2). The extremum of a dispersion curve approximated by Eq. (19) is situated at   1 A2 λ2 ln − , (22) dmax = λ1 − λ2 A1 λ1

M. Or-Guil et al. / Physica D 135 (2000) 154–174

167

Fig. 7. Eigenvalues (21) of the linearization of the set of ODEs (4) around the fixed point (0, 0, 0) as a function of the parameter . The single pulse velocity c was taken from the branch of pulse solutions HA1 , cf. Fig. 10. The positive (negative) eigenvalues give the spatial decay rates of the tail (front) of the pulses.

provided that the coefficients A1 and A2 have different signs. When |A1 λ1 | < |A2 λ2 | (|A1 λ1 | > |A2 λ2 |) holds, an extremum exists while λ1 > λ2 (λ2 > λ1 ). The value of dmax goes to infinity when, at a critical parameter value, either λ1 approaches λ2 or one of the Ai with i = 1, 2 changes sign. Here, we observe λ1 = λ2 at the value  = L (see Fig. 7). The extremum disappears below that critical value. Consequently, one expects that the maximum of the dispersion curve moves towards higher values of d for decreasing . A transition from anomalous to normal long range behavior is thus predicted to take place at  = L . The shift of the maximum of the dispersion curve to larger d is illustrated in Fig. 8. NED solutions and the emergence of NED branches out of an ED branch can also be discussed in this framework. When the maximum of the curve moves towards larger distances as  approaches L from above, the short distance db between the pulses of bound pulse solutions also gets larger (see Fig. 9). At the same time, their velocities approach the c∞ of the single pulse. The perturbation ansatz and the reduced dynamics are valid up to order of e−2λi db , where λi is the slower of the two relevant decay scales λ1 , λ2 .

Fig. 8. Dispersion curves c(d) predicted by perturbation theory with Eq. (42) (solid line) in comparison with curves calculated with the full equations (dashed line) for three different values of . When  is decreased, the maximum of the curves is shifted to higher values of d, until a transition from anomalous to normal behavior takes place at  = L = 0.069.

168

M. Or-Guil et al. / Physica D 135 (2000) 154–174

Fig. 9. Bound pulse pair solutions, approximated by non-equidistant two-humped pulse solution in a ring of length L = 70, for different values of the parameter . The “smaller” interpulse distance db increases as  decreases and goes to L/2 as  approaches the bifurcation point L

6. Pulse interaction and homoclinic bifurcations Stationary solutions of Eqs. (2) are also captured by the three-dimensional system Eqs. (4) of traveling wave ODEs, where the velocity c is considered as a parameter. There, periodic wavetrains correspond to limit cycles, and single resp. mult-pulse solutions to one-loop resp. multiloop homoclinic orbits at L → ∞. Thus, the existence range of bound states could be bordered by homoclinic bifurcations where the bound pair solutions (double loop homoclinics) are born from a branch of solitary pulses (single-loop homoclinics). In general, this will be a codimension-2 bifurcation of homoclinic orbits in the traveling wave ODEs. For a review on this topic, see Ref. [24]. Results on homoclinic bifurcations of higher codimensions have also been successfully applied in the study of the interaction of stationary, symmetric pulses in nonlinear optics [25]. For the model studied here, Zimmermann et al. [20] and Krishnan et al. [26] discuss existence and bifurcations of various solution branches in the ODE system (4). Fig. 10 shows some of these results for homoclinic branches and bifurcation points in the parameter plane (c, ). Here, the branch HA1 corresponds to an orbit homoclinic to the fixed point (u, w, v) = (0, 0, 0) (single pulses with rest state A in a system with infinite length). This branch transforms into the branch HB1 of homoclinic orbits to fixed point (b/a, 0, 0) (single pulses with unstable rest state B) upon increase of  through a heteroclinic bifurcation called 2 T-point [20]. Corresponding two-loop homoclinic orbits are denoted HA2 and HB2 . For comparison, the sites corresponding to solutions shown in the Figs. 1(a) and 6(b) for very large systems are denoted in Fig. 10 by triangles. Bound pulse pairs correspond to the branch HA2 . The bifurcating NED branch at  = 0.103 (Fig. 6(b)) ends up for large system size L in a solution asymptotically approaching the branch HB2 (two-loop homoclinic to the unstable fixed point B). A single pulse with rest state B (corresponding to the branch HB1 ) does not exist for these parameters. It is convenient, when working with the traveling wave ODEs, to use the speed c as a bifurcation parameter (and infer the ring-length L); this is exactly the converse of what we did in the pseudospectral PDE discretization case in Section 4. In this view, the maxima of the dispersion curves correspond to saddle-node bifurcations of limit cycles in the ODE system. The branching of NED pulse pairs from a site near these maxima corresponds then to period-doubling bifurcations of a limit cycle. The formation of a pulse triple is due to a degenerate bifurcation leading to a tripling of the spatial period, and so on. The sites of these saddle-node and period-doubling bifurcations are depicted likewise in Fig. 10. In the previous sections we have argued that bound states emerge due to a non-monotonic dispersion curve. The transition from monotonic to non-monotonic dispersion takes place when two relevant eigenvalues in Eq. (18) become equal. This situation corresponds to a bifurcation point (termed L-point) of pulses (resp. homoclinic orbits) in Eqs. (4) [20]. Our results suggest that at this L-point a variety of multiloop homoclinics bifurcate

M. Or-Guil et al. / Physica D 135 (2000) 154–174

169

Fig. 10. Homoclinic orbits of the traveling wave ODEs (4), corresponding to stationary solutions of Eqs. (2), in the parameter plane (c, ). SN (PD) denotes saddle-node (period-doubling) bifurcations of limit cycles. The other lines represent homoclinic orbits. The triangles denote the sites corresponding to the respective solutions for very large system lengths shown in the Figs. 1 and 6(b) ( = 0.1 and  = 0.103). The branches HA1 , HA2 and the bifurcation lines SN and PD do not represent the actually computed values but were slightly shifted in order to make them visible.

off the original single-loop homoclinic orbit branch corresponding to the single-humped pulse. The study of Zimmermann et al. [20] shows only the two-loop homoclinic orbit HA2 . The “newborn” solutions predicted by reduced dynamics correspond to multihumped pulses that represent bound states, where the number of pulses is equal to the number of loops of the homoclinic orbit. Since Eqs. (4) contain all possible solutions with any number of pulses, there is an infinity of homoclinic branches that emerge from the L-point. These branches can be labeled HA3,1 , HA3,2 , HA4,1 , HA4,2 , HA4,3 , HA5,1 . . . . So far, the main features of Fig. 10 conveniently summarize the main results obtained earlier from the reduced dynamics description of pulse interaction. A closer inspection reveals a number of important differences stemming from the approximative nature of the reduced dynamics. Besides the multiloop homoclinics (bound pulse states), a line of period doublings (origin of NED solutions) and a line of saddle-node bifurcations (maxima of the dispersion curve) emerge from the L-point. Reduced dynamics suggests that the velocity of bound states should be equal to the solitary pulse velocity c∞ . In the traveling wave ODE, Eq. (4), each of these bound states corresponds to a multiloop homoclinic orbit to the fixed point A. Since it is impossible to obtain multiple homoclinic orbits to the same fixed point for exactly the same parameters (here wave speed c), the predictions of reduced dynamics cannot be exact. In fact, perturbation theory here is accurate only up to an order O(e−2λi db ). Thus, we conjecture that the bound solutions branch HA2 (see Fig. 10) is situated at a distance |cb − c∞ | = O(e−2λi db ) from the single pulse branch HA1 . Similar arguments apply to the distances between other branches of homoclinic orbits (e.g., HA1 and HA3,1 , HA2 and HA3,2 and so on). A similar estimate holds for the distance between the period doubling and the saddle node branch; it should be of order O(e−2λi dmax ), corresponding to the aproximation error of reduced dynamics near the maximum of the dispersion curve. These conjectures cannot be easily checked, since the differences in velocity are very small, and the numerical results for very large systems may not be accurate enough.

7. Discussion We have studied the dynamics of interacting pulses in excitable media by a kinematic theory (derived from a perturbation ansatz) as well as numerical stability and bifurcation theory. We focused on media that exhibit a dispersion curve with a single maximum, connecting a region of anomalous dispersion at large wavelength with

170

M. Or-Guil et al. / Physica D 135 (2000) 154–174

one of normal dispersion at small wavelength. Kinematic theory predicts an instability of periodic trains for large d and bifurcations of new solutions at the maximum of the dispersion curve. Furthermore, it provides simple criteria for the existence and stability of non-equidistant wave trains, approaching bound pulse states for large total system lengths. The predictions of reduced dynamics have been tested and shown to be accurate for large enough interpulse distances through comparison with numerical computations of the full PDEs as well as with bifurcation results of the traveling wave ODEs. However, bound states often involve interpulse distances shorter than the validity range of kinematic theory. Then, the numerical tools of Section 4 have to be employed to accurately establish the existence and stability of bound states of pulses as well as other periodic and nonperiodic pulse trains. Similar methods have been previously used in the study of pulses for a modified Kuramoto–Sivashinsky equation (MKSE) that possesses an oscillatory dispersion curve ([27] and references therein). The MKSE, however, cannot form stable bound states of pulses, because it does not possess a stable homogeneous state (rest state) and consequently its solitary pulse solutions have eigenvalues with positive real parts in the continuous spectrum. Nevertheless, as long as the pulse density is large enough, an extended version of reduced dynamics assigning two degrees of freedom per pulse reproduces the dynamic evolution of waves in the MKSE satisfactorily [28]. Finally, we discussed the transition from normal monotonic dispersion, typically found in excitable media, to the case considered here: an anomalous dispersion curve with a single maximum. This transition can happen when two eigenvalues of the traveling wave ODEs, Eq. (4), become equal and the relevant decay scale of the pulse is no longer unique. The existence region of bound states should be bounded by homoclinic bifurcations where these solutions appear/disappear. In general, this will be a codimension-2 bifurcation of homoclinics in the traveling wave ODEs. Near this so-called L-point, dispersion curves are correctly reproduced by a two-scale perturbation theory, which we derived as an extension of the standard treatment. The analytical expressions for c(d) depend only on the solitary pulse solution in the medium. Anomalous long range dispersion is more likely to appear away from the typical excitable medium limit with fast activator and slow inhibitor (here   1). In that case, the decay dynamics is usually slaved to the slow inhibitor dynamics and arguments based on singular perturbation [4,5] theory will lead to normal dispersion for large d. In summary, our results suggest a simple strategy to search for the existence of bound pairs of traveling, strongly asymmetric pulses: (a) compute the dispersion curve c(d) and check if it is non-monotonic. (b) Trace the branches of solutions emerging from maxima of the dispersion curve in the traveling wave ODE (here Eqs. (4)). (c) If they bifurcate towards large d, check their behavior for d → ∞. If they approach the dispersion curve in that limit, they fulfill a necessary condition for existence of a bound state. (d) Check their stability by the described numerical method. In our experience, a combination of all methods is able to definitely establish the picture of the interaction dynamics of pulses as well as an answer to why and when bound states of pulses do form. The methodology used here in the study of reaction–diffusion systems might be helpful to investigate pulse dynamics and bound states in different fields of physics, e.g., nonlinear optics [29] or convection in fluids [30]. Acknowledgements IGK would like to acknowledge the hospitality of the Max Planck Institute for the Physics of Complex Systems in Dresden, as well as the partial support of NSF and of the Alexander von Humboldt Foundation. Appendix A Consider a system of the form: ∂t U = LU + N (U),

(23)

M. Or-Guil et al. / Physica D 135 (2000) 154–174

171

where U(x, t) ∈ Rm , x ∈] − ∞, ∞[, L is a linear operator which depends on ∂x and denotes the linear part of the dynamics; N (U) : Rm 7→ Rm denotes a purely nonlinear polynomial. We assume that Eq. (23) is invariant under translation, and that a perturbation corresponding to a translation of the solution is the only neutral mode. Reaction–diffusion systems are special cases of Eq. (23). We suppose that this system possesses a stable, single pulse solution H(z) with z := x − c∞ t, which moves with constant velocity c∞ . Now, we perform a multiple scale perturbation analysis to estimate the relevant pulse interactions. We consider a train of equidistant pulses (ED) going rightwards with large interpulse distances to be the unperturbed solution and show that reduced dynamics with merely nearest neighbor interactions follow from this ansatz. We follow the procedure of Elphick et al. [7], adopting their notation. However, the calculation is extended to pulses that possess different relevant spatial decay scales λi and λj . The inequality 1 λi < <2 (24) 2 λj is assumed to hold for all relevant pairs i, j . These scales are obtained from the real parts of the eigenvalues of the linearization around the rest state, as discussed in Section 5. We show that, if (24) holds, a perturbation analysis including these scales can be performed to low order. Then, influences of interactions on the pulse shape and on pulses other than next nearest neighbors can be neglected in the multiscale ansatz. Thus, we still obtain nearest neighbor interactions, as is the case with one spatial scale. We suppose, without loss of generality, that the rest state into which the pulse decays is U = 0. Furthermore, it is assumed that there are only two relevant spatial decay scales λ1 , λ2 with the property (24). Both λ1 and λ2 are assumed to be positive and real. Thus, we study a case where the only relevant eigenvalue pair λ1 and λ2 describes the back of the pulse. An extension to the case of oscillatory decay as well as to more than two decay scales of the back and/or the front of localized excitations with the property (24) is straightforward. We start with an ED pulse train as an ansatz for infinitely large interpulse distances and calculate weak perturbations for large but finite distances, showing then that to lower order, only nearest neighboring pulses contribute to the relative pulse dynamics. An ED N-pulse train can be approximated by a superposition of N single pulses. To account for weak, distance dependent interactions, we introduce the small parameter µ := e−λ1 d .

(25)

With α :=

λ2 , λ1

(26)

it follows that µα = e−λ2 d ,

(27)

where 1/2 < α < 2 holds. We introduce two slow time scales t1 := µt,

t2 := µα t,

(28)

so that ∂t → µ∂t1 + µα ∂t2 .

(29)

The ansatz in the comoving frame z = x − c∞ t is: U=

N X n=1

Hn + µR1 (z, t1 ) + µα R2 (z, t2 ),

Hn := H(z + nd − 8n (t1 , t2 )).

(30)

172

M. Or-Guil et al. / Physica D 135 (2000) 154–174

The functions Hn denote single pulses in different positions. The function 8n with 8n (0, 0) = 0 accounts for deviations from the original pulse positions in a periodic pulse sequence. The functions Ri account for small deviations from the shape of the solitary pulse. In order to avoid ambiguities, we require the Ri to be orthogonal to each other. Inserting Eq. (30) into Eq. (23) and using ∂t → ∂t − c∞ ∂z and ∂x → ∂z in a comoving frame, we obtain the following equation: N N X X 0 α ∂t2 8n Hn0 −µc∞ ∂z R1 − µ c∞ ∂z R2 − µ ∂t1 Φn Hn − µ α

n=1

= µLR1 + µα LR2 + N

n=1

N X R 1 + µα R 2 H n + µR

!

n=1 N X

−

N (Hn ) + O(µ2 ) + O(µ2α ).

(31)

n=1

In order to sort the terms of Eq. (31) according to their different orders in µ, it is necessary to focus on the underlined, nonlinear term. Consider its Taylor-series: ! ! N N X X α Hn + µR1 + µ R2 = N H n + µ∇Ned R1 + µα ∇Ned R2 + O(µ2 ) + O(µ2α ), (32) N n=1

n=1

where the Jacobian P ∇N isa m × m-matrix denoting the functional gradient with respect to the argument. The term N ∇Ned = ∇N n=1 Hn is calculated around the superposition of equally spaced pulses. This Jacobian can be estimated using the linearization around the single pulse: for a given site z near the nth pulse, we obtain, from a Taylor expansion, Ned (z) = Nn (z) + O(µ), where ∇Nn = ∇N (Hn ). As Nm (z) = O(µ) if z is far from the mth pulse, we can thus write for all z: ∇Ned =

N X

∇Nn + O(µ).

(33)

n=1

Consider now the underlined term of the r.h.s. of Eq. (32). For it, ! N N N N X X X X Hn = N (Hn ) + ∇Nn Hk + . . . N n=1

n=1

(34)

n=1k=1,n6=k

holds. We focus on the second term in the r.h.s. of this equation. Regarding the shape of a single pulse helps to pick the expressions which are of order lower than O(µ2 ) + O(µ2α ). Relative to the site zn of the nth pulse, the vector Hk can be approximated by Hk = a1 eλ1 (z−zn −|k−n|d) + a2 eλ2 (z−zn −|k−n|d) + O(µ2 ) + O(µ2α ), where ai =



ai0 , z < zn , 0, z ≥ zn

(35)

(36)

a10 and a20 are suitable, constant vectors in Rm which can be found by analyzing some point far away from the pulse position. Only when k = n − 1 holds, do we obtain in Eq. (35) terms of lower order than O(µ2 ) + O(µ2α ). When considering that, Eq. (34) can be written as

M. Or-Guil et al. / Physica D 135 (2000) 154–174

173

!

N

N N N X X X Hn = N (Hn ) + ∇Nn Hn−1 + O(µ2 ) + O(µ2α ) n=1

n=1

=

n=1

N X

N N X X ∇Nn a1 eλ1 (z−zn ) + µα ∇Nn a2 eλ2 (z−zn ) + O(µ2 ) + O(µ2α ).

n=1

n=1

N (Hn ) + µ

(37)

n=1

Inserting Eqs. (37) and (33) in Eq. (31) and separating the terms of orders O(µ) and O(µα ), we obtain for each n two independent linear equations of the form: (−c∞ ∂z − L − ∇Nn )Ri = ∇Nn ai eλi (z−zn ) + ∂ti 8n Hn0

(38)

with i = 1, 2. Because of the existence of a one-dimensional kernel spanned by the Goldstone mode Hn0 , the linear operator Hn := −c∞ ∂z − L − ∇Nn is non-invertible. Thus, the condition to avoid secular terms when solving † Eq. (38) requires its r.h.s. to be orthogonal to the kernel of the adjoint operator Hn . This operator is defined by the † equation hU|Hn Vi = hHn U|Vi, where h|i denotes spatial integration over the whole space. † Suppose the kernel of Hn is spanned by a vector Pn 6= 0. Then, the solvability condition demands orthogonality between the rhs of Eq. (38) and Pn up to an error of order O(µ). Carrying out the projection, we obtain, ∂ti 8n =

hPn |∇Nn ai eλi (z−zn ) i + O(µ) := Ai + O(µ) hPn |Hn0 i

(39)

˙ n = µ∂t1 8n + µα ∂t2 8n . This is equivalent for i = 1, 2. Transforming back to the original time variable t, we get 8 to: ˙ n = A1 e−λ1 d + A2 e−λ2 d + O(µ2 ) + O(µ2α ). 8

(40)

These results indicate that, at leading order, all deviations of pulse positions in the reference frame moving with the solitary pulse speed c∞ evolve at a constant rate that only depends on the spatial period d. This is usually described by the dispersion curve in the laboratory frame. Switching coordinates back to the original frame, we obtain Eq. (19) ˙ n + c∞ . c(d) = 8

(41)

We can use Eq. (40) in a reduced description that is based on the dispersion relation, as long as all di are large enough to make µi = e−λi di “small” parameters. Using Eq. (35), we can replace Eq. (40) by ˙n = 8

hPn |∇Nn Hn−1 i + O(µ2 ) + O(µ2α ), hPn |Hn0 i

(42)

thus avoiding the explicit computation of the vectors ai . Note that for the functions Ri to depend only on the slow time scales ti , one has to demand them to be orthogonal to the vector Pn . In that way, they do not affect the dynamics of pulse position. In our case, the pulse solution Hn , its linearization ∇Nn and the kernel Pn are calculated numerically in a discretized representation of Eqs. (2). The leading order term of the rhs of Eq. (38) is then compared to the dispersion curve (cf. Eq. (37)). The results are shown in Fig. 8. In case λ1 and λ2 have different signs, not only the decay of the tail but also that of the front of the localized structure have to be considered. Accordingly, the impact of the pulse Hn+1 on Hn has to be included to leading order in Eq. (34), yielding dynamics which depend on the distances to both nearest neighbor pulses.

174

M. Or-Guil et al. / Physica D 135 (2000) 154–174

References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30]

J. Ross, S.C. Müller, C. Vidal, Science 240 (1988) 466. R. Kapral, K. Showalter (Eds.), Chemical Waves and Patterns, Kluwer Academic Publishers, Dordrecht, 1995. S. Jakubith, A. von Oertzen, H.H. Rotermund, G. Ertl, Phys. Rev. Lett. 65 (1990) 3013. J. Keener, J. Tyson, Physica D 32 (1988) 327. E. Meron, Phys. Rep. 218 (1992) 1. C. Elphick, E. Meron, E.A. Spiegel, Phys. Rev. Lett. 61 (1988) 496. C. Elphick, E. Meron, E.A. Spiegel, SIAM J. Appl. Math. 50 (1990) 490. C.P. Schenk, P. Schütz, M. Bode, H.-G. Purwins, Phys. Rev. E 57 (1998) 6480. J. Rizel, K.J. Maginu, in: C. Vidal, A. Pacault (Eds.), Non-Equilibrium Dynamics in Chemical Systems, Springer, Berlin, 1984. C. Elphick, E. Meron, J. Rinzel, E.A. Spiegel, J. Theor. Biol. 146 (1990) 249. C. Elphick, G.R. Ierley, O. Regev, B.A. Spiegel, Phys. Rev. A 44 (1991) 1110. A.T. Winfree, Physica D 49 125 (1991). F. Siegert, C.I. Weijer, Physica D 49 (1991) 224. H. Linde, M.G. Velarde, A. Wierschem, W. Waldhelm, K. Loeschcke, A.Y. Rednikov, J. Coll. Int. Sci. 188 (1997) 16. J. Christoph, M. Eiswirth, N. Hartmann, R. Imbihl, I.G. Kevrekidis, M. Bär, Phys. Rev. Lett. 82 (1999) 1586. M. Or-Guil, M. Bode, C.P. Schenk, H.-G. Purwins, Phys. Rev. E 57 (1998) 6432. D. Barkley, Physica D 49 (1991) 61. M. Bär, N. Gottschalk, M. Eiswirth, G. Ertl, J. Chem. Phys. 100 (1994) 1202. J. Rinzel, G.B. Ermentrout, in: C. Koch, I. Segev (Eds.), Neuronal modeling: From Synapses to Networks, MIT Press, Cambridge, 1989. M.G. Zimmermann, S.O. Firle, M. Natiello, M. Hildebrand, M. Eiswirth, M. Bär, A. K. Bangia, I.G. Kevrekidis, Physica D 110 (1997) 92. M. Bär, M. Or-Guil, Phys. Rev. Lett. 82 (1999) 1160. E.J. Doedel, X.J. Wang, AUTO 94: Software for continuation and bifurcation problems in ordinary differential equations, Tech. rep., Center of Research on Parallel Computing, California Institute of Technology, Pasadena, CA 911125, CRPC 95-2. J.M. Flesselles, A. Belmonte, V. Gaspar, J. Chem. Soc. Faraday Trans. 94 (1998) 851. A.R. Champneys, Y.A. Kuznetsov, Int. J. Bifurcation and Chaos 4 (1994) 785. B. Sandstede, C.K.R.T. Jones, J.C. Alexander, Physica D 106 (1997) 167. J. Krishnan, I.G. Kevrekidis, M. Or-Guil, M.G. Zimmermaun, M. Bär, Comput. Meth. Appl. Mech. Engrg. 170(3)(4) (1999) 253. N.J. Balmforth, G.R. Ierley, R. Worthing, SIAM J. Appl. Math. 57 (1997) 205. H.C. Chang, E.A. Demekhin, E. Kalaidin, SIAM J. Appl. Math. 58 (1998) 1246. N. Akhmediev, A. Ankiewicz, J.M. Soto-Crespo, Phys. Rev. Lett. 79 (1997) 4047. P. Kolodner, Phys. Rev. A 44 (1991) 6466.