Traveling spiral wave chimeras in coupled oscillator systems: emergence, dynamics, and transitions

Systems of coupled nonlinear oscillators often exhibit states of partial synchrony in which some of the oscillators oscillate coherently while the rest remain incoherent. If such a state emerges spontaneously, in other words, if it cannot be associated with any heterogeneity in the system, it is generally referred to as a chimera state. In planar oscillator arrays, these chimera states can take the form of rotating spiral waves surrounding an incoherent core, resembling those observed in oscillatory or excitable media, and may display complex dynamical behavior. To understand this behavior we study stationary and moving chimera states in planar phase oscillator arrays using a combination of direct numerical simulations and numerical continuation of solutions of the corresponding continuum limit, focusing on the existence and properties of traveling spiral wave chimeras as a function of the system parameters. The oscillators are coupled nonlocally and their frequencies are drawn from a Lorentzian distribution. Two cases are discussed in detail, that of a top-hat coupling function and a two-parameter truncated Fourier approximation to this function in Cartesian coordinates. The latter allows semi-analytical progress, including determination of stability properties, leading to a classification of possible behaviors of both static and moving chimera states. The transition from stationary to moving chimeras is shown to be accompanied by the appearance of complex filamentary structures within the incoherent spiral wave core representing secondary coherence regions associated with temporal resonances. As the parameters are varied the number of such filaments may grow, a process reflected in a series of folds in the corresponding bifurcation diagram showing the drift speed s as a function of the phase-lag parameter α.


Introduction
Coupled oscillator systems have played a fundamental role in our understanding and modeling of physical systems since the seminal work of Huygens [1,2].More recently, coupled oscillator systems have been used to model certain aspects of neural activity in the brain [3][4][5], collective dynamics of cilia carpets coupled by hydrodynamic interaction [6,7], as well as oscillatory and excitable media such as those modeling chemical oscillations [8][9][10][11].In all these examples the key question of interest is the extent of synchronization, both in frequency and in phase, within the oscillator population and its dependence on the oscillator properties and the nature of their coupling.Systems of two coupled oscillators are well described by the Adler equation [12] for the phase difference between the oscillators, but substantial progress is possible in the limit of an infinite number of oscillators, both with all-to-all coupling and with sparse coupling corresponding to different types of oscillator networks.The theory distinguishes between a phase description, in which the amplitude dynamics are adiabatically eliminated, and an amplitude-phase description that is required for highly nonlinear systems.In either case one tends to find clusters of synchronized oscillators, that is, states of partial synchrony.The theory has numerous key applications, ranging from coupled laser arrays [13] used to maximize power output to the working of the electrical power grid [14].In each case the identification of the clusters, their number and extent, as well as their stability, represent a major challenge to our understanding of large arrays of coupled oscillators.
If the oscillators are identical and interact attractively, one may expect that all will ultimately oscillate in synchrony.That this is not inevitable was pointed out in 2002 by Kuramoto and Battogtokh [15], who discovered a remarkable state of partial synchrony in a system of coupled identical phase oscillators, in which a subset of the oscillators synchronizes but the remainder remains incoherent.Subsequently called chimera states [16], these states of partial synchrony have been studied by numerous authors in recent years [17][18][19][20][21][22].In the simplest case, the oscillators are on a ring and coupled nonlocally via a periodic coupling function.In the presence of a phase lag α in the coupling, the resulting coherent region may drift through the system [23], with oscillators at the leading edge kicked into synchrony while those at the trailing edge fall out of synchrony, thereby maintaining, at least approximately, a constant size of the traveling coherent region.The dynamics of these synchronization and desynchronization fronts thus provide a key to the understanding of the formation of traveling chimera states and their stability properties.
Chimera states are also found in systems of nonidentical oscillators.While less surprising, the resulting states of partial synchrony exhibit similar properties, and in particular motion, whenever the parameter α is nonzero.We distinguish this type of system, in which motion arises spontaneously, from situations in which the motion is forced, for example, via an asymmetrical coupling function [24].
While studies of one-dimensional coupled oscillator arrays are common, similar studies of two-or even three-dimensional oscillator arrays are less frequent.In planar arrays of identical oscillators, one tends to find incoherent cores, surrounded by a rotating coherent spiral wave [25][26][27] somewhat reminiscent of the rotating spiral waves familiar from oscillatory or excitable media [28,29].One may also find coherent cores embedded in an incoherent background [30,31], although stripe and spot patterns are also possible [32].In some cases, these may become destabilized via a Hopf bifurcation, leading either to standing oscillating chimera states or to traveling structures as predicted by abstract theory [33].Three-dimensional arrays support a greater variety of states, most of which are only known via direct numerical simulations [34].Remarkably, in the special case of a sinusoidal coupling function it is possible to establish the existence and stability properties of such states semi-analytically [31,35,36], and in particular to predict the onset of spontaneous motion and even the stability of the resulting moving structures.This is so even though the transition from stationary to time-dependent chimeras is accompanied in general by the appearance of complex filamentary structures within the incoherent core representing secondary coherence regions associated with temporal resonances between the spiral wave frequency and the spatial translation.These structures, first observed in [35], are a property of quasiperiodic states as explained in [31], and form regardless of whether this state is a standing oscillation as in [31,35] or, as shown here, a traveling chimera.We show here that the properties of these resonance structures are responsible for much of the remarkable complexity of the associated bifurcation diagram.
Most of the above results have been obtained using Kuramoto's model of phase-coupled oscillators, although some recent work has been devoted to more realistic oscillator systems, among which coupled Stuart-Landau oscillator systems are most popular [37].However, arrays of both coupled van der Pol oscillators [38] and coupled FitzHugh-Nagumo oscillators [10] have also been studied from this point of view.These models extend the work on the Kuramoto model to include amplitude dynamics in addition to the phases, and in the case of the van der Pol oscillators, to coupled relaxation oscillators, i.e. to oscillators with a strongly nonlinear phase evolution, as well as to excitable systems.
From an experimental perspective, there is a great deal of evidence for the existence of chimera states in one-dimensional arrays of coupled oscillators [8,11,[39][40][41][42], while two-dimensional oscillator arrays (not to mention three-dimensional arrays) remain poorly studied.To the best of our knowledge, the only example of their laboratory realization involves nonlocally coupled Belousov-Zhabotinsky chemical oscillators [8,9].In these experiments, spiral wave chimeras were indeed observed, but usually in the form of moving structures.In some cases, this motion resembles a two-dimensional random walk that may be associated with finite-size fluctuations by analogy with one-dimensional systems [43].On the other hand, spiral wave chimeras with persistent drift motion were also reported.Motivated by the latter observation, we seek the simplest model for studying such drifting structures, and one that allows a detailed study of the emergence and stability properties of uniformly drifting spiral wave chimeras and their parameter dependence, as well as their relationship to other synchronization patterns in the system.
Since the mathematics behind traveling partially coherent states in more realistic systems involving both amplitude and phase dynamics remains largely beyond current reach, we revisit here the Kuramoto model, within which traveling structures can be simulated and, depending on the coupling function adopted, computed semi-analytically, at least in the continuum limit described by the Ott-Antonsen ansatz [44].It is important to realize that for identical oscillators the Ott-Antonsen approach precludes a simple self-consistent description of traveling chimera states (see [45, lemma 2] and [46, section 4]).However, this is no longer the case when the oscillators are nonidentical, and in this paper we therefore assume that the oscillator frequencies are drawn from a Lorentzian frequency distribution of width γ.Thus γ becomes an additional (and key) parameter of the system that manifests itself as a damping term in the continuum description.
In recent work [47] we have studied a discrete two-dimensional oscillator system of this type with a nonlocal coupling function and showed that this system is able to support bound states of two counter-rotating spiral waves with incoherent cores that drift either rigidly or exhibit more complex meandering motion.The specific example considered in [47] was a two-dimensional array of phase oscillators {θ jk (t)} N j,k=1 evolving according to This equation implies that each oscillator (j, k) interacts only with its neighbors within the circular region where the distances m − j and n − k are considered mod N and σ ∈ (0, 1/2) is the relative coupling radius.The interaction is normalized by the number of points |B σ (j, k)| in the region B σ (j, k) and involves a phase lag parameter α ∈ [0, π/2).In addition, it is assumed that the oscillators are heterogeneous in the sense that their natural frequencies ω jk are drawn randomly and independently from a Lorentzian distribution of width γ > 0.
The system (1) was found to exhibit a large variety of different moving spiral wave chimeras associated with the complex spatial structure of their incoherent cores, as summarized in figure 1 for two different values of the phase lag parameter α.In particular, these states exhibit staggered coexistence as a function of α (see below), behavior that is associated with different numbers of crescent-shaped filaments in the core of a moving spiral, hereafter referred to as fingerprint patterns (see figure 1(b-d)).Similar slanted snaking bifurcation diagrams have been observed for spatially localized states in both fluid and optical systems [48][49][50][51][52] and are a consequence of the nonlocal nature of the system (1).As a result the relationship between the presence of stable, albeit moving chimera states and the system parameters is exceedingly intricate and remains to be elucidated.
In this paper we are able to identify, for the first time, the main prerequisites necessary for the emergence of moving spiral wave chimeras of different types and to establish their stability properties.To this end we formulate a version of problem (1) that is tractable semi-analytically, and employ extensive numerical continuation to follow distinct states through parameter space, together with their stability properties.This approach builds a picture of the parameter space of the problem, enables us to identify the different chimera states that are possible, and ultimately allows a detailed understanding of the system.Our results are corroborated using extensive direct numerical simulations of this system and provide a roadmap for understanding more realistic coupled oscillator systems.

Results
The model (1) is a particular case of a more general nonlocally coupled system G jk;mn sin θ jk − θ mn + α with where G(x, y) is a non-constant coupling function, which is 2π-periodic with respect to x and y and satisfies the symmetry conditions and the normalization condition Indeed, if we assume in the square domain (x, y) ∈ [−π, π] 2 , then equation (3) reduces to model (1).The resulting coupling function G(x, y) is shown in figure 2(a).It is well-known [19,54,55] that in the continuum limit N → ∞ the long-term dynamics of equation (3) settle down on the Ott-Antonsen manifold parametrized by a complex-valued function z(x, y, t).This function is the local order parameter at the position (x, y) at time t and quantifies the synchronization degree of oscillators θ jk (t) with (−π + 2π j/N, −π + 2π k/N) ≈ (x, y).Importantly, in the case of Lorentziandistributed natural frequencies ω jk the evolution of z(x, y, t) is described by the integro-differential equation where the damping parameter γ > 0 is determined by the width of the distribution in equation ( 2) and is a convolution-type integral operator with the coupling function G(x, y) from (3).The above observation allows us to use equation ( 6) as a mathematical tool for investigating the properties of moving spiral wave chimeras in the system (3).The discrete system (3) is invariant under discrete translations in two directions, while the continuum description ( 6) is invariant under continuous translations.Both are in addition invariant with respect to the group D 4 of rotations and reflections of a square inherited from the coupling function.The problem ( 6) is therefore invariant under the group T 2 +D 4 , the semidirect product of the two-torus of translations and the discrete group D 4 .The solutions of this problem may respect certain subgroups of this symmetry group or have no symmetry.This observation applies to both stationary and drifting states.In particular, every uniformly drifting state satisfying equation ( 6), including the simplest moving spiral wave chimeras, corresponds to a solution of the form where (s x , s y ) T ∈ R 2 is the velocity vector and Ω ∈ R is the collective phase frequency.In the case s x = s y = 0, the pattern is called motionless or stationary.If, in addition, we also have Ω = 0, the corresponding pattern is called static.From the symmetry conditions (4) it follows that among all possible moving solutions (7) there are two special solution types related to these symmetries: solutions of the form z (x, y, t) = a (x, y − st) e iΩt (8) with s ∈ R and a(−x, y) = a(x, y) (hereafter Z 2 symmetry) corresponding to uniform motion in the y direction, and solutions of the form with s ∈ R and a(y, x) = a(x, y) (hereafter Z2 symmetry), corresponding to uniform motion along a diagonal.We refer to these solutions as symmetric spiral waves but distinguish them by their symmetries Z 2 and Z2 under reflection.Note that rotations by 90 • rotate a Z 2 -symmetric state into another Z 2 -symmetric state, and similarly for Z2 -symmetric states.

Top-hat coupling function
The bifurcation diagram for Z 2 -symmetric spiral wave chimeras (8) with the top-hat coupling function (5) with σ = 0.4 is shown in figure 3. The diagram shows the speed s as a function of the phase lag parameter α, obtained by inserting ansatz (8) into equation ( 6), discretizing the resulting equation on a square grid with 128 nodes in x and y directions and using an arc-length continuation scheme and a standard Newton solver to follow the branch of equilibria of the resulting system of 128 × 64 nonlinear equations.(Note that owing to the reflection symmetry of a(x, y), the number of equations is reduced by half.)This procedure allows us to compute both stable and unstable solutions; in all cases, the predictions of this approach were confirmed using direct numerical simulations of the discrete system (1).cores as α increases.In both regions the bifurcation diagram resembles slanted snaking [48][49][50][51][52].As discussed below, these filamentary structures are a consequence of the two-frequency nature of the traveling chimera state: its global oscillation frequency Ω and its translation frequency determined by the speed s.
Figures 4 and 5 show the corresponding behavior in regions III and IV.The former shows the gradual approach of the two cores to one another along the branch, resulting in the presence of filaments that extend across the whole domain beyond the location indicated by the red star, while the latter shows the gradual onset of the translation motion from α = 0 until α = 0.4 where the speed s is large enough for prominent filaments to be present in the spiral cores.Finally, figure 6 shows the effects of decreasing the top-hat radius σ on the behavior in region I.We observe that as σ decreases so does the corresponding value of α resulting in smaller speed s and a smaller spiral core, all for a given number of core filaments.Overall, however, the behavior remains qualitatively unchanged.We remark that even with a relatively small number of discretization points the time needed to calculate the diagram shown in figure 3 turned out to be extremely long (ca.6 weeks on a dedicated computer with large RAM).In contrast, the results we are going to describe in the next section were obtained much more rapidly (ca.3-4 days on a laptop with a double number of discretization points in each direction).This remarkable computational speed-up was achieved thanks to the use of a special analytical technique for calculating periodic orbits in the Ott-Antonsen manifold proposed in [56] by one of the authors of this paper.

Trigonometric coupling function
The results for the top-hat coupling function ( 5) illustrate some of the complexities inherent in the present problem.Since our ultimate goal is to understand the properties of spiral wave chimeras in the case of general coupling functions G(x, y), we need a deeper analytical approach to this problem.This becomes possible if we limit ourselves to a narrower class of D 4 -symmetric functions of the form where A and B are real parameters such that A 2 + B 2 = 0.This expression arises as a truncation of a Fourier series representation of the coupling function G(x, y), assumed to be an absolutely integrable 2π-periodic function: where Moreover, if this function satisfies (4), then In addition, if its integral is normalized to the identity, ˆπ −π ˆπ −π G (x, y) dx dy = 1, then ĝ0,0 = 1/(2π) 2 .The simplest non-constant truncation of ( 11) is obtained on keeping terms with indices n, m = 0, ±1 only, leading to (10) with For example, if we calculate the Fourier coefficients ĝ1,0 and ĝ1,1 of the top-hat coupling function ( 5) and insert them in the above formulas, we obtain the fit σ → (A(σ), B(σ)) shown in figure 7.
In the following we study the resulting spiral chimeras as a function of the coefficients A and B in the relevant range revealed by figure 7, together with their stability properties, starting with stationary two-core chimeras for α = 0.The damping parameter γ is fixed at γ = 0.01.

Static patterns
For α = 0 equation ( 6) is variational and all chimera states are therefore time-independent and motionless, and correspond to equilibria of the form z (x, y, t) = a (x, y) .
In section 3.1 below, we construct the self-consistency equation determining these equilibria and explain how we solve it.In addition, we describe how the linear stability analysis of such equilibria can be performed.Using these mathematical tools, we calculate bifurcation diagrams for various values of the coupling parameters A and B. Specifically, we show that all static solutions of equation ( 6) with the coupling function (10) and α = 0, including two-core chimeras, have the form where w (x, y) ∈ Span {1, cos x, cos y, sin x, sin y, cos x cos y, cos x sin y, sin x cos y, sin x sin y} .
Patterns with symmetry under the reflection x → −x, i.e. with a(−x, y) = a(x, y), are Z 2 -symmetric, and are described by a subset of the admissible functions w(x, y), namely w (x, y) ∈ Span {1, cos x, cos y, sin y, cos x cos y, cos x sin y} , while patterns satisfying a(y, x) = a(x, y), i.e. with reflection symmetry in the diagonal, are Z2 -symmetric, and are described by functions w(x, y) of the form w (x, y) ∈ Span {1, cos x + cos y, sin x + sin y, cos x cos y, sin x sin y} .
Among these states we distinguish between fundamental states that are independent of the parameter B in (10), and compound states that depend on B. The former include the following: (c) Partially coherent splay state a (x, y) = qe iy , q ∈ (0, 1) .
The above states are all Z 2 -symmetric, but similar expressions can be written for fundamental Z2 -symmetric states.This is a consequence of the following.

Equivalence of the cases A = 0 and B = 0
From the identity it follows that every pattern observed for (A, B) = (A 0 , 0) has its counterpart rotated by the angle π/4 and observed for (A, B) = (0, 2A 0 ).Therefore, the dynamics of equation ( 6) in the case A = 0 and in the case B = 0 are identical modulo the above spatial rotation and rescaling.This result is independent of the value of the phase lag parameter α.

Bifurcation diagrams for static patterns
Thanks to the simplicity of the above formulas, the existence and stability properties of the completely incoherent state, the partially coherent uniform state and the partially coherent splay state can be described comprehensively.In particular, in section 3.1 we show that: (i) the completely incoherent state is linearly stable if γ ⩾ max(1/2, A/4, B/8) and unstable otherwise, (ii) the partially coherent uniform state exists only for γ < 1/2 and is linearly stable if A ⩽ 2 and B ⩽ 4 and unstable otherwise.
The stability analysis of the partially coherent splay state can be performed by generalizing the analytical scheme proposed in [35,57]; the stability analysis for four-core spiral patterns in the case B = 0 was performed in [31,36].
To display our results, including stability results, we employ the global order parameter However, a number of the states listed above have the same global order parameter, R = 0, as the completely incoherent state.To distinguish among these states, we introduce the quantity measuring the contribution of sin y in the Fourier expansion of a(x, y).(green curve, state (d)).For this state |a(x, y)| has period π in the y direction.These states lose stability at a secondary pitchfork bifurcation creating a branch of fundamental Z 2 -symmetric solutions with period 2π in |a(x, y)| in the y direction (pink curve, state (e)).This new arc-shaped branch in figure 8(a) is everywhere unstable and connects the green branch of period π states with the branch of partially coherent splay states with R = 0.
Along the pink branch of period 2π states there are two tertiary pitchfork bifurcation points, with a S-shaped branch of Z 2 -symmetric states modulated in both x and y directions in between (orange curve, state (h)).These states correspond to stationary compound Z 2 -symmetric two-core spiral chimeras and so depend on the parameter B. At B = 0 these states are all unstable (dashed orange curve).
On the dashed brown branch we find unstable compound D 2 = Z 2 2 -symmetric four-core chimeras (state (j)) described by the formula (12) with w(x, y) ∈ Span{1, cos x, sin y, cos x sin y}.This branch connects the period π fundamental states with the fundamental D 4 -symmetric four-core chimeras with R = 0.A tertiary branch of unstable compound states (dashed red curve, state (k)) bifurcates from the brown branch.(ii) B = 0: to distinguish between the different states with R = 0 we show in figure 9 the same bifurcation diagram but showing the quantity √ R 2 + S 2 as a function of A instead of R. In this figure, the yellow line shows the partially coherent splay states (state (c)), while the blue line corresponds to D 4 -symmetric four-core spiral chimeras (state (f)).We see that each of these states bifurcates from the completely incoherent state a(x, y) = 0 at A = 4γ and loses stability with decreasing A in subcritical bifurcations, generating unstable states with R > 0.Moreover, two other unstable fundamental branches, a branch of anti-phase states (light green curve, state (d 0 )) and a branch of generalized four-core patterns (dark purple curve, state (g)) also emerge at the same parameter value, A = 4γ, indicating the highly degenerate nature of this point.(iii) B > 0: figure 8 shows how the behavior of the compound states changes with the coefficient B in the coupling function.The figure shows that while the fundamental Z 2 -symmetric solution is independent of B its stability may change as B changes.A similar statement applies to fundamental Z2 -symmetric solutions.The stability of tertiary states may likewise change.In particular, at B = 2, there appears a narrow parameter range A ≈ (1.93, 1.95) with stable Z 2 -symmetric two-core spiral chimeras.This range corresponds to the segment of the orange solution branch with negative slope.(iv) At B = 4, the segment of the orange branch with negative slope is broader but the stability range of Z 2 -symmetric spiral chimeras is limited by a (subcritical) quaternary bifurcation to (unstable) asymmetric two-core states (purple branch) that connect the states with Z 2 symmetry to similar states with Z2 symmetry (blue branch, state (i)).By asymmetric we mean any state that is neither Z 2 -nor Z2 -symmetric.Stability calculation points to a narrow region of bistability between the Z 2 -and Z2 -symmetric states.Stable asymmetric states are present at larger values of B (purple branch, state (ii)).(v) At B = 5, the Z 2 -symmetric spiral chimeras are unstable, but Z2 -symmetric spiral chimeras can still be stable in a certain parameter range (blue branch, state (i)).A branch of Z2 -asymmetric spiral chimeras (light green branch, state (iii)) bifurcates from the Z2 -symmetric spiral chimeras and some of these may also be stable.We call the resulting states super-asymmetric to distinguish them from the other states previously called asymmetric.The reason for this terminology is the following.Every Z 2 -or Z2 -symmetric state as well as every asymmetric state is represented by expression (12) with an appropriate function w(x, y) in the nine-dimensional manifold determined by expressions ( 22) and ( 25) from section 3.1.In contrast, for super-asymmetric states such a representation is not possible.In this case, the corresponding function w(x, y) is still given by ( 22) but with fully complex coefficients, except for those in the pinning conditions (24).
For every static two-core spiral chimera shown in figure 8 we performed continuation of its stability boundaries.As a result, we identified five partly overlapping stability regions in the (A, B) plane shown in figure 10.Each of these regions is bounded by two fold bifurcation curves (solid lines) and one pitchfork bifurcation curve (dash-dotted line).3) with the same coupling function (10) and N = 512.

Moving spiral wave chimeras
Moving spiral wave chimeras are observed in the non-variational case, i.e. for α = 0. To continue such states in parameter space we used a numerical method based on the self-consistency equation from section 3.2 together with a non-iterative algorithm for calculating periodic orbits in the Ott-Antonsen manifold [56].To describe typical bifurcation diagrams as the phase lag α changes, we consider the case B = 3. Figure 11(a) shows a diagram analogous to figure 8 that indicates that for A = 2 and α = 0 there are three coexisting motionless Z 2 -symmetric two-core spiral wave chimeras (pink diamonds).Only one, the one with an intermediate value of R, is stable.When α increases each of these spiral patterns persists as a moving spiral wave chimera and its stability remains unchanged for small α (pink curves in figure 11(b)).For larger α the dependence of the speed s on α becomes nonlinear (figure 11(b)) and the bottom and intermediate branches annihilate in a fold point at α ≈ 0.74 while the upper, s > 0 branch continues beyond this point.
If we take A = 2.1 the upper, s > 0 branch continues to exist (not shown) while a loop composed of the two s < 0 branches detaches from α = 0 (black curves) but stable solutions continue to exist on the upper portion of the resulting isola.This loop shrinks rapidly with increasing A and disappears by A = 2.2, thereby eliminating stable two-core spiral chimera states.
A qualitatively different scenario occurs if we decrease the parameter A. For A = 1.9 the middle branch now corresponds to s > 0 forming part of a looped branch that grows in size and intersects the horizontal axis at some α * ≈ 0.645 (blue curves).Thus the spiral wave chimeras travel in opposite directions for α < α * and α > α * .
In contrast, for A = 1.8 and α = 0 the remaining two-core spiral wave chimera (figure 11(a)) results in moving states with s < 0 (red curve in figure 11(b)) with the two upper branches detached from α = 0 forming an S-shaped curve whose intermediate section corresponds to numerically stable two-core spiral wave chimeras.Snapshots of the solutions along this branch (figure 12) show that with increasing |s| the incoherent regions around the phase defects develop internal structure in the form of crescent-shaped filaments.Along the bottom part of the branch, the development of each new filament corresponds to an S-shaped fold in the solution branch in the (α, s) plane reminiscent of slanted snaking of spatially localized states [48][49][50][51][52].However, in other cases, especially when the number of filaments is large (more than 10 or so) or s is close to zero, the speed s evolves monotonically with α in a manner reminiscent of smooth snaking [51].Both behaviors are characteristic of nonlocal systems.
Note that the shape of the solution branch for A = 1.8 resembles qualitatively the behavior of the solution curve in figure 3 calculated for the top-hat coupling function.Moreover, the changes in the profiles a(x, y) along the curve are also reminiscent of those shown in figure 3.

Static patterns for equation (6) with α = 0
In the case α = 0, equation ( 6) is variational and its long-term dynamics therefore correspond to static patterns or equilibria of the form z (x, y, t) = a (x, y) . ( All two-core spiral wave solutions of equation (6) shown in figure 8 have this form.Below we describe the mathematical methods used to carry out continuation and stability analysis of these states.For this, we adapt the techniques from [36,45].

Self-consistency equation
Inserting ansatz (15) into equation ( 6) we obtain where Solving equation ( 16) for a and choosing the square root branch that ensures the inequality |a| ⩽ 1, we obtain Expressions ( 17) and ( 18) agrees with one another iff the function w(x, y) satisfies the self-consistency equation where

Reduced self-consistency equation
Let us consider equation (19) in the case of the trigonometric coupling function (10).First, we define an inner product on the space C([−π, π]; C), Then, it is easy to see that for the trigonometric coupling function (10) we have where  Owing to the finite-rank nature of the resulting integral operator G every solution of the self-consistency equation ( 19) can be written in the form w (x, y) = 9 k=1 ŵk ψ k (x, y) (22) with complex coefficients ŵk .Inserting (22) into equation ( 19) and equating similar terms on the left and right sides, we obtain the finite-dimensional system The system (23) inherits the three continuous symmetries of equation (19).Therefore a unique solution requires that we impose three pinning conditions, for example Note that owing to the first pinning condition the value of ŵ1 always coincides with the global order parameter R defined in equation (13).
It turns out that the vast majority of stable static solutions of equation ( 6) with the coupling function (10) and in particular almost all the patterns shown in figure 8 are represented by solutions of the system ( 23) and ( 24) on the invariant manifold determined by the identities Im ŵk = 0 for k = 1, 2, 3, 5, 9, Re ŵk = 0 for k = 4, 6, 7, 8. ( The invariance of this manifold follows from the fact that the basis functions ψ k with k = 1, 2, 3, 5, 9 have the reflection symmetry while the remaining four functions with k = 4, 6, 7, 8 satisfy the antisymmetric relation The restriction of the system ( 23), (24) to the manifold (25) almost halves the dimensionality of the system and allows a significant speed-up of its solution.
We may further reduce the dimensionality of the system ( 23), ( 24) by looking for Z 2 -symmetric solutions of equation (19).Such solutions satisfy the relation w(−x, y) = w(x, y), so we must assume ŵ7 = ŵ8 = ŵ9 = 0.The resulting six equations require only the first two pinning conditions in (24) in order to obtain a unique solution.A similar approach applies to Z2 -symmetric solutions as well.

Initial conditions for system (23)
To apply a Newton solver to system (23) we need to have a good initial guess.In addition, it is desirable that this guess satisfies the pinning conditions (24).We perform this task as follows.
We run the numerical simulations for the oscillator system (3) with the coupling function (10) until it approaches a stationary state of interest.Using the last snapshot of the phases θ jk , we calculate a discrete analog of formula (17), where the complex exponent e iθmn appears instead of z(−π + 2π m/N, −π + 2π n/N).Next, using the discrete Fourier transform, we calculate the necessary Fourier coefficients ŵk .These coefficients do not, in general, satisfy the pinning conditions (24).To impose the pinning condition on ŵ1 , we apply a transformation Next, we perform a transformation that ensures the pinning condition for ŵ4 ŵ3 → ŵ3 cos y 0 − ŵ4 sin y 0 , ŵ4 → ŵ3 sin y 0 + ŵ4 cos y 0 , ŵ5 → ŵ5 cos y 0 − ŵ6 sin y 0 , ŵ6 → ŵ5 sin y 0 + ŵ6 cos y 0 , ŵ8 → ŵ8 cos y 0 − ŵ9 sin y 0 , ŵ9 → ŵ8 sin y 0 + ŵ9 cos y 0 , where and from the two signs in the last formula, the one that makes the imaginary part of ŵ4 positive is chosen.Finally, we perform a third transformation, where The sign in this expression is chosen such that the imaginary part of ŵ7 positive.The resulting coefficients ŵk satisfy all the pinning conditions (24) exactly.

Stability analysis
To analyze the linear stability of the equilibria of equation ( 6), we proceed as follows.We insert the ansatz into the equation and linearize it with respect to the infinitesimal perturbation v. Thus where the subscript 0 on η 0 (x, y) indicating that we are considering a static pattern.Owing to ( 17) and ( 18), we have implying that η 0 (x, y) is real and satisfies |η 0 (x, y)| ⩾ γ.
The structure of equation (26) implies that the spectrum consists of two parts: the essential spectrum and the point spectrum σ pt consisting of isolated eigenvalues of finite multiplicity.To find the eigenvalues λ ∈ σ pt we employ the ansatz which yields a solution to equation ( 26) provided the eigenvalue λ and the components (v + , v − ) T of the eigenmode satisfy Applying the integral operator G to both sides of equation ( 28), we obtain the nonlocal eigenvalue problem for the components of the eigenmode of the local mean field V ≡ Gv.
Recall that for the coupling function (10), the integral operator G has a finite-dimensional range spanned by the functions ψ k (x, y), k = 1, . . ., 9. Therefore, every solution of equation ( 29) can be written in the form Vk ψ k (x, y) with some Vk ∈ C 2 .Inserting this ansatz into equation (29) yields a system of nonlinear equations for the nine pairs of complex coefficients Vk .Collecting these coefficients into a single vector V ∈ C 18 , we can rewrite these equations as an equivalent matrix equation, where we solve for the eigenvalue λ and the corresponding vector V ∈ C 18 .The matrix B(λ) has the structure and consists of 2 × 2 blocks: , where the ξ m are defined by (20) and we have replaced a(x, y) by expression (18).Thus the matrix B(λ) is completely determined by the solution w(x, y) of the self-consistency equation (19).
The eigenvalues λ can be found as solutions of the characteristic equation where I n denotes the n × n identity matrix.If all solutions λ = 0 of equation (30) lie in the left half-plane, Re λ < 0, then the corresponding equilibrium a(x, y) is linearly stable.In contrast, if equation ( 30) has at least one solution λ = λ * such that Re λ * > 0, then the equilibrium a(x, y) is unstable.

Stability of the completely incoherent state
The completely incoherent state corresponds to the zero solution w(x, y) = 0 of equation ( 19) and hence to the zero solution z(x, y, t) = 0 of equation (6).It is present for all values A, B ∈ R. In this case, η 0 (x, y) = γ from equation (27), and therefore Owing to the mutual orthogonality of the trigonometric functions ψ m and the relations (21), we easily see that equation (30) factorizes into three equations and which determine three eigenvalues each of double multiplicity.Consequently, the completely incoherent state z(x, y, t) = 0 is linearly stable if γ ⩾ max(1/2, A/4, B/8) and is unstable otherwise.

Stability of the partially coherent uniform state
This state corresponds to nonzero constant solutions of equations ( 19) and ( 6) of the form w (x, y) = p.
Inserting this ansatz into equation (19) we obtain and η 0 (x, y) = γ 2 + p 2 = 1 − γ; see equation (27).Moreover, Due to the mutual orthogonality of the trigonometric functions ψ m and the relations ( 21), we easily see that equation (30) factorizes into nine equations (m = 1, . . ., 9), and that each factor determines two eigenvalues, For γ < 1/2 and c m ⩾ 0, the largest eigenvalue is Thus the linear stability condition for partially coherent uniform states reads If either of the above two inequalities is violated, then the corresponding partially coherent uniform state is unstable.

Moving spiral wave chimeras for equation (6) with α ̸ = 0
In this section we describe the mathematical tools used to calculate moving spiral wave solutions of equation ( 6) and to determine their stability.We focus on Z 2 -symmetric solutions of the form z (x, y, t) = a (x, y − st) e iΩt (31) with speed s, collective frequency Ω and a reflection-symmetric profile a(x, y) = a(−x, y).

Self-consistency equation
Inserting ansatz (31) into equation ( 6), reordering the terms and dividing the resulting equation by s = 0, we obtain We look for 2π-periodic solutions in both x and y satisfying the inequality |a(x, y)| ⩽ 1.Therefore, if we happen to know s, Ω and Ga, then for each fixed x ∈ [−π, π] equation ( 32) can be read as a periodic boundary value problem for the complex Riccati equation where a(x, •) is an unknown function, x is a parameter, It is known [56] that, for every ζ / ∈ iR and every w(x, •) ∈ C per ([−π, π]; C), equation ( 33) has a unique 2π-periodic solution u(y) depending on x as a parameter such that |u(y)| < 1.The corresponding solution operator is denoted by U(w, ζ).Then equation (32) with the additional condition |a(x, y)| ⩽ 1 can be recast into the equivalent form We interpret equation ( 34) as a self-consistency equation analogous to equation (19).It is to be solved for the function w(x, y) and the two scalars s and Ω.This can be done if we equip equation (34) with two pinning conditions, and recall that w(x, y) must share the reflection symmetry of the function a(x, y).
In the case of the trigonometric coupling function (10) we can assume w (x, y) = 6 k=1 ŵk ψ k (x, y) .
The basis function ψ 7 , ψ 8 and ψ 9 do not appear in the sum because they do not satisfy the symmetry relation ψ(−x, y) = ψ(x, y).Equation ( 34) is thus equivalent to the six-dimensional complex system for six complex unknowns ŵj and two real unknowns Ω and s.The balance between the number of equations and the number of unknowns is ensured by the two pinning conditions ( 35) and (36), which are equivalent to the two scalar constraints Im ŵ1 = 0 and Im ŵ4 = 0.

Solution operator U(w, ζ)
There is no explicit expression for the operator U(w, ζ).But, as shown in [56], its value can be determined by solving only four initial value problems for equation (33).This possibility follows from the fact that the Poincaré map of equation ( 33) coincides with the classical Möbius transformation.

Stability analysis
The linear stability analysis of a uniformly drifting state (31) can be performed by analogy with the static case, see section 3.1.For this, we insert the ansatz z (x, y, t) = (a (x, y − st) + v (x, y − st, t)) e iΩt into equation ( 6) and linearize the resulting equation with respect to the infinitesimal perturbation v.We thereby obtain a linear partial integro-differential equation for v, where ξ ≡ y − st is the comoving variable and η (x, ξ) ≡ γ + iΩ + e iα a (x, ξ) Ga.
The spectral problem for the eigenvalues λ ∈ C and eigenmodes (v + , v − ) T can be derived from equation (37) on using the ansatz After substitution into equation (37), we obtain where it is assumed that v + and v − are 2π-periodic with respect to x and ξ.Rigorous analysis of the properties of the spectrum defined by equation ( 38) can be performed using the approach proposed in [24, section 4].On the other hand, a naive way to calculate the spectrum numerically is to discretize equation ( 38) on a uniform grid, approximate the derivatives by finite differences and the integrals by the trapezoid rule.This procedure leads to a matrix eigenvalue problem, which can be solved by standard numerical routines.However, this method has a significant drawback.Since the eigenmodes v + (x, ξ) and v − (x, ξ) depend on two variables, each of these functions must be approximated by arrays of minimal size 100 × 100 = 10 4 .The corresponding eigenvalue problem thus involves a huge, dense matrix of size 10 4 × 10 4 .As a consequence, the computation of the eigenvalues is extremely time-consuming and we do not perform it in this paper.As an alternative stability analysis scheme, we use the following phenomenological approach.Given a solution (31) of equation ( 6), we calculate an initial condition of the corresponding system (3) using the formula We then run simulations of the system (3) for 10 4 time units.If at the end we obtain a state resembling the expected solution (31), we consider it a stable solution of equation (6).Moreover, in this case, we use the last 5000 time units to calculate the mean drift speed of the corresponding pattern, which is plotted on top of the theoretically predicted curves in figure 11(b).

Discussion and conclusion
We have shown that both stationary and moving partially synchronized states of two-dimensional arrays of nonlocally coupled nonidentical phase oscillators are well described in the continuum limit by a self-consistency condition originally derived by Ott and Antonsen [44].In particular, we have shown that we can correctly compute the speed s of these structures as a function of the phase-lag parameter α, with s solving a nonlinear nonlocal eigenvalue problem.We have also shown that the same procedure can be used to determine the stability of these states and confirmed the results using direct numerical simulations of the discrete oscillator system.Our results explain the multiplicity of the different stable states in systems of this type, as well as the various bifurcations responsible for transitions between them.
Although the work was motivated by simulations performed with a top-hat coupling function, much of our progress was based on a two-parameter truncation of the Fourier expansion of this function.The advantage of using coupling functions of this type is apparent from earlier work in [31,35,57].It turns out, for reasons that we do not fully understand, that the values of the parameters A and B typical of different top-hat coupling functions in fact capture the most interesting regimes in the (A, B) plane, regimes where stable stationary and moving chimeras are present (figures 7 and 10).
Our approach enabled us to understand how stable two-core spiral chimeras are generated for different values of the parameters (A, B) as well as of α, albeit for a single value of the width of the Lorentzian frequency distribution, γ = 0.01.We have seen that when α = 0 two-core symmetric spiral chimeras, either Z 2 -symmetric or Z2 -symmetric, are found on tertiary branches, following three successive symmetry-breaking bifurcations, the first of which generates a π-periodic state |a(x, y)| from the partially coherent state, while the second creates a 2π-periodic state |a(x, y)|, after which the incoherent cores localize in the transverse direction.For stability these two-core states require nonzero values of both coupling coefficients A and B; the optimal conditions are such that the functions 1, A cos x, A cos y and B cos x cos y have comparable L 2 -norms, i.e. for A ≈ 2 and B ≈ 4. Four-core states are created by a similar process.These states are all static.
When α > 0, these states begin to drift, resulting in a quasiperiodic state.Such two-frequency states are associated with the appearance of curved filaments in the incoherent cores.The number of such filaments increases with α, and the addition of each new filament is associated with a fold in the bifurcation diagram, resulting in states with a great many filaments in the core such as state (viii) of figure 12.The filaments of moving symmetric chimera states are curved in the direction opposite to the direction of motion.Moreover, they typically appear in the front part of the moving structure, in contrast to the one-dimensional traveling chimera states, where similar spatial oscillations emerge in the wake of the moving structure.However, the speed s of such states is in general a highly nonmonotonic function of α and may vanish for nonzero α.It is remarkable that results of this complexity can be accessed semi-analytically and that their stability properties can likewise be established by similar techniques.Direct numerical simulations of the discrete system have confirmed these results with exquisite accuracy (figure 11(b)).
In this paper, we focused on the properties of symmetric moving spiral wave chimeras with a trigonometric coupling function.We expect that with some modifications, the computational scheme described in section 3.2 can also be applied to study asymmetric spiral wave chimeras.Our results for the variational case α = 0 show how such states are related to each other.Moreover, our methods readily extend to three-dimensional arrays of phase oscillators with trigonometric coupling functions [36] as well as to spatially-coupled arrays of phase oscillators with finite response times [58].
Although most of our results are for the trigonometric coupling function, the results qualitatively reproduce many of the direct numerical simulation results obtained for the top-hat coupling function.The trigonometric coupling approximation to this and related coupling functions has the advantage that both stable and unstable spiral chimera states can be readily computed.The latter are of inestimable value for establishing the sequence of bifurcations that are required to generate the observed stable states.We expect that the scenarios we have identified for the trigonometric coupling function carry over to other, more realistic coupling functions, where their presence can be established only with the greatest difficulty.
At a practical level, our findings can be used to explain the appearance of traveling spiral wave chimera-like patterns in phase oscillator models related to the dynamics of hydrodynamically coupled cilia carpets [6,7].On the other hand, they can potentially suggest new experimental designs for experiments on Belousov-Zhabotinsky chemical oscillators [8,9].In a broader context, we may expect that qualitatively similar spiral waves and corresponding bifurcation scenarios can be found in neural field models [59,60] or models of cardiac tissue electrophysiology [61], where moving spiral waves are common.Potential systems that also support spiral wave chimeras include excitable optical systems, such as semiconductors lasers with [62,63] or without absorber-saturable medium [64] or with two distinct time delays [65].We mention, finally, two recent papers describing novel realizations of chimera states, an experimental system of degenerate optical parametric oscillators mimicking neuronal dynamics in biological systems [66] and a two-component Bose-Einstein condensate of ultracold atoms, a Hamiltonian system exhibiting spiral chimera states that may also be realizable in experiments [67].

Figure 1 .
Figure 1.Moving spiral wave chimeras (two left columns) in the discrete oscillator array (1) and their continuum limit counterparts (two right columns) determined from equation (6).The arrows indicate the drift direction.Panels (a) and (e) focus on the behavior of the phase near a single spiral core while panels (b) and (f) present snapshots of the corresponding phase θ jk (t) of the full array of oscillators, revealing the presence of a pair of incoherent spiral cores surrounded by a coherent region in the form of a rotating spiral.The pattern moves in the y direction at a constant speed.Under appropriate conditions, the incoherent cores break up into nested filamentary structures (panels (b) and (c)).Panels (c) and (g) show snapshots of arg z(x, y, t) while panels (d) and (h) show snapshots of |z(x, y, t)|, a representation that eliminates the global spiral wave frequency Ω and focuses on the fingerprint structure of their spiral cores.The continuum limit evidently provides an accurate description of the discrete system, even for a moderate number of oscillators.The supplementary material[53] provides an animation of the state in panels (c) and (g).Parameters: N = 128, γ = 0.01, σ = 0.25 with α = 0.6 (top row) and α = 0.7 (bottom row).

Figure 2 .
Figure 2. Two examples of coupling functions considered in this paper.(a) Top-hat function with σ = 0.4.(b) The trigonometric function (10) with A = 1.63 and B = 2.62, corresponding to the leading order Fourier approximation of the top-hat function with σ = 0.4.

Figure 3
identifies four key regions in this diagram with different behavior.Enlargements of regions I and II are included in the figure together with solution snapshots at the locations labeled in the bifurcation diagram.These depict |z| in the left half of each panel and arg z in the right half and show that the sequence of folds in regions I and II is associated with an increasing number of filamentary structures in the spiral

Figure 3 .
Figure 3.The speed s of symmetric spiral wave chimeras moving in the y direction versus the phase lag parameter α in equation (6) for the top-hat coupling function with σ = 0.4 and γ = 0.01.Four regions of distinct behavior are indicated, two of which are enlarged in the figure and accompanied by snapshots of the solutions at the locations indicated.The snapshots show the whole solution, with |z| in the left half of each panel and arg z in the right half.The supplementary material [53] shows an animation indicating how the solution profile changes along the solution branch in the top panel starting from α = 0 and confirms that the folds in regions I and II are associated with the successive addition of new filamentary structures in the spiral cores in a behavior that resembles slanted snaking.Solid (dashed) lines indicate stability (instability).

Figure 4 .
Figure 4. Detail of region III of figure 3 showing the coalescence of two spiral cores and the creation of an extended filament.The bottom panels present the cross section of |z| at y = 0.

Figure 5 .
Figure 5. Detail of region IV of figure 3 showing the onset of motion for two-core spiral chimeras as α increases from zero and the subsequent development of filamentary structures in the spiral cores at larger speeds s.The black arrows indicate the direction of propagation.

Figure 6 .
Figure 6.Detail of region I of figure 3 for different values of the coupling radius σ.Increasing σ displaces the branch to larger α resulting in larger speed s and a larger spiral core.

Figure 7 .
Figure 7.The parameters A and B in equation (10) corresponding to the leading order Fourier approximation to the top-hat coupling function (5) with different σ.

Figure 8 .
Figure 8.The B = 0 panel (top left) shows the global order parameter R for different static solutions of equation (6) with the trigonometric coupling function (10) and α = 0 as a function of the parameter A. The four other panels with B > 0 show only two-core spiral chimeras, namely Z2-symmetric chimera states (orange curves), Z2-symmetric chimera states (blue curves), asymmetric chimera states (purple curves), and super-asymmetric chimera states (light green curve).Stable (unstable) states are shown in solid (dashed) lines.Snapshots at the top right and at the bottom show sample stable states in terms of arg z and |z|.

Figure 8 ,
top panel, shows the global order parameter R for static states computed from equation(6) for different values of the parameter A in equation (10) when B = 0, α = 0. We observe:(i) B = 0: the partially coherent uniform state R = √ 1 − 2γ (gray line, state (b)) becomes unstable through a pitchfork of revolution at A = 2.This bifurcation gives rise to fundamental Z 2 -symmetric states

Figure 9 .
Figure 9. Static patterns z = a(x, y) from equation (6) for different values of the parameter A in equation (10) when B = 0, α = 0, shown in terms of the quantity √ R 2 + S 2 with R and S defined in equations (13) and (14) for comparison with figure8, top panel.Corresponding branches are indicated by like colors, while stable (unstable) states are shown in solid (dashed) lines.The yellow branch shows the partially coherent splay state, while the blue line shows the D4-symmetric four-core spiral chimera.Snapshots at the right show sample states in terms of |z| and arg z.

Figure 10 .
Figure 10.The (A, B) parameter plane (panel (e)) showing regions of stable (shaded) static spiral wave chimeras from the Ott-Antonsen equation (6) with α = 0 together with solution snapshots corresponding to the color-coded diamond symbols, showing the phase θ and modulus |z| of the complex order parameter z.Filled circles indicate the parameters of stable spiral wave chimeras found in direct numerical simulations of equation (3) with N = 512.The solid and dash-dotted stability boundaries correspond to fold and pitchfork bifurcations, respectively.(a) Z2-symmetric chimera, (b) Z2-symmetric chimera, (c) Z2-asymmetric chimera, (d) Z2-symmetric chimera and (f) Z2-asymmetric (super-asymmetric) chimera.Profiles (b) and (d) correspond to different stability intervals on the same solution branch (see blue branch in figure 8).

Figure 11 .
Figure 11.(a) The global order parameter R of Z2-symmetric two-core spiral wave chimeras versus the parameter A of the trigonometric coupling function (10) when B = 3 and α = 0. (b) The speed s of symmetric spiral wave chimeras moving in the y direction versus the phase lag parameter α computed from equation (6) with the trigonometric coupling function (10) with four different values A (color-coded) and B = 3. Open circles indicate the speed of stable spiral wave chimeras found in direct numerical simulations of equation (3) with the same coupling function (10) and N = 512.

Figure 12 .
Figure 12.The solution branch for A = 1.8, B = 3 shown in red in figure 11(b) together with the corresponding fingerprint patterns at the locations labeled in the main panel.