Breakup of spiral waves caused by radial dynamics: Eckhaus and finite wavenumber instabilities

In this paper, we link numerical observations of spiral breakup to a stability analysis of simple rotating spirals. We review the phenomenology of spiral breakup, important applications in pattern formation and the state of the art in numerical stability analysis of spirals. A strategy for the latter procedure is suggested. Phenomenologically, spiral breakup can occur near the centre of rotation (‘core breakup’) or far away from it (‘far-field breakup’). It may be accompanied by instabilities of the spiral core in particular spiral meandering that affect also the stability of waves in the far-field, because an unstable core acts as a moving source and introduces a (nonlinear) Doppler effect. In general, breakup of non-meandering spirals is related to an absolute instability of the planar wavetrain with the same wave number. To simplify the stability problem, we consider a one-dimensional cut (‘1D spiral’) with a fixed core position in simulations and compare the results with a stability analysis of planar wavetrains. These 1D spirals approximate the radial dynamics of non-meandering 2D spirals. To fully account for instabilities of 1D spirals, it is not sufficient to compute the direction of propagation of the unstable modes of the wavetrains, one also needs to compute the so-called absolute spectrum of these wavetrains. This allows us to decide whether the instability is of convective or absolute nature. Only the latter case implies an instability of spirals in finite domains. We carry out this programme for the case of core breakup in an excitable reaction–diffusion system, the modified Barkley model. Our analysis yields that core breakup can result from the absolute variant of a novel finite wavenumber instability of the radial dynamics where the critical perturbations are transported towards the core. From these results, we can confirm that a simple spiral breaks up if the wavetrain in the far-field is absolutely unstable. The central new result is the discovery of respective convective and absolute instabilities of wavetrains connected with modes of finite wavenumber that propagate in the direction opposite to the wavetrain. Hence, in spirals, perturbations are moving inward and core breakup becomes possible even in the absence of meandering.


Introduction
Rotating spiral waves are frequently observed in homogeneously and heterogeneously catalysed chemical reactions [1,2] and various biological systems, namely slime mold aggregation [3], cardiac tissue [4,5] and calcium waves in frog eggs [6]. In dissipative (quasi two-dimensional) chemical systems the rotating spirals appear as iso-concentration contours of the spatial distribution of the reactants [7,8]. The temporal evolution of the concentration patterns has been modelled by partial differential equations of reaction-diffusion type. Such models include oscillatory, excitable or bistable systems with either none, one or two linearly stable homogeneous states [9,10].
Patterns form either due to the instability of a steady state in oscillatory media or due to local suprathreshold, finite-amplitude perturbations of homogeneous stable steady states in bistable and excitable media. In the latter case, they have a stable rest state and respond to a suprathreshold perturbation with an amplification followed by saturation and recovery. The existence of an excitation threshold is a generic feature. A popular model for oscillatory systems is the complex Ginzburg-Landau equation (CGLE) [11] that is suited to describe media near a supercritical Hopf bifurcation [12]. The so-called activator-inhibitor systems are widely used to model excitable media [13].
In this paper, we analyse instabilities of rotating spirals that lead to spatiotemporally chaotic dynamics [14] via spiral breakup. The problem has received considerable attention for various reasons. For instance, the observation of 'defect-mediated' turbulence in numerical simulations of the CGLE has been attributed to the breakup of spirals [15]- [20]. Until the early 1990s, it was unclear if spatiotemporal chaos or irregular activity is possible in homogeneous two-dimensional (2D) continuous excitable media. The first observation of chaotic wave patterns were reported in discrete models [21,22]. The availability of faster computers led to the discovery that various models exhibit spatiotemporal chaos. Among them are models of cardiac tissue [4], [23]- [26] and activator-inhibitor models of FitzHugh-Nagumo type designed to capture essential aspects of pattern formation [27]- [30].
An important motivation for the study of excitable media has been the quest for the cause of irregular electrical activity in cardiac muscle [31]. Experiments in thin sheets of heart tissue displayed only stable spirals in contrast with the irregular activity seen in experiments with whole hearts [32]. Consequently, it has been suggested that irregular activity in the heart might be a genuinely 3D phenomenon [33]. Thus more realistic three-dimensional, anisotropic models of the heart and excitable media have been investigated and revealed various sources of irregular activity on the surface including intricate dynamics of scroll waves [34,35] and the analogue of breakup in three dimensions [36,37]. The reason for the onset of ventricular fibrillation as well as possible treatments still remain subjects of intense experimental and theoretical research [38]. In pattern-forming chemical reactions, progress in the design of open reactors has finally also yielded experimental results that demonstrate a controlled transition to spatiotemporal chaos via spiral breakup in the Belousov-Zhabotinsky reaction [39]. Additional examples of transitions from spiral patterns to irregular spatial organization have been reported in catalytic surface reactions [40].
In what follows, we shall concentrate only on destabilization against modes in the radial direction.
It is important to note that in all mentioned examples of radial spiral breakup in models and experiments two different scenarios are observed-spirals may break first close to their core or alternatively far away from the core [45]. Breakup near the core is found in simulations in excitable media [24,27,29] and in experiments with a chemical reaction [50], whereas breakup far away from the core is typically seen under oscillatory conditions both in chemical experiments [39] and in simulations of the CGLE [18,19]. Recent studies of a model of intracellular calcium waves show breakup far away from the core under excitable conditions [51].
Here, we will present results on a simple activator-inhibitor model that exhibits both types of breakup as well as oscillatory and excitable kinetics depending on the chosen control parameters [29,45]. After we have discussed the phenomenology of breakup in this model, we will turn to a stability analysis of periodic waves in two variants of the model. The rationale of the treatment of periodic waves is as follows: stable rotating spirals do emit periodic waves, and far away from the core the concentration patterns resemble planar periodic wavetrains. If these wavetrains are unstable, then the spirals may also disappear (breakup).
For the CGLE it is possible to compute the stability of periodic wave solution analytically, because the waves and the corresponding eigenfunctions have the form of Fourier modes [17]. For given parameters, periodic waves in the CGLE become unstable below certain wavelength through an Eckhaus instability. If the spatially homogeneous oscillations becomes unstable to spatial perturbations, one uses the term Benjamin-Feir instability. Eckhaus instabilities first appear through long-wavelength perturbations that lead to spatially slowly varying modulations of the periodic waves. In the nonlinear stage, Eckhaus unstable wavetrains may evolve to modulated amplitude waves (MAWs). In the CGLE, the Eckhaus instability is of convective nature and several groups have pointed out that spiral breakup requires an absolute instability of the emitted wavetrain [17,18,19].
In excitable media, stability analysis has been originally restricted to the kinematic treatment of the phenomenon of alternans [46,47]. Alternans is a period-doubling instability of periodic waves, wherein the width of excitation pulses oscillates between two different values. It has been identified to play an important role in the process of breakup in models of cardiac tissue like the Beeler-Reuter model and some simplifications [4,24]. Analytical approaches are also successful in the case of very sharp interfaces, respectively large separation of time-scales of the activator and inhibitor [52]. When an analytic stability analysis is not feasible, then the numerical analysis of the 1D asymptotic wavetrain provides important insight into the basic mechanism of spiral breakup. It has been shown that the stability properties of 2D spirals related to their continuous spectrum in unbounded domains and periodic 1D waves are identical up to a transformation of the imaginary parts of the eigenvalues [53].
The continuous spectrum of spirals in large but bounded finite domains is however different from the continuous spectrum of spirals in unbounded domains. It is, nevertheless, still related to the stability of the corresponding 1D wavetrains. In mathematical terms, the relevant spectrum is the absolute spectrum of these wavetrains, whereas the usual stability analysis of wavetrains in an unbounded system yields the so-called essential spectrum. For more details, we refer to work by Sandstede and Scheel [53,54]. In practice, instabilities of the absolute spectrum often indicate an absolute instability of wavetrains, whereas instabilities of the essential spectra already capture convective instabilities of the wavetrains. It is worth noting that absolute and essential spectra of a given wavetrain are obtained from the same linearized equation. The differences between absolute and essential spectra are caused only by the different boundary conditions that in turn lead to different eigenfunctions; for a simple example, see [55].
Potential discrepancies may arise from nonlinear behaviour not captured by the linear stability analysis or instabilities of the discrete spectrum, like spiral meandering [10,48], that do not have a counterpart in one dimension. Numerical stability analysis has been crucial in the understanding of the meander instability of rotating spirals. It was identified as a Hopf bifurcation that introduces a second frequency apart from the rotation frequency in the spiral movement [49]. The spiral tip no longer follows a circular trajectory during spiral meandering, but instead gives rise to flower-like trajectories. However, numerical stability analysis in two dimensions is computationally expensive and in practice often restricted to rather small domain sizes (typically up to system sizes of a few spiral wavelength) as well as to the iterative computation of the largest eigenvalues.
The aim of this work is twofold: first, we like to review earlier work by us and other groups that dealt mostly with the analysis of the Eckhaus instability of the plane waves in the far-field and its implication for spiral breakup. Second, we present an new instability of the planar waves to finite wavenumbers (in contrast with the zero wavenumber of the critical Eckhaus mode) that appears almost simultaneously with an Eckhaus instability in the model under consideration and show that this instability becomes absolute earlier than the Eckhaus modes. Therefore, it is more relevant to spiral breakup and the instability of waves in 1D sources. Our approach is confined to spiral breakup as an instability of the continuous spectrum of a simple rotating spiral and does not give a full explanation of breakup phenomena in meandering spirals. In the latter case a linear stability analysis has not been achieved, neither analytically nor by numerical means. Breakup of meandering spirals does however often involve a subtle interaction between the modulation of waves caused by the Doppler effect [29] created by the spiral core motion and the stability of plane waves in the far-field. Therefore, an analysis of the latter property of a spiral (such as the ones presented here) offers at least a first step to a full understanding.
The paper is organized as follows. In section 2, we introduce a modified version of the Barkley model which exhibits spiral breakup and describe the different phenomenologies of spiral breakup arising therein. The main distinction between far-field and core breakup is also reproduced in a 1D source with stationary 'core' position. In section 3, we perform numerical stability analysis of planar wavetrains and compare the results with the observations in simulations of 1D sources. Two main reasons for the spiral instability are found: the absolute Eckhaus instability where the perturbations travel away from the spiral core and a novel finite-wavelength instability where perturbations travel towards the core. The latter instability causes spiral breakup near the core, whereas the absolute Eckhaus instability produces far-field breakup for sufficiently large group velocities of the outward propagating fastest growing modes. Finally, we summarize and discuss the results in section 4.

Phenomenologies of spiral breakup
To study spiral breakup, we use here a reaction-diffusion model that was originally introduced by Barkley [56] but with modified reaction kinetics [57]. It describes the interaction of a fast activator (u) and a slow inhibitor (v) variable: The form of h(u) describes a delayed production of the inhibitor and the equations have been used to model patterns in a catalytic surface reaction [57]. The change of h(u) from the standard choice h(u) = u [56] leads to the possibility of spatiotemporal chaos in 2D due to spiral breakup for > BU [29,58]. For smaller , the system settles into stable rotating, respectively travelling, waves. The parameter choice a < 1 yields excitable (oscillatory) behaviour for b > 0 (b < 0) and 0 < 1. In both cases, an unstable focus exists with (u, v) = (u 0 , v 0 ) in the local dynamics of equations (1). In the excitable case, two more fixed points appear: (u, v) = (0, 0) is the stable rest state and (u, v) = (b/a, 0) is a saddle that marks the threshold of the excitable medium. Throughout this study, a is fixed to 0.84 and b and are varied.
Zero-flux boundary conditions have been employed for simulations in 2D. We complement these by simulations of 1D sources where a fixed boundary condition (Dirichlet) is chosen such that it emits a wave of identical wavelength as that measured in 2D simulations. This intermediate step between simulations in 2D and stability analysis of periodic waves in 1D enables us to discriminate between contributions from the spiral's core (cause discrepancies between 1D and 2D simulations) and its radial dynamics (agreement of 1D and 2D simulations). Furthermore, the influence of secondary spirals is suppressed, which in 2D are created in the breakup process of the initial main spiral. The simulations of 1D sources may also reveal possible nonlinear effects as discrepancies to the results of the linear stability analysis. Moreover, the linear stability analysis is performed in 1D systems with periodic boundary conditions. We perform these three steps of our strategy for a large number of values of the control parameters.  In the following, results of numerical simulations for the model are described. We have investigated the dynamics of spiral waves as a function of the excitation threshold parameter b and the ratio of the timescales of the local dynamics. Almost independent of the choice of b, we find a transition from rotating spiral waves to spatiotemporal chaos via spiral breakup at values between 0.07 and 0.08 in the excitable as well as the oscillatory regime [29,45]. The simulations are usually started with a spiral initial condition that has been obtained at nearby values of . The first run has been performed at = 0.02 and then is increased in small steps ( = 0.0025) until destabilization occurs. In the vicinity of the spiral breakup we have used much finer 'resolution' in . The type of destabilization of spirals just above the critical , however, depends on the choice of b.
For the given model we find two different scenarios that are exemplified in figure 1. For excitable conditions, spiral breakup appears first close to the centre and the irregular pattern spreads then outward (cf figure 1(a)). For oscillatory conditions, we observe a different behaviour. Spirals first break up far away from the centre and eventually relatively large spiral fragments surrounded by a 'turbulent' bath remain (cf figure 1(b)). The size of the surviving part of the spiral shrinks if is further increased until finally no trace of the initial spiral persists in the long run.
At closer inspection we observe that, even inside the surviving part of the initial spiral, modulations in wavelength and period appear. The amplitude of these modulation grows with growing distance to the core. With increasing the modulation gets stronger and the critical modulation amplitude necessary to cause breakup is reached at distances closer to the centre. In addition, the breakup in the excitable regime is usually preceded by a meander instability of the rotating spirals, while for oscillatory conditions breakup bounds a region of stable spiral rotation [29]. The breakup far away from the core has been also observed in experiments in the Belousov-Zhabotinsky reaction [39] and in numerical simulations of the CGLE [18,19,39]. More recent studies provide also an experimental verification of the core breakup scenario; see [59].
As discussed in the introduction, further insight into the nature of the spiral instability may be obtained by a reduction to one spatial dimension. Here, we demonstrate this procedure for equations (2). Therein, a 1D wave source is studied as an analogue of a spiral in 2D. We will show that the essential properties of spiral breakup depend mostly on the asymptotic selected wavetrain and not on the detailed structure of the source (e.g. the spiral core geometry or the centre of a target pattern). Results for 1D sources and 2D spirals should allow for a quantitative comparison if their selected wavelengths are the same. A 1D source of periodic waves is given by the equations Dirichlet (u(0) = u 1 ) and zero-flux (∂u(r)/∂r| L = 0) boundary conditions are used at respective ends of the 1D system of length L. Equations (2) represent the radial part of equations (1) far from the core, where curvature effects can be neglected. The Dirichlet boundary at r = 0 constitutes a source of waves similar to the core of the spiral in 2D. By changing u 1 the wavelength of the emitted wavetrain can be varied. Here, we picked a value of u 1 that selects roughly the same wavelength as spirals in 2D just below the breakup. In the oscillatory case, u 1 = u 0 achieved that goal. Space-time plots from the integration of equations (2) are presented in figure 2 both for oscillatory and excitable conditions. In both cases, the wavetrains emitted from the left boundary at r = 0 exhibit an instability upon increase in at a critical value C . For b = −0.045, we find a scenario similar to the behaviour in 2D (cf figures 2(a) and (b)). In the excitable case with b = 0.07, the instability appears always very close to the source (see figures 2(c) and (d)).
Breakup occurs first far away from the Dirichlet boundary, i.e. the source of waves. Upon increase of the breakup of the wavetrain moves towards the source until it gets to a distance of about one wavelength. The wavelength λ dir and period τ dir of the waves selected by the Dirichlet boundary is for values near the breakup almost the same as the selected wavelength λ 0 and period τ 0 of the spirals in 2D in line with earlier results in the CGLE [17]. Nevertheless, the critical value for the 1D source C ≈ 0.081 > BU ≈ 0.074. The transition between the two scenarios occurs around b = 0.04. The values of C are not very sensitive to changes in b and are slightly larger than 0.08. Thus, we find the counterpart of breakup near In (a) a stable wavetrain is established in the long run; in (b) waves break roughly three wavelengths away from the boundary. and far away from the spiral core in the 1D model. In particular, the appearance of different qualitative scenarios does not depend on specific 2D ingredients like meandering and curvature.

Numerical stability analysis of spirals near breakup
Let us consider the basic arguments in spiral stability analysis. The spiral selects, for a given set of parameters, a particular wavelength λ 0 and a related temporal period τ 0 . At the same moment, there exists a one-parameter family of periodic wave solutions with varying speed c, temporal period τ and wavelength λ = cτ. The one-parameter family may be described by a dispersion curve c(λ). At small wavelength, the wavetrains either exhibit an instability (Eckhaus in the CGLE, alternans in cardiac tissue models) or cease to exist (saddle-node or drift pitchfork bifurcations in FitzHugh-Nagumo systems). These minimum stable wavelengths and periods shall be called λ min , respectively τ min , in the following. In a few cases (CGLE [17], Rinzel-Keller model [60]), it is possible to compute the instability or bifurcation analytically. In other cases, kinematic theories (alternans [46,47]) or approximations by singular perturbation theory (excitable media) may yield useful information on the nature of this instability. For the models treated here, none of these methods are sufficient to yield an understanding of the instability at λ min . This calls for alternative methods, namely numerical stability analysis. Here, we will mostly present quantitative results found by numerical stability computations. For excitable media given by reaction-diffusion equations of the form presented in equations (1), a large body of theoretical work has been done in the limit of very small . In particular, scaling laws originally predicted by Fife [61] have been explicitly derived for equations of the type used here by Karma [62] with the constraint 1/3 1. Applied to the model studied here, the scaling argument goes as follows: for small enough , the selected spiral wavelength scales like λ 0 ∝ α 1 , whereas the smallest wavelength that is stable behaves like λ min ∝ α 2 . Usually, it is found that α 2 > α 1 and for = C the condition λ 0 = λ min is fulfilled. For > C , the selected wavelength corresponds to an unstable wavetrain in the far-field and breakup should result. The analytically found scaling laws in the present case follow qualitatively the same trend as the numerically obtained data. They are, however, not very accurate quantitatively and cannot be used to estimate the value for C . Hence, one has to carry out the stability analysis numerically as we will show below. Altogether, it is not too surprising that analytical scaling arguments fail quantitatively, because the breakup occurs at quite large values of . Scaling arguments have been used with more success in the spiral breakup in cardiac tissue, where the instability appears upon decrease of and the scaling properties of the spiral waves are qualitatively different from the standard excitable media investigated here [24,25].
In the following, we will numerically compute the 1D wave solutions and their linear stability along the dispersion curve. Continuation software as AUTO has been used to obtain the spatial profiles of periodic waves as stationary solutions u(z), v(z) in a co-moving frame z = r−ct [63]. Consequently, we compute periodic orbits that correspond to a single wavelength with periodic boundary conditions. The evolution of small perturbations of a spatially periodic travelling wave U 0 (z) in the reaction diffusion equations is described by W jn (z)e ω jn t . Insertion of the perturbed wave solution into equation (2) allows us to expand the equation in powers of the small perturbation which in linear order poses an eigenvalue problem for the eigenvectors W jn and eigenvalues ω jn . The indices j, n will be specified below. There are infinitely many eigenvalues and eigenfunctions. For the system length L → ∞, ω jn are located on continous curves in the complex plane. The nature of the posed problem requires that L is an integer multiple of the wavelength of the periodic wavetrain, i.e. L = Nλ. Consequently, the wavetrain in the system of length L consists of N maxima.
The solution is stable when all eigenvalues ω jn have negative real parts. The wavelength λ is used as a bifurcation parameter. The method yields accurate information on the minimum stable period λ min for the given medium, the eigenvalue spectra and the eigenfunctions describing the dynamics of perturbations to the wavetrains. These results are then compared with the selected wavelength of the 1D sources λ dir near the breakup threshold.
If we consider periodic travelling waves with constant shape U 0 (z) and constant speed c in 1D, a few facts on the eigenvalues and the eigenfunctions are known a priori. Since our media are homogeneous and translation-invariant, any translation of a given periodic solution is also a valid solution with identical stability properties. This translational symmetry of the wavetrains is reflected by an eigenvector W 00 = dU 0 /dz with zero eigenvalue ω 00 = 0 (Goldstone mode). Symmetry arguments require the eigenfunctions of the periodic operator, obtained by linearization around a wavetrain with wavelength λ, to be Bloch functions W jn (z) = e i(2πn/L)z Φ jn (z) = e ikz Φ jn (z) with Φ jn (z) = Φ jn (z + λ) and n = 0, . . . , N − 1 [64]. The above Bloch form can be used in an ansatz that reduces the stability problem for the infinite domain to a stability problem in a domain of length λ with the wavenumber k = 2πn/L as an additional parameter.
So far, we have assumed a purely imaginary exponent ikz in the prefactor of the Bloch eigenfunctions. With this constraint, one obtains the so-called essential spectrum. To assess the stability of spirals in finite domains, one needs to compute also the absolute spectrum. Below, we will do that by employing an ansatz W jn (z) = e (µ+ik)z Φ jn (z) for the eigenfunction.
In general, these exponentially weighted eigenfunctions yield spectra depending on the exponential weight, i.e. ω jn = ω jn (µ). The group velocities of the associated modes are then ∂Im(ω jn (µ))/∂k = ∂Re(ω jn (µ))/∂µ. The sign of the group velocities tells us in which direction modes are travelling. In addition, we can compute the most relevant parts of the absolute spectra ω jn (µ * ) by requiring ∂Re(ω jn (µ))/∂µ| µ=µ * = 0 and verifying the pinching condition [53]. In this fashion, we obtain the growth rates of the modes with zero group velocity that determine the stability in a finite system. We found further contributions ω jn (µ * ) with non-zero group velocity to be located far inside the stable half-plane which are not relevant here.

Eckhaus instability and spiral breakup scenarios
We start by computing the spectra of periodic waves for equation (1); for a more detailed discussion see [45]. The eigenvectors corresponding to the eigenvalues with the largest real part are modulations of the Goldstone mode of the approximate form W 0n ≈ e i(2πn/L)z W 00 (z). The amplitudes of the eigenfunctions are largest in the fronts and backs of the pulses in the wavetrain in contrast with the Fourier eigenfunctions in the CGLE [17]. Thus, the fastest growing modes correspond to an alternating compression and expansion of subsequent pulses in line with the observations in the 2D simulations of figure 1 and experiments in the Belouzov-Zhabotinsky reaction [39]. The calculations for both models reveal an instability reminiscent of the Eckhaus instability in the CGLE. As the wavelength λ is shortened, the spectra of the wavetrains shift towards larger real parts and cross the imaginary axis at λ min . The leading part of the spectrum switches from a parabola opening to negative real parts and touching zero to a segment of a parabola opening to positive real parts. A symmetric curve results with two positive maxima of the real part of ω jn at non-zero imaginary values; for an illustration see figure 4 in [45].
The dispersion curves describe the dependency among speed c, wavelength λ and period τ. It is sufficient to plot two of these three quantities, since λ = cτ. Here, we plot c versus λ and display stable (unstable) wavetrain solutions along the dispersion curve with full (dashed) lines in figure 3. The dispersion curve ends towards small wavelength and periods at a non-zero value of the speed c and neither at a saddle-node bifurcation nor at zero velocity as known from other excitable media [10,60]. Closer inspection reveals that the wavetrains instead bifurcate from one of the additional unstable equilibria with zero amplitude. Such a scenario may be typical for excitable media with additional unstable fixed points and is of course impossible in models which just have the rest state of the medium as a fixed point. The Eckhaus-type instability, however, appears somewhere along the dispersion curve and cannot be deduced from the form of the dispersion curve alone (as could a saddle-node bifurcation). The instability often happens just before the speed has a local minimum along the dispersion curve (see figure 3). The goal of the previous analysis has been to provide a comparison between the minimum stable wavelength λ min and the selected wavelength λ 0 and λ dir found in simulations near spiral breakup and its one-dimensional analogue. Now, we can compare the results of the simulations with the Dirichlet source and the stability analysis of wavetrains. Figure 4 shows the comparison between λ min (red dashed line, triangles) and λ dir , the selected period of the Dirichlet source (blue full line, circles). The upper panel of figure 4 shows the data for the oscillatory case (b = −0.045). There is a narrow interval where stable 1D sources and spirals are found with λ 0 resp. λ dir < λ min . This is related to the convective nature of the Eckhaus instability. Using the above outlined method of exponential weights to compute the absolute instability line, Sandstede and Scheel [53] could indeed confirm  this picture. Perturbations in the regime of convective Eckhaus instability are simply advected out of the system boundaries. The radial dynamics of spirals in 2D are analogous to the Dirichlet case; the amplitude of the waves in the spiral approaches zero in the centre of rotation. This corresponds to a 'self-imposed' Dirichlet boundary condition. We also note that the spiral breakup appears at lower values of than the instability of the Dirichlet source despite the fact that in both cases practically the same period and wavelength are selected. This is most probably due to non-linear effects not captured by the linear stability analysis. Near breakup it becomes increasingly difficult to prepare the initial condition and arbitrarily small steps in would be required to remain in the basin of attraction of the stable spiral solution [54]. For the excitable case (breakup close to the centre, b = 0.07) shown in the lower panel of figure 4, one gets very close to λ dir = λ min at the instability of the 1D source. While this is in line with naive expectations, the result is somewhat surprising at second glance. It basically suggests that the instability does not strongly depend on the boundary conditions, which is typical for an absolute instability. The Eckhaus instability, however, is expected to be convective for generic reasons. To resolve this paradox, a plain linear stability analysis is no longer sufficient-one needs to compute the propagation properties of unstable modes as well as the absolute spectra using exponential weights on the eigenfunctions.

Core breakup and finite wavenumber instability
So far, we have established the existence of two different scenarios of spiral breakup and showed that core breakup practically coincides with the convective instability of the planar wavetrains far from the spiral centre, whereas far-field breakup requires an absolute instability of the planar wavetrain similar to the scenario originally proposed for the CGLE [17]. Recent efforts have yielded a more detailed mathematical analysis of convective and absolute instabilities and their relation to boundary conditions [54]. Application of these results to the example described above indicate that the CGLE arguments hold indeed also for the case of far-field breakup in the modified Barkley model [45,53]. In the same work, the results for core breakup turned out to be less clear. One difficulty is that new defects once created in the far-field tend to be quite robust and are usually not pushed towards the boundary, even if the spiral is still linearly stable. This leads typically to deviations between the predictions from linear stability analysis and two-dimensional simulations already noticed in extensive studies of the CGLE [18]. As we have indicated above, a reduction to 1D sources allows for a better comparison between linear stability prediction and numerical simulations.
Sandstede and Scheel [53] show that in the case of an Eckhaus instability perturbations are always travelling away from the core (convective instability). Following their method, we find that the window between the Eckhaus and the corresponding absolute instability is tiny and below the resolution of the numerical results presented above in figure 4. Similar to the results in [53], the group velocity of the perturbation is also much smaller than in the case of the apparent 'convective', far-field breakup. Altogether such a convective Eckhaus instability leads to an apparent core breakup. Quantitative differences in the group velocity at the convective Eckhaus instability lead to a qualitative change in the breakup scenario for the 1D source shown in figure 2(d). The corresponding breakup in 2D spirals shown in figure 1(a) is, however, strongly affected by meandering [29], which is an instability of the spiral core and thus may naturally lead to core breakup.
In the remainder of this section, we will consider a different model equation and present an interesting alternative breakup mechanism: the absolute Eckhaus instability is preceded by an absolute instability to short-wavelength modes. These modes propagate inward and therefore cause a core breakup. We consider the following variant of the Barkley model: Below we use equation (3) with h(u) = u 3 supplemented again with Dirichlet (u(0) = u 1 ) and zero-flux (∂u(r)/∂r| L = 0) boundary conditions at respective ends of the 1D system of length L. As before, these boundary conditions are suitable to realize a 1D source of waves.
2D simulations of the model corresponding to equation (3) reveal a spiral breakup near the core influenced by a small-amplitude meandering [65]. The breakup phenomenology in 2D near the parameters investigated here is similar to the one shown in figure 1(a). At this stage, it is difficult to decide if this core breakup is mainly caused by meandering or by the unstable radial modes, which we will analyse in detail below. In another work, we have shown that meandering can induce also far-field breakup [69]. There, small-amplitude meandering excites the Eckhaus instabilities of the far-field. Large-amplitude meandering breaks the spiral near the core regardless of the far-field stability properties.
In the remainder, we consider only the 1D source case. For fixed a = 0.6, b = 0.06 we compute wavetrain solutions u(r − ct), v(r − ct) that depend on the choice of and their spatial period. The velocity c is uniquely determined by a choice of the former two parameters. The comparison with the numerical stability analysis of wavetrains will confirm the picture described above and in [45].
For small the wavetrains are linearly stable. As increases, either a band of longwavelength perturbations starting with zero imaginary part (marked blue in the following figures) or a band of finite-wavelength perturbations with non-zero imaginary parts (marked red) becomes unstable first; the second band follows upon further increase of . In our terminology, we assume that there is a relation between the imaginary part of the eigenvalue and the wavenumber of the corresponding eigenfunction. Visual inspection of the eigenfunction (not shown here) verifies this assertion. We have determined the corresponding critical values of for a range of wavelength of the wavetrains. Figure 5 shows the locations where the two instabilities occur in the essential spectra of the wavetrains (dashed blue and red line) in parameter space. The leading part of the (essential) eigenvalue spectrum is shown in figure 6 for selected solutions indicated by coloured dots in figure 5. The spectra in figure 6 have the same colour as the corresponding dots and figures 6(a) and (b) belong to wavetrains with periods 7.7 and 8.0, respectively. For periods smaller than 7.9 the finite wavelength instability (red) occurs first. Similar spectra (not shown here) are obtained for the parameters that yield core breakup in equation (1); therein, however, the Eckhaus instability appears always slightly before the finite wavenumber band and the two instabilities do not cross as in figure 5.
Near the onset of the instabilities we have computed the eigenvalue spectra of eigenmodes weighted by a factor exp(−µr). Sandstede and Scheel [53] pointed out that this method reveals the direction of transport of small perturbations. If eigenmodes transport towards positive r then a (sufficiently large) positive weight µ will damp their growth at positive r and stabilize them, i.e. shift their eigenvalues to smaller real parts. On the other hand, negative µ will damp the growth of modes that transport to negative r. One example of these weighted spectra is shown in figure 7. The thick curve corresponds to the essential spectrum (µ = 0) and thin curves are shifted due to weights of different magnitude. We find that long-wavelength modes (marked blue) similar to the Goldstone mode transport in the same direction as the wavetrain travels but finite wavelength modes (marked red) transport in the opposite direction. The corresponding numerical simulations of the one-dimensional system are shown in figure 8. The emitted waves are stable roughly as long as the periodic wave is linearly stable ( figure 8(a)). In figure 8(b), perturbations propagate inward with a group velocity ≈ − 0.4. It is difficult to conclude from these 1D simulations, how the corresponding spiral instability should appear in 2D, because the initial conditions are quite different in the 1D simulations shown, where we let the source boundary condition simply select the emitted wavetrain, and 2D spiral simulation which have to be started from a spiral computed at slightly different parameters. As mentioned above, 2D simulations for parameters similar to the ones in figure 8 show core breakup similar to figure 1(a). Different appearances of 1D and 2D simulations due to different initial conditions of the simulations are also seen for the far-field breakup case-compare figure 1(b) and figures 2(a) and (b).
The inward propagation of the instability suggests that it is related to the finite-k modes rather than to the Eckhaus modes. By using exponential weights, we can also measure the group velocity of the unstable modes from the linear stability analysis. We find that the group velocities amount to +3.0 (Eckhaus) and to −0.3 (finite k) for the fastest growing modes of the   respective bands for the parameters of the simulation. In addition, figure 6 clearly indicates that the finite-k modes have a much larger maximum growth rate (real part of the eigenvalue) than the Eckhaus modes in the region where both instabilities are already present. These observations hint at an important role of the finite-k modes for the instability observed in the numerical simulations which shall be proven below.
To complete the analysis, we turn to the discussion of the so-called absolute spectra. So far, we have computed the essential and weighted spectra of the wavetrains. From earlier work [17], it is however known that the relevant threshold of instabilities in finite domains are often given by the onset of an absolute instability of the relevant structures (here the wavetrains emitted from the 1D sources). Sandstede and Scheel [54] have rephrased this consideration by arguing that the so-called absolute spectra are the relevant indicators of instability in finite domains. In their view, an instability of the absolute spectra either indicates an absolute instability or more generally a so-called remnant instability. The term remnant instability indicates that unstable modes with positive and negative group velocities are present. This definition includes also the case of absolute instability, but it is also fulfilled by the presence of two different convective instabilities with unstable modes propagating in opposite directions. In our example, here we have indeed two different instabilities of a wave, which are convective (i.e. the modes at onset propagate with non-zero group velocities). This is a quite unusual situation and we are not aware of previous findings of this nature for the stability of periodic wavetrains. A remnant instability is realized for a wavetrain of given wavelength λ, if one is to the right of the two (red and blue) dashed lines. While such a remnant instability is a necessary condition for the The thin blue and red lines are again weighted eigenvalue spectra.
relevant instability of the absolute spectrum, it may not be sufficient. We can decide the reason for the instability only if we compute the so-called absolute spectrum following the definition given in [53,54]. For our example, it turned out that the absolute spectrum with largest real parts contains simply the exponentially weighted modes with zero-group velocity. These modes can be computed by finding the minima of Re(ω jn (µ)) when changing µ from µ = 0; see also the explanation at the end of section 2. The absolute spectra for two examples are displayed in figure 9. The absolute spectra consist of a line emerging from the Eckhaus band and a parabola related to the finite-k band. The absolute Eckhaus instability is found when the tip of the line enters into the right half-plane, whereas the absolute instability of the finite-k band appears when the parabola moves into the right half-plane. Upon increase of the control parameter the absolute instability of the finite-k band always appears clearly before the Eckhaus instability (see figure 9). The two absolute instabilities are shown in figure 5 as full lines (again red shows the finite-k absolute instability, whereas the absolute Eckhaus instability is indicated by the blue line).
In addition, we have performed careful numerical simulations with different values of u 1 in the boundary conditions, to compare the breakup in numerical simulations similar to the ones shown in figure 8 with these stability results. In figure 5, filled squares correspond to stable 1D spirals (as in figure 8(a)), whereas open squares indicate the breakup of 1D spirals (as in figure 8(b)). The comparison shows that clearly the absolute instability of the finite-k band (red line in figure 5) is responsible for the breakup observed in the numerical simulations, whereas the absolute Eckhaus instability responsible for breakup in 1D sources and spirals in previous investigations appears clearly at considerably larger values of and can thus not be the reason for the instability seen in the simulations. These results thus provide another example for the criterion originally derived for the CGLE [17], which requires an absolute instability beyond the convective Eckhaus instability of the emitted wavetrains. The requirement of an absolute instability holds also for the case of breakup in the 1D source from the finite wavenumber instability analysed here. The difference is that the modes connected with the finite wavenumber instability travel towards the core, whereas the modes responsible for the Eckhaus instability travel away from the core. The direction of propagation of the finite wavenumber modes in the essential spectrum or alternatively the negative exponential weight found in the computation of the absolute spectrum offer a plausible scenario for core breakup originating from the far-field.

Summary
In this paper, we have studied the phenomenon of spiral breakup by comparing numerical simulations with a stability analysis of periodic wavetrains. We used a modification of the Barkley model, which is very similar to the FitzHugh-Nagumo model. Therein, one observes two distinct phenomenologies of spiral breakup: core breakup and far-field breakup (see figure 1). The core breakup is accompanied by the meander instability, which makes the analysis more difficult, because it introduces a Doppler effect into the waves emitted from the spiral core. It is possible to suppress the meandering in simulations of a 1D source with a fixed 'core' boundary condition of Dirichlet type. Such a source is in some cases a good approximation of the radial dynamics of a non-meandering spiral. For the systems studied here, we found that the phenomena of core and far-field breakup persist in the 1D source and thus should not depend crucially on the presence or absence of meandering. Altogether, one can distinguish four distinct phenomenologies of spiral breakup: • Far-field breakup without meandering as seen in figure 1(b), the CGLE and the Belousov-Zhabotinsky reaction. • Core breakup without meandering as reported in the modified Barkley model by Sandstede and Scheel [53]. • Far-field breakup with meandering as seen in experiments by Zhou and Ouyang [68,59], in simulations of calcium waves by Falcke et al [51]. For a recent theoretical analysis, see the work of Brusch et al [69]. • Core breakup with meandering as seen in figure 1(a), in many cardiac models [24,28,35] and in experiments with the Belousov-Zhabotinsky reaction [59]. This route is often termed Doppler-induced instability.
The second main topic treated in this paper is the numerical stability analysis. We suggest the following strategy: • Identify spiral breakup in two-dimensional simulations of a reaction-diffusion model.
• Try to simulate the radial dynamics for a 1D source and suppress potential meandering influences on breakup. • Find all periodic wavetrain solutions from continuation and compute their spectra numerically. • In case of instability, determine the direction of propagation of unstable modes and/or compute the absolute spectrum to decide if the instability is convective or absolute. • Check if the selected wavetrain in the 2D spiral, respectively 1D source, is linearly stable, convectively or absolutely unstable and try to establish the relation between the stability results and the simulations.
In the above examples, we found that the instabilities of the wavetrain are always of convective nature and precede the breakup of waves near the source (see section 3.2). To obtain breakup in stationary sources or non-meandering spirals (in the absence of transversely unstable modes), an absolute instability of the plane waves in the far-field is required. This observation is related to the fact that the absolute spectra determine the stability in finite domains.
The linear stability analysis is useful to answer two main questions: (i) Is spiral breakup caused by a linear instability of a rotating spiral? (ii) Can the phenomenology of breakup be deduced from a particular type of instability and its properties like wavenumber of the critical modes and their direction of propagation? The evidence presented in this paper and in the recent literature allows us to answer the first question in a positive way. This is not at all trivial, because there are alternative possibilities. Other potential scenarios include the non-normal dynamics responsible for turbulence in simple flow geometries [70], extremely long-lived chaotic transients towards a final non-chaotic or even periodic attractor [71,72] or a coexistence of simple periodic solutions like spirals or rotating spirals and complex spiral-defect-chaos patterns. The latter possibility has been investigated in detail for chaotic convection patterns in Rayleigh-Benard convection experiments and simulations [73]. These cases have been ruled out for the reactiondiffusion systems studied here. We found that spiral breakup in simulations is clearly correlated with linear instabilities. The second question is more difficult to answer. Far-field breakup of simple spirals can in some cases be linked with the convective nature of the Eckhaus instability whereas, as shown in section 3.2, core breakup may stem from a finite wavenumber instability with critical modes that propagate towards the core. Earlier work, however, points to the fact that the core breakup in a 1D source reported in figure 2(d) is also related to an Eckhaus instability [45]. Sandstede and Scheel [53] noted first that the critical modes can also propagate outward in the case of core breakup, although their group velocity is in this case much smaller than for far-field breakup. In section 3.2, we applied the method of spatially weighted perturbations as suggested by Sandstede and Scheel [53] to determine the sign of the group velocity of the unstable modes and to compute the absolute spectra of the wavetrains. For our slightly modified model, we found a new finite wavenumber instability that becomes absolute for parameter values near the convective Eckhaus instability and thereby much earlier than the absolute Eckhaus instability. This instability has been identified to be responsible for the numerically observed core breakup in figure 8(b). Application of the same method to the case shown in figures 2(c), (d) and 4 showed that core breakup results also from an Eckhaus instability with small group velocity of the critical mode. For small group velocity, the parameter window between the convective Eckhaus and the absolute instability becomes extremely narrow and for coarser numerical resolution the difference between the two points may be overlooked. It is important to note that (i) core breakup consequently can result from different types of instability and that (ii) the same instability scenario (here: absolute version of Eckhaus instability with dominating modes propagating away from the core) may lead to phenomenologically different appearances like the different scenarios in figures 1 and 4. This statement is further illustrated by the many cardiac models, in which core breakup is presumably caused by a period doubling known as alternans [25,28,35]. For far-field breakup, one may also find explanations different from the absolute Eckhaus instability unveiled in the CGLE and the modified Barkley model studied above.
A further complication in the analysis of spiral breakup is spiral meandering. A Dirichlet boundary condition suppresses the meandering instability which could therefore be neglected in the above analysis. Simulations in 2D reveal this meandering instability and yield breakup near the core at slightly smaller than the critical values for the instability of the emitted wavetrain. Meandering leads to a Doppler effect in the waves emitted from the spiral core. Due to the modulation of the spiral tip curvature a motion of the source of waves results in a modulation of the wavelength and frequency of the emitted wavetrains. Originally, this effect has been analysed in simulations with core breakup [29]. Core breakup is expected also from the shape of the linear eigenmodes related to the meandering, which have been found to be largest near the spiral core [49] and may decay exponentially in large systems [67]. Meandering spirals can, however, also display far-field breakup as has been demonstrated also with simulations of the Atri model for intracellular calcium waves [51]. The influence of meandering and the resulting Doppler effect in core breakup has recently been verified in experiments with the Belousov-Zhabotinsky reaction [50]. Far-field breakup of (presumably) meandering spirals with a non-decaying modulation of the wavetrain in the far-field ('superspirals' [66,67]) has been studied in the same reaction [68,59]. Here, an interpretation of the breakup in terms of the non-linear consequences of the Eckhaus instability is possible within the CGLE [69]. Altogether, we admit that breakup with meandering is barely analysed in terms of linear stability analysis and may be in fact very difficult to settle. This is true in particular for the case of core breakup with meandering. If the meandering amplitude gets too large, waves may break already close to the centre and new spirals can appear. A stability analysis would however require a consideration of time-periodic waves, which is technically more difficult than the stability analysis of travelling waves with constant profiles. The technique has, however, been proven successful in other circumstances [74,75]. From the four main phenomenologies of breakup, we have discussed in detail examples for the two scenarios without meandering. In other words, we have supposed that spiral breakup appears as a primary instability of the continuous spectrum of the spiral that is related to the continuous spectrum of the planar waves [53]. Spiral breakup from meandering is instead a secondary instability following the primary meandering instability of the discrete spectrum of the spiral.
In retrospect, it is not surprising that linear stability analysis of spirals leads to a much richer variety of breakup scenarios than the simple phenomenology, because it only requires an absolute instability of the planar wavetrain far from the spiral core or a sufficiently strong instability of the spiral core. This leaves a large number of possibilities, whereas the phenomenological distinction of far-field and core breakup is only a crude scheme. The good news, however, is that the spiral's stability analysis may in some cases be reduced to the analysis of radial dynamics, which can be done with standard numerical procedures for most reactiondiffusion models. The main new contribution of this paper is the discovery of a finite wavenumber instability associated with inwardly propagating modes in contrast with the previously analysed Eckhaus instability, where the modes propagate outward from the spiral core. The absolute finite wavenumber instability may thus explain a core breakup without the presence of a core instability itself (similar to meandering). Computations show, however, that meandering is quite difficult to suppress in reaction-diffusion systems and a verification of our scenario in 2D requires further work. As the above examples show, a comparison of stability analysis and simulations leads to advances in the understanding of the mechanisms of spiral breakup and allows for a mathematically precise classification of such a transition way beyond the simple phenomenology of core and far-field breakup.