# Helicoidal instability of a scroll vortex in three-dimensional reaction-diffusion systems

\[

## Abstract

We study the dynamics of scroll vortices in excitable reaction-diffusion systems analytically and numerically. We demonstrate that intrinsic three-dimensional instability of a straight scroll leads to the formation of helicoidal structures. This behavior originates from the competition between the scroll curvature and unstable core dynamics. We show that the obtained instability persists even beyond the meander core instability of two-dimensional spiral wave.

\]

Spiral waves arising in diverse physical, chemical, and biological systems are now one of the paradigms of non-equilibrium dynamical phenomena . Examples include the Belousov- Zhabotinsky (BZ) reaction  and catalytic oxidation of CO on Pt substrates , concentration waves in colonies of aggregating amoebae , waves in cardiac tissue , and many others. These seemingly unrelated phenomena share a common feature: they often allow for a description within the framework of two-component reaction-diffusion type systems, for which spiral solutions are generic. 

The three-dimensional analog of a spiral wave is a scroll vortex, which can be represented by translating the spiral wave along the third direction. Thus, the point singularity of two-dimensional spiral wave (tip) develops into a line singularity (vortex filament). The dynamics of the vortex filaments in reaction diffusion systems has attracted a great deal of attention  in connection with sudden heart fibrillation, where it is believed that scroll and ring vortices play a crucial role . Extensive numerical simulations  and recent experiments on gel-immobilized BZ reaction  show that in a wide range of parameters the scrolls are unstable and may assume helicoidal or even more complicated dynamic configurations. However, the theoretical analysis performed first by J. Keener, predicts an ultimate collapse of vortex rings and stability of straight vortex filaments . This conclusion was drawn on the basis of multiscale analysis, which shows that a vortex ring is similar to an elastic line with a positive line tension. The persistence of nontrivial vortex configurations and turbulence in reaction-diffusion systems was attributed to a negative line tension of the filament .

In this Letter we demonstrate, on the basis of numerical and analytical calculations, that the formation of helicoidal vortices can be related to the intrinsic three-dimensional instability of a straight scroll, caused by a nontrivial response of the filament core to a bending of the filament. We show that the limit of this instability, resulting in the formation of spontaneous helicoidal vortices, goes beyond the corresponding two-dimensional core meander instability.

The dynamics of a scroll vortex can be consistently described by the two-component reaction-diffusion system

$`_tu=ϵ^2u+{\displaystyle \frac{f(u,v)}{ϵ}},`$ (1)   
$`_tv=\delta ϵ^2v+g(u,v),`$ (2)   
where $`ϵ`$ is a small positive parameter, and $`\delta =D_v/D_u`$ is the ratio of diffusion coefficients of the variables $`v`$ and $`u`$. The functions $`f(u,v)`$ and $`g(u,v)`$ are chosen so as to make Eqs. (1)–(2) excitable . In a wide range of parameters Eqs. (1)–(2) have a spiral wave solution in two dimensions and a scroll vortex in three dimensions.

Let us derive the equation of motion for the core of a weakly-curved scroll vortex subject to the meandering instability . In two dimensions this instability was studied both via numerical simulations of the system (1)–(2), and using the numerical solution of the linearized problem . Recently it was proven that the spiral interface undergoes a core-meander instability via a supercritical Hopf bifurcation, as the diffusion coefficient of the slow field $`\delta `$ decreases . The mechanism of the meandering instability in a certain parameter limit has been elaborated on in a recent work .

The core meander was described by five equations for the frequency, coordinates, and the velocity of the spiral tip . These equations can be considerably simplified if, instead of the spiral tip, one considers the instant center of spiral rotation. At the threshold of instability one obtains a single complex Landau-type equation for the “complex” velocity of the rotation center $`\widehat{C}=c_x+ic_y`$, where $`c_{x,y}`$ are the components of the velocity:

$$_t\widehat{C}=\alpha \widehat{C}\beta |\widehat{C}|^2\widehat{C}.$$
(3)   
Here $`\alpha `$ and $`\beta `$ are complex coefficients, $`\alpha =\alpha _1+i\alpha _2,\beta =\beta _1+i\beta _2`$, that have to be determined numerically. For $`\alpha _1<0`$ the symmetry center is stable, and the spiral makes a pure rotation. When $`\alpha _1`$ becomes positive, the fixed point solution of Eq. (3) loses stability, and the rotation center itself performs a circular motion, which implies the meandering (composite rotation) of the spiral tip . As follows from Eq. (3), the absolute value of the velocity of the rotation center, in the saturated meandering regime, is $`C_0=\sqrt{\alpha _1/|\beta _1|}`$, and the corresponding rotation frequency is $`\omega _0=\alpha _2\alpha _1\beta _2/\beta _1`$.

Consider now a three-dimensional weakly curved scroll (the curvature $`\kappa `$ in the third dimension is much smaller then the local curvature of the spiral front). In this limit one expands the Laplacian in equation (1) as follows, $`^2_{2D}^2\kappa 𝐍`$, where $`_{2D}^2`$ is the Laplacian in the cross-section, and $`𝐍`$ is the unit vector pointing toward the center of the filament curvature. Therefore, the curvature in the third dimension plays the role of an advective field directed along $`𝐍`$ and causing the drift of the filament . To lowest order, the curvature $`\kappa `$ enters into the equation of motion (3) linearly:

$$_t\widehat{C}=\alpha \widehat{C}\alpha \gamma \kappa +\mathrm{}$$
(4)   
where $`\gamma =\gamma _1+i\gamma _2`$ is a (complex) constant which can be determined numerically from two-dimensional simulations  or analytically in a large core limit .

Equation (4) readily implies the short-wavelength instability of a straight filament. Indeed, for an almost straight filament parallel to the $`z`$-axis, the curvature vector $`\kappa 𝐍`$ is simply $`(x_{zz},y_{zz})`$. Using $`\widehat{C}=_t\widehat{x}`$, where $`\widehat{x}=x+iy`$, we obtain from Eq. (4), for a periodic perturbation of the filament $`\widehat{x}(z)\mathrm{exp}[\lambda (k)t+ikz]`$, that $`\lambda ^2=\alpha (\lambda +\gamma k^2)`$. The latter equation has the following unstable solution

$$\lambda (k)=\frac{\alpha }{2}+\sqrt{\frac{\alpha ^2}{4}+\alpha \gamma k^2}.$$
(5)   
Eq. (5) represents the eigenvalue of the linearized problem (4) and is valid only for slightly curved filaments, i.e. for small wavenumbers $`k`$. If this near-threshold instability saturates, the saturated structure corresponds to the most unstable mode, which implies the formation of stable helicoidal vortices.

By analogy with the vortex filaments in the CGLE , we can expect that, for $`k>0`$, the growth rate of the obtained three-dimensional instability, $`Re[\lambda (k)]`$, substantially exceeds that of the two-dimensional meandering instability, $`Re[\lambda (0)]=Re[\alpha ]`$. In that case the filament curvature plays a destabilizing role in the dynamics of the filament. To analyze the stability, one can find the coefficients $`\alpha `$ and $`\gamma `$ from the dynamics of the two-dimensional system (1)–(2), and then substitute them into Eq. (5).

To determine $`\alpha `$ and $`\gamma `$ numerically, we have simulated Eqs. (1)–(2) in two dimensions using the EZ-spiral code of D. Barkley . We have chosen Barkley’s model with the functions $`f(u,v)=u(u1)[uu_{th}(v)]`$ and $`g(u,v)=uv`$, where $`u_{th}(v)=(v+b)/a`$. To determine numerically $`\alpha ,\beta `$ and $`\gamma `$, one should start with an unstable rigidly rotating spiral, which is not available in numerical simulations. To overcome this difficulty, we used the following approach. We started with a spiral with already developed meander. Then we applied a localized control technique, developed in Ref. , to turn the spiral motion to a pure rotation around its symmetry center. Following Ref. , we applied a pinning source to equation (2), in the form of a localized inhibiting term $`\mu h(rr_0)`$, where $`r_0`$ is the coordinate of the domain center. Here $`\mu `$ plays the role of a control parameter and is governed by the equation

$$_t\mu =a_1\mu +b_1[v(r_0)v_0],$$
(6)   
where $`v_0`$ is the expected value of the slow field at the spiral center, and $`a_1`$, $`b_1`$ are coefficients that should be properly chosen. We have elaborated the controlling technique by allowing $`v_0`$ to vary adiabatically (slowly compared to $`\mu `$) as follows

$$_tv_0=a_2\mu +b_2[v(r_0)v_0],$$
(7)   
with another pair of coefficients $`a_2`$, $`b_2`$. In this way, we have obtained a purely rotating spiral with the parameter values corresponding to a meandering regime. The final value of the control parameter was as small as $`\mu =10^4`$, which made the source term added to equation (2) a small perturbation. After “switching off” the control (setting $`\mu =0`$), we allow the linear meander instability to develop.

In order to find the constant $`\alpha `$, we fitted the trajectory of the spiral tip to the developing meander, in the form $`x(t)=x_0+x_1\mathrm{sin}(\omega _1t+\varphi _1)+x_2\mathrm{exp}(\alpha _1t)\mathrm{sin}(\alpha _2t+\varphi _2)`$, where $`x_0`$ is the average position of the spiral tip, $`\omega _1`$ is the main spiral rotation frequency, and $`\alpha =\alpha _1+i\alpha _2`$. The same fitting was carried out for $`y(t)`$. Thus, we have retrieved the motion of the center of spiral rotational symmetry. The coordinates of the center, $`x_c`$ and $`y_c`$, as functions of time are given in Fig. 1. We have found that, for the set of parameters given in the caption to Fig. 1, $`\alpha _10.012`$ and $`\alpha _20.336`$. To find the constant $`\gamma `$ we have applied a homogeneous advective field to equation (1), along the $`x`$-axis, as we have done in Ref. . Then $`\gamma `$ is determined from fitting the spiral tip trajectory to the meandering in the above form with an additional term $`\gamma _1Et`$ ($`\gamma _2Et`$), in the expression for $`x(t)`$ ($`y(t)`$), due to the drift of the spiral. The result for $`\gamma `$ is $`\gamma _11.312`$ and $`\gamma _20.154`$. In principle, it is easy to obtain the value of the parameter $`\beta `$ by matching the spiral center trajectory to the saturated meander. However, we do not do it because it is not important for the studied mechanism of the instability.

We have performed numerical simulations of the three-dimensional Eqs. (1)–(2). We have studied the behavior of an almost straight scroll. We have prepared a two-dimensional purely rotating (unstable) spiral, with the parameter values in the meandering regime, using the elaborated controlling technique described above, and translated it along the third dimension to build an initial scroll. Then, a periodic perturbation $`\mathrm{exp}[ikz]`$ has been applied. The three-dimensional plot of the helix, obtained as the result of the instability saturation, is given in Fig. 2. In Fig. 3 we have plotted the real and imaginary parts of the eigenvalue $`\lambda `$ versus $`k`$, both for the theoretical prediction given by (5), and for the results of the simulations. As we have predicted, the observed instability of the filament with initial perturbation of a finite $`k`$ appears to be substantially stronger than the two-dimensional core meander instability ($`k=0`$). We have also tried initial conditions other than purely periodic, and have always found a saturated helix. Almost periodic initial conditions have been mainly used to cut CPU time (the transient from random initial conditions is very long).

As we can see from Fig. 3, the $`k0`$ asymptotics of $`Re[\lambda ]`$ and $`Im[\lambda ]`$ obtained numerically are well matched by our theoretical prediction. As $`k`$ increases, however, the small curvature approximation used for the derivation of (5) no longer works, and the prediction becomes invalid. This is why we see the divergence of the theoretical curves from the numerical ones for larger $`k`$.

We have also studied numerically the dynamics of filaments in the space of parameters $`a`$ and $`b`$ of Barkley’s model . We have performed the stability analysis of the harmonic with the growth rate near the maximal one (see Fig. 3(a) ). In Fig. 4 we plot the bifurcation line of the obtained three-dimensional instability together with that of the core meander instability. We see that the curvature of the filament enhances its instability, so that it is unstable even in a domain of the $`a`$$`b`$ plane, for which the two-dimensional spiral is stable.

In conclusion, we have demonstrated that the intrinsic three-dimensional instability of a straight scroll leads to the formation of stable helicoidal structures. We have shown that persistent curved vortex configurations are not necessarily related to a ”negative line tension” of the filament, but originate from the underdamped core dynamics. Our current analysis is restricted to untwisted scrolls. We expect that twisted scrolls are unstable even in a wider space of the parameters, and exhibit in general even more violent dynamics. Helicoidal structures with twist were observed recently in gel-immobilized BZ reaction . The twist was generated by an external temperature gradient along the scroll axis. The helicoidal instability was observed above specific value of the twist and disappeared when the temperature gradient (and, therefore, the twist) was removed. Presumably, by tuning the parameters of the BZ reaction, one may approach the limit of the intrinsic three-dimensional instability for untwisted scrolls. We can also speculate that the helicoidal instability of scroll vortices in large aspect ratio reaction-diffusion systems is one of the mechanisms that supports more complex vortex configurations and drives turbulent states.

We are grateful to R. Goldstein, A. Pertsov, and S.-Y. Chen for illuminating discussions. The work of IA was supported by the U.S. Department of Energy under contracts W-31-109-ENG-38 and by the NSF, Office of STC under contract No. DMR91-20000.