A response function framework for the dynamics of meandering or large-core spiral waves and modulated traveling waves

In many oscillatory or excitable systems, dynamical patterns emerge which are stationary or periodic up in a moving frame of reference. Examples include traveling waves or spiral waves in chemical systems or cardiac tissue. We present a unified theoretical framework for the drift of such patterns under small external perturbations, in terms of overlap integrals between the perturbation and the adjoint critical eigenfunctions of the linearised operator (i.e. `response functions'). For spiral waves, the finite radius of the spiral tip trajectory as well as spiral wave meander are taken into account. Different coordinates systems can be chosen, depending on whether one wants to predict the motion of the spiral wave tip, the time-averaged tip path, or the center of the meander flower. The framework is applied to analyse the drift of a meandering spiral wave in a constant external field in different regimes.


I. INTRODUCTION
Spiral waves are remarkable self-sustained patterns which arise in various extended systems, including oscillating chemical reactions [1], catalytic oxidation [2] and biological systems. In biology, spiral waves have been observed across vastly different spatial scales, as they organise intracellular calcium waves [3], slime mould aggregation [4], waves of spreading depression in the brain [5], and cardiac contraction during arrhythmic events [6][7][8][9]. In the context of cardiac arrhythmias, different drift regimes are thought to correspond to different heart rhythm disorders [7]. Moreover, tracking of rotation centers of cardiac rotors was shown to improve ablation therapies of atrial arrhythmias [10].
The development of asymptotic description of drift of spiral waves has been contingent on two important issues: localization of adjoint symmetry modes [11][12][13][14][15][16] and an appropirate choice of the drift system of reference [17][18][19][20]. On this pathway, the choice of an "optimal" system of reference naturally depends on the unpertubed spiral wave dynamics, i.e. whether it is (i) a rigidly rotating spiral wave with its tip following a perfect circle, (ii) a bi-periodic "meander" regime when the spiral wave periodically changes its shape and its tip describes a biperiodic flower-like trajectory, or (iii) more complicated cases dubbed 'hypermeander'. While the asymptotic theory of the drift of rigidly rotating spiral waves is the most advanced [19,20], the theory of hypermeander drift is the one least developed [21][22][23]. The current work aims to lay out a technical framework for the drift of a bi-periodic meander regime, and illustrate its application on a simple example of electroforetic drift.
The propagation of cardiac excitation, as well as many other systems in which spiral waves are observed, can be described by reaction-diffusion (RD) systems (see, e.g. [24]), ∂ t u = ∂ j P jk (⃗ r, t)∂ k u + F(u, ⃗ r, t).
Here, u = u(⃗ r, t) ∈ R n is the list (column-vector) of n state variables, u = u 1 , . . . , u n T , e.g. the concentrations of reacting species, ⃗ r = (x j ) ∈ D ⊂ R N is the vector of Cartesian coordinates, N = 1, 2 or 3, and t is time. The vector-function F ∶ R n × R N +1 → R n describes the local dynamics of state variables, e.g. the reaction rates. The space-time dependent matrix-tensor P jk ∶ R N +1 → R n×n × R N ×N contains the diffusivities of the state variables on its diagonal (of which some may be vanishing). Non-diagonal terms of P are present in systems with cross-diffusion; see e.g. Ref. [25] and Fig.  1(b).
The asymptotic description follows the paradigm that first an idealized mathematical situation is considered, which is an unbounded, perfectly stationary, homogeneous and isotropic medium governed by the isotropic reaction-diffusion equation (RDE) and then any deviations from the idealized picture are considered perturbatively. Here the diffusion matrix P is constant and the kinetic functions F(u) are taken uniform in space and time. In the context of pattern formation, a bistable, oscillatory or excitable point system is commonly used. For well-chosen sets of reaction kinetics, i.e. functions F(u), the system (1) may sustain waves of finite amplitude, which propagate through the medium, see Fig. 1(a-b). Although our methodology works for any reactiondiffusion system allowing pattern formation, we will use below examples with the reaction kinetics of FitzHugh-Nagumo [14,26] P = diag(1, 0), Barkley [27] P = diag(1, 0), and Fenton and Karma. The equations and parameters of the latter model (with 3 state variables) are found in [28].
In N ≥ 2 spatial dimensions, the reaction-diffusion equation (1) may sustain rotating spiral-shaped solutions, as illustrated in Fig. 1(c-f). Spiral waves are commonly characterised by their tip trajectory, which can be found in different ways [28,29], e.g. as the locus where where j, k are indices of two selected state variables, j ≠ k, and u * j , u * k are two appropriately chosen constants. In N = 2 spatial dimensions, Eq. (2) means the intersection of isolines of two selected field components, where one isoline may correspond to a threshold level of the activator field, and the second condition separates the "front" and "back" parts of the first isoline based on the values of the inhibitor field. E.g. for Barkley kinetics, we use j = 1, u * j = 0.5, k = 2, u * 2 = 0.5a − b. It is often also convenient to take ∂ t u j = 0 as the second condition which corresponds to a point where the velocity of the activator isoline vanishes. Example tip trajectories are shown in Fig. 1(c-f). Another way of analysing spiral waves is by computing their phase, e.g. defined as a polar angle χ in a suitably chosen plane in the phase space of local kinetics: This mapping will produce a phase singularity (PS) in the vicinity of the spiral tip [30]. If j, k, u * j , u * k are set to the same values in Eqs. (2) and (3), the PS coincides with the spiral wave tip definition. Otherwise, different observers will generally not exactly agree on the tip or phase singularity position [31].
In three spatial dimensions, the solution to the reaction-diffusion equation (1) consisting of a stack of spiral waves is known as a 'scroll wave'. When tracking the scroll wave's tip or PS, Eqs. (2) and (3) produce a curve that is known as the 'scroll wave filament'. The same remark as in 2D holds here: different algorithms and chosen thresholds will yield different filament curves, which generally lie in each other's vicinity, i.e. in the tubular scroll wave core region.
However, if the tip trajectory is circular, all tip paths and PS trajectories will in 2D describe circles around a unique rotation center C at the same angular frequency ω. In previous works, the dynamics of circular-core spiral waves have been analysed in terms of the motion of the rotation center C [32][33][34][35]. For circular tip paths, the region inside it never excites and is sometimes called the 'spiral wave core'. The dynamics of circular-core spiral waves is reasonably well understood, and their drift response to small external stimuli can be found by projecting the stimulus on the inertial manifold of the system using so-called 'response functions' [17,19,36,37].
Interestingly, in several experiments and numerical simulations, both in chemical [38,39] and biological [9,28,[40][41][42][43][44] systems, the spiral wave was found to perform not a rigid rotation, but a more complicated motion, a phenomenon known as 'spiral wave meander' [45][46][47]. In excitable models, meander can arise due to the wave front interacting with the wave back from the previous excitation. In the simpler cases, the meander is quasiperiodic, in the sense that the evolution is periodic up to a Euclidean transformation. This transformation can always be written as a pure rotation or a pure translation (see below), which is why one can refer to it as 'biperiodic' meander, even when the tip trajectory is not a superposition of two circular motions. In more complicated cases, often called "hypermeander" [47], the tip motion can have more than two periods [22] or be chaotic, resulting in deterministic Brownian motion of the tip position [21]; we will not consider such cases here. A significant step in understanding meander of spiral waves was made by Barkley et al. [48] and Karma [49], who showed that the transition to meander from the circular-core case exhibits typical features of a Hopf bifurcation. Barkley [50] then showed that the process can be described by a set of five coupled ordinary differential equations. As a result, a tip trajectory is formed that is a superposition of two circular rotations, making a flowerlike tip trajectory. This picture was put on constructive footing in [51][52][53][54], using "skew-product" representation of the reaction-diffusion system exploiting its symmetry with respect to Euclidean motions of the plane. These works described the equivariant Hopf bifurcation (i.e. the regime close to the transition to meander), while many chemical and cardiac models show a qualitatively different tip trajectory, i.e. zig-zagged, star-like path known also as the "linear core" case. It turned out, however, that skew-product representation does not need to be restricted to the vicinity of the transition zone, and the more complicated meander pattern can be described as relative periodic solutions not necessarily related to a Hopf bifurcation [55,56].
Correspondingly, evolution of biperiodically meandering spirals, both flower-like and star-like, in response to small perturbations, can be analysed by linearising the system around a relative periodic orbit [16]. To show the steps of this process in detail is the purpose of this manuscript; it extends the procedure for rigidly rotating spirals found in e.g. [14,19,36]. Since periodically deforming waves in one spatial dimension are obtained as Qualitatively different one-and two-dimensional patterns that arise as solutions to the reaction-diffusion equation (1), depending on the chosel reaction model F(u) and initial conditions. a) Stationary traveling wave in the FitzHugh-Nagumo model [26] (α = 0.3, β = 0.68, γ = 0.5). b) Modulated 1D traveling wave in the FitzHugh-Nagumo model with cross-diffusion, shown at t = 3000 and t = 5000; reproduced from [25]. c) Rigidly rotating spiral wave in Barkley's model [27] (a = 0.52, b = 0.05, ς = 0.02). d) Spiral wave with Barkley kinetics exhibiting flower-like meander with inward petals (a = 0.58, b = 0.05, ς = 0.02). e) Meandering spiral wave with a linear-core in the Fenton-Karma Guinea Pig cardiac tissue model [28]. Labels indicate the order in which petals are visited. f) Spontaneous spiral wave drift in the case of resonant meander, in Barkley's model (a = 0.625, b = 0.05, ς = 0.02).
a special case, we also describe their dynamics.
A significant part of this paper will be dealing with introducing different coordinate systems that are suitable to capture the drift dynamics of solutions to the RDE (1) because, depending on which definition of filament one adopts one can obtain simple or complex laws of motion. These laws can be further simplified if one averages in time over rotation cycles of the spiral or scroll waves [36]. For non-stationary spiral or scroll waves, a possible averaging method is to analyse the motion in a Fourier series and only keep the non-oscillating, secularly growing terms [57,58]. Some recent works have computed the leading order eigenmodes in the laboratory frame of reference [15,59]. At the end, different descriptions need to be compatible, and it is our aim to illustrate that they all describe the same dynamics. This paper is organised as follows. In Sec. II we discuss the concept of symmetry breaking by the formed patterns, which will determine the number of degrees of freedom required to uniquely describe such patterns. We first illustrate this concept in Sec. III on traveling waves in one dimension, treating the cases of rigidly moving and periodically deforming waves. In Sec. IV we consider spiral wave patterns in two spatial dimensions, in the regimes of rigid rotation, meander and resonant meander. We consider in detail multiple possible frames of reference, which are defined with respect to either instantaneous or time-averaged dynamics. In Sec. V we linearise around the solution and discuss the critical eigenmodes of the system. This knowledge is used in Sec. VI to derive the equations of motion for modulated wave patterns, in a manner that is valid for both traveling and spiral waves. In Sec. VII we re-interpret our approach in terms of the geometry of the phase space of the reactiondiffusion system considerd as a dynamical system. This part is optional, but we found it useful to include as few works make the connection between the 'physics' style of description and the more abstract dynamical systems viewpoint.
In Sec. VIII the drift response of a spiral wave in a constant external vector field is analysed. We show how the dynamics can be averaged over rotation and temporal phases in order to find a simpler, manageable equation of motion. In particular, we demonstrate that if a meandering spiral wave is phase-locked to a constant external field of magnitude E, its drift velocity can be larger than order E: it is close to the 'orbital velocity' which which the spiral wave circumscribes the meander flower.
We conclude with a discussion of our present work in Sec. IX. To assist the reader in remembering the notations used, we provide a list of notations in the Appendix.

II. SYMMETRY CONSIDERATIONS
It is customary to consider Eq. (1) in a bounded spatial domain with Neumann boundary conditions, which agrees with many physical situation where such models are relevant, and this is used in all numerical simulations of this work. However, to develop the theory of symmetry-breaking, we here consider the whole space R N , with N = 1 for wave fronts and N = 2 for spiral waves. The two sets of boundary conditions do not need to conflict, as numerical evidence in various reactiondiffusion patterns [11][12][13][14][15]58] shows that the sensitivity of patterns to external stimuli, including boundary conditions, is strongly localised. The mathematical condition for this is that the adjoint critical eigenfunctions ('response functions' [19]) W M defined in Eqs. (17) below decay exponentially far from the wave front (in 1 spatial dimension) or far from the spiral wave tip (in 2 spatial dimensions). Therefore, faraway boundaries are not felt at e.g. the wave front or spiral tip.
Since the system (1) involves the spatial variables only via the Laplacian ∆ in the whole space R N , it is equivariant (symmetric) with respect to the isometric trans- is also a solution, for any isometry γ. For our purposes it is sufficient to consider only the orientation-preserving transformations, that is, γ ∈ SE(N ). Besides, being autonomous, this system is equivariant with respect to the group of translations of the time axis, SE(1) ∼ R, that is if u(⃗ r, t) is a solution, thenũ(⃗ r, t) = u(⃗ r, t − T ) is also a solution for any T ∈ R. (Note that this is different from saying the solution is 'invariant' with respect to time translation, which would imply a solution that is constant in time.) Thus, we can say overall that (1) is symmetric with respect to the inhomogeneous Galilean group G = SE(N ) × R. The dimensionality of this Lie group is dim(G) = N (N + 1) 2 + 1.
Any given solution u(⃗ r, t) of Eq. (1) may be invariant with respect to none, or all, or part of the symmetries of the equation itself. Invariance of the solution with respect to a purely spatial symmetry γ ∈ SE(N ) can in general be reduced to a lower-dimensional solution and is not further considered. E.g. a plane wave in N = 2, which is invariant for translations along the wave front. This case can be reduced to the study of a traveling wave in N = 1 spatial dimension. An interesting case arises when a solution is invariant with respect to a spatiotemporal symmetry. E.g. for a rigidly translating traveling wave in N = 1 at constant speed c, u(⃗ r + cT , t + T ) = u(⃗ r, t) for any T , which implies that there are members of G that leave the solution invariant.
Let the group of symmetries leaving a solution invariant be H, H ≤ G. If this is also a Lie group, then application of all possible transformations γ ∈ G to u(⃗ r, t) will produce a manifold of solutions, of dimensionality N BS = dim G − dim H. We refer to this dimensionality as the number of broken symmetries (implying that the "unbroken ones" are in the subgroup H). By construction, the tangent vectors to this manifold of solutions are solutions of the (1) linearised on u(⃗ r, t) which do not exponentially grow or decay with time, so they make the center subspace of this solution. If the solution is an attractor to the system, this manifold of solutions will be called an 'inertial manifold', see details in Sec. VII.
In the following text, it will sometimes be convenient to indicate for every quantity, e.g. wave front or tip coordinates, the reference frame in which it is defined, using different math accents. We shall call the 'simplest' frame that makes the solution constant or periodic the center frame of reference. For any quantity f , its value pertaining to this frame of reference, will be denoted as ○ f . The frame which best follows the wave front or spiral wave tip will be called the 'tip frame', with the corresponding notationf . An list of all notations and symbols used throughout this manuscript is given in the appendix.

A. General framework
Let us first consider a generic choice for the frame of reference to describe a traveling wave U(x, t). We denote the origin of the lab frame with Cartesian coordinate x as O. In the lab frame, we can define the wave front position x = Z(t) by thresholding a state variable u j (Z(t), t) = u * j and demanding that ∂ t u j (Z(t), t) > 0 to distinguish it from the wave back. This is one possible definition for the wave front position as a moving point Z(t), with lab frame coordinate x = Z(t). However, alternative definitions are possible, e.g. for the wave shown in Fig. 1(b) the barycentre of u 2 2 was used to denote wave position.
Next, we introduce a (yet unspecified) moving frame of reference with origin X that has lab frame coordinate X(t) and denote the new spatial coordinate in the moving frame as ρ: The moving frame can be chosen in different manners; an overview is given in Tab. I.
The instantaneous wave front position is denoted Z(t), it is found at the position ρ = ζ(t) in the moving frame:

B. Symmetry breaking
In a one-dimensional homogenous spatial domain (N = 1), the RDE (1) is invariant with respect to translation in time and space, i.e. G ∼ R × R, and dim(G) = 2. A traveling-wave solution does not have translational symmetry, neither in time nor space, so we can say it breaks the translational symmetry. However, if the wave travels with constant speed c f without changing its shape, then we can say that u(x − X, t − T ) = u(x, t) as long as X and T are related by X = c f T . That means, that the solution is invariant with respect to a one-parametric group of transformations, so dim H = 1, and N BS = 1. This is the case of rigidly translating waves, to be treated in subsection III C. However, if there is no frame of reference moving with constant speed in which the solution is constant in time, then the subgroup H is trivial or at most discrete, so dim H = 0, and N BS = 2. The case where this dynamics is periodic, so H consists of translations The non-periodic case, when H is trivial, falls outside our present scope.

C. Rigidly translating waves
Firstly, if the wave is rigidly moving, Z(t) = Z(0) + c f t with c f the constant traveling wave speed. In this case, we simply take X ≡ Z and X ≡ Z and speak of the 'comoving frame'. In this frame, the traveling wave is a solution to the RDE (1) of the form: Note that Eq. (5a) describes the type of solution and (5b) gives the equation of motion for the collective coordinate X(t) (wave front position) that is introduced since the pattern breaks the translational symmetry along x. If a perturbation of order η is added to the system, Eqs. (5a) and (5b) will gain additional terms of O(η) in their right-hand sides, describing the wave profile change and wave speed correction, respectively.

D. Periodically modulated waves
Traveling wave may also possess a variable shape and variable propagation speed c f , in which case This situation has been encountered both in experiment and simulations [25,60], see Fig. 1b for an example. We only treat periodically modulated waves here, where c f (t) has temporal period T . In the presence of perturbations, it will be convenient to rescale the period to 2π, define the positive constant Ω = 2π T and call Ψ = Ωt the temporal phase of the solution, which labels how far the wave front shape has gone through its cycle. That is, Ψ is essentially a scaled time variable, and Ω is the scaling constant.
Throughout this work, we will use the bar notation to denote the average over one temporal cycle: for any 2π-periodic f (Ψ). We now have multiple options to choose the frame of reference.
The traveling wave solution is now periodic in a moving frame, i.e.
Due to the scaling of time, u * is 2π-periodic in Ψ.

Fully co-moving frame
In the 'front frame', we take again X ≡ Z and X ≡ Z, at the expense of a non-constant propagation speed c f (Ψ).

Constant-speed co-moving frame
Alternatively, one may choose to let the frame move at the time-averaged velocity Since c is constant here, this frame is simpler than the fully co-moving frame. Keep in mind, however, that X now describes the average progression of the wave front position; the true wave front position Z can be found from

A. General framework
Let x a , a ∈ {1, 2} be Cartesian coordinates in the lab frame with origin O. In the lab frame, the spiral wave tip Z as defined by (2) describes a path x a = Z a (t) that may be circular or flower-like, see Fig. 1.
To capture spiral wave dynamics, we will introduce a moving frame with origin X (t), whose lab frame coordinates are x a = X a (t). At all times, the orientation of the new frame of reference with respect to the lab frame is given by the rotation phase Φ(t), with rotation rate ∂ t Φ = ω that may be time-dependent. So, ω takes positive values for counterclockwise frame rotation and negative values for clockwise frame rotation. If ω is constant, Remark that for meandering spirals, ω coincides with the spiral rotation rate only if the 'tip frame' (defined below) is chosen. For the other frame choices in Tab. II, ω equals the (typically slow) precession rate of the meander pattern.
The spatial coordinates in the moving frame are denoted ρ A , and co-moving time coordinate is τ . The generic coordinate transformation between the frames is thus given by: The rotation matrices connecting lab frame coordinates x a to rotating coordinates ρ A are given by Since the generator of rotations is the Levi-Civita symbol A B ( 1 2 = − 2 1 = 1 and zero otherwise), it commutes with rotation matrices, and . We will write the polar angle in the chosen frame of reference as θ, while the spatial orientation of the reference frame will be the rotation phase Φ. Due to rotational invariance of the RDE (1), the solution will generally depend only on This relation is the angular equivalent of Eq. (4).
In the moving frame, the spiral wave tip position is found at position ρ A = ζ A (τ ). Therefore, the lab frame tip position Z a is given by: FIG. 2. Different coordinate systems suitable to analyse and circular-core spiral wave. a) Center frame: the origin is taken at the center of the circular tip trajectory. b) Tip frame: the origin is taken at the current position of the spiral wave tip, and the ρ 2 axis is tangent to the tip trajectory.

B. Symmetry breaking
The Euclidean plane R 2 is invariant under the Euclidean group SE(2) of translations in x and y and rotations. However, the spiral wave solutions are not invariant under SE(2), since shifting or rotating a solution will yield a solution to the RDE (1) that is distinct from the original one. So, at least 3 symmetries are broken by the spiral wave solution.
As in the case of 1D pulses, one may introduce a moving frame of reference, which in this case is also rotating; see below for different options. If in that frame, the spiral wave is stationary, one has N BS = 3 and the spiral wave is rigidly rotating; see Par. IV C. Else, the time symmetry is also broken and N BS = 4. Here, we consider only the case where the spiral wave solution in the co-moving frame of reference is periodic in time. Returning to the lab frame of reference, this case represents a meandering spiral wave, either in the non-resonant (see Par. IV D) or resonant case (see Par. IV E).

C. Rigidly rotating spiral waves
When going to a frame of reference rotating at ω = Kω s , (i.e. the spiral wave frequency), the solution becomes periodic: with c A and ω constant. We distinguish two frames of reference here:

Center frame
The simplest coordinate system is found when taking its origin in the middle of the circular tip trajectory, as Spiral type N BSS N BT S Frame name Origin of moving frame Choice Rotation Translation rigidly rotating 3 0 center frame middle of tip traj.
minimally rotating co-moving frame on a line Overview of reference frames used in the text to describe the motion of rotating spiral wave patterns. N BSS : number of broken continuous spatial symmetries, N BT S : number of broken continuous temporal symmetries, O: lab frame origin, C: middle of the tip trajectory, Z: spiral wave tip, X : chosen origin of the moving frame ("filament"), c a : velocity of X in the lab frame, c A : its velocity in the rotating frame, K = ±1 is the sign of the (period-averaged) frame rotation rate ω, introduced so that ωs > 0 or Ω > 0, respectively. N BSS + N BT S = N BS , the total number of broken symmetries.
shown in Fig. 2(a), so that This frame was used in the vast majority of previous works on spiral wave theory [17, 18, 32-35, 57, 61, 62]. This approach, however, is not optimal for describing the interaction of the spiral wave with localised heterogeneities, as dynamics is most influenced by stimuli close to the spiral wave tip rather than rotation centre, and the distance between the two can be significant.

Tip frame
Therefore, we introduce another frame of reference, with the origin of coordinates ρ A at the spiral wave tip, X = Z as shown in Fig. 2b: In this 'tip frame', the velocity of the spiral wave tip c A becomes constant. One can even choose the orientation of the frame such that c 1 = 0, c 2 = ωr with r the radius of the tip trajectory.
The tip frame closely relates to the kinematic approach of Zykov et al. [63] and is particularly suited to describe large-core spiral waves.
Finally, note that the freedom in the choice of the definition of the tip may be exploited by taking u * j = u j (C), u * k = u k (C), so that Z = C; then we have c A = 0 and recover the center frame as a special case.

D. Meandering spiral waves
The case of biperiodic meander studied in this work is characterized by the fact that the spiral wave solution can be made time-periodic in a moving frame of reference. As in the case of periodically modulated waves, we introduce a temporal phase Ψ, which in the absence of perturbations increases from 0 to 2π over the time interval T : Ψ(t) = Ωt, Ω = 2π T > 0. Without loss of generality, we may take Ψ(0) = 0. In the context of spiral waves, we refer to Ψ as the meander phase, or temporal phase, as opposed to the rotation phase (angle) Φ.
A meandering spiral solution is thus of the form: with u * 2π-periodic in Ψ. An overview of reference frames discussed below is given in Tab. II.
The Euclidean transformation mapping u(x, y, t) to u(x, y, t + T ) can, without loss of generality, be always thought of as either a pure rotation or a pure translation, since the composition of a rotation and a translation can always be written as a single rotation in 2D. Below we first treat the 'pure rotation' case, corresponding to nonresonant meander; the second case (pure translation), corresponding to resonant meander will be discussed in Sec. IV E.

Definition of the angle β between petals
In the case of non-resonant meander, the spiral pattern repeats itself after time T , but rotated around an angle β around a point C that will become the center of the spiral tip trajectory ('meander flower'); see Fig. 3.
Since we only impose that the solution be periodic, β can be chosen up to an integer multiple of 2π. Two values can play a special role, however. First, we can FIG. 3. Different coordinate systems suitable to analyse and predict the motion of a meandering spiral wave. Panels a-d show cycloidal meander in Barkley's model with parameters as in Fig. 1d; panels e-f illustrate linear-core meander in the FK-Guinea Pig cardiac model. Panels (a,e) show the 'center frame', where the origin is taken at the center of meander flower, with the frame rotating at the constant rate ω = β T , where β is the angle between consecutively visited petals. Panels (b,c) illustrate frames tailored describe the finite-core, which may either minimally rotate (over an angle β0 (see b) or co-rotate with the spiral, i.e. rotate over βs = β0 + 2πK where K = ±1 labels chirality. Panels (d,f) show a tip frame of reference, with the origin X of the coordinate system at the current tip position and the ρ 2 -axis tangent to the tip trajectory. The rotation rate of the tip frame is ω(Ψ), non-constant but periodic.
define this angle by following the movement of the spiral. If we attach a reference direction to the spiral wave tip (e.g. ⃗ ∇u 1 evaluated at the tip), this vector will turn over the angle β s during one spiral period, and for a faraway observer, the average spiral wave frequency (expressed as a positive number) will be ω s = Kβ s T , where K = ±1 denotes chirality: K = 1 for counterclockwise and K = −1 for clockwise rotation of both the spiral wave and the frame of reference.
Secondly, we may choose the β as the element of the set β s + 2πZ that has the minimal absolute value, and call this minimal value β 0 . For the cases studied here, we find

Center frame
We pick a frame rotating at constant frequency ω = β 0 T around C, i.e. we take X = C. In this 'center frame of reference', the meandering spiral solution becomes periodic. The solution is given by (8), with This set of equations describes the spatial position the spiral wave in terms of the position of the center of the meander flower, and was used in [16]. Since we chose X = C, Eq. (8b) expresses that the center of the meander flower does not move in the absence of perturbations to the RDE (1).
It is also possible to relax the convention to take the period T as the minimal time interval in which the solution is periodic modulo Euclidean group actions. E.g. when linear cores rotate over β ≈ 180 ○ [28], one could define the period to span two such cycles, bringing β ≈ 0 ○ instead. This formalism can then be used to study phase-locking to constant external fields, generalising the results in [16].
Another theoretical possibility is to use the same setting as described by system (8), but with the choice β = β s instead of β 0 . This choice would ensure that the spiral wave solution, considered in the rotating frame, oscillates but not rotates and remains approximately stationary far from the centre. Neither of these two theoret-ical alternatives is used later in this paper, so we mention them only for completeness.

Minimally rotating finite-core frame
It can be advantageous to allow the origin of the coordinate system to better follow the spiral wave tip. One possible choice is shown in Fig. 3b for reaction kinetics with cycloidal meander. We let the new origin move at constant angular velocity on a circular trajectory that approximates the exact tip trajectory, and refer to this situation as the 'minimally rotating finite-core frame':

Co-rotating finite-core frame
We can let the minimally rotating finite core frame rotate around its axis, at a rate KΩ, to better follow the spiral wave rotation, bringing the absolute rotation rate to β s ; see Fig. 3c. We will refer to this case as the 'corotating finite-core frame'. It is again described by Eqs.
In this frame, the spiral solution u * (ρ A , Ψ) will not rotate if Ψ is increased from 0 to 2π.

Tip frame
Lastly, we let the origin X of the co-moving frame of reference coincide with the spiral wave tip position, i.e. X = Z, X a (t) ≡ Z a (t); see Fig. 3c and f. Then c A and ω become non-constant, 2π-periodic functions of Ψ.
After one meander period of duration T , the frame will have turned over an angle If the tip trajectory is smooth, we can without loss of generality suppose that c 1 = 0, c 2 = c(Ψ). This can always be achieved by exploiting the freedom in the definition of frame orientation angle Φ.

Comparison with the classical theory of meander
In the classical theory of rigidly rotating spiral waves, the critical eigenmodes for translation have eigenvalues ±iω 1 , where ω 1 is the rotation frequency of the rigidly rotating spiral wave. If due to a parameter change in the medium, another complex eigenvalue pair crosses the imaginary axis at ±iω 2 , we have a Hopf bifurcation in the quotient system (rotating frame). If this bifurcation is supercritical, this results in an epicycloidal tip trajectory as shown in Fig. 3(a-c). In Barkley's notation [50], it is seen that beyond the Hopf bifurcation the absolute tip velocity will oscillate at ω 2 . Hence, we conclude that Barkley's ω 2 equals Ω used throughout this work. Secondly, in the circular-core regime just before the Hopf bifurcation, the spiral's rotation frequency is ω 1 , corresponding to ω s in our notation. We conclude that in the tip frame (where β = β s ), ω = ω 1 and in the minimally rotating frames, (where β = β 0 ), ω = ω 1 − KΩ.
The case of resonant meander happens in the classical theory when ω 2 → ω 1 . In our description, this happens when ω → 0 in the center or minimally co-rotating frame, or when ω → Ω in the tip frame or fully co-rotating frame.

E. Resonant meander
We now consider the case where the meandering spiral returns after the temporal period T to the same state, but translated over a vector d a = c a T ; see Fig. 1(f). Different authors call this sort of solutions cycloidal motion [46], 0 ○ isogon contours [47], modulated travelling waves [50,64], resonant linear motion [53] or linear drift [23]. We shall refer to it here as resonant meander.
Note that during dynamics under external perturbations, the direction of spontaneous drift may change, such that a rotation phase needs to be introduced nevertheless. This makes this case distinct from the 1D modulated traveling wave from Par. III D.
The solution is still given by (8), with c A and ω depending on the chosen reference frame: 1. Minimally rotating co-moving frame First, we let the origin of our coordinate system move on a straight line with velocity c 1 , c 2 , without frame rotation, parallel to the line along which the resonant spiral travels over time. This situation is sketched in Fig. 4a. Then, we find: Since the frame is not rotating, the function u * (ρ A , Ψ) shows a spiral that rotates a full turn when Ψ is raised from 0 to 2π.

Co-rotating, co-moving frame
We can also let the previous frame rotate at ω = KΩ, while its origin moves on a straight line at constant speed parallel to the displacement direction of the resonant spiral. We again obtain Eqs. (8), with FIG. 4. Reference frames for resonant meander. a) The 'minimally rotating co-moving frame' is rigidly translating such that the solution becomes periodic. During one period (i.e. between the frames shown), the frame rotates over an angle 2πK. b) The 'co-rotating co-moving frame' is obtained from the average frame by letting the moving coordinates rotate at constant frequency KΩ = 2πK T . c) Tip-based frame of reference, with the origin of the coordinate system at the current tip position and the ρ 2 -axis tangent to the tip trajectory. Rotation rate of the frame is ω(t), non-constant but periodic.
and c a constant. The spiral wave profile u * (ρ A , 0, Ψ) does not rotate a full turn when Ψ is raised from 0 to 2π; it only deforms ('breathes'), see e.g. [56].

Tip frame
As in the non-resonant case, it is also possible to find a tip-based frame of reference: Note that during the time interval T , the angle Φ needs to turn over ±2π, whence ω = KΩ.
Here, one can see that resonant meander is indeed a special case of Eqs. (8), in which the average value of ω tends to KΩ. In the tip frame, u * represents a spiral that performs no net rotation when Ψ is raised from 0 to 2π, it merely deforms around the same orientation.

F. General case
The foregoing patterns can all be captured by noting that each broken continuous symmetry delivers a collective coordinate X M : We shall denote this as x M → X M ; there is one pair for every broken continuous symmetry. We also convene that The use of indices M vs. m thus distinguishes between the moving and lab frames of reference, as did A vs. a for the translational modes only. Note that X a is tip position in the lab frame, while x a is a lab frame Cartesian coordinate that can label any point.
The parameters of the solution x M are chosen from {x, y, θ, t} and X M are the collective coordinates. Then we find that all aforementioned sets of equations (5), (6), (7), and (8) that capture the evolution of particular reaction-diffusion patterns can be written as given a time-dependent spatial coordinate transformation that has as many degrees of freedom as broken Galilean symmetries by the solution.
Note that the representation (10) is not unique. E.g. for a rigidly rotating spiral wave, one could choose U(x a , t) = u * (ρ A , Ψ) that is periodic in its third argument Ψ ≡ ϑ and ρ A identical to the lab frame coordinates x a . However, one can also choose to let u * in the right-hand side of Eq. (10) only depend on ρ A and Ψ.
For the actual calculations performed below, we work with Eqs. (8) for a meandering spiral in a tip-based frame of reference, since they contain all other patterns and reference frames as a special case, i.e. Eqs. (5), (6), (7)and (8). Those can be recovered from Eq. (8) by considering u * independent of a collective coordinate. Thereafter, the evolution equation for that coordinate becomes irrelevant and can be dropped. E.g. the theory of rigidly rotating spirals follows from the meander case by taking u * independent of Ψ and leaving out the evolution equation for Ψ. From this case, the 1D traveling wave case follows by considering u * to be also independent of Φ and ρ 2 , such that Eq. (5) follows.

V. PROPERTIES OF THE LINEARIZED PROBLEM A. Right-hand zero modes
Since the original RD Eq. (1) is invariant under spatial, rotational and temporal shifts, one has that, when a solution is rotated, or shifted in space or time, it is still a solution to Eq. (1). This can be made explicit by substituting u(x +ε, y, t) in Eq. (1), and similar for other coordinates, yielding in first order inε that As expected, shifted values of the solution are right-hand zero modes of the linearized operatorˆ associated to the RDE, in the lab frame.
In several previous works, a moving frame was chosen in which the solution was either stationary [18,32,34,61,65] or periodic [16]. This comes down to expressing the derivative in Eq. (11) in the moving frame using the chain rule. From Eqs. (8) then follows Here, θ is the polar angle around the origin of the chosen moving reference frame. Thus, u * obeys Here, ω and c A may depend on Ψ (which is the case in the tip frame of a meandering spiral wave). The linear operator associated to Eq. (14) iŝ Since the original RDE (1) is invariant in space, one has that in an infinite medium, if U(x a , t) is a solution, so should U(x a +ε a , t) for any constant displacementε a . We find in one spatial dimension that L∂ ρ u * = 0.
In two spatial dimensions, we remark that ∂ a U = ∂ A u * R A a (Φ), such that we find in first order inε a : By taking complex combinations V ± = −(1 2)(∂ x u * ± i∂ y u * ), one obtains true eigenmodes ofL, asLV ± = ±iωV ± [14]. Note that in some of the chosen reference frames, ω may depend on Ψ, and the 'eigenmodes' of the system are 2π-periodic functions of Ψ.
For 2D patterns, we can state that a solution that is rotated around an angleε around the current origin (i.e. spiral tip position or meander centre position) will still be a solution. With ∂ θ = A B ρ A ∂ B , one finds that, if U(x a , t) +ε∂ θ U(x a , t) is a solution, then L∂ θ u * = 0.
We will denote −∂ θ u * as V Φ and −∂ A u * as V A . Thirdly, expressing that a time-shifted solution U(x a , t +ε) also solves (1) yieldŝ where ∂ t u * , also denoted V Ψ , is given by Eq. (13). In the case of meandering spirals, V Φ is linearly independent of V Ψ , since otherwise a shift in time would be equivalent to simple rotation, and we would find ourselves in the non-meandering case.
To summarise, equation (12) defines the set of N BS zero modes for the linearized operatorˆ in the laboratory frame; in the comoving frame, according to (15), this produces a set of N BS eigenvalueŝ In quantum field theories, bosons appearing due to spontaneous breakdown of continuous symmetries are called Goldstone bosons; extending the analogy to the classical nonlinear field, the eigenfunctions corresponding to breakdown of continuous modes are sometimes referred to as 'Goldstone modes'.

B. Adjoint problem
Let us associate toL the operator Note that u * is 2π-periodic in Ψ, and in the space of functions 2π-periodic in Ψ,L † is the adjoint operator tô L, in the sense that where N = 1 for 1D waves and N = 2 for 2D wave patterns. For non-deforming solutions, we note that u * is independent of Ψ, so the inner products in Eqs. (16a) and (16b) coincide, for any f and g defined by u * . Given thatL † is the adjoint toL, we assume that it also has N BS critical eigenmodes W M that are 2πperiodic in Ψ: These 'adjoint critical modes' are known as sensitivity functions or response functions (RFs) [19], as will be explained in Sec. VI.

C. Instant orthogonality of left and right critical modes
The set of critical adjoint modes can be normalised as Moreover, the orthogonality of critical and adjoint critical eigenmodes holds instantaneously [16,59,66]: Here, we show where the proof in [16] needs to be adapted to accommodate for non-constant rotation rates ω(Ψ) and the case of resonance ( ω → Ω).

Let us suppose thatL
We note thatL † is adjoint toL with respect to the inner product ⟨⋅ ⋅⟩, which follows from integration by parts and the tempered nature of the adjoint eigenmodes W M . Denoting ⟨W M V N ⟩ as I M N , it follows that Hence Now, if we have λ M = λ N , then I M N is constant. Then, setting ⟪W M V N ⟫ = δ M N already yields (18). Else, if R(λ M ) ≠ R(λ N ), and I M N (0) ≠ 0, then I M N (Ψ) grows or decays exponentially in time, which cannot happen since it should be 2π-periodic. Therefore, I M N (Ψ) = 0 in this case. Finally, it is possible that R(λ M ) = R(λ N ) but I(λ M ) ≠ I(λ N ), e.g. when considering the inner product between critical eigenmodes. In that case, λ M − λ N = iω(Ψ) with ∈ {0, ±1, ±2}. Then, If not in the resonant case, β mod π ≠ 0, and the exponential factor in Eq. (21) cannot be periodic, whence I M N = 0. In the case of resonance, we are free to choose a non-rotating frame, where ω = 0, such that (18) still holds.
◻ In the laboratory frame of reference, the critical modes are true zero modes, whence all λ m = 0, and the preservation of the inner product immediately follows from Eq. (20).

A. Derivation of the drift equations
Using the ingredients defined above, it is possible to predict how a stable reaction-diffusion pattern (e.g. a plane wave or a spiral wave) reacts to a small perturbation. We mainly follow the derivation for rigidly rotating spiral waves [36] but extend it to the case of meander. In comparison to [16], we offer more flexibility in the frame of reference, such that also meandering spirals close to resonance can be treated.
We start from the perturbed RDE: where h = O(η) is a small perturbation. A more generic form of perturbation can include dependence on the solution itself, i.e. h(x, y, t, u, ∇) which however is reduced to (23) when the perturbation is evaluated at the unperturbed solution, i.e. h(x, y, t, u * (x, y, t), ∇).
If the initial state is close to a stable solution (i.e. travelling or spiral wave), the net effect of h will be to cause a spatiotemporal drift of that pattern, which can be inferred from the collective coordinates. As before, we will present the result from the most general case of a meandering spiral wave (Eqs. (8)), from which all other cases can be inferred by eliminating some of the collective coordinates.
Thus, we approximate the solution as whereũ = O(η). The coordinate transformation is now given by in which the collective coordinates' temporal evolution is perturbed by yet unknown drift terms v M , which are also O(η): We can make the decomposition (24) unique by imposing that This is possible by shifting the solution. E.g. if one approximates at a given instance of time the solution u * (ρ 1 , ρ 2 , Ψ) as u * (ρ 1 +ε, ρ 2 , Ψ) +ε∂ 1 u * (ρ 1 , ρ 2 , Ψ), then ⟨W 1 ũ⟩ =ε ≠ 0, and the proposed solution can be better shifted to match the true solution, i.e. until Eq. (25) holds. A second interpretation of Eq. (25) is that the deviation vector u should be orthogonal to the inertial manifold [36]. Next, plugging the ansatz (24) into the perturbed RDE (23) delivers: Note that by including the minus sign in the definition of the GMs, i.e.
a plus sign appears in front of v M V M . Projection onto the response function W M delivers, due to the meander lemma (18): The first term vanishes, due to the gauge condition (25) and the fact thatL is adjoint toL with respect to ⟨⋅ ⋅⟩: Hence, we find the simple result: It turns out that at all times, the drift induced by a small perturbation can be found by taking the overlap integral between the perturbation and the response function W M corresponding to that degree of freedom. This result justifies to call the critical adjoint eigenmodes 'response functions' [19]. In engineering terms, they are also the spatiotemporal 'impulse response' to a localised stimulus. This property can also be used to estimate RFs numerically or in future experiments [58].
In the above-discussed moving frames of reference, the equations of motion become which can also be summarised as [16]: Note the presence of a zeroth order motion c A (Ψ), which accounts for spiral tip motion even in the absence of perturbations. Moreover, the rotation rate c ϑ ≡ ω is allowed to depend on the meander phase Ψ.
To find the net drift effects, motion generally needs to be time-averaged in the laboratory frame of reference, starting from: We will provide some elementary examples in Sec. VIII.
At this stage, we can understand why different frames of reference may be used depending on the context. Let us suppose that we describe a process in which the perturbation varies only slightly across the essential support of the RFs of the system, say h(⃗ r) = E a (⃗ r)Q∂ a u, where Q ∈ R n×n is a constant matrix acting in the space of reactive components (state variables). Then, the overlap integrals will typically be approximated as [17,18,32,35] and for non-homogeneous E N the next-to-leading order will be small only if the response functions W M are well localised near the point ⃗ X. Although the extent of the RF is fixed by the parameters of the model, the observer can thus describe the reduced system more accurately by an appropriate frame and tip choice. E.g. the interaction of a three-dimensional scroll wave in a detailed cardiac geometry, the results using the tip frame are expected to be more accurate than the center frame description. A similar argument holds for wave fronts, whose critical adjoint eigenmodes are also localised around the wave front.

VII. GEOMETRIC INTERPRETATION: DYNAMICS ON THE INERTIAL MANIFOLD
The methods and results presented above can be represented geometrically in the language of dynamical systems on manifolds [55,67]. The formal application of the idealized scheme presented below, as used in this paper as well as, explicitly or implicitly, in some other previous studies, has well known technical difficulties related to the facts that our phase space is a functional space, i.e. is infinite-dimensional, and that SE(2) is a non-compact group. Some of the implications of these are discussed in the concluding section; for now, we proceed ignoring these technical difficulties, in order to describe interesting physical phenomena, leaving rigorous mathematical treatment for subsequent studies. The first part of this section is an elaboration on the short discussion given in Sec. II above.

A. Phase space orbits
Every possible state of the system (1) at a fixed time, e.g. U 0 (⃗ r), can be represented as a point of an infinitedimensional phase space V. The set U(⃗ r, t) obtained by evolution from the initial condition U(⃗ r, t 0 ) = U 0 (⃗ r) according to (1) is the orbit of U 0 . For the formal development of the theory, we take t ∈ R, but we will not use backward time evolution, which will be ill-posed for the systems considered. Of the examples shown in Fig. 1, only rigidly rotating spiral waves and meandering spiral waves with rational 2π β have closed orbits.

B. Spatial vs. spatiotemporal symmetries
Since time is treated as an parameter in phase space dynamics (i.e. it is not a 'direction' of the phase space), it will be useful to distinguish between different symmetry groups that include or not include time symmetry. As in Sec. II, we take G = SE(N ) × R and H the largest continuous subgroup of G which leaves a given orbit invariant. The largest continuous subgroup of SE(N ) leaving the given orbit invariant is denoted J. Then, we have The number of broken symmetries N BS , broken space symmetries N BSS and broken time symmetries N BT S by the solution are then given by:

C. Quotient space
The group orbit of a state U(⃗ r) under a group B is given by It will be useful to study the dynamics of the system, disregarding spatial Euclidean symmetries. Therefore, we consider equivalence classes: U 1 ∼ U 2 iff there exists γ ∈ SE(N ) such that U 1 = γU 2 . The space of equivalence classes is called the quotient space, denoted Q = V SE(N ) and contains as members the group orbits under SE(N ). The original dynamical system (1) induces a (reduced) dynamics on the states in Q.
The equilibrium points of the dynamics in Q are termed relative equilibria in V, as they correspond to equivalence classes under the action of the Galilean group in V. Hence, from the examples discussed above, the rigidly translating waves and rigidly rotating spirals are relative equilibria of the reaction-diffusion system. (For comparison, the resting state is a 'true' equilibrium).
The limit cycles of the dynamics in Q are termed relative periodic orbits in V. Of the examples above, modulated traveling waves in one spatial dimension, biperiodically meandering spiral (both flower-like and star-like) and spirals in resonant meander are relative periodic orbits.
We will below further describe solution to the RD system that are relative equilibria or relative periodic orbits. These have an attractor A ⊂ Q of dimension 0 or 1. The attractor is not necessarily unique to the medium. For instance, for the case of a 2D medium supporting spiral waves, one already finds 4 different attractors: the resting state, a traveling plane wave and spiral waves of two opposite chiralities.

D. Inertial manifold
The presence of stable persisting patterns in the system puts a special structure on the phase space. Consider the set of states in V that belong to the same relative equilibrium or relative periodic orbit: This set M (implicitly depending on a given choice for A) forms a manifold in the phase space V. From our previous assumption thatˆ has no eigenvalues with positive real part, it follows that nearby trajectories are attracted with exponential speed to M, whence it is an inertial manifold [68] [69]. Since unidirectional attraction is implied here, the emergent structure of an inertial manifold is typical for dissipative systems.
A given system (1) may have various inertial manifolds corresponding to different prototypical solutions U, such as the resting state, a propagating pulse, or a spiral wave with either chirality.
Inertial manifolds constitute a special case of slow manifolds [67]. When the Euclidean symmetry is not broken (h = 0), the dynamics on the inertial manifold is trivial and given by the orbit of that solution. In the presence of a perturbation, however, the solution will not follow the original orbits anymore, but slowly drift from one original orbit to the other. This regime has been the starting point for many analytical studies of non-linear wave dynamics, e.g. [35,65,[70][71][72][73], including extension to 3D scroll waves [16-18, 32, 33, 74].
A remarkable property is that in the high (infinite)dimensional phase space, the inertial manifold associated to the patterns above is finite-dimensional. Given that the orbit U(⃗ r, t) has dimension 1 and dim(SE(N )) = N (N + 1) 2, the maximal dimension of M would be N (N + 1) 2 + 1. However, a state U(⃗ r, t 0 ) may be invariant under a subgroup J of SE(N ), in which case J is called an isotropy subgroup. E.g. a plane wave in 2D is invariant under translations along its wavefront.
The elements of J correspond to the preserved symmetries of the solution, such that the number of broken spatial symmetries amounts to: For solutions that are relative equilibria, time evolution can be represented by a Euclidean group action: U(⃗ r, t) = γU(⃗ r, 0) with γ = γ(t) ∈ G for all t ∈ R. In this case, the time symmetry is not explicitly broken, and we say that N BT S = 0. For relative periodic orbits, U(⃗ r, t) = γ(t)U(⃗ r, 0) with γ(t) ∈ G only when t is an integer multiple of the period T of the limit cycle in Q.
In this case, we say that the number of broken time symmetries N BT S = 1.
From the above, it follows that the dimension of the inertial manifold equals the number of broken continuous symmetries, i.e.

E. Collective coordinates parameterize the inertial manifold
By definition, the inertial manifold is flow-invariant, so the dynamics on it can be described by a system of N BS coupled ordinary differential equations. Since it is a local attractor, this system will approximate the dynamics within the vicinity of that manifold.
The different frames shown in this work serve the same purpose: to uniquely characterize a dynamical state of the system. Different frame choices thus correspond to different parameterizations of the N BS -dimensional inertial manifold.
Since every spatial collective coordinate is related to a broken symmetry, one has that where ρ M is a co-moving frame coordinate and the prototypical solution U(x, t) is assumed to depend on the collective coordinates X M as parameters. E.g. since a traveling wave has the shape U(x, t) = u * (x − X(t)), one has Similarly, for spiral waves where θ is the polar angle in the chosen frame of reference and ϑ is the spiral's rotation phase.

F. Tangent spaces
One can introduce tangent vectors to the inertial manifold by differentiating the state with respect to a collective coordinate: Here, γ −1 (X) has the effect of working in the quotient space (or fully co-moving frame). If one differentiates with respect to the time-like collective coordinate (Ψ), one obtains a tangent vector to the orbit of the state under time evolution according to the RDE (1).
These tangent vectors to the inertial manifold thus correspond to the Goldstone modes (GM) of the theory, or right-hand critical modes of the system, also called V M in sections V-VI. The non-critical modes of the system then point towards the other directions in phase-space, off the inertial manifold.
In every point of the inertial manifold, one can also introduce a set of dual vectors W M , that correspond to the response functions defined and used above. In every point of the inertial manifold it is possible to choose

G. Effect of a small perturbation
Now, suppose one starts with a pattern state that exactly equals the prototypical solution U to the RDE. Hence the initial dynamics of this state is on the inertial manifold. Then, applying a small stimulus h at time t = 0 has the effect of shifting the state in phase space. Unless h is carefully tuned, the state after applying the perturbation does not lie on the inertial manifold anymore. However, we assume that the wave solution U is dynamically stable, such that the resulting state evolves towards the inertial manifold and converges to an orbit in it.
Then, the resulting orbit will generally not be the same state as would have been reached without the perturbation. That is, one retrieves the same solution (spiral or traveling wave), but with different collective coordinates X M . A practical way of finding this shift, is to measure which part of the perturbation was tangential to the inertial manifold, and that is precisely what Eqs. (26) mean. The part of the perturbation that is orthogonal to the inertial manifold will decay over time, and only represent transient dynamics (as long as the amplitude of the perturbation is small enough).

H. Coordinate change in the tangent manifold
Geometrically speaking, the different coordinate systems which we presented in Sections III and IV are alternative ways to define collective coordinates, and therefore alternative bases vectors for the tangent manifold. With that respect, we can state that such that where J is the Jacobian of the coordinate transformation and J H is its Hermitian transpose. Eq. (42) is an example of the adjoint representation, and will be exemplified by Eqs. (50), (55).

A. Motivation
A simple perturbation field that can be considered is a constant vector field that couples to the gradient in one or more state variables: In chemical systems, this term can model an applied electrical field ⃗ E that induces drift of different ions u j with respective mobilities M j that are found on the diagonal of the matrix M. The resulting phenomenon in a chemical context is known as 'electroforetic drift'.
In 3D systems supporting spiral waves, the motion induced by a small filament curvature κ also generates a term of the form (43), with ⃗ E = κ ⃗ N and M = P, where ⃗ N is the unit normal vector to the filament curve. For circular core spirals, the law of motion (27) with the stimulus (43) yields with The relations (44) hold both in the tip frame (with c A constant) and the center frame (with c A = 0). However, the RFs and Goldstone modes are different in both cases, such that in the end both descriptions are equivalent, as we now show.
Recasting the second equation for t(Φ) instead of Φ(t), we immediately get since the average of a rotation matrix over its angle vanishes. Next, dividing the first equation of (44) by the second, we obtain Integrating over Φ on [0, 2π] delivers the period-averaged drift velocity: with In the center frame, with c A = 0, this is the classical result from [18]. However, does the calculation in the tip frame yield the same result? Here, we can make the discussion of Par. VII H explicit. The center and tip frames are for a circular-core Call the polar angle in either frame ○ θ andθ; the unperturbed spiral solution is To find the RF transformation law, we consider a localized perturbation of the j-th state variable at a given position: and in the tip frame After substituting Eqs. (51) and (52) into Eq. (53), we conclude from the zeroth order term in η that and from the first order term in η that

⎛ ⎝W
Hence, using Eq. (54) we can verify thať and similar for γ 2 . This result resolves an apparent paradox: the expression for the filament tension γ 1 is different in the center frame and the tip frame, as shown by the first and second line of Eq. (56). However, both expressions are equal since the RFs and GMs appearing in these expression differ for both frames. The final result is thus independent of the chosen frame of reference.
C. Phase-locking of a meandering spiral under a constant external field.
In systems where spiral waves are found where the angle β between subsequent petals of the meander flower is close to a simple fraction of 2π, the rotation phase may lock to the external field if that is strong enough [16].
To capture this phenomenon quantitatively, we take the minimally rotating finite-core frame of Fig. 3b. Then, when E = ⃗ E = 0, the X b (t) describes the motion of a point moving on a circle with radius equal to the mean radius of the meander flower. To recover the precise tip motion from X b (t) one can take Plugging Eq. (43) into Eq. (27) then delivers: The matrix element M M B is formally defined as in Eq. (45), but now M ∈ {x, y, Φ, Ψ}. Note that functions M M B (Ψ) and R a A (Φ) are all 2π-periodic. Due to our choice of coordinate system, c A and ω are non-zero but constant.
In our analysis, we first treat the case where β ≪ 2π, i.e. near the 1 ∶ 1 resonance. Then, we assume that both ω and E = ⃗ E are O(η), such that a locked solution may exist. Under those conditions, Ψ(t) will remain a monotonous function, and we can trade the time coordinate t for Ψ. Expanding up to linear order in η delivers: If there is phase-locking of the rotation angle Φ to the meandering phase Ψ, then equation (59b) will give a solution for Φ that remains in a bounded vicinity of some mean value Φ * , We will be looking for the simplest case of 1:1 phase locking, such that From Eq. (59b), the rate-of-change of Φ is also small (i.e. O(η)), while the rate-of-change of Ψ is Ω = O(1). Therefore, we can assume thatΦ = O(η) during one meander cycle, and average over the cycle. Denoting Φ(2πn) as Φ n , we obtain the difference equation: Then, a phase-locked angle is a fixpoint of Eq. (60). Without loss of generality, we can take ⃗ E = E⃗ e x here, to show that this condition can be formulated as: with A necessary condition for a phase-locked solution to exist is then EA ω > 1, i.e. [16] From Eq. (61), the rotation angle around which the locked equation will equilibrate is The second solution to Eq. (61), namely α − arccos (ω AE), is unstable.
Since in the phase-locked state, Φ = Φ * + O(η), Eq. (59a) can be readily integrated to find the net spatial displacement during a phase-locked meander cycle: where T = 2π Ω. From ∂ t Ψ = Ω + O(η), one finds moreover that the time needed to make a full meander cycle is T + O(η), such that the net drift speed during phase locking is In words, it can be concluded that the drift speed of a meandering spiral phase-locked to a constant vector field ⃗ E is of the order of the 'orbital' velocity c, i.e. the net speed at which the tip traverses the rim of the meander flower. Importantly, within the locking region the drift velocity is not proportional to the perturbation strength E! In terms of the velocity parallel or perpendicular to the applied field, Eq. (63) reads: a result which was already quoted in [16]. The magnitude of drift velocity is therefore expected to be V = c + O(η). Fig. 5 confirms that this is the case for simulations in Barkley's model. Until here, we have described phase-locking of a single meander cycle to an external field, which happens when ω ≪ Ω. Since Ω = 2π T and ω = β T with β the angle between subsequently visited petal tips, phase-locking is found generally when In some cases, a higher-order resonance may occur, i.e.
In this case, q meander cycles can be re-interpreted as a single cycle to which the theory (27)- (63) can be applied, with ω replaced by qω − pΩ.

D. Regime close to phase-locking
It could happen that the applied field ⃗ E is not strong enough to enforce phase-locking, which occurs when E < E crit . In that case, the net drift can be found from the instant equations of motion (27) by averaging of one sort or another. In the case of circular core spirals, this was done by sliding time averaging [36] or by imposing condition of periodicity on the perturbed solution that in turn requiring a solvability condition [20]. Since for meandering spirals the unperturbed solution is not periodic but biperiodic, the situation here requires a bit more care. In [58] we considered the case below but not too close to phase-locking threshold, and used Fourier series. For the case E ≲ E crit , this does not work, so we resort to the more robust method, by sliding averages. First we do averaging over the meander phase: Then, if ω ≪ Ω, As above, we can write the second term in the right-hand side as EA cos(Φ − α). Let us denote the time needed to complete a full meander flower as Under a perturbation this value changes to T . Then, straightforward integration of Eq. (64) delivers: Therefore, if E approaches E crit from below, the time needed to complete a full meander flower diverges. Beyond E crit the spiral fails to complete a meander cycle and is phase-locked.
E. Average drift speed of a non-phase-locked meandering spiral We continue our analysis of Eq. (27) in the regime where phase-locking does not occur, i.e. when β 2π is For any of the reference frames shown in Fig. 3, we can in a first step integrate Eqs. (58) to find the net displacement in space, time and rotation angle during the single meander cycle, say the n-th. At the start of this cycle, taking place at t = t n the collective coordinates of the spiral are denoted X a n , Φ n , Ψ n . For t ∈ [t n , t n+1 ], we can write: X a (t) = X a n + X a * (t) +X a (t).
In these right-hand sides, the first term is the initial value at the start of the n-th cycle, the second term is the evolution in the unperturbed case, and the third term the correction to it, of O(E). Without loss of generality, we can set the n-th cycle to start at Ψ = Ψ n = 2πn.
The unperturbed evolution is given by Due to the group structure of rotations, we can write We will below use Note that these functions are not periodic in Ψ anymore.
With Eqs. (65)-(66), the evolution equation for Ψ becomes Since Φ n is constant, Eq. (67) yields the duration of the meander period in the presence of a constant external field ⃗ E: Note that the resulting period depends on Φ n , i.e. on the orientation of the trajectory relative to the applied field ⃗ E.
The rotation phase during one cycle can be found by integrating ∂ Ψ Φ = ∂ t Φ ∂ t Ψ: . Evaluation at Ψ n+1 delivers the change in rotation angle over one meander period: Hence, an applied field will change the angle between consecutively visited petals of the meander flower.
Finally, we can find the net spatial displacement during a meander cycle: Thus, the net drift of a meandering spiral in an external field depends on the angle Φ between the meander pattern and the applied field ⃗ E. The average drift speed, measured over one meander period (i.e. one petal of the flower) will be: V a n =d a n T n = R a α (Φ n )Q α β c α + R a α (Φ n )R β b (Φ n )E b + O(η 2 ). with Q α β = m α β + c α M Ψ β Ω If we are outside the phase-locking regime, the ratio β 2π = ω Ω will not be close to a fraction p q with p and q small integers. Therefore, after several rotations, the angle Φ n will have taken different values that are uniformly distributed over the unit circle. In that approximation, the time average of a quantity (net rotation, displacement or time lapse) will be that quantity averaged over all possible phases Φ n at the start of the meander cycle: For the average drift speed, we then find where Since the trace of a matrix is invariant and = α β is the generator of rotations, we furthermore find that for any 2 × 2 matrix P and rotation matrix R holds Tr(RPR T ) = Tr(R T RP) = Tr(P), Tr( RPR T ) = Tr(R T RP) = Tr(R T R P) = Tr( P).

Hence
Expression (71) was also computed in the center frame in [16] using a double Fourier series. We find the same result here using the discrete maps for X a n and Φ n . From Figs. 6 and 7 it can be seen that the net drift velocity is well described by Eq. (71) with coefficients (72).

IX. DISCUSSION
In this work, we have introduced different coordinate systems that allow to describe the dynamics of one-and two-dimensional wave patterns. The origin of the chosen coordinate system is interpreted as the wave front or spiral wave tip position, or as the filament position for 3D scroll waves.
For spiral waves, different coordinate systems produce the same result where they are applicable, but they can differ in their applicability areas. So, the "minimally rotating" frame of reference is not good for describing drift of meandering spirals for parameter values close to the equivariant Hopf bifurction (Winfree's ∂M boundary), and one needs to use the "co-rotating finite-core" frame instead. Similarly, the "center" frame is inadequate for the parameters in the vicinity of the "resonant meander" parameters, and one is much better with the "tip" frame there.
From the simple application of electroforetic drift already, it can be seen that the time-scale at which the motion is analysed is important. Here, we have demonstrated that different time-averaging strategies can be taken in addition to the Fourier approach of [16]: analysis of the return map and a sliding average are possible.
We have deliberately included a brief description of the geometrical viewpoint on phase space dynamics in Sec. VII, which shows that our technical computations are part of a simpler geometric theory.
We have given two applications of the response function framework here, both related to spiral wave drift in a constant external field. First, we have shown that the mean drift speed is equal in the center frame and finite-core frame, as both overlap integrals produce the same physical tension coefficients. Secondly, we have computed that for phase-locked meander, the drift speed is not proportional to the applied field strength, but close to the orbital velocity along the meander flower.
As always, the limitations of the presented work suggest the directions for future research. We mentioned in the beginning of Section VII that the mathematical foundation of our formal approach has known outstanding technical difficulties, e.g. to describe spiral waves the reaction-diffusion system has to be considered in a Banach space in which the action of the Euclidean group is not even continuous, and that the effective spatial localization of the response functions is only an empirical fact. A possible route through the first obstacle was shown by [54] which informally can be summarized thus: rather than considering the whole Banach space, one ought to focus on its part occupied by actual spiral-wave solutions, and in that part, the group acts "nicely". For the second difficulty, a promising rigorous result was obtained in [75]: that response functions of one-dimensional analogues of spiral waves are indeed exponentially localized in space.
From the physical viewpoint, more interesting are the limitations related to the types of perturbations considered. Of special interest is wave propagation in confined irregular geometries such as cardiac tissue. Previous works have shown that these can also to some extent be described by a perturbative approach [15,76], although the boundary perturbation may take the shape of a Dirac distribution and is therefore not small at all.
Note that phase-locking of meandering spirals caused by perturbations that are purely time-dependent have been considered before [77,78]. Our approach outlined in [16] and detailed in here is potentially by far more universal. The simple examples considered so far are admittedly only first steps, and consideration of more realistic problems, such as cardiac tissue with its anisotropy and spatial heterogeneity, promises many new discoveries.