Computation of the Drift Velocity of Spiral Waves using Response Functions

Rotating spiral waves are a form of self-organization observed in spatially extended systems of physical, chemical, and biological nature. In the presence of a small perturbation, the spiral wave's centre of rotation and fiducial phase may change over time, i.e. the spiral wave drifts. In linear approximation, the velocity of the drift is proportional to the convolution of the perturbation with the spiral's Response Functions (RFs), which are the eigenfunctions of the adjoint linearized operator corresponding to the critical eigenvalues $\lambda = 0, \pm i\omega$. Here we demonstrate that the response functions give quantitatively accurate prediction of the drift velocities due to a variety of perturbations: a time dependent, periodic perturbation (inducing resonant drift); a rotational symmetry breaking perturbation (inducing electrophoretic drift); and a translational symmetry breaking perturbation (inhomogeneity induced drift) including drift due to a gradient, step-wise and localised inhomogeneity. We predict the drift velocities using the response functions in FitzHugh-Nagumo (FHN) and Barkley models, and compare them with the velocities obtained in direct numerical simulations. In all cases good quantitative agreement is demonstrated.


I. INTRODUCTION
Spiral waves are types of self-organization observed in physical [1][2][3], chemical [4,5], and biological [6][7][8][9][10][11] systems, where wave propagation is supported by a source of energy stored in the medium. The interest in the dynamics of spiral waves has significantly broadened in the last decade as the development of experimental techniques has permitted them to be observed and studied in an ever increasing number of diverse systems such as magnetic films [12], liquid crystals [13], nonlinear optics [14,15], novel chemical systems [16], and in subcellular [17], tissue [18] and population biology [19].
In the ideal unperturbed medium, the core of a spiral wave may be anywhere, depending on initial conditions. However, real systems are always subject to a perturbation. A typical result of a symmetry-breaking perturbation is drift of the spiral waves, which has two components, temporal drift, which is shift of spiral wave rotation frequency, and spatial drift, that is slow movement of the spiral's rotation centre. Drift of spiral waves, particularly the spatial drift, is of great practical interest to applications. In cardiac tissue, drift of re-entry circuits may be caused by internal tissue inhomogeneities, or by external perturbations, such as electrical stimulation. The possibility of control of arrhythmias by weak electrical stimulation has been a subject of intensive research for decades.
Understandably, the drift of spiral waves was mostly studied in the BZ reaction, which is the easiest excitable system for experimental study, and in the heart tissues and tissue cultures, which represents the most important application area. Examples of drift observed in experiments and numer-ical simualtions include "resonant" drift caused by (approximately) periodic modulation of medium properties through external forcing [20], constant uniform electric field that causes electrophoresis of charged ions taking part in the chemical reactions [21], a spatial gradient of medium properties [22][23][24][25] and pinning (anchoring, trapping) to a localized inhomogeneity [26][27][28]. Interaction with a localized inhomogeneity can be considered to be a particular case of the general phenomenon of vortex pinning to material defects, ranging from convective microvortex filaments in nanosecond lasermatter interaction to magnetic flux strings in the Sun's penumbra [29].
A most intriguing property of spiral waves is that despite being propagating waves affecting all accessible space, they, or rather their cores, behave like point-like objects.
Correspondingly, three-dimensional extensions of spiral waves, known as scroll waves, act as string-like objects. There have been several ad hoc theories of drift of spiral and scroll waves exploiting incidental features in selected models, e.g. [30][31][32][33]. Our present study is based on an asymptotic theory applicable to any reaction-diffusion system of equations in which a rigidly rotating spiral wave solutions exist. The theory was first proposed for autonomous dynamics of scroll waves for the case of small curvatures and small twists [34,35] and then extended to the drift of spiral waves in response to small perturbations [36]. In this theory, the particle-like behaviour of spirals and string-like behaviour of scrolls corresponds to an effective localization of so called response functions (RFs, see exact defininition later in this paper). The localization of RFs is the crucial assumption, which underpins the entire analysis. Originally [37] this property was only a conjecture based on the phenomenology of spiral waves in experiments and numerical simulations [31,[38][39][40]. The analytical calculation of the response functions appears to be infeasible. Numerical calculations in the Barkley model [41] and the Complex Ginzburg-Landau equation [42] have confirmed that indeed they are essentially localized in the vicinity of the core of the spiral. The asymptotic theory based on the response functions has been successfully used to quantitatively predict drift of spirals, for resonant drift and drift due to parametric inhomogeneity in the CGLE [43][44][45] and for drift in response to a uniform electric field in Barkley model [46]. Despite this success, so far the asymptotic theory has not become a generally used tool for the prediction of spiral wave drift. This is partly due to difficulties in the numerical calculation of the response functions. In our recent publication [47] we have presented an efficient numerical method of calculating response functions in an arbitrary model with differentiable right-hand sides. The complexity of calculating response functions with this method is similar to the complexity of calculating spiral wave solutions themselves. In the present paper, we describe the application of the asymptotic theory using the response functions for the prediction of several types of drift and show how it works for two of the most popular generic excitable models, the FitzHugh-Nagumo system [48][49][50], and the Barkley system [51]. We demonstrate that predictions of the asymptotic theory are in good quantitative agreement with direct numerical simulations. In addition, we demonstrate that the response functions are capable of predicting nontrivial qualitative phenomena, such as attachment of spiral waves to stepwise inhomogeneity and orbital movement around a localized inhomogeneity.
The structure of the paper is as follows. In Section II, we briefly recapitulate the asymptotic theory of the drift of spiral waves in response to small perturbation and present explicit expressions for drift parameters in terms of the spiral wave's response functions for several sorts of drift. In Section III, we describe the numerical methods used for calculating the response functions, for direct numerical simulations, and for processing of the results. The results are described in Section IV. We conclude the paper by discussion of the results and their implications in Section V.

A. General
We consider reaction-diffusion partial differential equations, where u( r, t) = (u 1 , . . . u ℓ ) T is a column-vector of the reagent concentrations, f (u) = (f 1 , . . . f ℓ ) T is a columnvector of the reaction rates, D is the matrix of diffusion coefficients, and r ∈ R 2 is the vector of coordinates on the plane.
A rigidly clockwise rotating spiral wave solution to the system (1) has the form where R = (X, Y ) T is the center of rotation, Φ is the initial rotation phase, and ρ( r − R), ϑ( r − R) are polar coordinates centered at R. For a steady, rigidly rotating spiral, R and Φ are constants. The system of reference co-rotating with the spiral's initial phase and angular velocity ω around the spiral's center of rotation is called the system of reference of the spiral. In this system of reference, the polar angle is given by θ = ϑ + ωt − Φ, with R = 0 and Φ = 0. In this frame, the spiral wave solution U(ρ, θ) does not depend on time and satisfies the equation where the unknowns are the field U(ρ, θ) and the scalar ω.
In a slightly perturbed problem where ǫh(u, r, t) is some small perturbation, spiral waves may drift, i.e. change rotational phase and/or center location. Then, the center of rotation and the initial phase are no longer constants but become functions of time, R = R(t) and Φ = Φ(t). In the co-rotating system of reference, time dependence will take form of a phase depending on time Thus, we consider three systems of reference: 1. laboratory, ( r, t) ; is the polar coordinate system centered at R; 3. co-rotating, (ρ, θ, φ), where θ = ϑ( r − R) + φ(t) is the polar angle, and φ = ωt − Φ(t) is the rotational phase, replacing time.
We shall look for a solution to (4) in the form of a slightly perturbed steady spiral wave solutioñ Then, assuming that˙ at leading order in ǫ, the solution perturbation g will satisfy the linearized system where whereh(U, ρ, θ, φ) is the perturbation h(u, r, t), considered in the co-rotating frame of reference. The linearized operator has critical (Re (λ) = 0) eigenvalues which correspond to eigenfunctions related to equivariance of (1) with respect to translations and rotations, i.e. "Goldstone modes" (GMs) In this paper we do not consider perturbations h(u, r, t) that depend on t other than 2π/ω-periodically (for a more general version of the theory free from this assumption see [36,45]). Thenh(U, ρ, θ, φ) is a 2π-periodic function in φ, and we look for solutions g(ρ, θ, φ) to equation (5) with the same periodicity. A solvability condition leads to the following system of equations for the drift velocities, where R = X + i Y is the complex coordinate of the instant spiral centre, the inner product · , · stands for the scalar product in functional space and the kernels W (n) (ρ, θ), n = 0, ±1, are the response functions, that is the critical eigenfunctions L + W (n) = µ n W (n) , µ n = −inω, n = 0, ±1, (9) of the adjoint operator L + , chosen to be biorthogonal to the Goldstone modes (8).
The drift velocities can be written as (henceforth we shall drop the O(ǫ 2 ) terms) where the "forces" F 0 and F 1 = (Re (F 1 ) , Im (F 1 )) T are defined by F n ( R, Φ) = W (n) (ρ, θ) , α n (ρ, θ; R, Φ) , and In the above formulae, the dependence on ( R, Φ) is explicitly included to emphasize that the response functions depend on coordinates (ρ, θ) in the corotating frame of reference whereas the perturbations are typically defined in the laboratory frame of reference, and the two systems of references are related via R and Φ.
Below we show how the forces (13), determining the velocity of the drifting spiral wave subject to a variety of perturbations, can be calculated using the computed response functions W (n) . We also compare the quantitative analytical prediction of drift velocities with the results of direct simulations.

B. Resonant Drift
Let us consider a spiral wave drifting due to the perturbation where A ∈ R ℓ is a constant vector. In the co-rotating frame the perturbation (15) will bẽ Substitution of (16) into (14) gives and, by (13), Hence the speed of the resonant drift of the spiral is whereas its direction is constant and arbitrary, as it is determined by the inial phase of the spiral Φ, or, rather, by the phase difference between the spiral and the perturbation, (19) is only valid in the asymptotic sense, and a more accurate formulation isΦ Hence, at finite ǫ the resonance is expected to be imprecise, and a typical trajectory of a resonantly drifting spiral is a circle of radius R rd = |Ṙ|/|Φ| = O(ǫ −1 ).

C. Electrophoretic Drift
Here we consider an anisotropic perturbation which breaks rotational symmetry where B ∈ R ℓ×ℓ is a constant matrix. This perturbation corresponds to action of an external electric field on a chemical reaction where some of the species are electrically charged. In this case matrix B is diagonal and its nonzero elements represent motilities of the ions of the reaction species. The same sort of perturbation appears in the asymptotic dynamics of scroll waves [34,35], where B = D.
In the co-rotating system of reference, the perturbation (21) can be written using the Goldstone modes (8), as which, by substituting into (14), gives (1) , which following (12) and (13) gives the velocity of the electrophoretic drifṫ which remains constant in time.

General
We now consider the case when the reaction kinetics f in (1) depend on a parameter p, and the value of this parameter varies slightly in space, Substitution of (25) into (1) gives, to the first order in ǫ, with the perturbation in the laboratory frame of reference Substitution of (26) into (14) gives where andp 1 (ρ, ϑ) is the parameter perturbation considered in the co-moving frame of reference. The final equations for the drift velocities can then be written in the forṁ where for brevity we introduce

Linear Gradient
Let p 1 vary linearly in a sufficiently large region containing the spiral tip and its subsequent drift trajectory. Specifically we consider p 1 = x − x 0 , where the x-coordinate of the trajectory remains near x 0 . In the co-moving reference frame, the linear gradient perturbation will bẽ Substituting (32) into (28) gives Then, by (27), (12) and (13), the velocity of the drift due to gradient of a model parameter will bė An important feature of equations (33) is that the first of them depends on X while the second does not. The dependence on X means that the drift velocity changes during the drift, unless the drift proceeds precisely along the y-axis. As it happens, at first order in ǫ, only the temporal drift, that is the correction to the frequency, shows this dependence. Namely, the first of equations (33) shows that the instant rotation frequency corresponds to the parameter value at the current centre of rotation, p = p 0 + ǫp 1 = p 0 + ǫ(X − x 0 ). The spatial drift, described by the second of equations (33), does not depend on X. That means that while the drift proceeds, its speed and direction remain the same, at least at the asymptotic order considered. This is an important observation, firstly, because it allows us to treat linear gradient induced drift in the same way as the electrophoretic drift, i.e. expecting drift along a straight line, and secondly, that unlike electrophoretic drift, the assumption is inherently limited to such X that ǫ(X − x 0 ) remains sufficiently small.

Step inhomogeneity
Here we consider a step perturbation located at x = x s , where H() denotes the Heaviside unit step function. In the co-moving frame of reference we havẽ Substitution of (34) into (28) gives We consider three intervals for xs−X ρ .
(2) ρ < |x s − X| , x s < X. Then H cos(ϑ) − xs−X ρ = 1, Thus, Substituting the above K n for the three intervals into (13) and (12), we get the velocities of the drift due to a step-wise inhomogeneity of a model parameter in the forṁ Note that bothṘ andΦ are functions of the current xcoordinate of the spiral with respect to the step, d = X − x s , andṘ is an even function of this coordinate.

Disk-shaped inhomogeneity
We now consider an inhomogeneity which is unity within a disc of radius R in centered at (x d , y d ), and which is zero outside the disc. Thus we havẽ Then calculations, similar to those for a stepwise inhomogeneity, lead to where l and ϑ 0 designate the distance and the direction from the current centre of the spiral to the centre of the inhomogeneity, i.e. x d = X + l cos ϑ 0 , y d = Y + l sin ϑ 0 . This leads to the equations for the drift velocities in the form It is straightforward to verify that if ϑ 0 = 0, x d = x s + R in and R in → ∞, that is when the disk is so large it turns into a half-plane at x > x s , then expressions (39) and (40) tend to expressions (35) and (36) respectively, as should be expected.
in accordance with the case of a pointwise, δ-function inhomogeneity considered in [52].

A. Models
We have considered two different kinetic models, both twocomponent, ℓ = 2, with one nonzero diffusion coefficient, for convenience. The kinetics FitzHugh-Nagumo system was chosen in Winfree [50] notation, with parameter values α = 0.3, β = 0.68, γ = 0.5 as in [47]. The Barkley [51] kinetics is given by with parameter values a = 0.7, b = 0.01 and c = 0.025, as in [46]. Note that both α and c are called ǫ in [50] and [51] respectively; however we use ǫ for the small parameter in the perturbation theory.

B. Response functions computations
For both the FitzHugh-Nagumo and the Barkley models, the response functions and the Goldstone modes have been computed using the methods described in [47]. The discretization is on disks of radii from ρ max = 12 up to ρ max = 50, using N θ = 64 of discretization intervals in the angular direction and a varying number N ρ of discretization intervals in the radial direction, up to N ρ = 1280. The components of the spiral wave solution and its response functions for Barkley model are shown in fig. 1(a). Similar pictures for the FitzHugh-Nagumo model can be found in [47].

C. Perturbations
We considered similar types of perturbations ǫh(u, r, t) in both FitzHugh-Nagumo and Barkley models, both for theoretical predictions based on response functions and in numerical simulations. Specifically, the perturbations were taken to have the following forms.
a. Resonant drift where ω is the angular velocity of the unperturbed spiral obtained as part of the spiral wave solution for the equation (3).
c. Spatial parametric inhomogeneities. As set out in Section II, a spatial dependence of a parameter p of the kinetic terms in the form p( r) For each of the two models, we consider inhomogeneities in all three parameters, namely p ∈ {α, β, γ} for the FitzHugh-Nagumo model, and p ∈ {a, b, c} for the Barkley model. The "linear gradient" inhomogeneity is of the form where x 0 is chosen to be in the middle of the computation box and close to the initial centre rotation of the spiral wave. The "stepwise" inhomogeneity is of the form where x s is varied and chosen with respect to the initial centre of rotation of the spiral wave. The − 1 2 term is added to make the perturbation symmetric (odd) about x = x s , to minimize the inhomogeneity impact on the spiral properties while near the step. As it can be easily seen, within the asymptotic theory, this term only affects the frequency of the spiral but not its spatial drift.
The "disk-shape" inhomogeneity is of the form where the position of the centre of the disk r d = (x d , y d ) T is varied and chosen with respect to the initial centre of rotation of the spiral wave.

D. Drift simulations
Simulations have been performed using forward Euler timestepping on uniform Cartesian grids on square domains with non-flux boundary conditions and five-point approximation of the Lapacian. The space discretization step ∆x has been varied between ∆x = 0.03 and ∆x = 0.1, and time discretiation step ∆t maintained as ∆t = 1 5 ∆x 2 . The tip of the spiral is defined as the intersections of isolines u(x, y) = u * and v(x, y) = v * , and the angle of ∇u at the tip with respect to x axis is taken as its orientation. We use (u * , v * ) = (0, 0) for the FitzHugh-Nagumo model and (u * , v * ) = (1/2, a/2 − b) for the Barkley model.

E. Processing the results
For coarse comparison, we use the trajectories of the instantanous rotation centre of the spiral wave. They are directly predicted by the theory. In simulations, they are calculated by averaging the position of the tip during full rotation periods, defined as the intervals when the orientation makes the full circle (−π, π], see fig. 1(b).
For finer comparison, we fit the raw tip trajectories, i.e. we use theoretical predictions including the rotation of the spiral. That is, if the theory predicts a trajectory of the centre as R = R(t; A, B, . . . ) ∈ C (a circle for resonant drift and a straight line for electrophoretic or linear gradient inhomogeneity drifts) depending on parameters A, B, . . . to be identified, then the trajectory of the tip is assumed in the form R tip (t) = R(t; A, B, . . . ) + R core e i(ωt+Θ0) where R core ∈ R is the tip rotation radius, ω ∈ R is the spiral rotation frequency and Θ 0 ∈ R is the initial phase. The parameters R core , ω and Θ 0 are added to the list A, B, . . . of the fitting parameters. arbitrarily selected parameter. For the resonant drift, the motion equations given by (18), (19) and (20) can be summarized, in terms of complex coordinate R = X + iY , as

A. Simple drifts
where p = 1 2 ǫ W (1) , A is predicted by the theory at leading order, and q = O(ǫ 2 ) is not, and we only know its expected asymptotic order. The theoretical trajectory is a circle of radius p/q, and the spiral drifts along it with the speed p. In the simulations, we determined both the radius and the speed by fitting. The speed is used for comparison and the radius is ignored.
For the other two types, electrophoretic drift and linear gradient inhomogeneity drift, the theory predicts drift at a straight line, according to (24) and (33) respectively. In these cases, we measure and compare the x and y components of the drift velocities separately.
For numerical comparison in the case of linear gradient inhomogeneity, we chose a pieces of trajectories not too far from x = x 0 , selected empirically to achieve a satisfactory quality of fitting.
A common feature of all graphs is that at small enough ǫ, there is a good agreement between theory and simulations. As expected, differences appears for larger ǫ with the disagreement occurring sooner (for smaller values of drift speed) for the linear gradient inhomogeneity drift. This is related to an extra factor specific to the inhomogeneity-induced drift: the properties of the medium where they matter, i.e. around the core of the spiral, changes as the spiral drifts. Since we require a certain number of full rotations of the spiral for fitting, faster drift meant longer displacement along the x axis and more significant change of the spiral properties along that way, which in turn affects the accuracy of the fitting. Fig. 3 illustrates numerical convergence of results with discretization parameters. We consider the simple drift cases and focus on forces, defined as the drift speed/velocity per unit perturbation amplitude ǫ. The discretization parameter that primarily dictates the accuracy of solutions is a spatial discretization step: ∆x in the simulations and the radial discretization step ∆ρ in the response functions calculations.

B. Numerical convergence
In simulations, the forces are determined for values of ǫ well within the linear range as determined in Fig. 2. These are calculated for different values of the space discretization step ∆x, where the time discretization changed simultaneously so that the ratio ∆t/(∆x) 2 remaines constant.
In theoretical predictions, the forces are given by the values of the corresponding integrals of response functions as described by Section II, and we have calculated the response functions and the corresponding integrals with various values of the radius discretization steps ∆ρ. Our discretization in both the theoretical and stimulations cases is second order in ∆x and in ∆ρ, so one would expect to see linear dependence of the drift forces on the squares of these discretization steps, (∆x) 2 and (∆ρ) 2 , at least for the values of these steps small enough. This is indeed what is observed.
We have gone further and extrapolated the calculated theoretical and simulation values of forces to zero ∆ρ and ∆x respectively, based on the expected numerical convergence properties. Such extrapolation gives the values of the forces which differ from the exact value only due to other, smaller discretization errors, which are: angular discretization and restriction to the finite domain in the theoretical predictions, and second-order corrections in ǫ and the boundary effects in the simulations. Comparison of such extrapolated data shows a very good agreeement between theory and direct numerical simulations (DNS) which is illustrated in Table I. Note that the values for fig. 3(e) and fig. 3(f) are also in good agreement with the results of [46].
For the extrapolation, we fitted the numerical data with the expected numerical convergence dependencies, which were different for theoretical calculations and for the simulations. In simulations, the central difference approximation of the Laplacian means that the next term after (∆x) 2 is (∆x) 4 . The expected error due to time derivative discretization is a power series in ∆t ∼ (∆x) 2 , hence the next term there after (∆x) 2 is again (∆x) 4 . The situation is different in the response functions calculations as there is no symmetry in the approximation of ρ-derivatives, therefore we expect that in the theoretical convergence, the next term after (∆ρ) 2 is (∆ρ) 3 . We note, however, that approximation of both theoretical and simulation data with similar dependencies, be it with a cubic or a quartic third term, gave very similar results.

C. Drift near stepwise inhomogeneity
The theoretical predictions for stepwise inhomogeneity, (45,47) and disk-shaped inhomogeneity considered next, are more complicated than the simple forms of drift considered up to this point. Because now the medium is inhomogeneous in the presence of the perturbation, the velocity depends on the instant position of the spiral center and as a result the spiral trajectories can be quite complex. Qualitative comparisons between theory and simulations can be made in the general case, but for detailed quantitative comparison we focus on the cases where the theory predicts simple attractors, e.g. a straightforward drift along the step for the stepwise inhomogeneity.
Equations (37) give a system of two first-order autonomous differential equations for X = Re (R) and Y = Im (R),       Graph Theory DNS Discrepancy fig. 3(a) 3.2795 + 1.8183∆ρ 2 + 1.8928∆ρ 3 3.2844 + 0.7166∆x 2 + 3.0663∆x 4 0.15% fig. 3 fig. 3(c)  where The right-hand sides of system (50) depend only on X but not on Y , that is, the system is symmetric with respect to translations along the Y axis. For this reason, the roots Note that the stability of invariant lines reverses with a change of sign of ǫ and also that F x (d) is an even function. Hence if ǫ = 0 and d * = 0 then either {x s + d * , Y } or {x s − d * , Y } will be an attracting invariant set. Fig. 4(a-f) show the theoretical predictions for the drift forces, i.e. velocity components per unit perturbation magnitude ǫ, F x (d) and F y (d), on the distance d = X − x s from the instant spiral centre to the step. This is done for both FitzHugh-Nagumo and Barkley models, for steps in each of the three parameters in these models. The roots of F x (d) are specially indicated. One can see from the given six examples, that existence of roots of F x () is quite a typical, albeit not a universal, event.
The qualitative predictions of the theory about a stable invariant line are illustrated by fig. 4(g) where we present results of numerical integration of the ODE system (50) and the results of direct numerical simulation of the full system. In the example shown, the positive root d * ≈ 2.644 of F x is stable and the negative root −d * is unstable. Hence the theoretical prediction for different initial conditions are: for X(0) > x s − d * and not too big, the spiral wave will approach the line X = x s +d * and drift vertically along it with the speed ǫF y (d * ) ≈ 0.8468ǫ; for X(0) < d * , the spiral wave will drift to the left with ever decreasing speed, until its drift is no longer detectable; for big |X(0)|, the drift will not be detectable from the outset.
As seen in fig. 4(g) this is indeed what is observed, both for the theoretical and for the DNS trajectories, and the visual similarity between theoretical and DNS trajectories is an illustration of the validity of the qualitative predictions of the theory.
Since the generic drift is non-stationary, a quantitative comparison for typical trajectories is difficult. However, the drift along the stable manifold X = x s + d * is stationary with vertical velocity given by ǫF y (d * ) so a comparison is easily made using the same methods as in the case of "simple" drifts considered in the previous subsections. The results are illustrated in fig. 4(h). As expected, we see good agreement between the theory and the DNS for small ǫ.
The phenomenological predictions are different for the case when F x (d) has no roots, or when its roots are so large that |F y (d * )| is so small that the drift cannot be detected in simulations. In such cases, the theoretical predictions for different initial conditions are: for |X(0)| not too large, the spiral wave will move with varying vertical velocity component but always in the same horizontal direction (to the right if ǫF x (0) > 0), eventually with ever decreasing speed, until its drift is no longer detectable; for |X(0)| too large, the drift will not be detectable from the outset. Note also that to get the drift velocities, F x and F y should be multiplied by ǫ which should be much smaller than c 0 = 0.025. So when the spiral is further than |X − x s | ∼ 2 from the step, the drift is very slow and hardly noticeable, even though according to the theory, there should be stable vertical drift around X = x s + d * , which is too slow to be observed in normal simulations.

D. Drift near disk-shape localized inhomogeneity
The theoretical predictions for the disk-shaped inhomogeneity (45,48), are more complicated but also more interesting. The theoretical spiral motion equation (41) has a rotational, rather than the translational symmetry of the stepwise        inhomogeneity. In polar coordinates (l, ϑ 0 ) centered at the center of the inhomogeneity, so that R = x d + i y d + le iϑ0 , equation (41) can be rewritten in the forṁ l = −ǫF r (l), where F r and F a are the radial and azimuthal components of the drift force, given by The minus sign in the first equation of (53) comes from the fact that in (41), the origin was placed at the instant rotation centre of the spiral, and the position of inhomogeneity is determined with respect to it, where as now we do the other way round: the origin is at the centre of inhomogeneity and the current position of the spiral rotation centre is determined with respect to it.
In the system (53), the axial symmetry is manifested by the fact that the right-hand sides of (53) depend on l but not on ϑ 0 , and the equation for l is a closed one. Hence roots l * of F r (l * ) = 0 represent invariant sets, which in this case are circular orbits. The movement along those orbits will have a linear speed ǫF a (l * ) and angular speed Ω = ǫF a (l * )/l * . The stability of these orbits is determined by the sign of ǫF ′ r (l * ): stable for positive and unstable for negative. Unlike the case of the stepwise inhomogeneity, now we do not have any mirror symmetries, as only positive l make sense, therefore for a given root l * a stable circular orbit is guaranteed only for one sign of ǫ but not the other. Fig. 4(a-f) show the theoretical predictions for the drift forces F r (l) and F a (l). This is done for both FitzHugh-Nagumo and Barkley models, for inhomogeneities in each of the three parameters in these models. The roots of F r (l) are specially indicated. One can see from the six given examples that existence of roots of F r () is rather common and often there is more than one root, l j , such that 0 = l 0 < l 1 < l 2 < . . . .
The qualitative predictions of the theory are illustrated in fig. 5(h,i) where we present results of numerical integration of the ODE system (53) and the results of direct numerical simulations of the full system. In the example shown in fig. 5(h), the predictions are given by fig. 5(e), which say that the smallest orbit has radius l 1 ≈ 3.724, with the or-bital speed F a (l 1 ) ≈ 0.003938, which is rather small compared to max(|F r (l)| ≈ 3.458 and max(|F a (l)| ≈ 8.534 and hardly observable in numerical simulations. Hence in this case, the radial component of the drift speed F r (l) is effectively constant-sign, and for negative ǫ, one should observe repulsion of the spiral wave from the inhomogeneity until it is sufficiently far from it, l ∼ 3, to stop feeling it, and for positive ǫ, the spiral wave will be attracted towards the centre of inhomogeneity from any initial position l ≃ 3. This is indeed what is observed in simulations shown in fig. 5(h) where the case of ǫ > 0 is shown, and the centre of the inhomogeneity, l = l 0 , is attracting for the spiral wave.
In the example shown in fig. 5(i), the inhomogeneity centre l 0 = 0 is repelling. Instead, the first orbit of radius l 1 ≈ 1.7722 is attracting. The perturbation amplitude ǫ = 0.3 in this case is quite large and comparable with the value γ 0 = 0.5 of the perturbed parameter itself. We see that although the numerical correspondence between theory and DNS in this case is not very good (note the distances between the open circles and between the filled circles), the qualitative prediction of orbital movement remains impeccable. As expected, the numerical correspondence becomes good for smaller values of ǫ, see fig. 5(g).

V. DISCUSSION
We have considered symmetry breaking perturbations of three different kinds: time-translation symmetry breaking, that is homogeneous in space and periodic in time ("resonant drift"); rotational symmetry breaking through differential advective terms ("electrophoretic drift"); and spatial translation symmetry breaking through space-dependent inhomogeneities ("inhomogeneity induced drift"). The latter type includes three sub-cases cases: a linear parametric gradient, a stepwise parameter between to half plane, and a parameter inhomogeneity localized within a disk.
a. Quantitative: drift velocity. We have demonstrated that asymptotic theory gives accurate predictions for spiral drift: in some cases the discrepancy between the theory and the direct simulations was as low as 0.02%. The discrepancy is affected by the numerical discretization parameters, both for the direct simulations and for response function computations, and by the magnitude of the perturbation.
b. Qualitative: attachment and orbiting. In the more complicated cases of spatial inhomogeneity, the response functions allow us to predict qualitatively different regimes of spiral motion, which we have been able to confirm by direct simulations.
In the presence of a stepwise inhomogeneity, the centre of spiral wave rotation may either be attracted to one side of the step where it gradually "freezes", or it may get attached to the step and drift along it with the constant velocity. In the latter case, the speed of the drift is proportional to the inhomogeneity strength, whereas the distance at which the attachment happens, does not not depend on the inhomogeneity strength at leading order. If the sign of inhomogeneity is inverted, the attachment occurs on the opposite side of the step and pro-ceeds in the opposite direction.
In disk-shape inhomogeneity, the situation is somewhat similar but more interesting. The spiral wave may be attracted towards the centre of the disk, or repelled from it. It may also be attracted to or repelled from one or more circular orbits. The drift velocity along the orbits is proportional to the strength of the inhomogeneity, whereas the radii of orbits do not depend on it at leading oder. The repulsion changes to attraction and vice versa, with the change of the sign of the inhomogeneity.
c. Prevalence of attachment and orbiting The possibilities of attachment to the step inhomogeneity and orbital movement for the disk-shape inhomogeneity are both related to the change of sign of the integrals of the translational response functions, which in turn are possible due to changes of sign of the components of those response functions. Not suprisingly, there is a certain correlation between these phenomena. The graphs fig. 5(a-f) may be viewed as deformed versions of the corresponding graphs fig. 4(a-f). Respectively, positive roots of F x (d) in fig. 4(a,d,e,f) have corresponding roots of F r (l) in fig. 5(a,d,e,f). However, the integrals in equations (51) and (54) are only similar but not identical, and the above correspondence between the roots is not absolute: the roots of F r (l) in fig. 5(b,c) and the smaller roots of this function in fig. 5(a,f) have no correspondences in fig. 4. Overall, based on results considered, orbital motion around a localized inhomogeneity seems to be more prevalent than attachment to a stepwise inhomogeneity. Moreover, the typical situation seems to be that there are multiple stationary orbits around a disk inhomogeneity. We have already discussed this situation in our recent preliminary short communication [52] where we have also illustrated how for the initial conditions between two stationary orbits, the spiral wave launched into one orbit or the other depending on the sign of the inhomogeneity.
The possibility of orbital drift, related to a change of sign of an equivalent to the function F r (l), has been discussed at a speculative level in [53]. The sign change of translational response functions was observed in oscillatory media described by CGLE [54,55]. The examples we consider here suggest that this theoretical possibility is in fact quite often realized in excitable media, and even multiple orbits are quite typical. Theoretical reasons for this prevalence are not clear at present. As stated in [52], the prevalence of multiple orbits may be understood in terms of asymptotic theories involving further small parameters. So, the version of kinematic theory of spiral waves suggested in [56] produces an equivalent of response functions, which is not only quickly decaying, but also periodically changing sign at large radii, with an asymptotic period equal to the quarter of the asymptotic wavelength of the spiral wave. Other variants of the kinematic theory, e.g. [57,58] did not, to our knowledge, reveal any such features on a theoretical level. However, numerical simulations of kinematic equations presented in [57] showed attachment of spiral waves to non-flux boundaries, which in a sense is similar to attachment to stepwise inhomogeneity. On a phenomenogical level, such attachment is, of course well known since the earliest simulations of excitable media, e.g. [38].
d. Orbiting drift vs other spiral wave dynamics. Properties of the orbital drift resemble properties of resonant drift when the stimulation frequency is not fixed as in the examples above, but is controlled by feed-back [59]. In that case, the dynamics of the spiral wave is controlled by a closed autonomous system of two differential equations for the instant centre of rotation of the spiral, like (50) or (53). In particular, depending on the detail of the feedback, this planar system may have limit cycle attractors, dubbed "resonant attractors" in [60], which may have circular shape if the system with the feedback has an axial symmetry. Apart from this being a completely different type of drift, we also comment that the second order ODE system is an approximation subject to the assumption that the feedback is instant, and in the situations when the delay in the feedback is significant due to the system size and large distance between spiral core and feedback electrode, the behaviour becomes more complicated.
For some combination of parameters, the trajectory of an orbiting spiral may also resemble meandering and may be taken for this in simulations or experiments. So, it is possible that orbital movement was actually observed by Zou et al. [61, p.802] where they reported spiral "meandering" around a "partially excitable defect"; although it is difficult to be certain as no details are given. The difference is that spiral meandering, in the proper sense, is due to internal instabilities of a spiral wave, whereas orbital motion is due to inhomogeneity. E.g. in orbiting, the "meandering pattern" determined by Ω/ω will change depending on the inhomogeneity strength.
The phenomenon of "pinning" of spiral waves to localized inhomogeneities has important practical implications for the problem of low voltage defibrillation [62][63][64][65]. In terms of spiral wave dynamics this is usually understood as attraction of the spiral centre towards the inhomogeneity locus. Practically interesting cases of pinning are usually associated with inexcitable obstacles, which are not small perturbations and therefore not amenable to the asymptotic theory considered in this study. However, the possibility of orbital motion around a weak inhomogeneity suggests that a similar phenomenon may be observed in strong inhomogeneities as well. This offers an unexpected aspect on the problem of pinning. Instead of a simplistic "binary" viewpoint, that a local inhomogeneity can either be attractive, which is the case of pinning, or repelling, which is the case of unpinning, there is actually a third possibility, which can in fact be more prevalent than the first two, namely, that at some initial conditions the spiral may orbit around one of a number of circular orbits, regardless of whether or not it is attracted to the center, which can be considered just as one of the orbits that happens to have radius zero. That is, there is more than one way that a spiral may be bound to inhomogeneities. e. Conclusion. We have demonstrated that the asymptotic theory of spiral wave drift in response to small perturbations, presented in [36,45], works well for excitable media, described by FitzHugh-Nagumo and Barkley kinetics models and gives accurate quantitative prediction of the drift for a wide selection of perturbations.
The key objects of the asymptotic theory are the response functions, i.e. the critical eigenfunctions of the adjoint lin-earized operator. The RFs have been found to be localized in most models where they have been calculated; however there are counterexamples demonstrating that surprises are possible [66]. Physical intuition tells that for the response function to be localized, the spiral wave should be indifferent to distant perturbations, which will be the case if the core of the spiral is a "source" rather than a "sink" in the sense of the flow of causality, for example as defined by the group velocity. Indeed, this localization property has been proven for onedimensional analogues of spiral waves [67] and there is hope that this result can be extended to spiral and scroll waves.
The effective spatial localization of the RFs on the mathematical level guarantees convergence of the integrals involved in asymptotic theory, and on the physical level explains why wave-like objects like spiral and scroll waves, while stretching to infinity and synchronizing the whole medium, behave respectively as particle-like and string-like localized objects. This macroscopic dissipative wave-particle duality of the spiral waves has been previously demonstrated for the Complex Ginzburg-Landau equation [45] which is the archetypical oscillatory media model. Here we confirmed it for the most popular excitable media models important for many applications.