Basins of Attraction for Chimera States

Chimera states---curious symmetry-broken states in systems of identical coupled oscillators---typically occur only for certain initial conditions. Here we analyze their basins of attraction in a simple system comprised of two populations. Using perturbative analysis and numerical simulation we evaluate asymptotic states and associated destination maps, and demonstrate that basins form a complex twisting structure in phase space. Understanding the basins' precise nature may help in the development of control methods to switch between chimera patterns, with possible technological and neural system applications.

Chimeras are known to arise in systems with nonlocal coupling that decays with increasing distance between phase oscillators, thus bridging the gap between the extremes of local (nearest-neighbor) and global (all-to-all) coupling ‡. Such long-range coupling is characteristic of many real-world technological and biological [34][35][36] systems. In many systems, chimeras are steady-state solutions stably coexisting with the fully synchronized stateS, not emerging via spontaneous symmetry breaking, and are thus only attained via a certain class of initial conditions [14,24,29]. Figure 1 graphically demonstrates this puzzling aspect of basins of attraction for chimera states: apparently similar initial conditions (panel B) can evolve to completely different steady-states (panel C). Thus, a natural question arising in any practical situation is: given a random initial phase configuration, how likely is the system to converge to a chimera state? Even though this important question was raised in 2010 [39], basins of attraction for chimera states have not yet been investigated systematically.
One difficulty with the examination of basins of attractions is that they are computationally expensive to obtain, e.g. via Monte Carlo simulation [40]. Here, in contrast, we use primarily analytic methods to explain the structure of the phase space and provide a systematic study of the basins of attraction leading to chimeras in the thermodynamic limit.
Model. The simplest realization of nonlocal coupling is achieved with two populations, where each population is more strongly coupled to itself than to the neighboring population (see Figure 1 panel A). It has been used as a model for several investigations of chimera states [15,26,27,[41][42][43][44]; here chimeras manifest themselves ‡ For oscillators with non-constant amplitudes it appears that local and global coupling can both be sufficient [31][32][33]. S Known exceptions where chimera states and fully synchronized states are not necessarily bi-stable are certain generalizations of the present system with amplitude dynamics [33,37] or non-linear delay feedback [21,38]. Three superficially similar oscillator phase distributions taken as initial conditions. C: Oscillator phase distributions after long-time evolution of system-each corresponds to the initial condition shown directly above it. DS: "desync-sync" state; SS 0 : "sync-sync" state; "SD": "sync-desync" state.
as a state with one synchronous and one asynchronous population. Accordingly, we consider the Kuramoto-Sakaguchi model with n = 2 populations [41,42] each of size N σ ,θ where θ σ k is the phase of the kth oscillator k = 1, . . . , N σ in population σ ∈ {1, 2} and ω is the oscillator frequency. For consistency with previous work [41,42,45], we assume the coupling is symmetric with neighbor-coupling K σσ = K σ σ = ν and self-coupling K σσ = µ. Imposing without loss of generality µ + ν = 1, the coupling can be parameterized by the coupling disparity A = µ − ν. We redefine the phase lag parameter via β = π/2 − α as chimeras emerge in the limit of near-cosinecoupling (β → 0) for this type of system [41,45,46]. The mean field order parameter R σ e iΦσ = (N σ ) −1 N σ j=1 exp (iθ σ j ) describes the synchronization level of population σ with R σ = 1 for perfect and R σ ≤ 1 for partial synchronization. We consider the thermodynamic limit N σ → ∞, allowing us to express the ensemble dynamics in terms of the continuous oscillator density f σ (θ, ω). This facilitates a low-dimensional description of the dynamics via the Ott-Antonsen (OA) ansatz [47][48][49] in terms of the mean-field order parameter of each population, ρ σ (t)e iφσ(t) = e iθ f σ (θ, t)dθ with 0 < ρ σ ≤ 1, see Appendix A and Appendix B.
By virtue of the translational symmetry φ σ → φ σ + const., the resulting dynamics State variables (ρ 1 , ψ, ρ 2 ) are interpreted as cylindrical coordinates. Phase space is structured by (i) two invariant rays, R 0 and R π (dashed); and (ii) two invariant surfaces, S 1 and S 2 , forming the side and top surfaces of the cylinder slab. Except for a set of measure zero, all trajectories converge to one of three locations: SD chimera state on S 1 (red), DS chimera state on S 2 (blue), or fully synchronized state SS 0 (yellow).

Invariant Manifolds and fixed points
Analysis of equations (2)-(4) reveals the existence of two invariant surfaces defined by S σ = {(ρ 1 , ρ 2 , ψ)|ρ σ = 1} ⊂ D (the top (blue) and lateral (red) surfaces of the cylinder displayed in Figure 2). The dynamics on these manifolds were studied previously [41]: chimera states are born in a saddle-node bifurcation and undergo a Hopf bifurcation for larger coupling disparity A = µ − ν. The resulting stable limit cycle grows with A until eventually it is destroyed in a homoclinic bifurcation. Studying the basins of attraction, we generalize the previous analysis by considering the entire three-dimensional phase space D.
Outside of S 1 and S 2 , trajectories follow a complex winding motion, structured around the two invariant rays R 0 and R π defined by ρ 1 = ρ 2 with ψ = 0 and ψ = π, respectively (see Figure 2 and Appendix C). Other than the origin, which is a repeller, there are no fixed points for ρ 1 , ρ 2 < 1 (see Appendix D). Thus, limit cycles in the interior of the phase space are also absent. In principle, a chaotic attractor could appear inside D but is not observed.

Numerical investigations
First insights regarding basins of attraction for chimera states were gathered via simple Monte Carlo integration of uniformly distributed random initial conditions for ρ 1 , ρ 2 ∈ [0, 1] and ψ ∈ [−π, π]. These computations reveal that the probability p(A, β) of ending up in a chimera state depends primarily on β with a maximum value for β → 0, see Figure 3. This approach provides information about the sizes of the basins of attraction, but it reveals little about their structure. We therefore ask: how is the three-dimensional phase space structured?
To better reflect symmetries of the phase space, the dynamics may be re-expressed in terms of the sum and difference of the order parameters (see Figure 2 In order to characterize the structure of the basins of attraction, we compute the destination maps for a set of initial conditions (s, d, ψ). Figure 4A shows a typical crosssection of the destination map with fixed s: basins form a spiraling structure around R π (the ray (d, ψ) = (0, π)), with SD and DS basins always separated by the (often thin) basin for SS 0 . The thickness of the basin spiral arms increases away from R π , with maximum near R 0 (the ray (d, ψ) = (0, 0)).
The area ratio between basins for SD (or, by symmetry, DS) and SS 0 is related Probabilities to obtain chimera states via random sampling of initial conditions (ρ 1 , ρ 2 , ψ). Chimeras appear within the wedge defined by a saddle-node bifurcation (SN, solid) for small A and a homoclinic bifurcation (HC, dotted) for large A [41]. Phase portraits in the ρ 2 = 1 plane are shown (insets) with stable nodes (full circles), unstable chimera (white circle) and saddle chimera (half-filled circles), together with its stable (solid) and unstable (dashed) manifolds. For intermediate A, the asynchronous order parameter undergoes a Hopf bifurcation (HB, dashed). Probabilities for ending up in either SD/DS chimera were measured by realizing 1000 random initial conditions (ρ 1 , ρ 2 , ψ) for each parameter value set. Further details are in Appendix E.
to the probability that a random initial condition will lead to a chimera state, and depends on parameters A and β as follows. For β → 0, the SS 0 basin occupies an infinitesimal fraction of the area. As β increases, the SS 0 basin increases its area until it occupies the entire plane at β = β SN (A) when the chimera state is annihilated through a saddle-node bifurcation. For A < A SN (β) or A > A HC (β), no chimera state exists and the entire basin belongs to the SS 0 state. With increasing A > A SN , the (total) basin area of SS 0 gradually decreases from 100% approaching a constant near the homoclinic bifurcation, see Figure E2.
As s increases from zero, basins merge and pinch-off in an alternating fashion (see Figure 4, Sec. 4 and Supplementary Video 1) so that the basin boundaries rotate clockwise about R 0 ((d, ψ) = (0, 0) in Figure 4A). Once s reaches s c ≈ √ 1 − A, this rotation stops, demonstrating that knowledge of the trajectory position in the s = s c plane is sufficient for determining its final fate.
The basin density appears singular near R π , with a nested structure that allows even tiny perturbations of the initial condition to strongly influence the final state. More generally, the highly alternating basin structure is reflected in the times to reach steadystate attractors, which are displayed in Figure 4 B. Figures 5 and E1 show destination times T along a section of that figure between the origin and (d, ψ) = (0, π), revealing a power-law behavior that may result from this nested spiral arm basin structure. Figure 5 also reveals that destination times diverge on the basin boundaries (see   Figure 4). Average destination times T grow like a power law as ψ → ±π, thus basin structure is self-similar around (d, ψ) = (0, π). T diverges at the boundary between SD/DS and SS 0 basins (inset), since these trajectories lie on stable manifolds leading to saddle points on invariant manifolds S 1 , S 2 .
inset), which is explained by the fact that these boundaries form separatrix sheets: these are the two-dimensional stable manifolds emanating from the saddle chimeras on S 1 and S 2 , originating in the saddle-node bifurcation that gives birth to chimeras, see Figure 2 and [41]. Numerical continuation of those sheets (see Figs. 6 and E4, and Appendix E.7) displays the same twisting motion as seen in Figure 4.

Analysis
A complete analysis of the basins for the entire phase space is difficult to achieve, but the main features of the basin structure have their origin in the invariant rays R 0 and R π about which a perturbation analysis can be made for small A and β (asymptotic results remain qualitatively in fair agreement for parameter values further off the origin).

Perturbation analysis around the invariant ray R 0
We consider the coupling constants (µ, ν) to be perturbed from global coupling (A = 0) by setting µ = 1 2 (1 + A), ν = 1 2 (1 − A) as in [41]. We then make the perturbative approximation that ψ, d, β and A are all small and of the same order (while keeping in mind that A > 2β is required for the existence of a chimera state in this limit). This means that near R 0 , we make the ansatz [41]: After making a change of variables x = d 1 , y = 1 2 sψ 1 , and then a second change of variables x = r cos(θ) and y = r sin(θ) we find the following equations: Note that the derivatives of r and s are both order while the derivative of θ is order 1 (when s is not close to 1). Thus θ evolves on a fast time scale while r and s evolve slowly. We may therefore use the method of averaging on the higher order terms involving θ to simplify these equations to The change of variables is chosen so that that the spiraling cycles become circular in shape.
This reveals a couple of properties. First of all, as expected, solutions will be spirals around the R 0 manifold. The radius r goes to 0 as t → ∞ (for the truncated equations). Thus, trajectories slowly converge toward the R 0 manifold until the approximations break down when dθ dt = 0 and higher order terms become significant. In other words, the R 0 -manifold is weakly attracting. The frequency of rotation is This is here referred to as the s c -plane. In this plane, the trajectories cease to have a spiral character and instead begin to separate and evolve toward the fully synchronized state or the DS or SD chimeras.
Note that there is an alternative way that the "critical plane" could be defined. Setting dθ/dt = 0 in Eq. (6) yields a minimum s solution s c = 1 − A 1 + O( 2 ) (possible only for particular θ values), thus s c ≈ 1 − A is the smallest value of s for which rotation about ray R 0 may stop. The difference between the two expressions s c ≈ √ 1 − A ≈ 1 − A/2 and s c ≈ 1 − A comes from whether averaging has been applied or not. In the former case (averaged equations), rotation about R 0 stops on average over all θ; in the latter case, rotation about R 0 stops only for some particular θ.
By symmetry, if a trajectory originating at (s, d, ψ) converges to the SD state, the trajectory originating at (s, −d, −ψ) must converge to the DS state. Therefore, the position relative to a separating boundary in the s c plane determines the final state. Numerical integration confirms that trajectories converging to SD and DS chimeras form opposing sides of a positively oriented double helix centered on R 0 (red and blue in Figure 2 and Supplementary Video 2).
The "rotation" of the basin boundary as s increases along the R 0 manifold (indicated symbolically in Figure 4 A) can also be understood analytically in the perturbative limit close to R 0 . Equation (10) can be solved explicitly to get Takingθ = (1 − s 2 )/2 to lowest order (from (9)), we can substitute in for s(t) and then integrate from t = 0 to t = t crit to approximate the total angle change about the R 0 manifold over the course of the trajectory. Here t crit represents the time at which the trajectory s(t) reaches the critical plane s = s c : Integratingθ gives a total angle change (This can also be written as ∆θ = β −1 ln(s c /s 0 ), and then this expression is valid for either definition of s c .) Thus, the boundary angle is proportional to β −1 ln(1 − A/2) − β −1 ln(s), yielding a rotation rate of (βs) −1 as the section plane s varies uniformly.
Since the angle of a trajectory at the critical plane determines the basin the trajectory belongs to, the appearance of the basin boundary in a section plane orthogonal to the ray R 0 is just a line with angle proportional to ∆θ.

Perturbation analysis around the invariant ray R π
We can perform a similar analysis around the R π ray by a similar ansatz: This time we make an (analogous) change of variables 1+s 2 ψ 1 ¶, and again convert to polar coordinates. This yields Again we find that the derivative of θ is larger than the other derivatives, allowing us to reduce the system to dr dt Similar to the results near R 0 , this analysis reveals a spiraling motion around the R π manifold, but here the radius diverges exponentially. In contrast to the previous case, three distinct time scales are present: the derivative of θ is an order of magnitude larger than the derivative of r and two orders of magnitude larger than the derivative of s. This means that the rotation around the R π manifold and the radial divergence away from the manifold occur more quickly than the translation along the manifold. In other words, each trajectory (and consequently any basin boundary) winds around R π within a plane with approximately fixed s (see cross section in Figure 4 and Supplementary Video 3. Rotation around the manifold occurs at a faster rate than divergence away, and both occur faster than translation along the manifold. ¶ The change of variables is chosen so that that the spiraling cycles become circular in shape; the particular shape differs from the one near the R 0 -manifold.

Shape of the basin boundaries
We can understand two qualitative aspects of the basin boundaries seen in numerics: (1) the basin boundaries are linear in the s c -plane near R 0 , and (2) the basin boundaries have spiral shape near R π .
The invariant ray R 0 is surrounded by the SS 0 basin, and as we move further away form R 0 , we enter the SD and DS basins, respectively (Figs. 4 and E3, F1C). The boundaries (separatrices) between SS 0 /SD and SS 0 /DS basins are the two stable manifolds leading to the SD/DS chimera saddles. The relative width of the SS 0 and SD/DS basins varies with parameters A and β (Fig. E2), and the SS 0 basin is thin close to the origin, and for smaller (A, β).
To understand the basin shapes near R 0 , it is helpful to consider trajectories generated by points along a straight line orthogonal to R 0 . We will thereby use the asymptotic results (5)-(7) derived in the previous section valid for smaller A and β. Consider the set of points along a line segment parameterized by k with |k| ≤ 1 where r = kr 0 , θ = θ 0 and s = s c . The trajectories generated by integrating the equations with initial conditions along that line segment will intersect the plane perpendicular to R 0 at some later time at s = s c + δ where δ . By symmetry, if a point along the line with k > 0 evolves toward the DS chimera, then the corresponding point with k < 0 will be mapped to the SD chimera. Similarly, if a point with k > 0 is mapped to the synchronized state, the corresponding point with k < 0 will also be mapped to the synchronized state. Suppose r 0 δ, so that all points along the line segment are chosen to be arbitrarily close to the R 0 manifold. According to Eqs. (5)- (7), dθ is independent of r, so the images of these points in the plane s c + δ will remain collinear. Now dr = O(δ 2 ) and ds = δ, so not only will the points along this segment remain collinear; as s increases, they will also remain arbitrarily close to each other and to R 0 . For at least one particular choice of θ 0 , this line segment will reach the surface of the cylinder in a direction tangent to the intersection of the invariant surfaces S 1 and S 2 . As discussed above, this intersection is itself an invariant manifold, and thus the entire line segment will be mapped to the synchronized state, being the only attractor on the manifold. Thus there exists a line segment in the s c plane that lies in the basin of attraction for the SS 0 state and that separates the basins for SD and DS chimeras. This suggests that, at least sufficiently close to R 0 , the basin boundary between SD and DS chimeras must be linear.
Near R π the picture is different. Because there are three distinct time scales, with the evolution of both θ and r faster than the evolution of s, trajectories initially close to R π generate spirals within a fixed plane perpendicular to s -this is the origin of the spiral shape of the basin structure near R π .
These qualitative arguments can be made rigorous in the limit where the SS 0 basin becomes a set of measure zero (infinitesimal thickness). The basin boundaries (separatrices) are visualized in Fig. 6 (also see Fig. E3 for a close up of the basins near R 0 ).

Control strategies
Determination of the structure of the basins of attraction for this system naturally invites the question of whether we can "control" the system. This can mean several things, among them: (1) can we intervene during an initial transient so as to direct the system to a desired equilibrium; (2) can we perturb the system to move it from one stable equilibrium to another; (3) can we stabilize an unstable equilibrium. Each of these questions can also be examined with the goal of finding an "optimal" strategy of some kind, where optimality is usually defined as minimizing some aspect of the intervention. To answer questions (1) and (2), knowledge of the basin structure in the thermodynamic limit N → ∞ is clearly useful, at least for sufficiently large N σ .
A full exploration of control and intervention strategies is beyond the scope of this paper, but as a demonstration of the power of our approach, we have performed a simple experiment. We wish to take a system at equilibrium in the DS chimera state, and to perturb it sufficiently that it goes to a different equilibrium (either SD or SS 0 ). We restrict ourselves to finite perturbations of the form θ where Q quantifies a uniform phase shift to all oscillators in the synchronous group (population 2). Since the perturbation leaves the system state on the invariant DS manifold, the final state may change from DS to SS 0 .
In the thermodynamic limit, this is equivalent to holding ρ 1 and ρ 2 = 1 constant while perturbing ψ via ψ → ψ − Q. Thus we expect the minimal required perturbation Q min to be determined by the size of the restricted DS basin of attraction (restricted to the surface ρ 2 = 1, the top surface of the cylinder shown in Figure 2).   Figure 7 shows the expected asymptotic value of Q min determined from our study of the basins of attraction presented here, as well as the Q min values determined via numerical experiment for finite values of N σ . The precise threshold varies slightly depending on the initial phases but asymptotes to the value observed in the continuum limit. This good agreement confirms that (1) our knowledge of the basins of attraction in the thermodynamic limit can indeed inform control strategies for chimera states, and (2) insight into the finite system is indeed gained by analysis of the thermodynamic limit.
Our analysis of the thermodynamic limit suggests that switching from the SS 0 state to a DS (or SD) chimera requires a perturbation to ρ 1 (or ρ 2 ), because ρ 1 = ρ 2 = 1 is an invariant manifold. Hence, while a uniform phase shift (θ i + Q) will perturb ψ, it will not be sufficient to desynchronize one of the two populations. Instead, a nonuniform phase shift that decreases the value of ρ 1 (or ρ 2 ) can accomplish the desired switching behavior.
Finally, we note that the control strategy presented here is quite naive. It is likely that more "optimal" strategies exist in the sense that the perturbation magnitudes could be reduced.

Discussion
The probability that a random initial condition evolves to a chimera state, while important for real-world applications, has not been a frequent topic of investigation [17,25]. Here, we have provided a detailed mathematical analysis unveiling the basin structure for a very simple system with two populations, allowing for insight into the chimera's relative rarity. It remains to be seen whether similar efforts applied to other models such as neuronal or pulse-coupled oscillators (where reduction methods have become available only very recently [50,51]) will bear fruit, and how basin structure in those systems will compare.
Oscillators on a ring with finite-range coupling exhibit chimera states that are (very long) transients with chaotic dynamics [52]. However, extensive computational analysis of the finite system (1) displays no such transient behavior while chaotic behavior is absent [44,53]; this difference in dynamic behavior (along with others) seems to be related to differing coupling topology [44]. Moreover, very small oscillator systems with N σ = 2, 3, 4 are shown to display asymptotic stability of chimera states [54].
Sampling the immense initial state space associated with the case of N σ < ∞ would be a burdensome task. Our analysis was facilitated by considering N σ → ∞, allowing us to focus on the low-dimensional order parameter dynamics on the OA-manifold [47]. While the higher dimensional dynamics off this manifold poses a challenge in its own right [30], the continuum theory allows us to gain useful insight by mapping the discrete to the continuous order parameter, R σ e iΦσ ≈ ρ σ e iφσ (identity for N σ → ∞). Though bifurcation boundaries may blur for very low N σ and the finely filigreed basin boundaries near R π may break down, general basin structures will look similar even for moderate N σ .
The stable manifolds of the saddles near SS 0 divide phase space into simply connected basins of attraction (see Figs. 6, E4 and Supplementary Video 4). Basins near the fully synchronized (SS 0 ) and chimera (SD/DS) states are simple in structure and relatively large + , resulting in robustness to perturbations (see also Fig. 7). The twisting motion around the invariant rays R 0,π , however, yields a complex basin structure which we explain analytically. As one approaches R π , the basin density diverges, and basins become locally intermingled [55]: perturbations in that region affect the fate of a trajectory drastically.
Continuum theory allows us to construct initial phase densities leading to chimera states via f σ (θ, t) = 1 2π 1 + ∞ n=1 ρ n σ e in(θ−φσ) + c.c. (see Appendix A). If instead initial phases θ k are sampled uniformly on [−π, π], the probability distribution for R σ is unimodal with mean ∼ 1/ √ N σ and variance ∼ 1/N σ . The probability distribution for Φ σ is invariant and uniform; this observation combined with the intermingled basin structure near R π explains why in practice random initial phases can lead to both chimera and fully synchronized states. Thus, this implies in particular that chimera + Local basin volumes of chimeras presumably scale like |x SD,DS − x SADDLE |, where x ∈ R 3 denotes coordinates of stable (SD/DS) and unstable (SADDLE) chimera states in the reduced phase space. states, given random initial conditions, are not always rare occurrences (depending on parameter values, see also Fig. 3).
The results we have presented focus on the case of two populations. While two populations allow for multi-stability between the fully synchronous state (SS 0 ) and two symmetrically equivalent chimera states (SD and DS), generalizations of such hierarchical structure to n > 2 populations [45,56] make accessible larger configuration spaces of size 2 n by variation of the synchronization-desynchronization patterns. One may wonder if any of the structures studied here retains relevance in cases of more than two populations [46,56]. It can be shown that the invariant hyper-ray corresponding to R 0 defined by ρ σ = ρ and φ σ = 0 exists for n > 2, and that the flow on this ray is ρ → 1. This suggests that the phase space is skeletonized similarly and a similar analysis may be feasible-a task left for a future study.
Biologically, systems with multiple coexisting chimera attractors have been proposed to describe 'metastable dynamics' required to modulate neural activity patterns [57], or to encode memory. Indeed, localized dynamical states are directly related to function in neural networks [58,59]; localized synchrony has been widely studied in neural field models as bump states [51,60,61], and are phenomenologically similar to chimera states. It is worth noting that chimeras occur in models of neural activity [62][63][64][65][66].
Implementations of chimera states could likely be achieved in micro-(opto)-electromechanical oscillators [67,68] where synchronization patterns may have technological applications. Conversely, as power grid network topologies evolve to incorporate growing sources of renewable power, the resulting decentralized, hierarchical networks [69,70] may be threatened by chimera states, which could lead to large scale partial blackouts and unexpected behavior.
The potential for applications-or threats-makes the dynamic re-configuration and switching between chimera configurations (possibly modulating functional properties of the underlying oscillator network) particularly relevant [71]; applications that modulate functional properties can only be achieved using detailed knowledge of the basin structure. As a test of principle, we successfully implemented a simple algorithm to move the finite oscillator system between different equilibria, demonstrating that an understanding of the basins of attraction for the N σ → ∞ system has value. We hope that future work will explore the construction of efficient control strategies to stabilize or prevent chimera states, with applications across many fields.

Acknowledgments
Research supported by the Dynamical Systems Interdisciplinary Network, University of Copenhagen (EAM). We thank C. Bick, R. Mirollo, S. Strogatz and O. Omel'chenko for valuable conversations and comments, and anonymous referees for helpful suggestions to improve the manuscript.
To make further progress, we consider the thermodynamic limit, i.e., the case of N σ → ∞ oscillators per population. This allows for a description of the dynamics in terms of the mean-field order parameter [47][48][49]. Eqs. (A.1) then give rise to the continuity equation where f σ (θ, t) is the probability density of oscillators in population σ, and v σ (θ, t) is their velocity, given by Here we have dropped the superscripts to simplify notation: θ means θ σ , and θ means θ σ . Following Ott and Antonsen [47,48], we consider probability densities along a manifold given by where * denotes complex conjugation and a σ (t) is given by Defining a σ (t) = ρ σ (t)e iφσ(t) where ρ σ and φ σ represent mean-field order parameters, the governing equation can be reduced to the system [28,41,45] The Ott/Antonsen manifold, in which the Fourier coefficients f n (t) of the probability density f satisfy f n (t) = a(t) n , is globally attracting for a frequency distribution with non-zero width ∆ [48]. For identical oscillators (∆ = 0), the dynamics for the problem (with n = 2 populations) can be described by reduced equations using the Watanabe/Strogatz ansatz [72], as shown in Pikovsky and Rosenblum [30]; the authors showed that Eqs. (1) may also be subject to more complicated dynamics than those described by the Ott/Antonsen ansatz. Studies by Laing [26,73] investigated the dynamics using the Ott/Antonsen ansatz for n = 2 populations for the case of nonidentical frequencies and found that the dynamics for sufficiently small ∆ is qualitatively equivalent to the dynamics obtained for ∆ = 0. It is therefore justified to discuss the dynamics for ∆ → 0 representing the case of nearly identical oscillators using the Ott/Antonsen reduction.

Appendix B. Governing equations for two populations
We restrict our attention to the case of n = 2 populations. Accordingly, we define the coupling parameters K 11 = K 22 = µ and K 12 = K 21 = ν; by rescaling time we can eliminate one parameter so that 1 = µ + ν without loss of generality. The remaining parameter is redefined via A = µ − ν, expressing the disparity of coupling between the two neighboring populations. By virtue of the translational symmetry, φ σ → φ σ + const., the dynamics of the system is effectively three dimensional. We introduce the angular phase difference ψ = φ 1 − φ 2 of the order parameter, and the resulting governing equations becomė where the state variables lie in the domain To investigate the basins of attraction, it proves useful to express the dynamics in terms of the sums and difference of the order parameters, i.e., we define and ψ as above. These variables belong to the domain defined by ψ ∈ [−π, π], s ∈ [0, 1] and d ∈ [−a, a] with a(s) = 1 2 − | 1 2 − s| (the back-transformation is ρ 1 = s + d and ρ 2 = s − d, without the factor of 2.). The governing equations are then expressed aṡ where σ = 1, 2 refers to the SD, DS manifolds, respectively. The dynamics in one manifold is identical to the other via the symmetry operation defined by operator Σ, see main text. The dynamics on these IMs is analyzed in [41].
One-dimensional invariant manifolds. Our numerical investigations indicate the presence of an invariant manifold at ψ = 0, π with ρ 1 = ρ 2 . Substituting these values into Eqs. (B.1)-(B.3), we geṫ 3) The first equation implies that any initial point on the rays with d = 0 and ψ = 0, π remains there for all times; if 0 < β < π, the trajectory moves towards the SS 0 attractor according to (C.2). Thus, two invariant rays exist, defined via Note that another one-dimensional invariant manifold S 12 is defined as the intersection S 1 ∩ S 2 , and any initial point with s = 1 on S 12 will therefore always end up in the SS 0 state.

Appendix D. Fixed points.
Fixed points on S 1,2 . The fixed points in the S 1,2 manifolds are the SD, DS chimera states and fully synchronized SS 0 states that are discussed in detail in [41]; note that since S 12 is an invariant manifold, there must be another fixed point in addition to SS 0 contained in it, with opposite stability: this source is found at (ρ 1 , ρ 2 , ψ) = (1, 1, π), which we refer to as SS π . Figure 2 illustrates how trajectories nearby are repelled from the ray R π . On S 1,2 , stable chimera states are born through a saddle node bifurcation, and undergo a Hopf bifurcation for sufficiently large disparity values A so that ρ σ < 1 is oscillatory; the associated limit cycle is destroyed in a homoclinic bifurcation with even larger A.
Chimera states. In addition to the in-phase (ρ 1 = ρ 2 = 1 and ψ = 0) and anti-phase (ρ 1 = ρ 2 = 1 and ψ = π) equilibrium points, there are also three equilibrium points with ρ 2 = 1 and ρ 1 = 1 (and three analogous fixed points with ρ 1 = 1 and ρ 2 = 1) [28]. These equilibrium points represent chimera states. Numerics suggest that two of these equilibrium points occur near ψ = 0 and one occurs near ψ = π. Using an ansatz motivated by these numerical results, we find that these equilibrium points satisfy the following scaling relationships (where A = µ − ν and µ + ν = 1: i) Stable Chimera near ψ = 0 (DS): ii) Unstable Saddle Chimera near ψ = 0 (UC): iii) Unstable Chimera near ψ = π (UC): These relationships are useful when trying to solve for the precise fixed point locations numerically and for approximating their stable and unstable manifolds in order to deduce the basin boundaries. We note that chimera states asymptotically approach either SS 0 or SS π as (A, β) → (0, 0).
Origin. The origin is an unstable fixed point, as can be seen by linearizing for small ρ 1 and ρ 2 in Eqs. (B.1)-(B.3).
For all other cases, let us consider the equations by introducing K = µ/ν > 1 and ρ rel = ρ 2 /ρ 1 : We note now that ρ rel > 0 and sin β > 0 by assumption and we can eliminate these expressions as follows Equating (D.7) and (D.8), it follows that a fixed point with 0 < ρ 1 , ρ 2 < 1 and 0 < β < π can only exist if which has the solutions However, by assumption, we have K > 1 and the solutions are complex. Therefore, even if we find real values ψ, ρ 1 that satisfy the third fixed point equation (D.6), there will be no real solutions for ρ rel , and thus also not for ρ 2 . Therefore fixed points in the interior of the domain can be excluded when β = nπ where n is an integer. It should be noted that the numerical experiment above assumes that ρ 1 , ρ 2 , and ψ are uniformly distributed. For systems with a finite number of oscillators N σ in each population, the expected value of the order parameter value ρ σ is O(1/ √ N σ ). Hence, the probabilities computed using the above scheme should not be interpreted as the probability that a state with randomly selected initial phases θ (σ) k would evolve toward a chimera state. Instead, they represent the size of the basins of attraction of the chimera states relative to the size of the basin of attraction of the fully synchronized state in the continuum limit N σ → ∞.

Appendix E.2. Destination Maps
Simulations for a given initial condition were carried out until trajectories to a fully synchronized (limit point, LP) or a stable (LP) or breathing chimera (limit cycle LC) occurred. The detection of these three types of states was carried out in two steps, described below. Integration of Eqs. (B.1)-(B.3) or (B.6)-(B.8) were carried out in Matlab TM using the ode45 solver routine with event detection (see below) on a high performance computation cluster, with a relative error tolerance of 10 −8 . Algorithms below are outlined for (ρ 1 , ρ 2 , ψ)-coordinates; analogous detections for (s, d, ψ)-coordinates are carried out by applying the related coordinate transformations. Below, dρ 1,2 ≈ρ 1,2 dt and dψ ≈ψ dt denote the approximate differential values evaluated by the o.d.e. integrator at discrete time steps.

Appendix E.4. Estimating time to attractor
The following algorithm is adopted for obtaining estimates for the time to reach the attractor, T , i.e., the traveling time from initial to end condition. When these times are not of interest, the previous scheme is preferred due to significant gains in computational speed.
(ii) If the integration fails to detect a fixed point, the algorithm enters a loop of max. 100 iterations, where: i.) Integration is carried out to detect k = 1, . . . , 10 events of type Event A, limit point and limit cycles). Periods of limit cycles and event states (ρ 1 , ρ 2 , ψ)| t=t k are stored.
(iii) If LC or LP is detected, the final state is detected as explained above. Otherwise, failed convergence is stored as a failed end state.  which regions of the state space will evolve toward this chimera state, as shown in Figure 6 and Figure E4. The manifold can be approximated as follows: Step 1: Compute the two stable eigenvectors of the saddle chimera to obtain a local approximation to the stable manifold near the saddle chimera. (There are two stable eigenvectors v 1 and v 2 and one unstable eigenvector v 3 for the saddle chimera located at p. The stable eigenvectors define a plane tangent to the stable manifold.) Step 2: Obtain a family of starting points x 0 (θ) for continuation by making small perturbations off of the saddle chimera in every direction within the stable manifold. (Define a vector of angles θ and a magnitude . The family of starting points are defined by x 0 (θ) = p + [cos(θ)v 1 + sin(θ)v 2 ].) In Figure 6, we used 23 angles θ between 0 and π with the vectors v 1 and v 2 chosen so that all of these perturbations led to relevant parameter values and a perturbation magnitude of = 10 −6 . In Figure E4, 94 trajectories were used.
Step 3: Integrate backward in time from each point until the trajectories reach x 1 with a predetermined distance ||x 1 − x 0 || from the start point, and plot the surface containing these trajectories. (We used ode45 to integrate the equations and then interpolated to determine when the trajectories had reached the desired length.) In Figure 6, the predetermined distance was set at 0.01 to obtain a high resolution near the manifolds. In Figure E4, the predetermined distance was set between 1 and 20 (and no additional refinement was performed) in order to reduce data points for quick rendering, and in particular to enhance the visibility of the separatrices while reducing the total number of points displayed. This way the point families are equidistant in space, rendering an accurate picture of the separatrix surface in all regions.
Step 4: The endpoints of the trajectories define a curve. Use evenly spaced points along the curve as new starting points and return to step 3 until enough of the stable manifold has been computed.
In Figure 6, we used a spacing of 0.01 near the saddle chimera and 0.05 once the trajectories had reached a distance of 0.2 from the saddle chimera. While this method yields satisfactory results for the problem at hand, we mention that more advanced and accurate continuation methods are available for the computation of the manifold, for an overview of such methods see [74].

Appendix F. Alternative coordinate representation
In the main text, we chose to use the parametrization with ρ 1 , ρ 2 and ψ, because this allows for visualization in a familiar cylindrical coordinate system, and because these coordinates have natural interpretations in terms of the distributions of phases in the finite oscillator system: ρ 1 and ρ 2 indicate the degree of synchrony in each population, and ψ defines the mean phase difference between the populations. However, an alternative coordinate representation is possible that better reflects the symmetries inherent to the system, as is discussed here. The equations describing the thermodynamic limit (B.1)-(B.3), before being transformed into polar coordinates, can be rewritten in terms of two complex amplitudes, z k = ρ k e iφ k , k = 1, 2, taking the form and the corresponding equation for z 2 with interchanged indices. This system exhibits a rotational symmetry according to (z 1 , z 2 ) → (z 1 e iφ , z 2 e iφ ).
This symmetry motivates a reduced coordinate system with 0 ≤ |γ| ≤ 1, and δ ∈ [−1, 1], from which we recover the original variables with the intuitive meaning provided that |z 1 | and |z 2 | are non-zero. This system is singular only at z 1 = z 2 = 0 and its geometry can be presented so that the symmetry of exchanging z 1 and z 2 is maintained, i.e., the reflection symmetry δ → − δ γ →γ.
For this parameterization, the fully synchronized states SS 0 and SS π are located at (γ, δ) = (1, 0) and (γ, δ) = (−1, 0), respectively. The invariant rays R 0 and R π are located on the same straight line given by δ = 0 with Im γ = 0. The invariant manifolds S 1 and S 2 are the two paraboloids defined via ±δ = 1 − |γ| 2 . States of interest are then in the region enclosed by these two paraboloids. Sample trajectories are shown in Fig. F1.  Figure F1. Trajectories with initial conditions near R 0 (A, close up in C) and near R π manifold (B), leading to the states SD (red), DS (blue) and SS 0 (yellow). Parameters are A = 0.1 and β = 0.025. The invariant surface manifolds S 1 and S 2 are colored red and blue, respectively.