A quantitative theory for phase-locking of meandering spiral waves in a rotating external field

Rotating spiral waves are among the most fascinating examples of self-organizing spatiotemporal patterns that have been observed in diverse systems outside of equilibrium [1, 2]. They occur, for example, as concentration waves in the chemical Belousov–Zhabotinsky (BZ) reaction [3, 4], as coverage patterns of CO oxidation on Platinum [5, 6], or as convection patterns in Rayleigh–Bénard systems [7]. Moreover, spiral waves seem to be ubiquitous in living systems including aggregating slime-mould cells [8, 9], the retina [10] and cardiac tissue [11–13], where they are thought to lead to life-threatening cardiac arrhythmias.

The dynamics of a spiral wave is often characterized and to some extent even determined by the motion of the tip [14, 15]. This point is also known as a 'topological defect' or 'phase singularity' as the excitation phase is ill-defined there. The simplest motion of a spiral wave is rigid rotation, and the corresponding orbit followed by the tip is a circle. Quasi-periodic motions of the spiral tip are also possible and this phenomenon was called 'meander' [16–18]. Such meandering of spiral waves has been experimentally observed, e.g. in BZ systems [19, 20] as well as in cardiac tissue [21]. Via a phenomenological model of 5 ordinary differential equations, the origin of meander with flower-like tip trajectories could be shown to arise via a Hopf bifurcation from the rigidly rotating state [18]. Close to the bifurcation, meander occurs generally in two types called outward meandering and inward meandering, dependent of whether their tip trajectory forms a flower-like orbit with petals pointing outward or inward.

The occurrence and control of spiral wave meander is important to applications. Firstly, under certain conditions, spiral waves can break up due to the strong meandering of the spiral tip [22–24]. This process results in the production of a large number of new interacting vortices evolving into spatiotemporal irregularity which is analogous to the patterns during fibrillation events in the heart. Secondly, it has been hypothesized that the lethal cardiac arrhythmia 'Torsade de Pointes' is caused by a meandering spiral wave [25, 26].

Therefore, how to control or suppress meandering spiral waves is continually attracting attention, and several control schemes have been proposed [20, 27–34]. For example, it has been shown that introducing a feedback stimulus could steer the meandering spiral waves to a desired location [29]. Furthermore, meandering spiral waves can be forced to rigidly rotating spiral waves by applying a proportional control algorithm as well as time delay autosynchronization [33].

Electric fields are also widely employed to control the behaviors of spiral waves. Taking into account the rotation symmetry of the spiral wave around the center, a circularly polarized electric field (CPEF) has been proposed [35] to control the drift of spiral waves. It has been implemented in the BZ experiment [36], which allows us to study the response of spiral waves to a rotating electric field in reaction-diffusion systems. Further, it was shown that the CPEF has some pronounced effects on spiral waves that were not observed in RD systems in the presence only of a dc or ac electric field [34, 37–39].

One particular numerical observation [34] is that applying a rotating electric field to a meandering spiral may convert it to a circular-core spiral that is phase-locked to the external field. The previous work [34] was purely phenomenological and lacks a quantitative explanation. Our present aim is to provide a quantitative theory and to elucidate the mechanism behind the phenomenon of stabilization of meandering spiral waves and synchronous motion of the tip along with the applied rotating electric field. (See our numerical results in figure 1 where we show the initial meandering spiral wave, the phase-locked spiral and the Arnold tongue for phase-locking (PL).) Specifically, in the present work we for the first time provide a theoretical prediction of the asymmetric Arnold tongue region (including the left, right boundary and a finite critical field strength) for the stable phase-locked spiral waves, which matches the outcome from direct simulations very well.

Zoom In Zoom Out Reset image size

To obtain quantitative results, we use response function (RF) theory. Knowing how spiral waves respond to small perturbations is an important issue as it allows us to develop more precise control schemes for practical applications. As demonstrated by Biktasheva and Biktashev, spiral waves behave like a wave-particle dualism in the sense that although spiral waves affect the medium as a delocalized, wave-like entity, they are only influenced by applied forcing when it is delivered close to the spiral tip; in that sense, spirals also act as localized particle-like objects [40]. This is mathematically expressed in terms of the spiral wave RFs, which are essentially nonzero only nearby the spiral tip [40–42]. It was first used by Keener in the context of spiral waves in [43]. Knowledge of RFs allows to quantitatively predict some motions, e.g. spatial drift of spiral waves caused by small perturbation of any nature [40] and PL behavior between the rigidly rotating spiral wave and a rotating electric field [37]. However, such previous theoretical works using RFs are only limited to rigidly rotating spiral waves rather than meandering spiral waves [40–42]. Only very recently, the RF approach was extended such that it is accessible to the case of meandering spiral waves [44, 45]. For instance, a drift law of meandering spiral waves was derived and applied to different external perturbations [45]. Note, however, that for the present application, we will use RFs of a circular-core spiral only.

This paper is organized as follows. In the next section, we present numerical models and methods and as an immediate result the PL region in (ω, E)-space that requires clarification. In our results section we first present theory in E for this which overall predicts the shape of the PL-region (Arnold tongue), first in linear and then in quadratic order in E. Thereafter, we present numerical results on RF calculations and direct numerical simulations of PL, and show that we get an excellent agreement between second order theory and numerical experiments. We conclude this paper by a short discussion, outlook and conclusions.

2.1. Reaction-diffusion model

Our results hold for a class of reaction-diffusion systems, to which a homogeneous external field $\vec{E}$ is applied that rotates at constant frequency ω_E:

$$\begin{eqnarray}&&{\partial }_{t}{\bf{u}}=\hat{{\bf{P}}}{\rm{\Delta }}{\bf{u}}+{\bf{F}}({\bf{u}})-\hat{{\bf{M}}}\vec{E}(t)\cdot \vec{{\rm{\nabla }}}{\bf{u}}.\end{eqnarray} \tag{ 1 }$$

Here ${\bf{u}}(\vec{r},t):{{\mathbb{R}}}^{2}\times {{\mathbb{R}}}^{+}\to {{\mathbb{R}}}^{N}$ describes the temporal evolution of N state variables. $\hat{{\bf{P}}},\hat{{\bf{M}}}\in {{\mathbb{R}}}^{N\times N};$ they describe which variables undergo diffusion or are sensitive to the external field. The nonlinear function ${\bf{F}}:{{\mathbb{R}}}^{N}\to {{\mathbb{R}}}^{N}$ describes local excitation kinetics of the system.

In the numerical examples for this study, we work with Barkley kinetics [46], which is a simple model (N = 2) for excitation-inhibition dynamics that is able to produce both circular-core and meandering spiral waves. It has ${\bf{u}}={\left[u,v\right]}^{T}$, $\hat{{\bf{P}}}=\mathrm{diag}(1,0)$, ${\bf{F}}({\bf{u}})={\left[f(u,v),g(u,v)\right]}^{T}$ with

$$\begin{eqnarray}\begin{array}{rcl}f(u,v) & = & {\epsilon }^{-1}u(1-u)\left(u-\displaystyle \frac{v+b}{a}\right),\\ g(u,v) & = & u-v.\end{array}\end{eqnarray} \tag{ 2 }$$

We furthermore model that only the first variable couples to the external field, i.e. $\hat{{\bf{M}}}=\mathrm{diag}(1,0)$.

Equation (1) has been used before to model chemical systems subjected to an electrical field $\vec{E}$. Constant $\vec{E}$ was shown experimentally and numerically to cause so-called electrophoretic spiral wave drift [47–49].

2.2. Direct numerical simulations of PL

3.1. Theory

3.1.1. Proximity of an unstable circular-core solution

Presently, we study cases where equation (1) with $\vec{E}=\vec{0}$ has a meandering spiral ${{\bf{u}}}_{{\rm{m}}}$ as a solution. However, under rotating external fields $\vec{E}(t)$ with well chosen frequency ${\omega }_{E}$, we see a rigidly rotating solution ${{\bf{u}}}_{E}(\vec{x},t)$ that is stable and phase-locked to the external field.

In the Barkley model that we presently investigate, the flower-like tip trajectory of the spiral indicates that the system is still close to the Hopf bifurcation, e.g. in contrast ionic cardiac models which have meander due to a dynamical wave block. In our case, it means that for the same model parameters as the meandering spiral, there still exists a circular-core spiral wave, which is essentially an unstable limit cycle of the system. Our main idea is that this circular-core solution is still there in the meander regime, and its stability can be restored by the applied external field.

Figure 2(a) shows the numerically obtained bifurcation diagram for varying Barkley parameter a. Let us denote c_min and c_max the minimal and maximal displacement speed of the spiral wave tip during several spiral wave rotations. Then, $({c}_{\max }-{c}_{\min })/2$ plays the role of a bifurcation amplitude.

Zoom In Zoom Out Reset image size

We now hypothesize that under a rotating external field, the previously unstable circular-core spiral regains stability. The underlying reasoning is the following. The meandering spiral in Barkley's model undergoes biperiodic motion and therefore has two phases: rotation phase Φ, which governs the spatial orientation of the spiral and temporal phase Ψ, describing the periodic deformation of the wave in a co-moving frame [45]. Now, suppose that under the applied field only one of the phases is locked. Then, the residual motion relative to the external field would be not constant but periodic. As a result, the resulting motion would still be flower-like in the lab frame of reference, in contrast to the observed circular motion in figure 1(b). Thus, we conclude that in the regime of our interest, both phases are locked to the external field, and the phase-locked spiral is rigidly rotating. Therefore, it makes sense to look at the nearby unstable circular-core spiral wave solution (at E = 0) and investigate whether its properties can explain the observed PL.

A geometric view on this process is offered in figure 2(b). A meandering spiral can be seen as a solution performing two composite motions: performing a full spiral rotation, together with a slow precession of the orientation of the pattern. In general the two frequencies are incommensurate, such that the path of the system in phase space (blue trajectory) does not close on itself but instead traces out a torus shape that is an invariant manifold of the system. If the amplitude of the wave modulation is decreased, the motion reduces to rotation without modulation, and the invariant manifold becomes a single closed orbit (red trajectory). In the present context, the rotating external field can restabilize the rigidly rotating solution (red) under certain conditions.

3.1.2. Critical eigenmodes of the unstable circular-core spiral

Our analysis is performed in a frame rotating at constant frequency ${\omega }_{E}$, equal to the rotation frequency of the externally applied field. In this frame, equation (1) becomes

$$\begin{eqnarray}&&{\partial }_{t}{\bf{u}}=\hat{{\bf{P}}}{\rm{\Delta }}{\bf{u}}+{\omega }_{E}{\partial }_{\theta }{\bf{u}}+{\bf{F}}({\bf{u}})-\hat{{\bf{M}}}{E}^{A}{\partial }_{A}{\bf{u}}\end{eqnarray} \tag{ 3 }$$

with E^A constant. Here and below, summation over repeated translation indices $A,B,\ldots$ is implied. (We never sum over the label E.)

Supposing that this equation has a stationary solution ${{\bf{u}}}_{E}$, i.e.

$$\begin{eqnarray}&&{\bf{0}}=\hat{{\bf{P}}}{\rm{\Delta }}{{\bf{u}}}_{E}+{\omega }_{E}{\partial }_{\theta }{{\bf{u}}}_{E}+{\bf{F}}({{\bf{u}}}_{E})-\hat{{\bf{M}}}{E}^{A}{\partial }_{A}{{\bf{u}}}_{E}.\end{eqnarray} \tag{ 4 }$$

We now verify when it will be linearly stable, such that it can be the phase-locked spiral observed in our numerical simulations.

The associated linear operator to equation (3) is

$$\begin{eqnarray}&&\hat{{\bf{L}}}=\hat{{\bf{P}}}{\rm{\Delta }}+{\omega }_{E}{\partial }_{\theta }+{\bf{F}}^{\prime} ({{\bf{u}}}_{E})-\hat{{\bf{M}}}{E}^{A}{\partial }_{A}.\end{eqnarray} \tag{ 5 }$$

Its adjoint with respect to $\langle {\bf{f}}| {\bf{g}}\rangle ={\iint }_{{{\mathbb{R}}}^{2}}{{\bf{f}}}^{H}{\bf{g}}{\rm{d}}{S}$ is given by

$$\begin{eqnarray}&&{\hat{{\bf{L}}}}^{\dagger }={\hat{{\bf{P}}}}^{T}{\rm{\Delta }}-{\omega }_{E}{\partial }_{\theta }+{\bf{F}}{{\prime} }^{T}({{\bf{u}}}_{E})+{\hat{{\bf{M}}}}^{T}{E}^{A}{\partial }_{A}.\end{eqnarray} \tag{ 6 }$$

In what follows, we consider this as a continuation problem, i.e. the solution ${{\bf{u}}}_{E}$ for general $\vec{E}$ follows from the case $\vec{E}=\vec{0}$, ${\omega }_{E}={\omega }_{0}$, denoted ${{\bf{u}}}_{0}$ with linear operator ${\hat{{\bf{L}}}}_{0}=\hat{{\bf{P}}}{\rm{\Delta }}+{\omega }_{0}{\partial }_{\theta }+{\bf{F}}^{\prime} ({{\bf{u}}}_{0})$. Furthermore, we denote $\delta ={\omega }_{E}-{\omega }_{0}$.

Since ${{\bf{u}}}_{0}$ is a circular-core solution that spontaneously evolves towards a meandering spiral, we expect that ${\hat{{\bf{L}}}}_{0}$ has no eigenvalues with positive real part, except for a pair ${\lambda }_{0}$, ${\lambda }_{0}^{* }$ that crossed the imaginary axis at the Hopf bifurcation at the onset of the meander regime:

$$\begin{eqnarray}&&{\hat{{\bf{L}}}}_{0}{{\bf{V}}}_{0}={\lambda }_{0}{{\bf{V}}}_{0},\quad {\hat{{\bf{L}}}}_{0}^{\dagger }{{\bf{W}}}^{0}={\lambda }_{0}^{* }{{\bf{W}}}^{0},\quad \langle {{\bf{W}}}^{0}| {{\bf{V}}}_{0}\rangle =1.\end{eqnarray} \tag{ 7 }$$

For finite $\vec{E}$ and ${\omega }_{E}\ne {\omega }_{0}$, this eigenvalue will be written λ, with left and right-hand eigenvectors ${\bf{V}}$, ${\bf{W}}$:

$$\begin{eqnarray}&&\hat{{\bf{L}}}{\bf{V}}=\lambda {\bf{V}},\quad {\hat{{\bf{L}}}}^{\dagger }{\bf{W}}={\lambda }^{* }{\bf{W}},\quad \langle {\bf{W}}| {\bf{V}}\rangle =1.\end{eqnarray} \tag{ 8 }$$

Now, the phase-locked solution will be stable if

• (1)  
  the phase-locked solution is stable with respect to rotation
• (2)  
  there exist combinations $(E,{\omega }_{E})$ for which ${\mathfrak{R}}(\lambda )\lt 0$.

In addition to the meander mode (8), ${\hat{{\bf{L}}}}_{0}$ has also a zero mode associated to rotational degrees of freedom: by differentiating equation (4) with E = 0, $\delta =0$ with respect to the polar angle θ, one finds, with ${{\bf{V}}}_{\theta }=-{\partial }_{\theta }{{\bf{u}}}_{0}$:

$$\begin{eqnarray}&&{\hat{{\bf{L}}}}_{0}{{\bf{V}}}_{\theta }={\bf{0}},\,{\hat{{\bf{L}}}}_{0}^{\dagger }{{\bf{W}}}^{\theta }={\bf{0}},\, \langle {{\bf{W}}}^{\theta }| {{\bf{V}}}_{\theta } \rangle =1.\end{eqnarray} \tag{ 9 }$$

3.1.3. PL condition

Note that although the control parameters in the experiment are $E,{\omega }_{E}$, it will be convenient to work with the pair E^A instead. Let us now expand in linear order in E:

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{{\bf{u}}}_{E} & = & {{\bf{u}}}_{0}+{{\bf{u}}}_{A}{E}^{A}+\widetilde{{\bf{u}}},\,\widetilde{{\bf{u}}}={O}({E}^{2}),\\ \lambda & = & {\lambda }_{0}+{\lambda }_{A}{E}^{A}+{O}({E}^{2}).\end{array}\end{array}\end{eqnarray} \tag{ 10 }$$

Our numerical simulations show that the phase-locked state is rigidly rotating, whence ${{\bf{u}}}_{E},{{\bf{u}}}_{A}$ and $\tilde{{\bf{u}}}$ are constant over time in the regime of interest. Therefore, putting equation (10) into (3) delivers

$$\begin{eqnarray}&&{\bf{0}}={\hat{{\bf{L}}}}_{0}({{\bf{u}}}_{A}{E}^{A}+\tilde{{\bf{u}}})+\delta {\partial }_{\theta }{{\bf{u}}}_{0}-\hat{{\bf{M}}}{E}^{A}{\partial }_{A}{{\bf{u}}}_{0}.\end{eqnarray} \tag{ 11 }$$

Now, projection on the rotational RF ${{\bf{W}}}^{\theta }$ yields as a necessary condition for PL:

$$\begin{eqnarray}&&\delta +{P}_{A}{E}^{A}=0,\end{eqnarray} \tag{ 12 }$$

where

$$\begin{eqnarray}&&{P}_{A}=\langle {{\bf{W}}}^{\theta }| \hat{{\bf{M}}}| {\partial }_{A}{{\bf{u}}}_{0}\rangle .\end{eqnarray} \tag{ 13 }$$

This relation can also be written as $-\delta =\vec{E}\cdot \vec{P}$, where $\vec{P}$ is a vector that co-rotates with the spiral. We write $| | \vec{E}| | =E$ and $| | \vec{P}| | =P$. It states that PL can only occur within a frequency band where $\delta \lt {EP};$ within this band the spiral orientation with respect to the applied field $\vec{E}$ is used as a degree of freedom to adapt to different field strengths E. Equation (12) is the sole criterion for PL of a circular-core spiral that is stable itself for E = 0; see [37].

Note that for given $\delta ,E$, equation (12) admits two solutions, only one of which is stable. This can be seen by considering the time-dependent problem: if one lets the spiral wave enclose an angle $\phi (t)$ with $\vec{E}$, one finds from equation (3):

$$\begin{eqnarray}&&{\partial }_{t}\phi =\delta +{P}_{A}{E}^{A}.\end{eqnarray} \tag{ 14 }$$

Linearizing around the phase-locked angle shows that the stable one is characterized by ${P}_{x}{E}_{y}-{P}_{y}{E}_{x}=(\vec{P}\times \vec{E})\cdot {\vec{e}}_{z}\gt 0$.

3.1.4. Stability with respect to the meander mode

To phase-lock a meandering spiral, we should also investigate the sign of ${\mathfrak{R}}(\lambda (E))$. The leading order coefficient in $\lambda ={\lambda }_{0}+{\lambda }_{A}{E}^{A}+{O}({E}^{2})$ can be found using the Feynman–Hellmann relation [50]:

$$\begin{eqnarray}\begin{array}{rcl}{\lambda }_{A} & = & {\left.\displaystyle \frac{\partial \lambda }{\partial {E}^{A}}\right|}_{E=0,\delta =0}\\ & = & {\left.\displaystyle \frac{\partial }{\partial {E}^{A}}\lt {\bf{W}}| \hat{{\bf{L}}}| {\bf{V}}\gt \right|}_{E=0,\delta =0}\\ & = & \langle {{\bf{W}}}^{0}| {\left.\displaystyle \frac{\partial \hat{{\bf{L}}}}{\partial {E}^{A}}\right|}_{E=0,\delta =0}| {{\bf{V}}}_{0}\rangle \\ & = & \langle {{\bf{W}}}^{0}| \left(-{P}_{A}{\partial }_{\theta }+{\bf{F}}^{\prime\prime} ({{\bf{u}}}_{0}){{\bf{u}}}_{A}-\hat{{\bf{M}}}{\partial }_{A}\right)| {{\bf{V}}}_{0}\rangle .\end{array}\end{eqnarray} \tag{ 15 }$$

To evaluate this coefficient, knowledge of ${{\bf{u}}}_{A}$, i.e. the wave profile deformation under an external rotating field is required. By identifying ${O}(E)$ terms in equation (11), we find

$$\begin{eqnarray}&&{\hat{{\bf{L}}}}_{0}{{\bf{u}}}_{A}{E}^{A}=\delta {{\bf{V}}}_{\theta }+\hat{{\bf{M}}}{E}^{A}{\partial }_{A}{{\bf{u}}}_{0}.\end{eqnarray} \tag{ 16 }$$

Let us now project into the space perpendicular to the direction of ${{\bf{V}}}_{\theta }$ using $\hat{{\boldsymbol{\Pi }}}={\bf{1}}-| {{\bf{V}}}_{\theta }\rangle \langle {{\bf{W}}}^{\theta }|$:

$$\begin{eqnarray}&&{\hat{{\bf{L}}}}_{0}{{\bf{u}}}_{A}{E}^{A}=\hat{{\boldsymbol{\Pi }}}\hat{{\bf{M}}}{E}^{A}{\partial }_{A}{{\bf{u}}}_{0}.\end{eqnarray} \tag{ 17 }$$

In the resulting subspace ${\hat{{\bf{L}}}}_{0}$ is invertible, yielding

$$\begin{eqnarray}&&{{\bf{u}}}_{A}={\hat{{\bf{L}}}}_{0}^{-1}\hat{{\boldsymbol{\Pi }}}\hat{{\bf{M}}}{\partial }_{A}{{\bf{u}}}_{0}.\end{eqnarray} \tag{ 18 }$$

Denoting ${\mathfrak{R}}({\lambda }_{0})={\sigma }_{0}$ and ${\mathfrak{R}}({\lambda }_{A})={M}_{A}$, we find that the real part of the growth rate of the meander mode is negative if

$$\begin{eqnarray}&&{\sigma }_{0}+{M}_{A}{E}^{A}\lt 0.\end{eqnarray} \tag{ 19 }$$

Hence, a finite field strength

$$\begin{eqnarray}&&{E}_{{\rm{c}}}={\sigma }_{0}/M\end{eqnarray} \tag{ 20 }$$

is required to convert a meandering spiral into a circular-core one (in contrast to the PL condition, which can be fulfilled for arbitrarily small E).

3.1.5. PL regions in the $({E}^{x},{E}^{y})$ and $({\omega }_{E},E)$-planes

To understand how the linear relations, i.e. conditions (12), (19) in E give rise to the wedge-shaped region shown in figure 1(c), it is instructive to study the coordinate transformation $({E}^{x},{E}^{y})\leftrightarrow ({\omega }_{E},E)$. As an intermediate step, we work in coordinates relative to $\vec{P}$: say $\vec{E}\cdot \vec{P}={E}_{| | },\vec{E}\times \vec{P}={E}_{\perp }{\vec{e}}_{z}$. Then,

$$\begin{eqnarray}&&\left\{\begin{array}{l}E=\sqrt{{E}_{| | }^{2}+{E}_{\perp }^{2}}\\ {\omega }_{E}={\omega }_{0}+{{PE}}_{| | }\end{array},\right.\qquad \left\{\begin{array}{l}{E}_{| | }=({\omega }_{E}-{\omega }_{0})/P\\ {E}_{\perp }=\sqrt{{E}^{2}-{E}_{| | }^{2}}\end{array}.\right.\end{eqnarray} \tag{ 21 }$$

To choose the sign of ${E}_{\perp }$, we note that of the two solutions to equation (12), only the one with ${E}_{\perp }\geqslant 0$ is stable with respect to rotation.

The relation between ${\omega }_{E}$ and ${E}_{| | }$ is linear, such that vertical lines are mapped to vertical lines. Figure 3 shows moreover that the entire upper half plane in $({E}_{| | },{E}_{\perp })$ is mapped to the triangular wedge region in $({\omega }_{E},E)$-space. Conversely, horizontal lines in $({\omega }_{E},E)$ become circles centered at the origin in $({E}_{| | },{E}_{\perp })$ space.

Zoom In Zoom Out Reset image size

The first condition (12) thus bounds the phase-locked region in the $({\omega }_{E},E)$ plane as a triangular wedge centered at ${\omega }_{0}$.

The second condition places a linear cut to this zone in the $({E}_{| | },{E}_{\perp })$-domain. Therefore, the border of the phase-locked region is the union of two curves, respectively produced by conditions (12) and (19). The PL condition (12) yields a linear boundary also in $({\omega }_{E},E)$-space. The second portion of the boundary of the phase-locked zone is related to stability of the phase-locked solution; it is shown to be a hyperbola in $({\omega }_{E},E)$ as follows.

We write the conditions (12) and (19) as

$$\begin{eqnarray}&&\delta +{PE}\cos {\beta }_{{PE}}=0,\end{eqnarray} \tag{ 22a }$$

$$\begin{eqnarray}&&{\sigma }_{0}+{ME}\cos ({\beta }_{{MP}}+{\beta }_{{PE}})=0,\end{eqnarray} \tag{ 22b }$$

taking into account the angles between $\vec{M}$, $\vec{P}$ and $\vec{E}$ during PL: ${\beta }_{{PE}}=\mathrm{atan}2({E}_{y},{E}_{x})-\mathrm{atan}2({P}_{y},{P}_{x})$ and similar for ${\beta }_{{MP}},{\beta }_{{ME}}$.

Eliminating ${\beta }_{{PE}}$ delivers after some arithmetic:

$$\begin{eqnarray}&&{E}^{2}=\displaystyle \frac{{\left(\delta /P-{E}_{{\rm{c}}}\cos {\beta }_{{MP}}\right)}^{2}}{{\sin }^{2}{\beta }_{{MP}}}+{E}_{{\rm{c}}}^{2}.\end{eqnarray} \tag{ 23 }$$

Thus, one side of the boundary of the PL region in $(\delta ,E)$-space is a hyperbola; the lowest field strength for synchronization, ${E}_{{\rm{c}}}$, is given by equation (20) and reached at

$$\begin{eqnarray}&&{\delta }_{c}=-{{PE}}_{{\rm{c}}}\cos {\beta }_{{MP}}.\end{eqnarray} \tag{ 24 }$$

Note that the first condition yields a linear bound with slope $1/P$, while the hyperbolic part has an asymptote with steeper slope $1/(P\sin {\beta }_{{MP}})$. Both curves touch in the point

$$\begin{eqnarray}&&{E}_{t}=\displaystyle \frac{{E}_{{\rm{c}}}}{\cos {\beta }_{{PM}}},\quad {\omega }_{t}={\omega }_{0}-{{PE}}_{t}.\end{eqnarray} \tag{ 25 }$$

The union of the two boundaries delineates an asymmetric wedge, confirming our numerical observation in figure 1(c).

3.2. Quadratic theory

In section 3.4 and figure 6(a) we will show that the linear theory in E exposed above predicts the overall wedge shape of the PL region, but the necessary parameters (M_A, P_A) do not seem to match those obtained by direct calculation of overlap integrals using equations (13) and (15). The reason for this is that in our theory, we linearize around ω_E = ω₀, E = 0, which lies at finite distance from the PL region. Therefore, we also compute quadratic corrections in the $\omega (\vec{E})$ relation, which will show an excellent agreement between theory and direct numerical simulations.

To this aim, we generalize equations (10) to

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}{{\bf{u}}}_{E} & = & {{\bf{u}}}_{0}+{{\bf{u}}}_{A}{E}^{A}+\displaystyle \frac{1}{2}{{\bf{u}}}_{AB}{E}^{A}{E}^{B}+O({E}^{3}),\\ \lambda & = & {\lambda }_{0}+{\lambda }_{A}{E}^{A}+\displaystyle \frac{1}{2}{\lambda }_{AB}{E}^{A}{E}^{B}+O({E}^{3}).\end{array}\end{array}\end{eqnarray} \tag{ 26 }$$

Putting this in the reaction-diffusion equation and projecting onto the rotational RF yields

$$\begin{eqnarray}&&{\omega }_{E}={\omega }_{0}-{P}_{A}{E}^{A}+\displaystyle \frac{1}{2}{H}_{{AB}}{E}^{A}{E}^{B}+{O}({\lambda }^{2}),\end{eqnarray} \tag{ 27 }$$

where

$$\begin{eqnarray}&&\displaystyle \frac{1}{2}{H}_{{AB}}=\langle {{\bf{W}}}^{\theta }| \left[-{P}_{A}{\partial }_{\theta }+\displaystyle \frac{1}{2}{\bf{F}}^{\prime\prime} ({{\bf{u}}}_{0}){{\bf{u}}}_{A}-\hat{{\bf{M}}}{\partial }_{A}\right]| {{\bf{u}}}_{B}\rangle .\end{eqnarray} \tag{ 28 }$$

Now, equation (27) can be used to convert the PL region in the (E_x, E_y) domain (see figure 6(c)) to the (ω_E, E) domain.

From figure 6(c), it can be seen that the point of the phase-locked region with minimal $E=\sqrt{{E}_{x}^{2}+{E}_{y}^{2}}$ is determined by the red line, i.e. ${\sigma }_{0}+{M}_{A}{E}^{A}=0$. Therefore, the minimal field strength for PL is still given by equation (20).

The explicit expression for the boundary of the Arnold tongue can be found as follows. Let us call the eigenvalues of H_AB h₁, h₂ and denote the angle between the first eigenvector and $\vec{P}$ as β_HP. Then, equation (27) becomes

$$\begin{eqnarray}\begin{array}{rcl}0 & = & {\omega }_{E}-{\omega }_{0}+{PE}\cos {\beta }_{{PE}}\\ & & +{E}^{2}(\tilde{a}+\tilde{b}\cos 2{\beta }_{{PE}}+\tilde{c}\sin 2{\beta }_{{PE}})+{O}({\lambda }^{3})\end{array}\end{eqnarray} \tag{ 29 }$$

with $\tilde{a}=-({h}_{1}+{h}_{2})/4$, $\tilde{b}=-\sin (2{\beta }_{{HP}})({h}_{1}-{h}_{2})/2$, $\tilde{c}=-\cos (2{\beta }_{{HP}})({h}_{1}-{h}_{2})/2$. At fixed E, the only degree of freedom here is β_PE, i.e. the orientation of the external field relative to the spiral wave. At the boundary of PL, in the regime of small E, this equation will have a double zero for a given β_PE. Setting the derivative of the right-hand side to zero, and expanding around the leading order solution β_PE = 0 delivers for the critical PL angle:

$$\begin{eqnarray}&&{\beta }_{{PE}}=\displaystyle \frac{2\tilde{c}}{P}E+{O}({\lambda }^{2}).\end{eqnarray} \tag{ 30 }$$

Substitution ${\beta }_{{PE}}=0+{O}(\lambda )$ into equation (29) then delivers

$$\begin{eqnarray}&&{\omega }_{E}={\omega }_{0}-{PE}-(\tilde{a}+\tilde{b}){E}^{2}+{O}({\lambda }^{3})\end{eqnarray} \tag{ 31 }$$

for the left-hand boundary of the PL region in the (ω_E, E) plane. (See the green line in figure 6(b).)

For the rightmost boundary of the Arnold tongue, we eliminate β_PE between equations (22b) and (29). With ${\beta }_{{HM}}={\beta }_{{HP}}+{\beta }_{{PM}}$, we define $B=\sin {\beta }_{{HM}}({h}_{1}-{h}_{2})/2$, $C=\cos {\beta }_{{HM}}({h}_{1}-{h}_{2})/2$, $A=\tilde{a}-B$. After introducing the relative frequency variable

$$\begin{eqnarray}&&X=\delta -{{PE}}_{{\rm{c}}}\cos {\beta }_{{PM}}-A-2{{BE}}_{{\rm{c}}}^{2}\end{eqnarray} \tag{ 32 }$$

one finds

$$\begin{eqnarray}&&X=\pm F\sqrt{{E}^{2}-{E}_{{\rm{c}}}^{2}}+A({E}^{2}-{E}_{{\rm{c}}}^{2})\end{eqnarray} \tag{ 33 }$$

with $F=(P\sin {\beta }_{{PM}}+2{{CE}}_{{\rm{c}}})$. The different signs denote the left and right branch of the red curve in figure 6(b). Equation (33) can be written similarly to (23) as

$$\begin{eqnarray}&&{E}^{2}=\displaystyle \frac{{\left(X-A({E}^{2}-{E}_{{\rm{c}}}^{2})\right)}^{2}}{{F}^{2}}+{E}_{{\rm{c}}}^{2}.\end{eqnarray} \tag{ 34 }$$

The red curve in figure 6(b) shows this improved bound of the PL region.

3.3. Numerical computation of RFs and overlap integrals

We used the open-source package DXSpiral [42] to compute the unstable circular-core spiral solution in Barkley's model for the selected parameters (a = 0.58, b = 0.05, $\epsilon =0.02$). Hereto, we started from circular-core parameters (a = 0.50) and numerically continued the solution into the meander regime.

Computations were performed on a disc of radius R = 10, with N_θ = 128 and N_r = 400 grid cells in the circumferential and radial directions.

Details of numerical methods can be found in [42]. The unstable meander modes ${{\bf{V}}}_{0}$, ${{\bf{W}}}^{0}$ were computed using the Cayley transform. We used $\lambda =0.1+1.3i$ as a first guess of the eigenvalue after which the Arnoldi method in DXSpiral converged to the desired eigenfunction.

The inversion of ${\hat{{\bf{L}}}}_{0}$ in (18) was performed using LU-decomposition, as in [51].

The unstable circular-core solution in Barkley's model is shown in figure 4, together with its critical eigenmodes and adjoint critical eigenmodes (RFs). It can be computed using Newton's method in spite of being unstable with respect to the meander instability. Figure 5 shows the radius of the tip trajectory for the phase-locked spiral waves, the unstable circular-core solution and the meandering spiral with E = 0. Table 1 summarizes the results obtained for the overlap integrals needed in the theory.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 1.  Eigenvalues and overlap integrals computed with DXSpiral. (Barkley model, a = 0.58, b = 0.05,  = 0.02.)

ω0	1.3871	σ0	0.0889
Px	0.2321	Py	0.7073
Mx	0.4543	My	0.3624
P	0.7444	M	0.5811
Hxx	−0.4938	Hyy	1.0310
(Hxy + Hyx)/2	−0.4298

3.4. PL region in Barkley's model

Our theory exposed above predicts the shape of region where a stable phase-locked spiral is observed when the external field is switched on.

Prior to DXSpiral calculations, we manually fitted the coefficients to describe the PL-boundary, see light blue lines in figure 6(a). The asymmetric wedge is well represented by this fitting. However, in this work we go one step further, to predict the actual position of the phase-locked region from RF computations. In this case, figure 6(a) shows that linear theory does not perform very well, as the border corresponding to condition (12) is poorly approximated. However, going to second order in E allows an excellent agreement, as shown in figure 6(b). In figure 6(c), we show that this shape corresponds to two linear cuts in the $({E}^{x},{E}^{y})$-plane, as explained in section 3.1.5.

Zoom In Zoom Out Reset image size

In this work, we have explained a mechanism by which the meandering spiral wave may be suppressed by an externally applied rotating field. Alternative approaches could consist of linearization around the meandering spiral solution, or around the critical solution at the Hopf bifurcation point. These three approaches will be in a sense similar, since both the circular-core solution and the meander solution ought be related to the spiral wave properties at the Hopf bifurcation point. We think therefore that it would be interesting to study RF properties around the Hopf bifurcation point, extending the theories of Barkley [18] and Biktashev [52] to the tangent space of the dynamics.

The detailed calculations which provide excellent agreement between theory and numerical results demonstrate that the observed PL can be explained as a forcing of the spontaneously meandering solution to a nearby rigidly rotating spiral wave solution. This may also explain why we have not found examples of this reduction in parameter regimes or reaction models far from the onset of meander (e.g. with star-like or 'linear' spiral cores). Nevertheless, we have shown that the conversion of meander to circular-core spiral wave is robust as it was observed numerically in 3 different reaction models. The work here only deals with the Barkley model. However, similar numerical results were previously reported in the Oregonator and Bär model [34]. What is more, additional computations in a Bär model with oscillatory kinetics suggest that PL occurs also in the oscillatory regime, that an asymmetric Arnold tongue is found which can be described with our theory using manual fitting of the parameters ω₀, P, E_c and β_PM. In fact, the derivations of our theory are quite general and do not involve the specific model. We only use proximity of a Hopf bifurcation and the existence of rotating spiral wave solutions which seem also valid in the oscillatory spatially distributed systems. Therefore, we expect our results to hold not only in excitable but also in oscillatory systems.

To obtain the first of two conditions (12), we have projected the dynamical system onto the zero mode for rotation. The second condition (19) can be obtained by projecting onto the meander mode, but here we chose to use the Feynman–Hellmann theorem instead as the derivation is shorter.

To see until which field strengths our theory ought be valid, let us estimate the magnitude of the perturbation term in equation (1). First, note that the perturbation term $\vec{E}\cdot {\bf{M}}\vec{{\rm{\nabla }}}{\bf{u}}$ is non-zero where the gradient of state variables is non-vanishing, i.e. near the wave front and wave back away from the core region. There, one is close to the traveling wave regime, such that ${\partial }_{t}{\bf{u}}\approx {c}_{0}| | \vec{{\rm{\nabla }}}{\bf{u}}| |$ with c₀ the plane-wave propagation velocity in the system. Therefore, our theory should be valid when

$$\begin{eqnarray}&&| | {\bf{M}}| | E\ll {c}_{0}.\end{eqnarray} \tag{ 35 }$$

Since we work with $| | {\bf{M}}| | =1$, it means that $E\ll {c}_{0}\approx 3.7$. Therefore, we expect that for $E\approx 0.4$, there will be about 10% error in the quantities computed from linearized theory. This order-of-magnitude calculation is in agreement with the errors found in the predicted width of the Arnold tongue in figure 6(a).

Since we were not satisfied with this magnitude of error, we went to a second-order calculation for the left-hand boundary of the PL region. This step was necessary to get a good match for the boundary for the PL condition at finite distance (δ, E) from the circular-core spiral solution around which our expansion is developed. With this, it was not necessary to perform a similar quadratic expansion for $\lambda (\vec{E})$, which would involve the computation of both ${{\bf{u}}}_{{AB}}={\partial }^{2}{\bf{u}}/\partial {E}^{A}\partial {E}^{B}$ and ${\partial }_{A}{\bf{V}}$. We believe the best evidence for the validity of our theoretical framework is the excellent match between the observed and predicted PL region in (ω_E, E)-space, which is improved when going from linear to quadratic theory.

The reaction-diffusion system (1) with $\vec{E}=\vec{0}$ has been used in a much wider context, e.g. for biological morphogenesis, catalytic oxidation waves and cardiac arrhythmia modeling. It can be expected that in some of these disciplines, other internal or external effects may take up the role of $\vec{E}$, as the simplest coupling of a vector field to the state variables of a reaction-diffusion system is through a gradient term as in equation (1).

The present findings are highly likely to be observed and tested in experiments with the following considerations: (1) The rotating electric field has been implanted in the BZ reactions as shown by Ji et al [36] and the existence of the meandering of spiral waves in BZ systems has also been verified [19, 20]; (2) The work here only deals with the Barkley model, however, similar numerical results were previously reported in the Oregonator and Bär model [34]. That similar results are found in quite different models indicates that the phenomenon is quite robust.

In summary, we have investigated the dynamics of the meandering spiral wave in a reaction-diffusion system under a rotating electric field. In line with previous numerical simulations, we observed that under certain conditions, the meandering spiral can be stabilized to the rigidly rotating spiral that also rotates synchronously along with the applied rotating electric field. To uncover the underlying mechanism, we developed a quantitative theory based on the RFs approach. From the theory, we can predict the stabilization regime and synchronization conditions which are in agreement excellent with the numerical results directly from the reaction diffusion systems. Considering the recent realization of rotating electric fields in BZ systems and the robust phenomena through numerical results, our findings presented here are highly likely to be observed and tested in the laboratory.

This work was supported by the National Natural Science Foundation of China under Grants Nos. 11847212, 11775186 and 11675141. BWL would also thank the support from Natural Science Foundation of Zhejiang Province under Grant No. LY16A050003. The authors thank the team that created DXSpiral and shared it with the world.