Travelling chimera states in systems of phase oscillators with asymmetric nonlocal coupling

We study travelling chimera states in a ring of nonlocally coupled heterogeneous (with Lorentzian distribution of natural frequencies) phase oscillators. These states are coherence-incoherence patterns moving in the lateral direction because of the broken reflection symmetry of the coupling topology. To explain the results of direct numerical simulations we consider the continuum limit of the system. In this case travelling chimera states correspond to smooth travelling wave solutions of some integro-differential equation, called the Ott–Antonsen equation, which describes the long time coarse-grained dynamics of the oscillators. Using the Lyapunov–Schmidt reduction technique we suggest a numerical approach for the continuation of these travelling waves. Moreover, we perform their linear stability analysis and show that travelling chimera states can lose their stability via fold and Hopf bifurcations. Some of the Hopf bifurcations turn out to be supercritical resulting in the observation of modulated travelling chimera states.


Introduction and main results
Many living systems, including neurons, cardiac pacemaker cells, and fireflies are capable to produce rhythmic outputs and thus behave as self-sustained oscillators [1]. If such systems are coupled together, their rhythms start to interact with each other leading to the appearance of spatiotemporal patterns with different degree of synchrony between the interacting units. These patterns of synchrony often are relevant to certain physiological states of an organism or perform specific functions in the population evolution. For example, it is known that some neurological diseases (e.g. Parkinson's disease or epileptic seizures) are characterized by pathologically strong synchronization of neural activity [2,3]. On the other hand, it is also known that the information in the mammalian's working memory is encoded by spatially localized patches of neural activity, called bump states [4]. These and many other synchronization phenomena found in physics, chemistry and biology [1,[5][6][7] were the motivation for the development of mathematical theory dealing with synchronization transitions on networks [8,9] and pattern formation in oscillatory media [10][11][12][13].
One of the most spectacular phenomena studied intensively for oscillatory networks in recent time, are so-called chimera states, or patterns of coexisting coherence and incoherence [14,15]. These patterns are interesting, because they appear as a result of spontaneous symmetry breaking in homogeneous and hence highly symmetric networks. On the other hand, chimera states have many remarkable dynamical features, e.g. random wandering in space or sudden collapse, analogous to those of more complicated real-world phenomena such as bump states in networks of spiking neurons [16] or epileptic seizures [17]. For realistic oscillators, e.g. oscillators with a two-or higher-dimensional individual dynamics, chimera states usually are investigated by means of direct numerical simulations, see [11,13] and references therein, with application of different statistical measures [18,19]. Much deeper analysis, including rigorous stability analysis and continuation techniques, is possible if the high-dimensional oscillator's dynamics can be reduced to a one-dimensional phase model [6,[20][21][22]. The resulting phase model may be a very rough approximation of the original model, but in many cases it gives a good qualitative insight into the mechanism of the underlying synchronization phenomenon. In particular, in the continuum limit of infinitely many phase oscillators many chimera states have a relatively simple mathematical representation allowing their detailed bifurcation analysis [12,23].
In this paper we consider the prototype model of chimera states suggested in [24]. This model describes the dynamics of a heterogeneous network of nonlocally coupled phase oscillators. The state of each oscillator is represented by its phase θ k (t) ∈ R, which evolves according to the equation The natural frequencies ω k are randomly and independently drawn from the Lorentzian distribution g γ (ω) = 1 π γ ω 2 + γ 2 with the width γ > 0. The interaction between phases is described by the Kuramoto-Sakaguchi sinusoidal function with the phase lag parameter α ∈ (−π/2, π/2). The coupling func tion G(x) is supposed to be a non-constant piecewise continuous 2π-periodic function. This implies that the coupling strength between the kth and the j th oscillators depends nontrivially on the distance between them and the system is periodic with respect to the oscillator indices.
In [24] it was shown that for a symmetric coupling function G, i.e. for a function G such that G(−x) = G(x), the model (1) supports stationary chimera states, figure 1(a). The position of such a chimera state is pre-determined by initial condition and remains fixed in time.
If the reflection symmetry of the coupling function is broken, then the chimera state starts to drift with a constant speed [25,26] and thus becomes a travelling chimera state, figure 1(b). Importantly, the latter state is the result of forced symmetry breaking and must be distinguished from travelling chimera states arising from spontaneous symmetry breaking [27][28][29]. Apart from the above examples, travelling chimera states were also found in more complicated models with two-and three-dimensional oscillator dynamics [30,31]. However, to the best of our knowledge, the systematic stability analysis of such states has not been carried out yet.
The coupling function G(x) in equation (1) usually is chosen to mimic a specific real-world oscillatory system. For example, the spiking neurons model of bump states in [16] involves a symmetric Mexican hat interaction kernel, which also appears in other neural field models with a nonlocal interaction [32][33][34]. In contrast, in the neural field models of direction selectivity [35][36][37] the Mexican hat kernel is usually perturbed by a suitably chosen asymmetric (non-even) function. In all these cases the corresponding 2π-periodic coupling function G can be written as a Fourier series with real coefficients A k and B k . Retaining in the sum (2) only the leading order non-constant terms, we obtain the trigonometric coupling function which can be used to carry out a qualitative analysis of model (1) for both symmetric (B = 0) and asymmetric ( B = 0) coupling function cases.
In this paper we are going to consider the model (1) with the coupling function (3) and study the behaviour of travelling chimera states for different values of the asymmetry param eter B. It turns out that for varying B the chimera states undergo a variety of qualitative transformations mediated by a non-trivial bifurcation scenario. To explain them let us define a complex-valued function where x ∈ [−π, π] is a parameter, x k = −π + 2π(k − 1)/N denotes the position of the kth oscillator and #{·} denotes the number of indices satisfying the condition in the curly brackets. By analogy with the Kuramoto order parameter [6] we call the quantity z(x, t) the local order parameter. Notice that for each time moment t the modulus |z(x, t)| characterizes the coherence of the oscillators with the positions x k ≈ x, while arg z(x, t) yields the most probable phase value of these oscillators. In particular, perfect coherence (synchronization) of oscillators with x k ≈ x corresponds to |z(x, t)| = 1, while perfect incoherence (disorder) of the oscillators is represented by |z(x, t)| = 0. The phenomena found in the numerical simulations of the model (1) and (3) can be summarized as follows.
(i) For small values of B we observe a pinning of the chimera's position, figure 2. More precisely, there exists a critical value B cr > 0 such that system (1)   region of the chimera state breaks into a bunch of tilted more synchronized filaments separated by a less coherent background, figure 1(b). (iii) In some cases, forward-backward scans with respect to B reveal a hysteretic behaviour indicating the coexistence of several travelling chimera states for the same asymmetry value, see figure 2. (iv) When the asymmetry B grows further, a travelling chimera state becomes modulated such that its lateral speed and spatial profile begin to oscillate in time, figure 1(c). (v) Moreover, for relatively large asymmetries B the coherent region of a travelling chimera state becomes 'twisted', see the snapshot in figure 1(c), and shows a similarity with the partially coherent twisted states described in [38].
Remarkably, most of the above numerical observations can be explained using the continuum limit analysis of model (1). Suppose that for N → ∞ the parameter δ in (4) tends to zero in such a way that the number of oscillators satisfying the inequality |x k − x| < δ tends to infinity. Then, as shown in [24], the dynamics of the local order parameter z(x, t) obeys the nonlinear integro-differential equation -π 0 π Oscillator, x 0 300 Oscillator, x |z| -π 0 π Oscillator, x is an integral operator of convolution type, the parameters α, γ as well as the function G(x) in formula (6) are the same as in equation (1), and z denotes the complex conjugate of z. We call equation (5) the Ott-Antonsen equation, because its derivation relies on the invariant manifold reduction technique suggested in [39,40]. Equation (5) can be discretized and solved numerically. For example, using the parameters from figure 1 we obtain solutions shown in figure 3. The solutions in columns (a) and (b), obviously, have the form where Ω ∈ R, s ∈ R and a(ξ) is a smooth 2π-periodic complex-valued function. The former solution is a stationary wave with s = 0, while the latter is a travelling wave with s = 0. Notice that the modulated travelling wave in figure 3(c) has a more complicated dynamics and therefore is not described by the ansatz (7). The main part of this paper is devoted to the analysis of travelling wave solutions (7) in equation (5). Recall that in the author's previous work [26] it was shown that for small asymmetries B, i.e. |B| << min(1, |A|), the profile of travelling waves is almost identical to the profile of the corresponding stationary wave (B = 0), and the lateral speed s grows proportionally to B. In this paper we go beyond this asymptotic analysis scheme. In sections 3 and 4 we describe a continuation and stability analysis algorithm allowing us to compute the bifurcation diagram shown in figure 4. According to this diagram the lateral speed s and the collective frequency Ω turn out to be non-monotone functions of the asymmetry parameter B. The diagram also indicates that the pinning of the chimera's position is a finite size effect that disappears in the continuum limit N → ∞. The bistability of travelling waves is naturally expected along the folded parts of the graph in figure 4(b). However, this diagram does not contain any information about the modulated travelling waves, which apparently are born at the Hopf bifurcation points separating the black and the red parts of the solution curve in figure 4(b). Figure 5 shows the spatial profiles z(x, t) of the travelling waves obtained for several points in figure 4. The travelling wave for B = 0.01 looks very similar to the stationary wave for B = 0. In contrast, for larger B the profiles of moduli |z(x, t)| become spatially modulated. When the asymmetry B grows the amplitude of this modulation first increases, while its wavelength decreases. But for B > 0.06 the modulation amplitude almost stabilizes and one observes a gradual increase of the modulation wavelength. For B = 0.085 and B = 0.097 figure 5 shows three coexisting travelling waves. In the former case two waves are stable and one unstable, while in the latter case only one wave is stable and two other unstable. The spectra of the travelling waves for B = 0.097 are shown in the insert panels.
Another remarkable transformation scenario occurs for B > 0.11, see figure 6. Notice that in all panels of figure 4 the total variation of the argument of z(x, t) for x varying from −π to π equals zero. Following the solution curve in the top right corner of figure 4 we find several points where the modulus |z(x, t)| touches zero. After each of these points the argument of z(x, t) becomes more twisted than it was before. Respectively, its total variation first jumps from zero to −2π, figure 6(a), then to −4π, figure 6(b), and finally to −6π, figure 6(c). This observation agrees with the appearance of the 'twisted' travelling chimera state shown in figure 1(c).
The paper is organized as follows. In section 2 we describe numerical simulations of the discrete oscillator system (1) and of the corresponding Ott-Antonsen equation (5). In both cases we explain how to extract the collective frequency Ω and the lateral speed s from the observed travelling chimera trajectory. In section 3 we consider a nonlinear equation expressing the dependence between the amplitude a(x), the frequency Ω and the speed s of a travelling wave (7) and the system parameters, i.e. the coupling function G(x), the phase lag α and the distribution width γ. This equation cannot be solved via the implicit function theorem because of two continuous symmetries, therefore we apply the Lyapunov-Schmidt reduction technique and derive a new system of equations suitable for the continuation of travelling wave solutions to equation (5). The section 4 deals with the stability of the obtained travelling wave solutions. There we derive the characteristic equation for the spectrum of the corresponding linearized operator and analyze it. In section 5 we summarize obtained results and outline some remaining open problems. Several auxiliary mathematical results are collected in the appendix.

Numerical simulation of system (1)
A primary chimera state was obtained in the system (1) with parameters N = 8192 and α = π/2 − 0.1. We used the symmetric trigonometric coupling function (3) with A = 0.9 and B = 0. The natural frequencies ω k were generated according to the formula ω k = γ tan ζ k , where γ = 0.01 and ζ k were random numbers from the interval (−π/2, π/2). The initial phases θ k with k N/2 were chosen randomly and independently from the interval [−π, π] while all other initial phases were set to zero.
To produce figure 2 we used the dynamical continuation approach where the asymmetry parameter B first varied with the step ∆B = 0.001 from zero to 0.13 and then in the backward direction. For each value B we analyzed a chimera trajectory of the length T = 10 4 discarding the preceding transient of the same length.
To determine the collective frequency Ω and the lateral speed s of an observed travelling chimera state we used the following approach. For each time moment t we calculated the local effective forces and the global order parameter Taking into account that the points (x k , |W k (t)|) tend to align along a sinusoidal-like curve we calculated Fourier coefficients and used them to extract the leading order non-constant spatial Fourier mode Then the value y(t) was identified with the instantaneous position of the chimera state and the lateral speed s was calculated by the formula Figure 6. Development of phase twists in travelling wave solutions to equation (5). The snapshots correspond to three crosses in figure 4.
We used the other formula to calculate the collective frequency Ω.

Numerical simulation of equation (5)
To integrate numerically equation (5) we discretized it on a uniform spatial grid consisting of 8192 points. All relevant parameters were chosen as in section 2.1. A stationary wave for the symmetric (B = 0) trigonometric coupling function (3) was obtained using the initial condition Then, stable travelling wave solutions to equation (5) were obtained via dynamical continuation for B = 0. For every such solution z(x, t) the lateral speed s and the collective frequency Ω were calculated using formulas (8) and (9). However, in (8) instead of y(t) we used the position of the maximum of the modulus profile |z(x, t)|, while in (9) we used the different definition of the global order parameter

Continuation of travelling waves
Inserting ansatz (7) into equation (5) we obtain the self-consistency equation connecting the travelling wave parameters a(x), Ω and s and the system (1) parameters γ, α and G(x). In the following we suppose that the parameters γ and α are fixed and the coupling function G(x) is chosen in the form (3) where the parameter A is fixed too. Thus, the asymmetry parameter B remains the only varying system parameter. Notice that if a triplet (a 0 (x), Ω 0 , s 0 ) solves equation (10) for some B 0 , then the triplet (a 0 (−x), Ω 0 , −s 0 ) solves equation (10) for B = B 0 , therefore without loss of generality we may consider non-negative values of B only.
In this section we are going to outline the continuation algorithm allowing one to reveal the dependence of the wave parameters a(x), Ω and s on the parameter B. The difficulty of this task originates from the fact that equation (10) has two continuous symmetries. Namely, it is invariant with respect to the phase-shift a(x) → e iφ a(x) with arbitrary φ ∈ R as well as with respect to the lateral shift a(x) → a(x − x 0 ) with arbitrary x 0 ∈ R. Therefore equation (10) cannot be solved by means of the classical implicit function theorem. However, the solution can be found using the Lyapunov-Schmidt reduction technique described in section 3.2.

Preliminaries
Throughout the paper we use the following convention: (1) normal letters (Ω, s, B etc) denote scalars, (2) bold letters (L, M, v etc) denote vectors and matrices, in particular, for any positive integer n, the symbol I n denotes the n-dimensional square unit matrix, (3) calligraphic letters (G , H, K etc) denote operators in vector spaces.
The symbol Θ(x) is used to denote the Heaviside step function such that Θ(x) = 0 for x < 0 and Θ(x) = 1 for x 0.
The space of all 2π-periodic continuous functions is denoted with the symbol C per ([−π, π]; C), while the space of all 2π-periodic C 1 -smooth functions is denoted with the symbol C 1 per ([−π, π]; C). For every pair w 1 , w 2 ∈ C per ([−π, π]; C) the inner product is defined by the formula This definition, obviously, implies that for every constant c ∈ C we have Given a bounded linear operator A acting on C per ([−π, π]; C) we call its adjoint operator the operator A † satisfying the identity According to this definition, for every pair of bounded linear operators A 1 and A 2 we have (13) with the piecewise-continuous kernel K(x, y), then the adjoint operator K † is also an integral operator of the form (13) but with the adjoint kernel K(y, x) . Notice that this rule yields Let us define the vector consisting of eight functions ψ k (x) (the terms u 1 (x) and u 2 (x) will be defined later by formulas (23)) . Then in the case of the trigonometric coupling function (3) for These formulas will be useful in the following analysis.

The Lyapunov-Schmidt reduction
Suppose that a triplet (a 0 , Ω 0 , s 0 ) satisfies equation (10) for some asymmetry parameter B 0 and we want to carry out the continuation of this solution for B ≈ B 0 . Let L denote the derivative of the nonlinear operator F with respect to its first argument It is easy to verify that the homogeneous equation La = 0 has two non-trivial solutions a = ia 0 and a = ∂ x a 0 . The former solution corresponds to the phase-shift symmetry of equation (10), while the latter solution corresponds to the translation symmetry of equation (10). Because of these solutions the operator L is not invertible and equation (10) cannot be solved using the classical implicit function theorem, therefore we need to use the Lyapunov-Schmidt reduction technique. To describe it let us make two assumptions: (A 1 ) Suppose that the functions ia 0 and ∂ x a 0 are linearly independent and i.e. the equation La = 0 has no other linearly independent solutions except of ia 0 and ∂ x a 0 .
Remark 1. According to the geometric interpretation of the inner product ·, · the functions ia 0 and ∂ x a 0 are linearly independent if and only if Assumption (A 2 ) and proposition 1 ensure the existence of the inverse operator Obviously, the product operator is a Fredholm operator on C per ([−π, π]; C) and ker K 2 L = ker L. Therefore, because of assumption (A 1 ), there exists a pair of linearly independent functions u 1 and u 2 such that ker K 2 L = span (u 1 , u 2 ). According to the Fredholm alternative, the cokernel of the operator K 2 L coincides with the kernel of the adjoint operator (K 2 L) † . Moreover, the dimension of ker (K 2 L) † is equal to the dimension of ker K 2 L. Along with the assumption (A 1 ) this implies that there exists a pair of linearly independent functions v 1 and v 2 such that ker (K 2 L) † = span (v 1 , v 2 ). Now the Lyapunov-Schmidt method tells that in contrast to the non-invertible operator K 2 L the modified operator is invertible and hence it is an isomorphism from C per ([−π, π]; C) onto itself.
Using this observation we define the new operator and consider the system Obviously, every solution to system (21) yields a solution to equation (10). Moreover, because of the identities ia 0 , a 0 = ∂ x a 0 , a 0 = 0 the triplet (a 0 , Ω 0 , s 0 ) is the solution to the system (21) with B = B 0 . The system (21) has a big advantage compared to equation (10) because the derivative operator is an isomorphism from C per ([−π, π]; C) onto C 1 per ([−π, π]; C). This fact will help us to formulate the continuation algorithm for travelling waves of the form (7) in section 3.5.
The next two sections play an auxiliary role. There we explain how in the case of the trigonometric coupling function (3) one can choose the functions u 1 , u 2 , v 1 and v 2 (in section 3.3) and compute the value of the inverse operator M −1 (in section 3.4).

Construction of the basis functions u 1 , u 2 , v 1 and v 2
Assumption (A 1 ) allows us to write the orthonormal basis of ker K 1 L = ker L explicitly. Since ker L = span (ia 0 , ∂ x a 0 ), following the Gram-Schmidt algorithm we obtain where the real constants C 1 and C 2 are chosen according to the normalization conditions u 1 , u 1 = u 2 , u 2 = 1.
In order to construct the orthonormal basis of the cokernel of the operator K 2 L we need to consider the adjoint equation (K 2 L) † v = 0 and find its two linearly independent solutions. This problem is solved in propositions 1 and 2 Proposition 1. The adjoint of the operator K 2 L reads where G † is given by formula (14), and Proof. Obviously, we have K † 2 = K 1 (the adjoint of an integral operator). Moreover, for Then, using the formulas (11), (12) and (14) we arrive at the identity that justifies the formula (24). ■

Proposition 2. Let
be a square matrix of the form where ψ j (x) are defined by formula (15), and K 1 denotes the integral operator (25). Then there exist two linearly independent vectors (w k 1 , w k 2 , . . . , w k 6 ) T ∈ C 6 , k = 1, 2, span ning the null space of the matrix L − I 6 . Moreover, ker (

Proof. Suppose that v is a solution to the homogeneous equation
where Inserting expression (27) into each of the equation (28) we obtain a linear six-dimensional system of the form where w = (w 1 , w 2 , . . . , w 6 ) T ∈ C 6 . Thus, there exists a one-to-one correspondence between the non-trivial solutions to equation (K 2 L) † v = 0 and the non-trivial solutions to system (29). Since ker (K 2 L) † is two-dimensional, the null space of the matrix L − I 6 is two-dimensional too, hence the homogeneous equation (L − I 6 )w = 0 must have two linearly independent solutions w 1 and w 2 . Inserting them into formula (27) we obtain two linearly independent functions ṽ k (x) which can be used to construct the orthonormal basis of ker (K 2 L) † following the Gram-Schmidt orthogonalization process. ■ Remark 3. If the entries of the matrix L are known only approximately, then one has to apply a special approximate procedure to identify the null space of the matrix L − I 6 . For this, one uses the singular value decomposition where U and V are real unitary matrices and S is a diagonal matrix with real non-negative elements. Then one chooses the two smallest diagonal elements of the matrix S and uses the corresponding columns of the matrix V as the approximate values of the vectors w 1 and w 2 .

Computation of the inverse operator M −1
Recall that M denotes the derivative operator ∂ a H(a 0 , Ω 0 , s 0 , B 0 ) relevant to the system (21).
In the following proposition we show how one computes the value of the inverse operator M −1 .

Proposition 3. Let
where ψ j (x) are defined by formula (15), and K 2 denotes the integral operator (19). Suppose that the matrix M − I 8 is nonsingular, then for every f ∈ C per ([−π, π]; C) the equation Mu = f has a unique solution given by the formula Thus, the operator M is an isomorphism from C 1 per ([−π, π]; C) onto C per ([−π, π]; C).

Proof. The equation Mu = f can be written as follows
Then using formulas (16) and (17) we obtain Acting on both sides of equation (30) with the inverse operator K 2 = (−s 0 ∂ x + η 0 ) −1 and then using the inner product operation ψ j , · with different j = 1, . . . , 8 we obtain where w j = ψ j , u . The latter is a system of linear equations, which can be solved because of the nonsingularity assumption about the matrix M − I 8 . Inserting the found w j instead of the inner products ψ j , u into equation (30)

Continuation algorithm
Let us show that system (21) does help to solve equation (10) H(a, Ω, s, B) = 0 we determine the function a =ã(Ω, s, B) such that ã(Ω 0 , s 0 , B 0 ) = a 0 . Moreover, using the formulas (10) and (20) we compute Therefore the partial derivatives of the function ã are given by the formulas Inserting the above function ã(Ω, s, B) into the second and the third equations of the system (21) we obtain This system can be solved with respect to the variables Ω and B if its Jacobian matrix is non-singular, then it delivers two functions Ω =Ω(s) and B =B(s) such that Ω (s 0 ) = Ω 0 and B (s 0 ) = B 0 . Along with the superposition formula a =ã(Ω(s), s,B(s)) this yields the s-parameterized solution to equation (10). The above consideration proves the existence of the solution to equation (10) for s ≈ s 0 , but it does not tell how to find this solution. In practice, this can be done using Newton's iterations. Recall that system (21) can be abbreviated as an operator equation of the form E s (a, Ω, B) = 0 in the vector space (a, Ω, B) ∈ C per ([−π, π]; C) × R 2 with s ∈ R being the parameter. Then Newton's iteration formula reads (37) where E s0 (a 0 , Ω 0 , B 0 ) −1 is the inverse of the derivative operator E s0 (a 0 , Ω 0 , B 0 ). To define this inverse operator explicitly we need to solve the system corresponding to the derivative E s0 (a 0 , Ω 0 , B 0 ). (38), which can be computed by the form ulas

Proposition 4. Suppose that the operator M is invertible and the matrix D is non singular, then for every
and Proof. If the operator M is invertible then the first equation of the system (38) can be written as follows see formulas (32)- (34). Inserting this expression into the second and the third equations of the system (38) we obtain a two-dimensional system of the form For non-singular matrix D this system is solved by the formula (39). Inserting the found values of Ω 1 and B 1 into (41) we obtain (40). ■

Remark 5.
The expression E s0 (a 0 , Ω 0 , B 0 ) −1 E s (a n , Ω n , B n ) in the formula (37) can be computed using the formulas (39) and (40) with f 0 = H(a n , Ω n , s, B n ), f 1 = u 1 , a n and f 2 = u 2 , a n . Now, we can summarize the developed continuation algorithm for travelling waves. If the triplet (a 0 , Ω 0 , B 0 ) solves equation (10) for s = s 0 , then the triplet (a(s), Ω(s), B(s)) satisfying equation (10) for some other s ≈ s 0 can be found using the following steps.
Step 1. Using remark 1 we check that the functions ia 0 and ∂ x a 0 are linearly independent.
Then we compute the basis functions u 1 and u 2 from the formulas (23).
Step 2. We check that the inequality (18) is satisfied. Then we use proposition 2 and remark 3 to find the basis functions v 1 and v 2 .
Step 5. Using the triplet (a 0 , Ω 0 , B 0 ) as an initial condition we start Newton's iterations = Ω n B n − D −1 u 1 , a n − b n u 2 , a n − b n , H(a n , Ω n , s, B n ). If the new value s lies close enough to s 0 , then the Banach fixed point theorem guarantees that these iterations are convergent. In practice, we stop them when the desired computational accuracy ε is achieved (in our simulations ε = 10 −6 ), i.e. when a n+1 − a n ∞ ε, |Ω n+1 − Ω n | ε and |B n+1 − B n | ε.

Remark 6. If the matrix D in
Step 4 is singular, this can be the indication of a fold bifurcation. Then the above algorithm does not work, but there may be other ways to find the solution to equation (10). For example, if the matrix is non-singular, then the system (35) can be solved with respect to the variables Ω and s. In this case we can seek the solution to equation (10) in the form (a(B), Ω(B), s(B)) where B is the new independent variable. Respectively, the Newton's iteration formula from Step 5 must be replaced with where b n = M −1 H(a n , Ω n , s n , B).

Remark 7.
In our numerical simulations, every function f ∈ C per ([−π, π]; C) was replaced with its discretization on the uniform grid with N = 8192 points. Respectively, we replaced all integrals over the interval [−π, π] according to the rule

Stability of travelling waves
Suppose that we know a travelling wave solution to equation (5), which can be written in the form (7). To analyze its stability we insert the ansatz into equation (5). Linearizing the resulting equation with respect to the small perturbation v we obtain where ξ = x − st is the wave variable and η(ξ) = γ + iΩ + e iα a(ξ)Ga. Notice that the function v is smooth with respect to its both arguments and 2π-periodic with respect to the variable ξ.
Below we analyze the stability of the zero solution to equation (42). For this we consider perturbations of the form Inserting this ansatz into equation (42) and equating the terms proportional to e λt and e λt separately, we obtain Using the notation v = (v + , v − ) T system (43) can be written in the operator form where is a two-component integral operator in C per ([−π, π]; C 2 ), and I is the identity-operator.
We are going to prove the following statements regarding equation (43).
(i) In proposition 5 we will show that for every λ ∈ C the operator D − λI + N is a Fredholm operator of index zero from C 1 per ([−π, π]; C 2 ) into C per ([−π, π]; C 2 ), therefore all nontrivial solutions of the spectral problem (43) correspond to isolated eigenvalues λ of finite multiplicity. (ii) In proposition 6 we will show that all eigenvalues of the problem (43) lie in a specific region of the complex plane. (iii) Finally, in proposition 7 we will show that in the case of the trigonometric coupling function (3) all eigenvalues of the problem (43) can be found solving an explicitly known characteristic equation.
Proof. Because of remark 2, there exists λ 0 ∈ C such that D − λ 0 I is an isomorphism from C 1 per ([−π, π]; C 2 ) onto C per ([−π, π]; C 2 ). Choosing this λ 0 we obtain The both operators (λ 0 − λ)I and N are compact because of the compact imbedding of the space of smooth functions C 1 per ([−π, π]; C 2 ) into the space of continuous functions C per ([−π, π]; C 2 ). Hence, the operator D − λI + N can be decomposed into the sum of an invertible operator and a compact operator. This ends the proof. ■ In order to prove the next proposition we need to require that the functions a and η are not only continuous but also smooth. An additional smoothening requirement is also imposed on the operator G .

Proof. Remark 2 implies that the operators
where c 1 > 0 is independent of λ, since in the formula (A.8) the constant c 0 depends on the difference ν − ν 0 only. Let us assume that the norm of Then for Re(λ + η m ) = 0 the operator D − λI is invertible, and where the constant c 2 > 0 may differ from the constant c 1 , but still it does not depend on λ.
Because of the assumptions made about a, η and G the operator N is a bounded linear operator on C per ([−π, π]; C 2 ), therefore there exists a constant c 3 > 0 such that In this case, due to [41,theorem IV.1.16], the operator D − λI + N is invertible and hence the equation (44) has no non-trivial solutions. Let us consider the first equation of the two-component system (44) written in the form Because of the smoothness assumptions imposed on the functions a and η as well as on the operator G there exists a constant c 4 > 0 such that Suppose that λ ∈ C satisfies two inequalities |λ| > 2 η ∞ and Re(λ + η m ) = 0, then we can apply remark 3 to equation (46) and obtain where c > 0 is the constant from the formula (A.10), which in our case is independent of λ.
Similarly, we consider the second equation of the system (44) and obtain an analogous inequality for v − ∞ . Altogether this implies that there exist two constants c 5 , c 6 > 0 independent of λ such that every solution v to equation (44) provided |λ| is large enough and Re(λ + η m ) = 0. This means that equation (44) does not have non-trivial solutions if Thus, all eigenvalues of the problem (44) lie in the region | Re(λ + η m )| c 7 /|λ| with some c 7 > 0. Assuming c * = c 3 and c * * = c 7 we obtain the above formulated spectral region estimate. ■ Next, we consider the case of the trigonometric coupling function (3) and derive the characteristic equation for eigenvalues λ determined by equation (43). For this we rewrite equation (43) in the form If (v + , v − ) T is a solution to equation (47), then there exist numbers v k + ,v k − ∈ C, k = 1, 2, 3, such that where (ϕ 1 (ξ), ϕ 2 (ξ), ϕ 3 (ξ)) T = (1, cos ξ, sin ξ) T are functions spanning the range of the operator G with the trigonometric coupling function (3). Inserting ansatz (48) into equation (47) we obtain Let us define two functions and and two integral operators Considering equation (49) as a two-dimensional ODE system with the initial condition v ± (−π) =v 0 ± we write its general solution in the form v This formula yields a 2π-periodic function if and only if v 0 where Along with equation (53) the coefficients v k ± in formulas (48) should also satisfy some self-consistency relations following from the definition of the operator G . Indeed, for every v ∈ C per ([−π, π]; C) the coupling function formula (3) yields Comparing this identity with ansatz (48) we find (Aϕ 3 (y) + Bϕ 2 (y))v ± (y)dy.
Inserting here v + and v − from formula (52) we obtain v 1 where for every j, k = 1, 2, 3 we denote Four equations (53) and (54)-(56) correspond to the eight-dimensional system Obviously, this system has non-trivial solutions if and only if λ ∈ C satisfies the characteristic equation Taking into account the constructive way of the derivation of equation (58), we obtain the following proposition.

Proposition 7.
In the case of the trigonometric coupling function (3) every eigenvalue of the spectral problem (43) corresponds to a zero of equation (58) and vice versa.

Remark 8.
Because of the phase-shift and translation symmetries of the Ott-Antonsen equation (5), the characteristic equation (58) has a double zero at λ = 0. If equation (58) has no other solutions λ = 0 in the right half-plane Re λ 0, then the corresponding travelling wave (7) is stable. In contrast, if equation (58) has at least one solution λ with Re λ > 0, then the corresponding travelling wave (7) is unstable. Notice that proposition 6 indicates that the operator D + N is of hyperbolic type and the linear stability principle may fail for such operators [42]. However, the spectral problem (43) is one dimensional in space (ξ ∈ R), therefore the relation between the stability of travelling wave (7) and the position of the rightmost roots of equation (58) follows from [43].
According to the definitions (50) and (51) the functions Φ + (ξ, λ) and Φ − (ξ, λ) are analytic with respect to λ in the whole complex plane, therefore the determinant det (I 8 − J(λ)) is also analytic for all λ ∈ C. This implies that the characteristic equation (58) has only isolated zeros of finite multiplicity. Moreover, the zeros have no accumulation points except of the point at infinity. Notice that proposition 6 tells that for sufficiently large |λ| the eigenvalues (if they exist) tend to the line Re λ = − Re η m , where η m is given by formula (45). Therefore if Re η m > 0, then the stability of a travelling wave is determined by finitely many eigenvalues lying in a bounded region around zero.

Conclusion
We considered a prototype model of nonlocally coupled heterogeneous phase oscillators and showed that in the case of a coupling topology with broken reflection symmetry this model can support travelling chimera states. As the coupling asymmetry grows these chimera states undergo a sequence of transformations, which can be adequately explained using the continuum limit Ott-Antonsen equation (5). For this equation we carried out a detailed analysis of travelling wave solutions of the form (7), including their numerical continuation and stability analysis. The continuation algorithm is based on the Lyapunov-Schmidt method and, in general, involves the same steps as similar continuation algorithms developed for the neural field models described by integro-differential equations [33,34]. The stability analysis of travelling waves uses standard PDE techniques [44,45], which allow us to prove the following assertions: (i) The spectrum of a travelling wave is purely discrete and is localized in a specific region of the complex plane, see propositions 5 and 6. (ii) Fold and Hopf bifurcations are the main destabilization mechanisms of these travelling waves. We expect that the bifurcation diagram in figure 4 shows typical features of the relation between the nonlocal coupling asymmetry and the corresponding travelling chimera state. In particular, it indicates the limitation of the chimera's position control suggested in [25,46]. Indeed, usually this control technique relies on a one-to-one correspondence between the chimera's lateral speed and the coupling function asymmetry. On the other hand, figure 4 shows that this requirement is satisfied for a relatively small range of parameter B only.
Notice that many formulas in sections 3 and 4 are written explicitly for the case of the trigonometric coupling function (3). However, both the continuation algorithm and the stability analysis scheme can be generalized for more complicated coupling functions G, e.g. for any truncation of the sum (2). Moreover, the numerical algorithm from section 3 can be also modified to carry out the continuation of travelling waves with respect to other parameters in equation (5), e.g. α and A. In a broader context, the approach developed in this paper allows one to study other non-stationary coherence-incoherence patterns, for example, travelling chimera states appearing due to drift instabilities [27][28][29] or travelling bumps in coupled thetaneuron models [47].
We want to emphasize that our work does not answer all possible questions related to travelling chimera states. For example, we did not explore the scaling behaviour of the pinning phenomenon in figure 2. We also did not carry out any theoretical analysis of modulated travelling waves, figure 1(c). Finally, we did not address a challenging mathematical problem concerning the continuum limit description of travelling chimera states in the case of identical oscillators [26]. All these questions remind us again how complicated can be the dynamics of such seemingly simple models as the oscillator system (1) or the Ott-Antonsen equation (5).
with a continuous 2π-periodic coefficient ν and a continuous 2π-periodic inhomogeneity f . Below we formulate sufficient conditions for the solvability of equation (A.1) in the space of smooth 2π-periodic functions and provide an explicit solution formula.