{'fraction_non_alphanumeric': 0.05734935812125059, 'fraction_numerical': 0.01809954751131222, 'mean_word_length': 4.273944498434294, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mathematical framework for breathing chimera states

About two decades ago it was discovered that systems of nonlocally coupled oscillators can exhibit unusual symmetry-breaking patterns composed of coherent and incoherent regions. Since then such patterns, called chimera states, have been the subject of intensive study but mostly in the stationary case when the coarse-grained system dynamics remains unchanged over time. Nonstationary coherence-incoherence patterns, in particular periodically breathing chimera states, were also reported, however not investigated systematically because of their complexity. In this paper we suggest a semi-analytic solution to the above problem providing a mathematical framework for the analysis of breathing chimera states in a ring of nonlocally coupled phase oscillators. Our approach relies on the consideration of an integro-differential equation describing the long-term coarse-grained dynamics of the oscillator system. For this equation we specify a class of solutions relevant to breathing chimera states. We derive a self-consistency equation for these solutions and carry out their stability analysis. We show that our approach correctly predicts macroscopic features of breathing chimera states. Moreover, we point out its potential application to other models which can be studied using the Ott-Antonsen reduction technique.


Introduction
Many living organisms, chemical and physical systems can behave as self-sustained oscillators [1]. Spiking neurons in the brain, flashing fireflies, Belousov-Zhabotinsky (BZ) reaction and swinging pendula of metronomes are just a few examples of such kind. When two or more self-sustained oscillators interact with each other, their rhythms tend to adjust in a certain order resulting in their partial or complete synchronization [2,3]. Such phenomena have been observed in many real-world systems and laboratory experiments, including Josephson junction arrays [4], populations of fireflies [5] and yeast cells [6], chemical [7] and electrochemical oscillators [8]. Moreover, synchronization turns out to underlie many physiological processes [9,10] and is, in some cases, associated with certain brain disorders, such as schizophrenia, epilepsy, Alzheimer's and Parkinson's diseases [11,12].
Mathematical modelling of synchronization processes often relies on different versions of the Kuramoto model [13,14,15,16], where the state of each oscillator is described by a single scalar variable, its phase. These models are derived as normal forms for general weakly coupled oscillator networks [17,18,19]. Moreover, the resulting equations are simplified additionally to reduce their complexity. This paper is concerned with a Kuramoto model of the form Here G(x) is a nonconstant continuous even function called coupling kernel and α ∈ (−π/2, π/2) is a phase lag parameter. System (1) describes dynamics of N identical nonlocally coupled phase oscillators θ k . Moreover, if G(x) is 2π-periodic, then the connections between oscillators have circular symmetry and system (1) is in fact a ring of coupled oscillators. Model (1) was first suggested by Kuramoto and Battogtokh in [20] and since then has been intensively studied in the context of chimera states. By chimera states one denotes specific dynamical regimes in system (1) where a part of oscillators get synchronized, while the others keep oscillating asynchronously. Such states are irrelevant to the circular symmetry of system (1), therefore their emergence seems to be counterintuitive, what explains the origin of their name [21]. Being first reported for nonlocally coupled phase oscillators, later chimera states have been found in a variety of other coupled oscillator systems, both in experiments [22,23,24,25,26,27,28,29], and in numerical simulations, see [30,31,32,33,34] and references therein. A typical example of a chimera state in the system (1) with a cosine coupling kernel is shown in Figure 1. Note that, for convenience, in the panels (a)-(c) we use the oscillator  (4). Parameters: N = 8192, A = 0.9 and α = π/2 − 0. 16. positions x k = −π + 2πk/N instead of the discrete indices k. The main difference between synchronized and asynchronous oscillators can be seen in the space-time plot (c). The former oscillators rotate at almost the same constant speed, while the latter oscillators have a more complicated dynamics. As a result, in the snapshot (a) there appear two regions: a coherent region, where the phases θ k lie on a smooth curve, and an incoherent region, where the phases are randomly distributed. Moreover, if we compute the effective frequencies Ω eff,k defined by then we obtain the arc-shaped graph (b) with a plateau corresponding to synchronized oscillators. Note that although the oscillator dynamics behind the chimera state is quite complicated and chaotic [35], in the statistical sense, this state is stationary as can be seen from the graph (d) of the global order parameter In this paper, we consider a more complicated and therefore less explored type of chimera states, called breathing chimera states, which is characterized by nonstationary macroscopic dynamics. An example of such state is shown in Figure 2. Roughly speaking, the typical features of breathing chimera states are: (i) multiple synchronized and asynchronous regions in the snapshot of θ k , (ii) multiple equidistant plateaus in the graph of the effective frequencies Ω eff,k , (iii) periodically varying (breathing) oscillator dynamics, (iv) oscillating modulus |Z N (t)| of the global order parameter. Apart from the above example breathing chimera states have been found in systems (1) with exponential [36] and top-hat [37,38] coupling kernels G(x). They were also observed in two-dimensional lattices of phase oscillators as breathing spirals with randomized cores [39,40]. Moreover, it was shown that breathing chimera states persist for slightly heterogeneous phase oscillators [41,42] and can emerge in systems of coupled limit cycle oscillators [43,44]. In spite of the number of these results the true nature of breathing chimera states still remains unclear. Even for relatively simple system (1) no mathematical methods are known to analyze their stability or predict their properties for given coupling kernel G(x) and phase lag α. This is the problem we are going to address below.
Our approach is based on the consideration of an integro-differential equation where z(x, t) is an unknown complex-valued function 2π-periodic with respect to x, z denotes the complex conjugate of z and symbol G denotes an integral operator with the coupling kernel G(x) identical to that in system (1). It is known [33,45] that in the limit of infinitely many oscillators N → ∞, Eq. (5) describes the long-term coarse-grained dynamics of system (1). If one assumes that x k = −π + 2πk/N is the physical position of the kth oscillator, then z(x, t) yields the local order parameter of the oscillators positioned around x ∈ [−π, π]. More precisely, this means where 0 < δ << 1 and #{·} denotes the number of indices k that satisfy the condition in curly brackets. If |z(x, t)| = 1, then the oscillators with x k ≈ x behave synchronously, in other words, they are coherent. In contrast, if |z(x, t)| < 1, then the oscillators with x k ≈ x behave asynchronously. This interpretation implies that only functions z(x, t) satisfying the inequality |z(x, t)| ≤ 1 are relevant to the oscillator system (1). Moreover, every chimera state is represented by a spatially structured solution z(x, t) composed of both coherent (|z| = 1) and incoherent (|z| < 1) regions.
In most of previous studies, chimera states were identified with the rotating wave solutions of Eq. (5) given by where a ∈ C per ([−π, π]; C) and Ω ∈ R. However, for breathing chimera states ansatz (7) does not work. Numerical results in [38,40] suggest that, in this case, one has to look for more complicated solutions where a(x, t) is a function 2π-periodic with respect to x and T -periodic with respect to t for some T > 0. To distinguish the cyclic frequencies Ω and ω = 2π/T in the following we call them the primary frequency and the secondary frequency, respectively. In general these frequencies are different and incommensurable, therefore the dynamics of the solution (8) appears to be quasiperiodic, thus breathing chimera states can also be called quasiperiodic chimera states.
Remark 1.1 Note that the product ansatz (8) with a T -periodic function a(x, t), in general, is not uniquely determined. Indeed, for every nonzero integer m we can rewrite it in the equivalent form z = a m (x, t)e iΩmt with a m (x, t) = a(x, t)e imωt and Ω m = Ω − mω.
To avoid this ambiguity, throughout the paper we assume that the function a(x, t) and the constant Ω in (8) are chosen so that Roughly speaking, we request that the variation of the complex argument of Y (t) remains bounded for all t ≥ 0. Importantly, the above calibration condition is well-defined only if Y (t) = 0 for all t ≥ 0, therefore we checked carefully that this requirement is satisfied for all examples of breathing chimera states shown in the paper.

Remark 1.2
The term quasiperiodic chimera state was previously used in [46] to denote some partially synchronized states in a two population model of globally coupled identical phase oscillators. In contrast to breathing chimera states considered here, the quasiperiodic chimera states from [46] are characterized by quasiperiodically (not periodically!) oscillating modulus |Z N (t)| of the global order parameter. Moreover, their mathematical description goes beyond the Ott-Antonsen theory [47] used in this paper and requires the application of a more general approach suggested by Watanabe and Strogatz [48,49]. The chimera state in Figure 2 demonstrates a simpler, periodic dynamics of |Z N (t)|, therefore, in the context of two population models, it should be compared with the breathing chimera states reported in [50].

Remark 1.3
It is easy to verify that Eq. (5) is equivariant with respect to the one-parameter Lie group of complex phase shifts z → ze iφ with φ ∈ R\2πZ. Therefore, in the mathematical literature, solutions of Eq. (5) of the form (7) and (8) are called relative equilibria and relative periodic orbits (with respect to this group action), respectively. These names also indicate that solution (7) becomes an equilibrium and solution (8) becomes a periodic orbit in the appropriate corotating frame. The existence and stability of relative equilibria of Eq. (5) were in the focus of many papers, in particular, those concerned with stationary chimera states. In contrast, relative periodic orbits of Eq. (5) were hardly considered. The main purpose of the present paper is to fill this gap.
In this paper, we demonstrate that ansatz (8) does describe breathing chimera states and can be used for their continuation as well as for their stability analysis. The paper is organized as follows. In Section 2 we consider an auxiliary complex Riccati equation where s is a real coefficient and w(t) is a continuous complex-valued function. We show that, in general, this equation has a unique stable periodic solution satisfying the inequality |u(t)| ≤ 1.
The corresponding solution operator is denoted by Its properties are discussed in Sections 2.1 and 2.2. Although the operator U is defined implicitly, it turns out that its value can be computed by solving only three initial value problems for the above complex Riccati equation.
In Section 3, we show that if Eq. (5) has a stable solution of the form (8), then its amplitude a(x, t) and its primary and secondary frequencies Ω and ω satisfy a self-consistency equation where A modified version of Eq. (9) is obtained in Section 3.1. Then in Section 3.2 we suggest a continuation algorithm allowing to compute the solution branches of Eq. (9) and thus to predict the properties of breathing chimera states. In particular, in Section 3.3 we explain how the continuum limit analogs of the global order parameter Z N (t) and of the effective frequencies Ω eff,k can be computed. Moreover, for the sake of completeness, in Section 3.4 we also explain how to extract the primary and secondary frequencies of a breathing chimera state from the corresponding trajectory of coupled oscillator system (1).
In Section 4, we carry out the linear stability analysis of a general solution (8) to Eq. (5). It relies on the consideration of the monodromy operator, which describes the evolution of small perturbations in the linearized Eq. (5). More precisely, the stability of the solution (8) is determined by the spectrum of the monodromy operator. The spectrum consists of two parts: essential and discrete spectra. The former part is known explicitly and has no influence on the stability of breathing chimera states, while the latter part is crucial for their stability but can be computed only numerically as explained in Section 4.1.
In Section 5, we illustrate the performance of the developed methods considering specific examples of breathing chimera states in system (1). In particular, using the continuation method, we find a branch of breathing chimera states starting from the chimera state in Fig. 2. Moreover, we also compute theoretical predictions of the graphs of |Z N (t)| and Ω eff,k shown in Fig. 2(b),(d). Furthermore, our analysis reveals that solutions (8) can lose their stability only via nonclassical bifurcations when one or two unstable eigenvalues emerge from the essential spectrum. Similar results were known for partially synchronized states in the classical Kuramoto model [51] and for stationary chimera states (7) [45]. But for breathing chimera states they are reported for the first time here. Finally, in Section 6 we summarize the obtained results and point out other potential applications of our methods.
Notations. Throughout this paper we use the following notations. We let C per ([−π, π]; C) denote the space of all 2π-periodic continuous complex-valued functions. A similar notation C per ([−π, π] × [0, 2π]; C) is used to denote the space of all continuous double-periodic functions on the square domain [−π, π] × [0, 2π]. Moreover, the capital calligraphic letters such as G or U are used to denote operators on appropriate Banach spaces. Finally, the symbol Θ(t) denotes the Heaviside step function such that Θ(t) = 0 for t < 0 and Θ(t) = 1 for t ≥ 0.

Periodic complex Riccati equation
Let us consider a complex Riccati equation of the form where s ∈ R and w(t) is a continuous complex-valued function. Let D = {v ∈ C : |v| < 1} denote the open unit disc of the complex plane and D = D ∪ ∂D be its closure. In this section we will show that for every (w(t), s) ∈ C per ([0, 2π]; C) × R such that |s| + max |w(x)| = 0, in general, there exists a unique stable solution to Eq. (11) lying entirely in the unit disc D. The nonlinear operator yielding this solution will be denoted by U(w(t), s).
Proof: Suppose that u(t) is a solution to equation (11), then According to Eq. (12), if |u(0)| = 1 then |u(t)| = 1 for all other t = 0, therefore the solution u(t) cannot blow up in finite time and hence it can be extended for all t ∈ R. On the other hand, Eq. (12) implies that every solution u(t) satisfying |u(0)| < 1 remains trapped inside the disc D, therefore it also can be extended for all t ∈ R.
Remark 2.2 Every solution u(t) to Eq. (11) satisfying the identity |u(t)| = 1 can be written in the form u(t) = e iψ(t) where ψ(t) is a solution to the equation Now we consider Eq. (11) with a 2π-periodic coefficient w(t). It is well-known [52,53] that the Poincaré map of such equation coincides with the Möbius transformation. Because of Proposition 2.1 this Möbius transformation maps unit disc D onto itself, therefore it can be written in the form

Remark 2.3
The fact that the Poincaré map of the periodic complex Riccati equation (11) coincides with the Möbius transformation (13) can also be justified in a different way, using the Lie theory, see [54,Sec. III].
The next proposition shows that the parameters θ and b in formula (13) can be uniquely determined using two solutions to Eq. (11) starting from the initial conditions u = 0 and u = 1.
Proof: The definition of the Poincaré map and Remark 2.2 imply The former equation yields b = −ζ. Inserting this into the latter equation we obtain a formula for e iθ . Notice that because of Proposition 2.1 we always have |ζ| < 1, and hence |b| < 1 too. Every 2π-periodic solution to Eq. (11) corresponds to a fixed point of the Poincaré map, or equivalently to a solution of the equation The periodic solution is stable or unstable, if the corresponding fixed point u * is stable or unstable with respect to the map M(u), and the latter condition can be easily verified by estimating the derivative M (u * ). Indeed, if |M (u * )| < 1, then the fixed point u * is stable.
On the other hand, if |M (u * )| > 1, then u * is unstable. Moreover, the special properties of the map M(u), see Remark 2.6, allow us to conclude that a fixed point u * with |M (u * )| = 1 is also stable provided it is non-degenerate.
In the next proposition, we show that every Poincaré map (13) with 0 < |b| < 1 has either a unique stable fixed point in the closed unit disc D, or a unique fixed point at all (in this case, the fixed point is degenerate and lies on the unit disc boundary ∂D).
Proposition 2.5 Suppose θ ∈ (−π, π] and 0 < |b| < 1, then Eq. (17) has a unique solution u 0 ∈ D such that |M (u 0 )| ≤ 1. This solution is given by the formulas Proof: For every |b| < 1 and u ∈ D equation (17) can be rewritten in the form Since b = 0 this is a quadratic equation which generically has two complex roots. We are going to check which of these roots lie in the unit disc D and what are their stability properties. To address the latter question we compute the derivative of the Möbius transformation (13) and evaluate its modulus (keeping in mind that |b| < 1) Depending on the sign of the difference |b| 2 − sin 2 (θ/2) we distinguish two cases. Case 1. Suppose |b| > | sin(θ/2)|, then two solutions to Eq. (19) read It is easy to verify that in this case |u + | = |u − | = 1. Moreover, we also obtain Obviously, for every θ ∈ (−π, π] and | sin(θ/2)| < |b| < 1 we have and hence Case 2. The other case is determined by the inequality |b| ≤ | sin(θ/2)|. If |b| < | sin(θ/2)|, then Eq. (19) has two solutions while for |b| = | sin(θ/2)| the values u + and u − given by (23) coincide. To estimate the moduli |u + | and |u − | we compute the difference Then for sin(θ/2) > 0 we obtain and hence |u + | > 1 and |u − | < 1. Similarly, for sin(θ/2) < 0 we obtain |u + | < 1 and |u − | > 1. Finally, we compute a difference relevant to formula (21) On the other hand, in the limiting case |b| = | sin(θ/2)|, formulas (20) and (23) yield Remark 2.6 Let us consider the formula (13) with 0 < |b| < | sin(θ/2)| in more detail. In this case, M(u) determines an elliptic Möbius transformation, see [55,Ch. 3.VII]. This means that it has two distinct fixed points that are neither attractive nor repulsive but indifferent. (Recall the equation |M (u ± )| = 1 from the Case 2 in the proof of Proposition 2.5.) Moreover, the transformation moves all other points of the complex plane in circles around the two fixed points. Therefore, according to the Lyapunov stability classification, both fixed points are stable, but not asymptotically stable.
Similarly one can verify that the other two cases | sin(θ/2)| < |b| < 1 and |b| = | sin(θ/2)| = 0 considered in Proposition 2.5 correspond to the Möbius transformations M(u) of hyperbolic and parabolic types, respectively. This is in accordance with the fact that M(u) has a pair of attracting and repulsive fixed points in the former case and a degenerate fixed point in the latter case. Note that the degenerate fixed point of a parabolic Möbius transformation is always unstable in the sense of Lyapunov [55,Ch. 3.VII].
Hence the assertion follows from formula (20) and from the first of two inequalities (22).

Solution operator U
In the previous section we showed that for every w ∈ C per ([0, 2π]; C) and s ∈ R the complex Riccati equation (11) has a uniquely determined 2π-periodic solution u(t) ∈ D that is stable in the sense of Lyapunov (or at least linearly stable in the degenerate case). Let us denote the corresponding solution operator . The definition of U is constructive and relies on the following steps: 1) Given w(t) and s one solves two initial value problems (14) and (15) and obtains coefficients b and θ of the Möbius transformation (13), see Proposition 2.4.
2) Using Proposition 2.5 one computes the initial value u 0 of the periodic solution u(t) that lies entirely in the unit disc D and, moreover, is stable provided |b| = | sin(θ/2)|. In the case b = 0, one assumes u 0 = 0, see Remark 2.7.
3) One integrates Eq. (11) with the initial condition u(0) = u 0 and obtains 2π-periodic solution u(t). Importantly, Propositions 2.4 and 2.5 ensure that the steps 1-3 can always be realized. Therefore, the mapping U : (w(t), s) → u(t) is well-defined. Remark 2.9 Note that the minimal period of the function u(t) = U(w(t), s) does not have to be 2π. In general, it can assume any value 2π/k with k ∈ N. Moreover, for certain values of the arguments (w(t), s) the operator U can also return a constant function u(t).

Proof:
We need only to show that M * = M (u 0 ) where u 0 = u * (0). Then the assertion follows from Propositions 2.1 and 2.5 and from Remarks 2.7 and 2.8.
Let us consider Eq. (11) for w(t) = w * (t) and s = s * . Inserting there ansatz u(t) = u * (t) + v(t) and linearizing the resulting equation with respect to small perturbations v(t) we Obviously, formula (24)

Derivatives of the operator U
In this section we show how to compute partial derivatives of the operator U.
Proof: This assertion has been proved in [ Proof: Differentiating this identity with respect to ε and inserting ε = 0 and s = s * , we obtain Now, using Proposition 2.12 we solve Eq. (29) with respect to v ε (t, 0, s * ) and obtain formula (26). Formula (27) is justified similarly. We differentiate (28) with respect to s and solve the resulting equation using Proposition 2.12.

Self-consistency equation
Suppose that Eq. (5) has a solution of the form (8) where a(x, t + T ) = a(x, t) for some T > 0. Let us define then the new function u(x, t) is 2π-periodic with respect to t and satisfies Dividing Eq. (31) by ω and introducing the notations we rewrite Eq. (31) in the form In Section 2 we showed that every stable solution to Eq. (33) that lies entirely in the unit disc D is given by the formula u(x, t) = U(w(x, t), s). Inserting this result into the definition of w(x, t) we arrive at a self-consistency equation which can be written in the equivalent form (9). Eq. (31) has several continuous symmetries. More precisely, the set of its solutions is invariant with respect to the following transformations: 3) time shifts u(x, t) → u(x, t + τ ) for τ ∈ R. All these symmetries are inherited by the self-consistency equation (9), therefore to select its unique solution w(x, t) we need to provide three pinning conditions. In practice, this number can be reduced by one if we restrict Eq. (9) to the space of even functions Indeed, for symmetric coupling kernels G(x) equation (9) is reflection symmetric with respect to x, therefore we can look for solutions w(x, t) satisfying w(−x, t) = w(x, t) only. In this case the spatial translation symmetry is eliminated automatically. Then two pinning conditions relevant to the complex phase shift and the time shift can be chosen in the form In the next sections we will show that the augmented system consisting of Eqs. (9), (34) and (35) is well-defined. This means that for fixed phase lag α and kernel G(x) it correctly determines the unknown even function w(x, t) and two scalar parameters ω and s.

Modified self-consistency equation
In this section we show that the phase shift symmetry can also be eliminated from Eq. (9). Then we decrease the number of equations and unknowns in the augmented system described above.
Let us define a linear operator which gives a constant part of the function w(x, t). Using this operator and the identity operator I, we rewrite Eq. (9) in the equivalent form 2ωe iα (I − P)w = (I − P)GU(w, s).
If we assume w(x, t) = p + v(x, t), where p ∈ (0, ∞) and v(x, t) ∈ {u ∈ X e : Pu = 0} , then pinning condition (34) is fulfilled automatically and can be discarded. Moreover, inserting the ansatz (38) into Eq. (37) and into the pinning condition (35) we obtain and Now instead of solving the system of equations (9), (34) and (35), we can look for solutions of the system comprising Eqs. (39) and (40). In this case p > 0 must be given, then the system of equations (39) and (40) has to be solved with respect to two unknowns: scalar parameter s and even function v(x, t) satisfying Pv = 0. As soon as such solution is found, one can compute the corresponding values of ω and α from Eq. (36) written in the form 2ωe iα = 1 p PGU(p + v, s).

Modified self-consistency equation for cosine kernel (2)
In this section we consider a specific example of integral operator G and show how system (39), (40) can be solved approximately using Galerkin's method. For this we assume that G(x) is the cosine kernel (2). Given a positive integer F let us define a set of 8F + 2 functions ψ k (x, t) √ 2 cos x, i √ 2 cos x, e imt , ie imt , e imt √ 2 cos x, ie imt √ 2 cos x, m = −F, . . . , −1, 1, . . . , F.
Note that m = 0, therefore constant functions √ 2 and i √ 2 are not included in the set. The order of ψ k (x, t) is irrelevant apart from the only place below where we will assume ψ 8 (x, t) = ie −it . It is easy to verify that ψ k (x, t) satisfy the orthonormality condition ψ k , ψ n = δ kn with respect to the scalar product where δ kn is the Kronecker delta. Therefore, functions ψ k (x, t) span a finite-dimensional subspace of {u ∈ X e : Pu = 0}. We look for approximate solution to Eq. (39) where c k ∈ R are unknown coefficients. Inserting (41) into Eq. (39) we write 8F +2 orthogonality conditions ψ n , (PGU(p + v, s))v = p ψ n , (I − P)GU(p + v, s) , n = 1, . . . , 8F + 2.
Since for the cosine kernel G(x) it holds PG = P, the above system can be written as follows To account for the pinning condition (40) we assume ψ 8 (x, t) = ie −it , then This means c 8 = 0. Inserting this identity into Eq. (42) we end up with a system of 8F + 2 nonlinear equations with respect to 8F + 2 real unknowns (these are 8F + 1 coefficients c k with k = 8 and the parameter s). The system (42) can be solved by Newton's method, using a semianalytic Jacobian expression involving the derivative representations obtained in Section 2.2. Note that breathing chimera states typically have a very fine spatial structure, therefore to approximate the integrals in (42) with the same accuracy, one needs to use either a nonuniform grid with a moderate number of nodes in the x-direction, or a uniform grid with a much larger number of nodes. For example, all numerical results reported in Section 5 were obtained using a nonuniform grid with ca. 10 3 discretization points in the x-direction (the distribution of points, in this case, was 10 to 100 times denser in the vicinity of the coherence-incoherence boundaries than in the other regions of the chimera state). On a uniform grid, the same accuracy would be acheived only with at least 10 5 discretization points, what would lead to extremely large computational times.

Formulas for global order parameter and effective frequencies
Every solution z(x, t) to Eq. (5) corresponds to a probability density f (θ, x, t) representing a specific statistical equilibrium of the large system of coupled oscillators (1), see [45]. More precisely, the function f (θ, x, t) yields the probability to find oscillator with phase θ at position x at time t and has the characteristic property 2π 0 f (θ, x, t)e iθ dθ = z(x, t).
Property (43) can be used to derive formulas for the global order parameter (4) as well as for the effective frequencies of oscillators (3). Indeed, formula (4) determines the average of e iθ k (t) over all oscillator indices k. For large enough N this average can also be computed as follows To obtain a similar formula for Ω eff,k we write Eq. (1) in the form Recall that x k = −π + 2πk/N is the spatial position of the kth oscillator, therefore for infinitely large N points x k densely fill the interval [−π, π]. In this case instead of the discrete set of functions W k (t) we can consider a function W (x, t) depending on the continuous argument x ∈ [−π, π]. Because of the property (43) we have W (x, t) = (Gz)(x, t) and Notice that formulas (44) and (45) are the continuum limit counterparts of formulas (4) and (3) for any solution z(x, t) to Eq. (5). If the solution z(x, t) is taken in the form (8), then we obtain the following proposition.
Proposition 3.1 Let the triple (w(x, t), ω, s) be a solution to the self-consistency equation (9) and let u(x, t) = U(w(x, t), s) and a(x, t) = u(x, ωt). Then the following formulas hold and

Extraction of the parameters of breathing chimera states
Suppose that we observe a breathing chimera state in system (1) and we wish to extract from the observation its primary Ω and secondary ω frequencies as well as the amplitude a(x, t) of the corresponding solution (8) to Eq. (5). This can be done in the following way. Let R min and R max be the minimal and the maximal values of |Z N (t)|, see (4), over a sufficiently long observation time interval for a breathing chimera state. If t k are consecutive time moments such that the graph of |Z N (t)| crosses the mid-level (R min + R max )/2 from below, then the period T of breathing chimera state can be computed as an average of all differences t k − t k−1 . On the other hand, if we compute the variation of the complex argument of Z N (t) from the time moment t k−1 till t k , then the quotient Remark 3.2 Note that the above method for determining the breathing period T relies on the assumption that the |Z N (t)|-graph crosses the mid-level (R min +R max )/2 only twice on the period (once from below and the other from above). If this is not the case, then one needs to select another mid-level value in the interval (R min , R max ), which guarantees the two intersections condition.
As soon as the period T and the primary frequency Ω are known we can find approximate values of the function a(x, t) in (8) by where the indices j are taken modulo N , x k = −π + 2πk/N is the scaled position of the kth oscillator and M = [ √ N /2] is the largest integer that does not exceed √ N /2. Finally, using formulas (10) we can compute the secondary cyclic frequency ω as well as the ratio s and the function w(x, t) appearing in the self-consistency equation (9). Remark 3.3 Note that using the method described above, we automatically obtain the function a(x, t) and cyclic frequency Ω which satisfy the calibration condition from Remark 1.1.

Stability analysis
Suppose that Eq. (5) has a solution of the form where a(x, t) is T -periodic with respect to its second argument. To analyze the stability of this solution we insert the ansatz z = (a(x, t) + v(x, t))e iΩt into Eq. (5) and linearize it with respect to the small perturbation v(x, t). In the result we obtain a linear equation with T -periodic coefficients where η(x, t) = iΩ + e iα a(x, t)Ga.
Along with Eq. (48) it is convenient to consider its complex-conjugate version These two equations can be written in the operator form where V (t) = (v 1 (t), v 2 (t)) T is a function with values in C per ([−π, π]; C 2 ), and For every fixed t the operators A(t) and B(t) are linear operators from C per ([−π, π]; C 2 ) into itself. Moreover, they both depend continuously on t and thus their norms are uniformly bounded for all t ∈ [0, T ]. Our further consideration is concerned with the stability of the zero solution to Eq. (50). Therefore, we are dealing only with the linear stability of the solution (47). We apply the methods of qualitative analysis of differential equations in Banach spaces [57]. Since A(t) and B(t) are uniformly bounded operators, we can define an operator exponent Then the solution of Eq. (50) with the initial condition V (0) = V 0 is given by the formula V (t) = E(t)V 0 . Recalling that A(t) and B(t) are T -periodic, we conclude [57, Chapter V] that the stability of the zero solution to Eq. (50) is determined by the spectrum of the monodromy operator E(T ). Roughly speaking, the necessary condition for the stability of the zero solution to Eq. (50) is that the spectrum of the operator E(T ) lies entirely in the unit circle of the complex plane. Otherwise, this solution is unstable.
The main problem in the application of the above stability criterion is concerned with the fact that the monodromy operator E(T ) acts in an infinite-dimensional functional space. Therefore, its spectrum consists of infinitely many points, which can be arbitrarily distributed in the complex plane. Below we use the explict form of the operators A(t) and B(t) and show the following properties of the monodromy operator E(T ): (i) The spectrum of the operator E(T ) is bounded and symmetric with respect to the real axis of the complex plane. It consists of two qualitatively different parts: essential spectrum σ ess and discrete spectrum σ disc .
(ii) The essential spectrum σ ess is given by the formula (iii) The discrete spectrum σ disc comprises finitely many isolated eigenvalues µ, which can be found using the formula µ = e λT where λ are roots of a characteristic equation specified below.

Proposition 4.1 The monodromy operator E(T ) can be written as a sum
where E 0 (T ) is a multiplication operator of the form and K is a compact operator from C per ([−π, π]; C 2 ) into itself.
Proof: Every function V (t) satisfying Eq. (50) and the initial condition V (0) = V 0 solves also integral equation where On the other hand, every solution to Eq. (53) can be decomposed into a sum where W (t) is a solution to the integral equation The Volterra integral equation (55) has unique solution W (t) that continuously depends on the initial value V 0 . Moreover, the mapping V 0 → W (T ) is a compact operator from C per ([−π, π]; C 2 ) into itself (recall the compactness of the operator G involved in the definition of the operator B(t)). This fact along with the formula (54) implies that the monodromy operator E(T ) is the sum of the multiplication operator E 0 (T ) and a compact operator abbreviated by W (T ). The spectrum of monodromy operator E(T ) consists of all numbers µ ∈ C such that the difference operator E(T ) − µI is not invertible. Because of the definition of A(t) and B(t) this spectrum is symmetric with respect to the real axis. Moreover, since A(t) and B(t) are uniformly bounded for t ∈ [0, T ], the monodromy operator E(T ) is bounded too, and hence its spectrum lies in a circle of finite radius of the complex plane. Other spectral properties of E(T ) follow from the decomposition formula (52).
Indeed, formula (52) implies [58] that the essential spectrum of monodromy operator E(T ) coincides with the spectrum of multiplication operator E 0 (T ). Using the definition of A(t) we obtain This allows us to calculate the spectrum of E 0 (T ) explicitly and obtain formula (51) for σ ess .
Remark 4.2 Suppose that we consider a solution (47) to Eq. (5) with the amplitude a(x, t), primary frequency Ω and secondary frequency ω, which satisfy the self-consistency equation (9) with w(x, t) and s defined by (10). Then inserting (30) and (32) into (49)  (is + 2w(x, t)u(x, t))dt and therefore formula (51) and Proposition 2.10 imply that every µ ∈ σ ess lies either on the boundary of the unit circle |µ| = 1 or on the interval (0, 1] of the real axis. Hence the essential spectrum σ ess cannot be relevant to any instability of the solution (47) obtained from the selfconsistency equation (9). Note that, in general, Eq. (5) may have solutions of the form (47) with unstable essential spectrum, but such solutions do not satisfy the self-consistency equation (9).

Formula (52) also implies
For every µ / ∈ σ ess the right-hand side of this formula is a sum of the invertible operator E 0 (T ) − µI and the compact operator K, hence it defines a Fredholm operator of index zero. This means that apart from the essential spectrum σ ess the monodromy operator E(T ) can have a discrete spectrum σ disc consisting of eigenvalues of finite multiplicity. Since the entire spectrum σ ess ∪ σ disc is confined in a bounded region of the complex plane, there can be at most finitely many such eigenvalues. These eigenvalues can be computed only numerically and in the following we outline the way how this can be done.
has a nontrivial T -periodic solution, then the number µ = e λT is an eigenvalue of the monodromy operator E(T ). Conversely, for every nonzero µ ∈ σ disc there exists a number λ ∈ C such that Eq. (56)   This is equivalent to the identity E(T ) = e λT I that ends the proof.

Remark 4.4
Notice that formula µ = e λT is not a one-to-one relation between λ and µ. For every λ ∈ C it yields one value µ. In contrast, given a nonzero µ one obtains infinitely many corresponding values λ, namely λ = (log µ + 2πki)/T with k ∈ Z.
Proposition 4.5 Let λ be a complex number such that e λT / ∈ σ ess , then for every continuous T -periodic function F (t) there exists a unique T -periodic solution of equation which is given by Proof: This assertion can be proved by analogy with [56,Proposition A.1]. Proposition 4.5 implies that for every λ ∈ C such that e λT / ∈ σ ess all T -periodic solutions V (t) of Eq. (56) satisfy also the integral equation This fact can be used to compute the discrete spectrum σ disc numerically. To this end let us choose C > 0 and consider Eq. (57) in the rectangular region If we find all λ ∈ Π such that Eq. (57) has a bounded nontrivial solution V (t), then according to Proposition 4.3 and Remark 4.4 we also find all eigenvalues µ = e λT of the monodromy operator E(T ) lying in the circular region e −CT ≤ |µ| ≤ e CT . Since the spectrum of the monodromy operator E(T ) is bounded, this ensures that for sufficiently large C we determine all eigenvalues µ ∈ σ disc relevant to the stability of the solution (47). Indeed, considering Eq. (57) for λ ∈ Π we may overlook eigenvalues µ in the circle |µ| ≤ e −CT . However, all these eigenvalues satisfy |µ| < 1 and therefore have no impact on the stability of the solution (47).

Computation of the discrete spectrum
Recalling the definitions of B(t) and D λ (t, t ) it is easy to see that Eq. (57) is a homogeneous Fredholm integral equation. In general, it cannot be solved explicitly, but its solutions can be found approximately using Galerkin's method. For this one needs to choose a set of linearly independent functions ϕ k (x, t), k = 1, . . . , K, which are 2π-periodic with respect to x and Tperiodic with respect to t. Without loss of generality it can be assumed that these functions are orthonormalized with respect to the scalar product such that ϕ j , ϕ k = δ jk . Then one looks for an approximate solution to Eq. (57) in the form Inserting this ansatz into Eq. (57) and writing the projected problem, we obtain This is a system of linear algebraic equations for the K two-dimensional coefficients V k . Obviously, it has a nontrivial solution if and only if λ satisfies the characteristic equation where Solving Eq. (58) one obtains approximate eigenvalues λ of Eq. (57) and hence the corresponding approximate eigenvalues µ = e λT of the monodromy operator E(T ). Taking into account that the functions ϕ k (x, t) appearing in the definition of matrix M 2K (λ) must be 2π-periodic with respect to x and T -periodic with respect to t it is convenient to choose them in the form of spatiotemporal Fourier modes. More precisely, let K x and K t be two positive integers, then we assume ϕ nm (x, t) = e inx+2πimt/T , n = −K x , . . . , K x , m = −K t , . . . , K t .
Notice that the functions ϕ nm (x, t) have an important property. If coupling kernel G(x) has a Fourier series representation then for all integer indices n and m we have This implies therefore in the case of functions ϕ nm (x, t) formula (59) can be written The main advantage of the latter expression is that it does not contain any operators, but only explicitly known functions. More precisely, the term D λ (t, t ) is a (2 × 2)-matrix with entries depending on x, t, t and λ, while B 0 (t ) is a (2 × 2)-matrix with entries depending on x and t . Importantly, the triple integration in (60) must be carried out for each entry of the resulting product matrix separately. As soon as all matrices B nmn m (λ) are determined they must be combined into the block matrix M 2K (λ) and then the left-hand side of the characteristic equation (58) can be computed.

Remark 4.6
In Section 3 we mentioned that Eq. (31) has three continuous symmetries. This implies that the monodromy operator E(T ) has three linearly independent eigenfunctions corresponding to the unit eigenvalue µ = 1. Respectively the characteristic equation (58) has the triple root λ = 0, which is embedded in the essential spectrum.

Example
Let us illustrate the performance of the methods developed in Sections 3 and 4. For this we consider the breathing chimera state shown in Figure 2. The primary frequency Ω and the secondary frequency ω of this state can be extracted from the time trajectory of the global order parameter Z N (t), see Section 3.4. Then using formula (46) we compute an approximate amplitude a(x, t) of the corresponding solution (8) to Eq. (5). Inserting these results into formulas (30) and (32) we obtain approximate values of the parameter s and function w(x, t). Finally using continuous symmetries of Eq. (31) we ensure that w(x, t) is even with respect to x and satisfies the pinning conditions (34) and (35). The obtained function w(x, t) can be represented as a Fourier series Then the leading coefficientsŵ 0,k ,ŵ 1,k with indices k = −10, . . . , 10 can be used as an initial guess in the Galerkin's system (42) with F = 10. The latter system was solved using Newton's method up to the accuracy 10 −9 . The obtained set of coefficients c k was transformed into function w(x, t) using formulas (38) and (41). Then the corresponding solution u(x, t) to Eq. (31) was computed using the operator U defined in Section 2.1.  shows an approximate amplitude a(x, t) computed using formula (46) for the breathing chimera state in Fig. 2. We also solved the self-consistency equation (9) for the same parameters A and α using Galerkin's system (42) and found a time-periodic solution u(x, t) to Eq. (31), see Figure 3(c),(d). As expected, the graphs of a(x, t) and u(x, t) agree with each other on a large scale, but have some fine structure differences which can be attributed to finite size effects. The assertion is confirmed by simulations of system (1) with more oscillators (not shown). In particular, several darker filaments protruding into the coherent region (yellow/bright) in Fig. 3(a) become thinner for growing system size and disappear in the limit N → ∞, see Fig. 3(c), in accordance with the coherence/incoherence invariance property described in Proposition 2.1.
The self-consistency equation (9) allows us to predict almost perfectly the graphs of the global order parameter Z N (t), see (4), and of the effective frequencies Ω eff,k , see (3). In Figure 4 these quantities, computed for a breathing chimera state in system (1), are compared with their continuum limit counterparts Z(t) and Ω eff (x) computed by the formulas from Proposition 3.1 where we inserted the functions w(x, t) and u(x, t) obtained from the Galerkin's system (42). Figure 5 illustrates another application of the self-consistency equation (9). We used it for computation of a branch of breathing chimera states in Eq. (5). The theoretically predicted primary and secondary frequencies are compared with the corresponding values of Ω and ω observed in breathing chimera states in the coupled oscillator system (1) with N = 8192. Again, the agreement between the theoretical curve and the numerical points is very good. A slightly recognizable mismatch can be attributed to finite size effects or to the small number of Fourier modes (F = 10) in the Galerkin's approximation. Note that the curves in Figure 5 fold for α ≈ π/2 − 0.145. This fact explains a sudden collapse of breathing chimera states to the completely synchronous state, which we observed in system (1) for α > π/2 − 0.145. For all breathing chimera states on the solution branch in Figure 5 we also computed the corresponding essential spectra σ ess . These spectra look identically, see Figure 6(a), and have the maximal possible size, see Remark 4.2 for more detail. The computation of the discrete spectra σ disc turned out to be a time demanding task, therefore at present we were unable to carry it out. However, because of Remark 4.6 we assert that σ disc includes a triple eigenvalue µ = 1 embedded into the essential spectrum. Since for every breathing chimera state the unit circle |µ| = 1 is a subset of its essential spectrum σ ess , the destabilization of such chimera state cannot occur via a classical bifurcation of finite codimension. Indeed, any unstable eigenvalue can emerge only from the essential spectrum and therefore this eigenvalue cannot be isolated at the bifurcation point. By analogy with other dynamical systems we may expect that breathing chimera states in general can lose their stability via a fold, symmetry breaking, period-doubling or torus bifurcation, which are associated with the appearance of one or two unstable eigenvalues µ from the essential spectrum on the unit circle as shown in Figure 6(b)-(d). Note that proper consideration of such bifurcations requires the use of generalized spectral methods [59,60,61], which, however, must be adapted to a situation where the reference solution is a relative periodic orbit rather than a simple equilibrium.

Conclusions
In summary, we showed that breathing chimera states observed in large size systems (1) are properly represented by solutions of the form (8) in Eq. (5). The self-consistency equation (9) and Proposition 3.1 allow one to predict the most important features of such chimera states. A general approach for stability analysis of breathing chimera states is formulated in Section 4. It relies on the consideration of the monodromy operator which describes the evolution of small perturbations in the system. We found that the spectrum of this operator consists of two qualitatively different parts: essential and discrete spectra. The former part was completely characterized in this paper, while for the latter part we suggested a numerical algorithm for its approximate evaluation. Although we were unable to compute discrete spectrum of any breathing chimera state we obtained a theoretical indication of the fact that such chimera states lose their stability in nonclassical bifurcation scenarios where one or two unstable eigenvalues appear from the essential spectrum. Note that this statement is still speculative and therefore needs to be confirmed by computable examples. We emphasize that the consideration scheme suggested in this paper can be applied to systems (1) with arbitrary coupling kernels G(x), including exponential [36] and top-hat [37,38] coupling. In particular, using the self-consistency equation (9) one can carry out a more rigorous asymptotic analysis of breathing chimera states reminiscent of that in [38]. Moreover, all above methods can be extended to the study of breathing spiral chimera states in two-dimensional  lattices of phase oscillators [39,40]. Furthermore, our results can also be applied to explore the appearance of pulsing and alternating coherence-incoherence patterns [41] and modulated travelling chimera states [62] in systems of heterogeneous nonlocally coupled phase oscillators, though in this case one needs to modify the definition of the solution operator U.
In a more general context, we believe that many theoretical constructions of this paper such as the solution operator of the periodic Riccati equation (11), the periodic self-consistency equation (9) and the essential and discrete spectra of the monodromy operator E(T ) can be adapted to a broader class of problems concerned with the application of the Ott-Antonsen ansatz [47]. These are, for example, the Kuramoto model with a bimodal frequency distribution, where one finds periodically breathing partially synchronized states [63], modular Kuramoto networks [64] as well as networks of theta neurons [65,66,67] (or equivalently, quadratic integrate-and-fire neurons [68,69]).