Attracting Poisson Chimeras in Two-population Networks

Chimera states, i.e., dynamical states composed of coexisting synchronous and asynchronous oscillations, have been reported to exist in diverse topologies of oscillators in simulations and experiments. Two-population networks with distinct intra - and inter-population coupling have served as simple model systems for chimera states since they possess an invariant synchronized manifold, in contrast to networks on a spatial structure. Here, we study dynamical and spectral properties of finite-sized chimeras on two-population networks. First, we elucidate how the Kuramoto order parameter of the finite sized globally coupled two-population network of phase oscillators is connected to that of the continuum limit. These findings suggest that it is suitable to classify the chimera states according to their order parameter dynamics, and therefore we define Poisson and Non-Poisson chimera states. We then perform a Lyapunov analysis of these two types of chimera states which yields insight into the full stability properties of the chimera trajectories as well as of collective modes. In particular, our analysis also confirms that Poisson chimeras are neutrally stable. We then introduce two types of"perturbation"that act as small heterogeneities and render Poisson chimeras attracting: A topological variation via the simplest nonlocal intra-population coupling that keeps the network symmetries, and the allowance of amplitude variations in the globally coupled two-population network, i.e., we replace the phase oscillators by Stuart-Landau oscillators. The Lyapunov spectral properties of chimera states in the two modified networks are investigated, exploiting an ansatz based on the network symmetry-induced cluster pattern dynamics of the finite size network.

Chimera states are a peculiar type of synchronization patterns in homogeneous oscillatory systems 1,2 where regions of synchrony and asynchrony form spontaneously 3 . They were observed in diverse experiments [4][5][6][7][8][9][10] and are believed to be important for certain biological manifestations, such as unihemispheric sleep of some animals or socalled bump-states of neural activity 11,12 . Also from a theoretical point of view, an understanding of chimera states plays an important role, as they mediate between order and disorder [13][14][15] . A detailed analysis of their dynamics is much facilitated with a simple topology, the simplest one consisting of two coupled populations [16][17][18][19][20][21][22][23][24][25][26][27] . For this minimal model, analytical results about the stability and bifurcations of chimera states could be obtained in the continuum limit 16 , and for the case of small populations it was shown the same type of bifurcations exist 17 . Yet, there are still many open questions, some of which we answer in this paper.
The incoherent dynamics of the two-population network depends sensitively on the initial conditions, and on the ensemble size 17,19,21 . In particular, when the initial conditions are obtained from the Poisson kernel, the incoherent motion is simpler than for general initial conditions 27,28 . Our paper is centered around the questions how the chimera states can be classified according to the initial condition and how the dynamics of large and small size populations are linked. Another question we address is how to make the special chimera state with the simpler dynamics of the incoherent oscillators attracting in more realistic situations. Our analysis suggests the definition a) Electronic mail: seungjae.lee@tum.de b) Electronic mail: krischer@tum.de of a Poisson chimera which gives a natural way to classify the chimera states arising from different initial conditions. The main methods employed is Lyapunov analysis [29][30][31][32][33][34] and network symmetry [35][36][37] .

I. INTRODUCTION
Chimera states were first discovered for non-locally coupled phase oscillators on a spatially one-dimensional ring 3 . To obtain a deeper understanding of the dynamics of chimera states, several mathematically more easily tractable models that still exhibit the primary dynamical properties of chimera states have been proposed [13][14][15] . The simplest of them is a network consisting of two populations of identical oscillators. All oscillators within one population are globally coupled to each other with a given intra-population coupling strength, which is the same for both populations. The coupling of the oscillators of different populations is all-to-all as well, but the inter-population coupling strength is different from the intracoupling strength. In a chimera state of such a two-population topology one population oscillates fully synchronously while the other one exhibits incoherent oscillations. The network topology makes sure that the synchronized oscillators live on an invariant sync-manifold, which causes the simpler mathematical accessibility of these chimera states compared to those in other networks, e.g., on the spatially one-dimensional ring 38,39 .
Studies with finite sized populations revealed a strong dependence of the chimera states on initial conditions (ICs) [19][20][21]27,28,42 . The simplest chimera dynamics was obtained when the ICs of the incoherent population were distributed according to the Poisson kernel. However, the chimera states in the identical phase oscillator model were shown to be neutrally stable in many directions 27 . In contrast, when heterogeneous populations were considered, the asymptotic dynamics even for slightly off Poisson kernel ICs was found to be attracting in the long time limit 25,28,43 .
In the following we will term such ICs Poisson initial conditions and abbreviate them with PIC, whereas all other initial conditions are referred to as non-Poisson ICs and abbreviated by n-PIC. In the case of PICs, the chimera states of small-sized populations exhibited pronouncedly different order parameter dynamics from large-sized populations, which has been attributed to finite-size fluctuations 17 . Moreover, for large populations, the numerical simulation suggested that the order parameter becomes indistinguishable from the one predicted by the continuum limit 16,17,19 .
In this paper, we elucidate the origin of both the impact of the initial conditions and of the population size on the chimera dynamics in two-population networks. In particular, we present evidence that there is a continuous change from the small to the large size populations up to the continuum limit. First, we consider the classical two-population network topology with identical Kuramoto-Sakaguchi phase oscillators and global intra-and inter-population coupling ( Fig. 1  (a)). We demonstrate that finite-sized chimeras emerging from PIC live in the neutrally stable Poisson submanifold, which corresponds to the Ott-Antonsen (OA) manifold in the continuum limit and on which the incoherent phase degrees of freedom (DOFs) are distributed according to the Poisson kernel 25,27,44 . To underline the different dynamical characteristics of chimera states arising from PICs and n-PICs, we introduce the concept of a Poisson chimera trajectory and illustrate that what has been so far considered as finite-size fluctuations of small-size chimeras is of fundamentally different nature in the case of Poisson chimeras and of chimeras resulting from n-PICs.
As the next step, we introduce two simple ways that render such Poisson chimera states stable in the sense that they attract nearby trajectories that start from n-PIC or at least evolve towards a close vicinity of the Poisson submanifold 28,42,43,45,46 . The first approach introduces a small topological perturbation of the network structure which leads to the simplest nonlocal intra-population coupling that is represented by a specific adjacency matrix that preserves the network symmetry as the system size increases ( Fig. 1 (b)). Then, we allow for amplitude degrees of freedom (DOFs) by coupling Stuart-Landau oscillators instead of phase oscillators 25,26 . Here, both the global and the nonlocal intra-population network topologies are used.
Our main method to access the properties of the various chimera trajectories is Lyapunov spectral analysis, which yields the spectra of the Lyapunov exponents (LEs) and the covariant Lyapunov vectors (CLVs) [29][30][31][32][33] . The analysis reveals FIG. 1. Schematics of the two-population network topologies considered in this paper. (a) Global intra-and inter-population topology, and (b) global inter-and nonlocal intra-population coupling. Here, only the connections from the first oscillator are fully depicted. The solid connections indicate the intra-population coupling with strength µ, and the dashed one the inter-population connections with strength ν. Note that in the nonlocal intra-population topology, each oscillator is connected to all the other oscillators except of the opposite one.
whether the incoherent oscillator population is attractive or not, as well as the full stability information of the synchronized population. In order to analytically address and approximate the Lyapunov exponents, an approach is introduced that is based on the network symmetry-induced cluster pattern analysis [35][36][37] . Here, we exploit the fact that the finite-sized two-population topology can be viewed as one network that possesses the inherent network symmetries represented by the automorphism group [47][48][49] . The details of the background theories are compiled in Appendices A-B.
The rest of this paper is organized as follows. In Sec. II, we investigate the properties of chimera states of phase oscillators according to the initial conditions and define Poisson chimeras as opposed to non-Poisson chimeras. Furthermore, we discuss the Lyapunov spectral properties of these chimeras. In Sec. III, we consider two ways that render Poisson chimeras attractive; nonlocal topology and amplitude variables. Finally, we summarize the results in Sec. IV.

A. Model and Observable Dynamics
In this section, we consider a set of identical Kuramoto-Sakaguchi (KS) phase oscillators arranged in the twopopulation network topology with global inter-and intrapopulation coupling of different strengths as depicted in Fig. 1 (a). This system is considered to be the simplest model that exhibits chimera states coexisting with a stable complete synchronization state. 16,17 .
Each of the two interacting populations is composed of N phase oscillators. The state of each oscillator is fully described by its phase φ i ∈ T = [−π, π) for i = 1, ..., 2N. The governing equations of the oscillators in the first population are with i = 1, ..., N, and those of the second population are with i = 1, ..., N. Notice that all the oscillators are identical, i.e., they have the same intrinsic frequency ω = 0 and the same Sakaguchi phase-lag parameter α = π/2 − β where β is small enough such that chimera states exist 13,15 . ν and µ are the inter-and intra-population coupling strengths (see Fig. 1).
We rescale time such that µ + ν = 1 and define A = µ − ν. Throughout this work, we set β = 0.08 and A either 0.2 or 0.35. This choice of parameters yields chimera states that are representative of so-called stationary and breathing chimeras, respectively, which are characterized by a stationary and oscillatory behavior of the magnitude of the Kuramoto order parameter with time for large populations 16 . The Kuramoto order parameters for the two populations are defined by r 1 (t)e iΘ 1 (t) = 1 N ∑ N j=1 e iφ j (t) and r 2 (t)e iΘ 2 (t) = 1 N ∑ N j=1 e iφ j+N (t) . Chimera states in a two-population network have one population consisting of perfectly synchronized oscillators with r sync (t) = 1 and the other one being composed of incoherent oscillators with 0 < r incoh (t) < 1. 17 .
Numerical solutions of Eqs. (1)(2) suggest that for each parameter set A and β the chimera trajectories can be divided into two groups, depending on the initial conditions. If the trajectory starts from PICs (the detailed description of ICs will be given in the next section), a chimera trajectory shows a simple, regular motion of the magnitude of the order parameter as depicted in Fig. 2. For large population numbers N as in Fig. 2 (b,d), the magnitude of the order parameter r incoh (t) of chimera states emerging from PICs is either stationary in time ( Fig. 2 (b)) or exhibits simple periodic oscillations ( Fig. 2 (d)), depending on the value of A. These dynamics were termed stationary and breathing chimeras, respectively 16 , and r incoh (t) is virtually indistinguishable from the one of the OA solution in the continuum limit. For small population sizes N, as in Fig. 2 (a,c), r incoh (t) is composed of two contributions: the motion it shows in the case of large N and a superposed, in the case of breathing chimeras secondary, oscillation. Note that throughout this paper, we name each chimera state according to its classification in the continuum limit at the given parameter set for the sake of simplicity. When the chimera trajectory starts from n-PICs, in contrast, r incoh (t) shows a more complicated motion, strongly depending on the given initial conditions (Fig. 3). This initial condition dependence of r incoh (t) has been pointed out previously 16,19,27,42 , and it has led many authors to use rather special initial conditions for their chimera studies. In this work, we will address the initial condition dependence in some detail, and introduce the concept of Poisson and non-Poisson chimeras in the next section. Furthermore, we explaining the stability of both synchronized and incoherent populations with a Lyapunov analysis.

B. Poisson and Non-Poisson Chimeras
As mentioned above, in order to obtain the simple motion of the magnitude of the order parameter as depicted in Fig. 2 and also in Refs. 16 and 17, a specific initial condition has to be used. We coin this initial condition Poisson initial condition (PIC) since the initial incoherent phases are generated from the Poisson kernel that corresponds to the OA manifold in the continuum limit 21,43,44,50 . To obtain PICs, one first has to solve the 2-dimensional Ott-Antonsen reduced equations for the incoherent population, which for the stable stationary chimera state with the parameter set A = 0.2 and β = 0.08 results in ρ 0 = 0.69998 and ϕ 0 = 6.11918, where ϕ 0 = ϕ 1 − ϕ 2 and ϕ i for i = 1, 2 is the OA phase variable for each population, respectively 16 . Then, consider the Poisson kernel a 0 e iφ n + c.c.
where a 0 = ρ 0 e −iϕ 0 , and its inverse cumulative distribution function (inverse CDF). For our finite-size chimeras, we want the initial incoherent phase distribution {φ i+N (0)} N i=1 to be as close as possible to Eq. (3). To obtain such ICs, equally spaced probabilities are used as arguments of the inverse CDF of the Poisson kernel, i.e., N initial phases of the incoherent population are numerically obtained from for i = 1, ..., N. For the synchronized population, the initial phases {φ i (0)} N i=1 are picked from the delta distribution f (1) (φ ) = δ (φ − φ 0 ) which manifests that this population consists of the perfectly synchronized oscillators.
Simulations of the governing Eqs. (1)(2) can also be initiated from an n-PIC. In this work, n-PIC consists of initial phases i=1 that are randomly and independently from each other picked from the uniform distribution within [−π, π) . Note that such initial conditions do not cover the entire manifold of the incoherent oscillator population off the Poisson submanifold but rather only correspond to some subset of the entire manifold corresponding to the incoherent population.
As we have pointed out above, starting from PICs, the magnitude of the order parameter exhibits one of two behaviors, depending on the population size N. For large N, r incoh (t) is virtually indistinguishable from the one of the continuum limit which is a solution of the OA reduced dynamics 16 . For small N, the motion of r incoh (t) is comprised of the main motion close to the OA dynamics superimposed by a regular secondary oscillation. Its clear and regular behavior suggests that the small-size behavior is not just a finite-size fluctuation 17 but rather has a deterministic origin. In the following, we disclose the source of the secondary motion of r incoh (t) of small-size chimeras that start from PIC.
To address the dynamical behavior of the small-size chimeras, we first focus on the stationary chimera states with A = 0.2. Numerical integration of Eqs. (1)(2) with PICs reveals that the instantaneous velocity of each incoherent oscillator {φ i+N (t)} N i=1 is in fact a periodic function, and, furthermore, all instantaneous velocities of the incoherent oscillators have the same functional form and share the same period T . On the level of the instantaneous velocities this behavior is reminiscent of the behavior of the instantaneous phases in a splay state 51-53 (see Fig. 4 (e)). Numerically, the period of the instantaneous velocity T has a value T ≈ 23.48, irrespective of the population size N. Hence, we assume that the instantaneous frequencies of the incoherent oscillators have the form of a splay state such thatφ in the definition of the order parameter, we obtain for all t ∈ R. Thus, r incoh (t) in Eq. (5) is indeed a periodic function and its period τ = T /N is continuously decreasing as N increases. In Fig. 4 (b), the numerical calculations of the period of the order parameter are plotted as a function of N together with the values predicted by Eq. (5). The nearly perfect agreement of both values confirms that the period τ(N) of the order parameter oscillations are indeed decreasing with N according to T /N. Next, we investigate the amplitude of the periodic order parameter of a small-size stationary chimera. As obvious from Fig. 4 (a), the amplitude of r incoh (t) also decreases with increasing N. To explain this, we here consider the Watanabe-Strogatz reduced dynamics ρ 2 (t), Φ 2 (t), and Ψ 2 (t) for the incoherent population 17,54 . These quantities are related to the Kuramoto order parameter according to 17,27,28 r incoh (t)e iΘ incoh (t) = ρ 2 (t)e iΦ 2 (t) γ 2 (ρ 2 , Ψ 2 ) where and {ψ (2) k } N k=1 are the constants of motion, which are determined by the given initial conditions and satisfy three appropriate constraints 28 . For PICs, the constants of motion comply with the uniform distribution ψ Numerical calculations suggest that the values of the radial variable ρ 2 (t) are consistent with the stationary OA radial variable, while exhibiting a very small finite-size oscillation that in this context we can ignore, even for the smallest chimera. Hence, we can assume ρ 2 (t) = ρ 0 in Eq. (7). Then, the Kuramoto order parameter can be rewritten as . The second term in Eq. (8) represents the secondary oscillation of the small-size stationary chimeras. In Fig. 4 (d), the amplitude of the secondary oscillation is plotted as a function of N. It decreases monotonically with N and approaches zero as N → ∞. Thus, the periodic behavior of r incoh (t) gradually disappears with increasing N, such that r incoh (t) → ρ 0 as N → ∞. Regarding the small-size breathing chimera state, r incoh (t) shows the main breathing motion while having the small secondary oscillation along it. It depends on the system size in a similar manner as the stationary chimeras do, namely according to where ρ 2 (t) is no longer a fixed constant but exhibits the main breathing motion (see Fig. 2 (c)). As in the case of the stationary chimeras, the secondary oscillation vanishes for sufficiently large system sizes since ρ(t) < 1 for ∀t ≥ 0, which makes (1 − ρ 2 (t))(−ρ 2 (t)) N → 0 as N → ∞ and the dynamics of the chimera states approach the one of the continuum limit. Our analysis has revealed that both period and amplitude of the secondary oscillation of r incoh (t) continuously decrease as the system size increases. From approx. N 24 on, the secondary oscillation is not discernible anymore. Rather, r incoh (t) displays a motion indistinguishable from the one of the OA dynamics in the continuum limit. We therefore classify chimeras with population sizes N 24 as large-size chimeras, those with N < 24 as small size chimeras. Yet, we would like to point out that there is a continuous change from the small-size to the large size chimeras and eventually up to the OA dynamics in the continuum limit as N → ∞.
On the other hand, when the chimeras started from n-PIC, a non-Poisson initial condition determines nonuniform constants of motion in the WS reduced dynamics. Then the stationary chimera states obtained from a given n-PIC with the same parameter set (A = 0.2 and β =0.08) show incoherent motion that is qualitatively different from the Poisson chimeras and depend on the specific initial conditions used, i.e., on the nonuniform constants of motion. Fig. 3 shows the temporal evolution of the magnitude of the order parameter for n-PICs and otherwise identical parameter values and system sizes as Fig. 2 does for PICs. Clearly, the behavior of r incoh (t) is more complicated in all four cases. In particular, the fluctuations of r incoh (t) do not disappear for the largesize chimeras and the overall motion of r incoh (t) of small-size chimeras is not composed of a superposition of the OA dynamics and the secondary oscillation. This is in line with the observation that the instantaneous velocities of the incoherent oscillators {φ i+N (t)} N i=1 do not form a splay state-like behavior but rather their shapes differ from oscillator to oscillator and the maxima are time-shifted by different amounts (Fig. 4 (f)). Notice that the quasiperiodic chimera states observed in Refs. 27, 28, and 42 are specific examples of non-Poisson chimera trajectories using a specific non-Poisson initial condition, or corresponding nonuniform constants of motion.
Finally, the red distribution in Fig. 4 (g) illustrates that if the chimera trajectory starts from PIC, then the incoherent phase distribution of this chimera state remains in the Poisson kernel as defined in Eq. (3) within an appropriate rotating reference frame. This is confirmed by the observation that the incoherent phases sorted by their magnitude and plotted against its index (normalized to the total number of oscillators) coincide with the inverse CDF of Eq. (3) (Fig. 4 (c), black dots). This observation is consistent with the fact that the OA manifold is invariant under the dynamics in the continuum limit 43,44,50 . For the finite-sized chimeras initially starting from PIC, we can deduce from the splay form ofφ i (t − τ) =φ i+1 (t) that at least at t = nτ for n ∈ N, the phases of the incoherent population are distributed according to the inverse CDF of the Poisson kernel since the splayed phase velocities result in the same constant shift for all the incoherent phases φ i (t − τ) = φ i+1 (t) + W . Beyond that, the numerical results indicate that the finite-sized Poisson submanifold along the chimera state starting from PIC is invariant under the dynamics. For example, let us define E(t) = e iφ (t) 2 − e 2iφ (t) where · is the ensemble average, then for large enough N, E(t) of the chimera trajectory starting from PIC is numerically found to be close to zero (more precisely, E(t) ∼ O(10 −5 )) revealing that the incoherent phases of such chimeras remain in the Poisson kernel. However, the large-size chimeras initiated from n-PICs do not have the incoherent phase distribution that satisfies the Poisson kernel (see Fig. 4 (c,g)), and after a long enough transient time E(t) ∼ O(10 −1 ). Thus, such chimera states initiated from n-PIC should definitely be distinguished from the Poisson chimeras. Notice that the incoherent motion of the breathing chimera with A = 0.35 starting from n-PIC is different from the incoherent motion of the Poisson chimeras, and also depends on the given n-PIC (see Fig. 3 (c-d)).
According to the above results we define a Poisson chimera trajectory in the two-population network topology as follows: A chimera trajectory is a Poisson chimera if the phase DOFs i=1 of a given ensemble of oscillators satisfy the following three dynamical characteristics: The sync-population is perfectly synchronized and invariant. Condition 3. Large-size Poisson chimeras are characterized by an incoherent order parameter being close to the one of the continuum limit, and the small-size Poisson chimeras by an incoherent order parameter whose motion is a superposition of the one of large-size Poisson chimeras and a secondary oscillation that continuously disappears through an increasing frequency and vanishing amplitude as N → ∞.
Chimera states in the two-population network topology that do not fulfill Conditions 1 -3 are termed a non-Poisson chimera trajectory. Note that the stationary Poisson chimera, whether small or large, has the additional property that the instantaneous frequencies of the incoherent oscillators are splayed within its period T such thatφ i (t − j N T ) =φ i+ j (t) for an arbitrary j ∈ {1, ..., N} and for i = N + 1, ..., 2N with φ 2N+1 (t) ≡ φ N+1 (t); however, the breathing chimeras do not.
For each parameter set, one can consider the manifold of the incoherent oscillator population. A state in this manifold can be characterized by (N − 3)-parameter family of invariant subspaces determined by N − 3 constants of motion, based on the WS framework (see Fig. 8 in Ref. 54). The incoherent oscillators of the Poisson chimeras remain in the Poisson kernel, which corresponds to the Poisson submanifold (OA manifold in the continuum limit) in the following denoted by M Poisson and the uniformly distributed constants of motion. However, the non-Poisson chimeras do not have such a property, corresponding to the invariant manifold outside of the Poisson submanifold denoted by M incoh and general non-uniform constants of motion. In Ref. 54, due to the constants of motion, the state for the identical oscillators described by the WS theory is neutrally stable in many directions. In the following, we will show the Lyapunov spectra in order to confirm the neutral stability of chimera states and then give two perturbations that render such chimera states attracting in the following sections.

C. Lyapunov Stability of Poisson and Non-Poisson Chimeras
In this subsection, we investigate the stability of Poisson and non-Poisson chimeras. Therefore, we consider each chimera state as a reference trajectory in phase space and first numerically determine the Lyapunov exponents and then the corresponding covariant Lyapunov vectors. The properties of the resulting Lyapunov spectra are then elucidated using an ansatz based on network symmetry-induced cluster patterns 35,36 . In particular, this method allows us to obtain approximate analytical expressions for the Lyapunov exponents associated with the synchronized population. Further insight into the Lyapunov exponents associated with the incoherent population is obtained from a Watanabe-Strogatz reduction of the dynamics. Finally, we present evidence of the existence of two collective modes. The detailed calculation for the synchronized population based on the network symmetryinduced cluster pattern dynamics is compiled in Appendix. C.
perturb , and the number of zero exponents has increased by 1 to N. These two type of partitions were characteristic for stationary and breathing Poisson chimeras, respectively, and independent of the system size N. trans . The approximate analytical expressions of them are given as is treated as an external forcing field, and {s m } N m=0 are the (coarse-grained) quotient dynamics of the chimera states according to the network cluster patterns discussed in Appendix. C. The transverse Lyapunov exponents in Eq. (9) are all negative and all degenerate, which confirms that the chimera state is stable in all directions transverse to the sync-manifold. Notice that the numerics ensures that −µcosα − ν N Z < 0. It also follows from numerical calculations that the covariant Lyapunov vectors corresponding to the LEs in Eq.
for κ = 2, ..., N where φ ch (t) ∈ T 2N stands for the given chimera trajectory and T φ ch (t) (T 2N ) is the tangent space at the point along such a chimera trajectory. These numerical CLVs = 0 which ascertains that these LEs correspond indeed to LEs transverse to the sync-manifold of the synchronized population.
We also discover in Fig. 5 (a-b) another negative exponent perturb for the synchronized population. The approximated value of it is given as where Z is again considered as an external forcing field. This Lyapunov exponent in fact corresponds to the perturbation along the sync-manifold (compare Eq. (C9)). Note that this LE mainly depends on the collective behavior of the incoherent oscillators {φ i+N (t) = s m (t)|i = m = 1, ..., N} (see Fig. 5 (f)) via the summation term in Eq. (11), i.e., the motion of the incoherent order parameter, and is much closer to zero than the transverse exponents in Eq. (9). The CLV corresponding to Λ Hence, we conclude that all the Lyapunov modes (CLVs) in the synchronized population, both transverse and parallel to it, are stable, and therefore the synchronized manifold is invariant under the evolution of Eqs. (1)(2). Note that the Lyapunov exponents corresponding to the syncpopulation obtained here in Eq. (9) and Eq. (11) are consistent with previous results in Ref. 17. Therein, the authors considered the Jacobian matrix of the synchronized oscillator dynamics by treating the incoherent oscillators as external forcing functions, and then calculated the eigenvalues of the Jacobian matrix for the synchronized oscillators.
All the chimera states in a global two-population network, regardless of the parameters, i.e., also regardless of whether they are of the stationary or breathing type, have the (N − 1)fold degenerate Λ (0) trans,κ for κ = 2, ..., N and Λ (0) perturb since it is dictated by the symmetries of the global network topology and the perfectly synchronized oscillators. Thus, in Fig. 5 (b), the same classes of the sync LEs for the breathing chimera state can be detected.
Next, we turn to the (N − 1)-fold degenerate zero Lyapunov exponents Λ (incoh) zero = 0 and the negative exponent Λ (incoh) ρ < 0 of the stationary Poisson chimeras (Fig. 5 (a)) that are associated with the incoherent oscillators. To better understand their origin, we consider the reduced dynamics according to the Watanabe-Strogatz transformation 27,44,54 .
where a = 1, 2 denotes the population index and ψ (a) i are the constants of motion determined by the initial condition. This transformation leads to the 6-dimensional reduced set of equations 17 for a = 1, 2. The mean-field forcing H a is given by where γ a is defined by the same way in Eq. (6) for each population. The 6-dimensional reduced dynamics in Eq. (13) with the tangent space dynamics along the corresponding chimera reference trajectory (ρ 1 (t) = 1 and ρ 2 (t) < 1) is associated with six Lyapunov exponents which can be determined numerically. In Fig. 5 (c-d) their values are shown versus the index for the same parameters which were used in the calculations of the full Lyaponov spectra depicted in Fig. 5 (a-b).
The results give further insight on the LEs of the incoherent population: the incoherent WS reduced dynamics resides in an invariant subspace of the phase space of the incoherent population that is determined by the N − 3 constants of motion, i.e., by the initial condition 54 (here, PICs and the uniform distribution of the constants of motion consistent with the Poisson submanifold), which yield N − 3 neutral directions, i.e., N − 3 zero LEs. In addition, there are two further zero exponents associated with the incoherent population that come from the two angular variables (Φ 2 , Ψ 2 ) in the reduced dynamics 55 . Hence, we obtain in total N − 1 zero exponents. Apart from these zero LEs, there exists one negative LE that corresponds to the stable fixed point of the radial variable ρ 2 (t) ∼ ρ 0 = const. whose value is determined by the parameter set. (Note that the remaining exponents in the WS reduced dynamics arise from the sync-group and the continuous timeshift symmetry.) Regarding the breathing chimera states, we find N-fold degenerate zero exponents in the incoherent population; an additional zero Lyapunov exponent results from the oscillating nature of the WS radial variable, i.e., the breathing motion of the order parameter of the incoherent population above the Hopf bifurcation 16,17 .

Collective Modes in Poisson Chimeras
As a last step of our analysis of the dynamics of Poisson chimeras, we investigate whether some of the CLVs correspond to collective perturbations, or modes. Therefore, we calculate the time-averaged inverse participation ratio (IPR) for various system sizes according to 32,33 where q = 2 and is close to 1 if the given vector is well localized but close to 1 2N if the vector components spread out through all the oscillators. Therefore, a CLV is a collective mode if IPR (i) (N) ∼ 1 N as N increases, whereas a CLV is localized when IPR (i) (N) ∼ const. as N increases 32,33 .
In Fig. 5 (e-f), the numerically obtained IPRs of the CLVs corresponding to Λ (0) perturb and Λ (incoh) ρ of the stationary Poisson chimera are plotted versus the system size. The proportionality of IPR(N) ∼ 1 N for large N strongly suggests that the corresponding CLVs are indeed Lyapunov collective modes. As discussed above, these modes are related to the incoherent oscillators and affected by the incoherent order parameter motion. This observation is confirmed by our Lyapunov analysis, i.e., by measuring the localization of the covariant Lyapunov vector. We stress that these Lyapunov modes (Λ  obtained from different n-PIC but otherwise identical parameters in the governing equations are depicted in Fig. 6 (a-b) together with the corresponding numerically determined Lyapunov spectra (c-d). In line with our discussion above in Sec. II B, non-Poisson chimera trajectories show different incoherent motions of the order parameter depending on a given n-PIC. In spite of this, since a non-Poisson chimera also lives on the two-population network, there are also (N − 1)-fold degenerate Λ (0) trans,κ for κ = 2, ..., N of the synchronized population given by Eq. (9). Likewise, the numerical CLV analysis confirms these are indeed transverse to the sync-manifold as in Eq. (10). What is different from Poisson chimeras, particularly in the synchronized population, is that the LE arising from the perturbation along the sync-manifold ( Eq. (11)) takes a different value than in Poisson chimeras. This is because Λ (0) perturb strongly depends on the motion of the incoherent oscillators through Z in Eq. (11), which is determined by the initial condition.
Concerning the LEs in the incoherent population, (N − 1)fold degenerate Λ (incoh) zero = 0 are also found from the WS reduced dynamics. However, since Λ (incoh) ρ strongly depends on the constants of motion determined by the non-Poisson initial condition, it also attains a value different from that of a Poisson chimera trajectory, See Fig. 6 (c-d).

III. TWO WAYS TO ATTRACTING POISSON CHIMERA
So far, many authors have observed that a small heterogeneity, e.g., nonidentical natural frequencies or noisy oscillators, makes the dynamics evolve towards at least a close neighborhood of the OA manifold and Poisson submanifold for the continuum limit and finite size system, respectively, and this stabilizing effect has been reported to be a generic consequence of the heterogeneity of the dynamics 21,27,[41][42][43]45,46 . In this section, we study two simple systems with identical oscillator populations that, according to the Lyapunov analysis, possess attracting Poisson chimeras. In the first system, we consider a nonlocal intra-population coupling, in the second one amplitude degrees of freedom of the oscillators, i.e., we employ Stuart-Landau amplitude oscillators rather than phase oscillators.

A. Topological variation: nonlocal intra-population network
While previous studies on nonlocal intra-population networks focused on randomly but systematically constructed topologies and on chimera states in the continuum limit 22,24 , we consider here the simplest regular and finite-sized nonlocal network. This allows us to take advantage of the symmetry of the network. As depicted in Fig. 1(b) the oscillators of each population are arranged on a ring. Compared to the globally coupled intra-population network, each oscillator has one intra-population connection less: it is not connected to the opposite oscillator. For this purpose, we only consider even numbers of the oscillators in each population here. The adjacency matrix of this nonlocal intra-population but global inter-population network is defined as where the i-th oscillator is disconnected to (i+ N 2 )-th oscillator of the same population.
The governing equations of the Kuramoto-Sakaguchi phase oscillators in the nonlocal intra-population topology are for the first population for i = 1, ..., N, and for i = 1, ..., N for the second one. A i j is the adjacency matrix that describes the nonlocal intra-population coupling of each population defined in Eq. (15). Although the nonlocally coupled system does not have a corresponding OA dynamics, we found that as long as we started from PIC ∈ M Poisson , chimera trajectories satisfy the dynamical characteristics of Poisson chimeras as defined in Sec. II B. For this nonlocal Poisson chimera state, the distribution of the incoherent phases remains in a close vicinity of the Poisson submanifold defined by Eq. (3) as the Poisson chimera distributions shown in Fig. 4 (c,g). Additionally, the nonlocal Poisson chimera also show the splay form of the instantaneous frequencies of the incoherent oscillators if A = 0.2. In Fig. 7, the simple motion of the magnitude of the order parameter of nonlocal stationary and breathing Poisson chimera states are depicted. For the parameter A = 0.2, the magnitude of the order parameter has a practically constant value for large size chimeras (slightly different from the global topology), and the small-size chimera displays the clear periodic motion that arises from the splayed instantaneous velocities. For the parameter A = 0.35, the order parameter of the large-size chimera state exhibits a main breathing motion as expected; however, the one of small-size chimeras does not show the main breathing motion superimposed by a secondary oscillation but rather it looks like that of the stationary Poisson chimera state.This might be interpreted as a hint that the nonlocality on the two-population network topology changes the Hopf bifurcation point for the small-size Poisson chimera as described for different non-complete networks in Ref. 22.
provided that Z is treated as an external forcing field. The numerical calculation of the CLVs confirms that the LEs in Eq. (18) are indeed transverse to the sync-manifold since the corresponding CLVs have the form v = 0 for κ = 2, ..., N. Note that as N increases, the gap between the transverse Lyapunov exponents in Eq. (18) is decreasing, and the numerical results in Fig. 8 reflect this fact.
Also, as can be seen in Fig. 8 there is another LE of the synchronized population, which arises from a perturbation along the sync-manifold. This perturbation brings forth the very negative exponent Λ (0) perturb = − ν N Z < 0 that strongly depends on the motion of the incoherent oscillators. Hence, we conclude that in the nonolocal intra-population topology the synchronized population of Poisson chimera states is also stable in both the directions transverse and parallel to the syncmanifold.

Incoherent Population : Paris of Two Near-degenerate Lyapunov Exponents
Next, we focus on the Lyapunov exponents corresponding to the incoherent oscillators. Although we cannot apply directly the Watanabe-Strogatz reduction ( Eq. (13)) in case of the nonlocally coupled oscillators, the classification of the incoherent LEs can be addressed as follows. The quotient dynamics for the incoherent population in Eq. (C13) contains discrete symmetries due to the topology of the nonlocal network (see Fig. 8 (a)). Since each oscillator is disconnected only from the opposite one, two oscillators s m (t) and s m+N/2 (t) are characterized by the same evolution equation. It is also known that such discrete symmetries cause neardegeneracy in the Lyapunov spectrum 34 . Thus, N/2 pairs of two nearly degenerate exponents occur in the incoherent population (see Fig. 8 (b,d)). Therefore, unlike the globally coupled Poisson chimeras, which are neutrally stable, the incoherent population of nonlocal Poisson chimeras is stable, as suggested by the fact that all pairs of the incoherent Lyapunov exponents are definitely negative, except for the two zero exponents which are connected to the continuous symmetries: the phase shift (v ps = (δ φ 0 , ..., δ φ 0 ) where |δ φ 0 | 1) and the time shift (v ts ∝φ ch = f(φ ch )), respectively, which in fact do not affect the stability of the trajectory 33 . For large-size Poisson chimeras in Fig. 8 (c,e), the near-degenerate pairs in the incoherent population are getting closer and closer to one another, until eventually, due to the nonlocal network symmetry, they tend to form two different continuous distributions, one of which consists of obviously negative LEs, whereas the other one consists of two (or some) zero and very slightly negative LEs, corresponding to slow but stable Lyapunov exponents (within our numerical ability).
On the other hand, we can also think of this attractiveness of nonlocal Poisson chimeras due to the heterogeneity of the system. Our nonlocal topology is generated by the least change from the global topology, and hence if we make global the summation term in Eq. (C13), then the disconnecting term due to the nonlocal topology between the two oscillators s m and s m+N/2 should be included in the uncoupled term outside withω m (t) = − µ N sin(s m+N/2 − s m − α) ∼ O(N −1 ). Thus we can interpretω m (t) as a small heterogeneity for the globally coupled incoherent oscillator population. Such a heterogeneity is known to confine the chimeras in a vicinity of the Poisson submanifold. 28,42,43,45,46

C. Dynamical variation: Stuart-Landau oscillators
As the second way to obtain attracting Poisson chimeras, we consider Stuart-Landau (SL) planar oscillators. This twopopulation network of SL oscillators has been studied recently in the continuum limit 25,26 , in which attracting chimeras states have been reported. Here, we consider the finite-sized ensemble and give a full Lyapunov stability analysis which gives further evidence that amplitude DOFs render Poisson chimeras attracting. The amplitude degrees of freedom introduce a small heterogeneity, which is, however, this time self-organized 25,41 .
In an ensemble of Stuart-Landau (SL) oscillators, each oscillator has a phase φ i (t) ∈ [−π, π) and an amplitude r i (t) ∈ (0, ∞) variable. The governing equations are for i = 1, ..., N, which depicts the evolution of the amplitude variables of the oscillators in the first oscillator population, and for i = 1, ..., N, describing the phase dynamics of the SL oscillators in the same population. The governing equations for the second population can also be easily obtained in the same way. In our further study, we fix some parameters: σ = 0.2 and ω = 0. Notice that as ε → 0, the system approaches the evolution equations (1-2) of the phase-only oscillators whose amplitude r i → 1 for all i = 1, ..., 2N 25 .
To study Poisson chimeras of the SL ensemble, we start from the PIC on the phase variables in Eq. (21) in one population, and set the phases of the second population to the same value and all the initial amplitudes in Eq. (20) to r i (0) = 1 for i = 1, ..., 2N (Note that the definition of Poisson chimeras involves only the phase DOFs). The states evolving from such a PIC satisfy all the dynamical properties in the definition of Poisson chimeras for the phase DOFs: one population remains perfectly synchronized, the incoherent phase distribution remains in the Poisson kernel, and finally large-and small-size behavior emerges according to the system size. In particular, the stationary chimera states show the splay form of the instantaneous incoherent frequencies that yield the periodic order parameter for the small-size stationary chimeras. Regarding the amplitude variables, all synchronized oscillators have an amplitude r i (t) = 1 for i = 1, ..., N and the amplitudes of the oscillators in the other population show some distribution with the degree of variation depending on the parameter ε.
For the SL oscillators, the coupling strength ε acts as a bifurcation parameter. For weak coupling strength, i.e., sufficiently small ε (here, we use ε = 0.01) the dynamics are close to the phase-reduced behavior. Hence, the evolution of the order parameter is very close to the one depicted in Fig. 2 for the phase reduced system; r incoh (t) is stationary for A = 0.2 and exhibits a breathing motion for A = 0.35. However, when increasing ε at constant A = 0.2 the stationary chimera undergoes eventually a Hopf bifurcation, giving rise to breathing chimeras, which are observed, e.g., for ε = 0.15, which is in line with findings reported in Ref. 25.

D. Lyapunov analysis on Poisson chimeras of Stuart-Landau oscillators
To study the Lyapunov exponents numerically, we exploit the real-valued coordinates of each Stuart-Landau oscillator 33 . The variables of an SL oscillator can be represented by for k = 1, ..., 2N where a k and b k are real-valued functions of time. Thus, the perturbation vectors in the tangent space are written in the form v (i) = (a 1 , ..., , a N , a N+1 , ..., a 2N , b 1 , ..., b N , b N+1 , ..., b 2N ) ∈ T x ch (t) (R 4N ). This coordinate transformation is a unitary transformation; hence, it can uphold the information on Lyapunov exponents.

Amplitude Degrees of Freedom
In Fig. 10, the numerically obtained Lyapunov spectra of chimera states for strong (ε = 0.1 (a,b) and 0.15 (c,d)) coupling are displayed. The spectra are composed of two parts, which correspond to the phase and amplitude degrees of freedom, respectively. The former are shown in the left column, the latter in the middle one.
First, consider the amplitude DOFs of the synchronized oscillators. They have (N − 1)-fold degenerate strongly negative Lyapunov exponents, which are transverse to the syncmanifold. The approximate values of these Lyapunov exponents are for κ = 2, ..., N (see Eq. (C18)). The numerically obtained CLVs confirm that these Lyapunov exponents are indeed transverse to the sync-manifold as they have the following form (amp,0) κi = 0 for κ = 2, ..., N. In Fig. 10 (b,d), we observe another negative exponent in the synchronized population of the amplitude DOFs caused by the perturbation along the sync-manifold as in Eqs. (C19 -C20). For the analytical value of it one obtains which is a slightly greater Lyapunov exponent than the transverse ones Λ perturb , in line with the numerical observations in Fig. 10.
The numerical CLV analysis reveals that this LE has the form v (inc) j = 0. Hence, we conclude that there is no perturbation direction in the amplitude DOFs, which corresponds to an unstable direction of the synchronized manifold, i.e., the sync-manifold remains invariant under the dynamics since for the sync-population in the amplitude DOFs the CLV modes both transverse and parallel to the sync-manifold are stable. For the other Lyapunov exponents in the amplitude DOFs, we guess that these stable Lyapunov exponents of the amplitude DOFs are linked to the incoherent oscillators through their quotient dynamics in Eq. (C21). Therefore, all the amplitude Lyapunov exponents are strongly negative, and the Poisson chimeras are strongly attracting in all the amplitude DOFs. Note that the amplitude DOFs of the SL ensemble for the weak coupling (ε = 0.01), depicted in Fig. 11, show the same behavior.

Phase Degrees of Freedom
In the phase degrees of freedom, the synchronized oscillators also have the (N − 1)-fold degenerate transverse Lyapunov exponents in Fig. 11 and 10, whose analytical approximate expressions are What makes Poisson chimeras of SL oscillators attractive are the incoherent LEs Λ (incoh) in Fig. 10 (see inset) and Fig. 11 (a,c). In an appropriate rotating reference frame, the quotient governing equations for the incoherent phase DOFs in Eq. (C23) are the same as for the phase-only oscillators in Eqs. (C4-C5) except for the amplitude variables that can be considered as a small self-organized heterogeneityΩ m (t) in the phase governing equations, Eq. (C23) 25 . For the strongly coupled systems with ε = 0.1 and 0.15 as in Fig. 10 (a,c), there are clearly negative Lyapunov exponents in the incoherent phase DOFs. For the stationary chimera (ε = 0.1), we have, in addition, two zero exponents, for the breathing chimera (ε = 0.15), besides the negative exponents, there are three zero exponents, one of which arises from the oscillatory nature of the breathing chimeras. The stable Lyapunov exponents arise due to the amplitude variables in the phase governing equations which present a heterogeneity that, in turn, renders the chimeras attractive 21,25,28,[41][42][43]45,46 . For the weak coupling case ε = 0.01 in Fig. 11 (a,c), the amplitude fluctuations are not that strong (R m ≈ R 0 = 1 for m = 1, ..., N) and as a result, Eq. (C23) can be approximated by Eqs. (C4-C5) like the phase-reduced model, and the Poisson chimeras and their Lyapunov exponents follow patterns similar to the ones obtained for the KS oscillators (see Fig. 11 (a,c) compared to Fig. 5). However, there is still a heterogeneity of the amplitude DOFs in the phase governing equations, and therefore we can expect the LEs to be still slightly negative (stable Lyapunov exponents) in the incoherent phase DOFs. Even for cases where these exponents are very close to zero in our numerical ability, compared to the KS phase-only LEs in Fig. 5, they are slightly decreasing to negative values in the index order which does not occur in the KS phase-only system. Hence, we tentatively conclude that also for weak coupling the stationary chimeras have only two zero LEs, all other exponents are weakly stable. Hence, in all cases, the Poisson chimeras are either at least weakly stable or clearly attracting, compared to the KS phase-only Poisson chimera states because the amplitude variables introduce a self-organized heterogeneity in the phase governing equations. As a consequence, even if starting from n-PIC, the chimera trajectories eventually approach the Poisson submanifold, as we could confirm with numerical simulations.
More than this, we also exploited weakly coupled (ε = 0.01) SL oscillators in the nonlocal intra-population network in order to see whether the Poisson chimeras are also attracting or not. The detailed results on the Lyapunov analysis are compiled in Appendix. D. In this case, the nonlocal topology leads to a stronger negative Lyapunov exponents than the globally coupled SL oscillators rather similar to the phaseonly system in Fig. 8. Therefore, the simultaneous perturbations also cause the Poisson chimeras to evolve towards a close neighbourhood of the Poisson submanifold, i.e. the Poisson chimera trajectory.
Finally, we also investigate whether the system has a Lyapunov collective mode or not by numerically evaluating the IPR function defined in Eq. (14), especially for the case of the breathing chimera ε = 0.15. As can be seen in Fig. 10 (e-g), for at least six Lyapunov modes the IPR shows the tendency to decrease according to IPR (i) (N) ∼ 1 N as the system size N increases. This strongly suggests that these modes, which correspond to the negative exponents in PART 1 in Fig. 10, are collective modes (Note that the stationary chimera ε = 0.1 shows the same collective modes, not shown here). In PART 2 and PART 3 in Fig. 10 for the amplitude DOFs, within our numerically tractable system sizes, at least one Lyapunov mode satisfies the inverse-proportional behavior of the IPR as a function of the system size in each PART, respectively. Consequently, these Lyapunov modes, Λ 2N+1 in PART 2 and Λ 4N in PART 3 are not localized but affect all the oscillators collectively (not restricted only on the incoherent group, but spread out over all the oscillators) and are strongly related to their collective motion in the state space, i.e., they are also Lyapunov collective modes.

IV. CONCLUSION
In this work, we have dealt with chimera states in twopopulation networks of identical oscillators. For the identical Kuramoto-Sakaguchi phase oscillators, the order parameter dynamics of the incoherent oscillator population strongly depends on the initial condition and the population size 17,27 . Once chimeras started from a special initial condition where all the initial phases of one population are in the Poisson kernel, i.e., the Poisson submanifold 41,44,50,54 , the phases remain in the Poisson kernel for all times, and we called this chimera a Poisson chimera. Poisson chimeras show a rather simple motion of the incoherent oscillator population that is virtually indistinguishable from the continuum limit OA solution for sufficiently large population sizes 16 . In contrast, the incoherent motion of a Poisson chimera with a small population size is drastically different from the simple OA dynamics 17 . This difference is not due to finite-size fluctuations, but has a deterministic origin: The magnitude of the order parameter of the incoherent oscillator population shows not only the main motion close to the OA dyanmics but also a superimposed secondary oscillation along the main motion. We demonstrated that this superposed oscillation is a consequence of the fact that the instantaneous frequencies of stationary Poisson chimeras exhibit a splay-form. Furthermore, the splayed distribution of the instantaneous frequencies bring about that the period of the superposed oscillation tends to zero with increasing N while the consideration of the WS global variables revealed how the amplitude of the secondary oscillation disappears with increasing N 17,27,28 . Consequently, our investigations have revealed that and how the order parameter changes continuously from small-size chimeras to large-size chimeras up to the continuum limit, eventually showing the same dynamics as the OA dynamics in the continuum boundary.
In contrast to such Poisson chimeras, the chimeras initialized outside the Poisson submanifold, called in this work non-Poisson chimeras, do not show such a simple order parameter dynamics, regardless of the system size, nor splay-formed instantaneous frequencies of the stationary chimeras, nor does the phase distribution stay in the Poisson kernel. Rather, they show complicated fluctuations along the main motion close to the OA dynamics. This complex, superposed trajectory exists for stationary as well as breathing chimeras, it does not disappear for the large population sizes and in the long time limit, and it depends on the particular initial condition outside the Poisson submanifold, i.e., a set of nonuniform constants of motion 27,28,42 .
In our numerical Lyapunov analysis and also in other previous results 54,55 , the stationary chimera states in twopopulation network with global intra-and inter-population coupling topology, whether it is a Poisson or non-Poisson chimera, are neutrally stable in N − 1 directions. Note that the other negative LE corresponds to the degree of the coherence, i.e., the global WS radial variable. Based on the WS theory, the neutral stability mainly originates from the constants of motion of the system. Any particular parameter set determines the value of the OA radial variable in the continuum limit. The phase DOFs of the neutrally stable Poisson chimeras are dictated by the Poisson initial condition, i.e., uniform constants of motion, and remain in the Poisson kernel 44 .
In contrast, the initial conditions for non-Poisson chimeras correspond to a non-uniform set of constants of motion that cause the different irregular motions of the incoherent oscillators outside the Poisson submanifold according to the different set of motion constants, i.e. the non-Poisson initial condition.
In the next step, we have considered two possibilities that make Poisson chimeras attractive or at least remain in a close vicinity of the Poisson submanifold. We have introduced two 'perturbations' to the Kuramoto-Sakaguchi phase oscillators on the global two-population network: a nonlocal intrapopulation topology and an amplitude degree of freedom, i.e. Stuart-Landau planar oscillators. Previously, many authors showed that the OA manifold in the continuum limit is attracting in the long time limit if the system exhibits some type of heterogeneity 42,43,45 . Considering the WS transformation, it was also shown that a finite-sized system is evolving towards at least a close vicinity of the Poisson submanifold when the system has an suitable heterogeneity, such as nonidentical natural frequencies, noisy oscillators, or experiences a heterogeneous mean-field forcing 28,46 . We have demonstrated that our two perturbations can be thought of as such a small heterogeneity for the incoherent oscillator population. Correspondingly, the Lyapunov analysis has revealed that the systems of nonlocally coupled phase oscillators and globally coupled Stuart-Landau amplitude oscillators have (slightly) negative Lyapunov exponents associated with the incoherent population of phase DOFs, and thus an attracting Poisson chimera trajectory 25,26 : Even when starting from non-Poisson ICs, the chimera trajectory evolved towards the Poisson chimera or to a close neighborhood of it in the long time limit.
As a concluding remark, we note that in real world systems heterogeneities of some type will naturally be present so that the Poisson submanifold becomes at least weakly attracting, which underlines the importance of Poisson chimera states.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are available from the corresponding author upon reasonable request.

ACKNOWLEDGMENTS
The authors would like to thank Maximilian Patzauer, Sindre W. Haugland and Felix P. Kemeth for fruitful discussions. This work has been supported by the Deutsche Forschungsgemeinschaft (project KR1189/18 'Chimera States and Beyond')

Appendix A: Lyapunov Analysis
To study the spectral properties of a chimera trajectory in state space, we perform a Lyapunov analysis. In this appendix, we review some basic concepts; the detailed descriptions can be found in Refs. 29-31, and 34.
First, our governing equations are represented by a set of autonomous ordinary differential equations. In the general vectorial notation, we considerẋ(t) = f(x(t)) with an initial condition x(0) = x 0 where x(t) ∈ R n is the dynamical variable, f is the vector field, and n is the dimension of the state space. A reference trajectory x ref (t) is a solution of the initial value problem, along which we want to study the spectral properties. In our context, it therefore should be a chimera state trajectory. Now we consider the tangent space at each state point along the reference trajectory, wherein the perturbation vector δ x(t) resides, i.e., δ x(t) ∈ T x ref (t) (R n ). Those perturbation vectors are governed by the Jacobian matrix of the vector field, evaluated along the reference trajectory, which can be represented as δẋ(t) = J(t; x ref (t))δ x(t) where the Jacobian matrix is defined by Oseledets' theorem 29,56 tells us that the limits (A1) exist and share the same real positive eigenvalues denoted by µ 1 > µ 2 > ... > µ n (Here, we only consider the nondegenerate case). The forward and backward Oseledets matrices are respectively defined by where stands for the transpose of a matrix and − for transpose and inverse of it. The forward/backward Oseledets matrix probes the future/past dynamics along the given reference trajectory. Those matrices have the eigenspaces spanned by the so-called forward/backward Lyapunov vectors d (i) ± (t). However, these vectors are not co-variant under the dynamics, i.e., it does not bear any information on the local expansion/contraction of the perturbation vectors. Nevertheless, we can construct the Oseledets' splitting that decomposes the tangent space according to the local expansion/contraction behavior along the reference trajectory. We define nested subspaces which construct the Oseledets' splitting in following (Γ (i) (t)) + = n j=i (U ( j) (t)) + and (Γ (i) (t)) − = i j=1 (U ( j) (t)) − where (U ( j) (t)) ± are the eigenspaces of the forward/backward Oseledets matrices spanned by {d 29 . Therefore, we have the decomposition of the tangent space such that where Ω (i) (t) = (Γ (i) (t)) + ∩ (Γ (i) (t)) − is called the Oseledets' splitting. This Oseledets' splitting is covariant under the given dynamics in the sense that Ω (i) (t) = M(t 0 ,t)Ω (i) (t 0 ). The spanning vectors {v (i) (t)} n i=1 of such Oseledets' splittings are called the Covariant Lyapunov Vectors (CLVs) 33 that hold the information on the local expansion/contraction direction of the perturbation vectors since they are norm-independent and also covariant under the dynamics. The exponential rate of such local expansion/contraction along the direction of the CLVs is called Lyapunov Exponents (LEs) and defined by Hence, the Lyapunov exponents characterize the exponential asymptotic growth rate ||M(t 0 ,t)v (i) (t 0 )|| ∼ ||v (i) (t 0 )||exp(Λ i t) and the covariant Lyapunov vectors indicate the stable/unstable directions of the perturbation vectors in the state space 34 .

Appendix B: Network Symmetry Analysis
The two-population topology we consider in the main text can, in fact, be seen as a finite-sized network with 2N nodes. This holds for both the global and nonlocal intrapopulation cases. Furthermore, the discrete network symmetries are represented by the automorphism group of a given network 35,36,47,48 . Recently, many authors have focused on such network symmetries to investigate the dynamics of various kinds of coupled oscillators on a given finite-sized network with abundant discrete symmetries 37,49,51,[57][58][59] . In the following, we exploit the same approach to study the spectral properties of the synchronized population of the chimera states both for the Kuramoto-Sakaguchi phase oscillators and Stuart-Landau amplitude oscillators. In this section, we introduce some important background theories introduced in Refs. 35 and 36. The automorphism group denoted by Aut(G ) of a given network G is a mathematical group consisting of all the automorphisms. An automorphism is a permutation σ of the set of nodes that preserve the adjacency relation among the nodes in the way that A i j = A σ (i)σ ( j) 47 . Consider the group action under a subgroup G ≤ Aut(G ). Then, an orbit partition of a given network G under the subgroup G is a set of orbits defined by ϕ(G, i) = {σ (i)|σ ∈ G} which defines a mathematical partition such that ϕ(G, i) = ϕ(G, j) for all j ∈ ϕ(G, i), and ϕ(G, i) ∩ ϕ(G, j) = / 0 if j / ∈ ϕ(G, i). This partition of a graph can be a candidate of a cluster synchronization (CS) pattern of a given dynamics on the network [35][36][37]49 since each oscillator in the same orbit should receive the same input from the others.
Let us now consider two different types of governing equations, one of which is called here the Pecora-type equation 36,37,49 and the other one the Kuramoto-type equation, which describes diffusively coupled oscillators 35,51,60 : for i = 1, ..., N where x i (t) ∈ R n denotes the dynamical variable, F(x) governs the uncoupled dynamics, H(x) the coupling function, and finally K denotes the coupling constant. For a given candidate of CS pattern, we consider the set of all clusters (orbits) where M is the number of clusters, including trivial clusters that have only one oscillator in it. An associated CS dynamics is described by the coarse-grained variables {s m (t) = x i (t)|i ∈ C m , 1 ≤ m ≤ M} under the quotient adjacency matrixÃ mm = ∑ j∈C m A i j for an arbitrary node i ∈ C m , which is nothing but the number of links from an arbitrary node in C m to all the nodes in C m . Hence, the quotient dynamics of the CS pattern is given bẏ    1 . The cluster-based coordinate transformation can blockdiagonalize a relevant matrix such as an adjacency matrix according to the given cluster pattern, which therefore reveals the spectral properties of the dynamics on each cluster 35,49 .
To study the spectral properties of the dynamics on each cluster, for the moment, we only consider the Pecora-type equation in Eqs. (B1-B2) and use the given CS pattern as a reference trajectory on which we inflict a small deviation. This, then, yields the coupled variational equations for all the clus- for i ∈ C m and DF and DH indicate the Jacobian matrices of the given dynamical functions. Notice that each variational equation in Eq. (B3) is coupled to all the others through the given adjacency matrix. However, if we see this in the cluster-based coordinates by following η is the eigenvalue of the adjacency matrix corresponding to the cluster C m with the eigenvector u (m) κ of the adjacency matrix. From those transversal variational equations, we can investigate the spectral information on the transverse direction of each cluster along our chimera states.
Appendix C: Lyapunov Exponents based on Network Symmetry-induced Cluster Patterns

Kuramoto-Sakaguchi Phase Oscillators
As a first step, we identify the cluster-synchronization (CS) pattern corresponding to the chimera state on the twopopulation network by assigning one of the two populations to the synchronized oscillators, and the other one to the incoherent oscillators. The population of the N perfectly synchronized oscillators can be thought of as just one giant cluster, which we denote by C 0 , whereas each incoherent oscillator in the other population is treated as a trivial cluster denoted by where the first column u √ N , 0, · · · , 0] indicates the sync-manifold direction, D = diag(1, ..., 1) ∈ R N×N indicating the incoherent trivial clusters, each O is a zeromatrix, and P ∈ R N×N−1 representing the directions transverse to C 0 , can be chosen to satisfy orthonormality and transversality such as The cluster-based coordinates decouple the variational equations according to the given CS pattern, as demonstrated in Appendix. B and in the Supplemental Material of Ref. 35. Our case is rather simple since our chimera state has only one nontrivial cluster for the synchronized oscillators. Considering the two-population topology as one large network consisting of 2N nodes with appropriately defined coupling weights, and describing a chimera state by a CS pattern defined above {C m } N m=0 , the Lyapunov exponents corresponding to the synchronized cluster C 0 can be analytically estimated. According to this approach, the governing equation can be written as for i = 1, ..., 2N where the uncoupled dynamics is F(φ ) = − µ N sinα (here, just a constant) and the coupling function is H(x) = sin(x − α). This is nothing but the Kuramototype equation discussed in Eq. (B1). The adjacency matrix B (c) i j ∈ R 2N×2N stands for the complete graph with 2N nodes, and the coupling weights are defined by K i j = µ N if i, j belong to the same population, and K i j = ν N if i, j belong to different populations, respectively, for i, j = 1, ..., 2N. From the CS pattern {C m } N m=0 , the quotient adjacency matrix is given as where A (c) ∈ R N×N is the adjacency matrix of the complete graph with N nodes that describes the global intra-population coupling. Note that the quotient adjacency matrix in Eq. (C3) is an R (N+1)×(N+1) matrix and the index is taken from 0 to N for the sake of simplicity:Ã mm for m, m = 0, 1, ..., N. Therefore, we obtain the (coarse-grained) quotient dynamics corresponding to our chimera pattern from Eq. (B2) with the CS variables denoted by s 0 (t) = φ i (t) (sync., C 0 ) and s m (t) = φ i+N (t) (incoh., C m ) for m = i = 1, ..., N: for each i ∈ C 0 where the deviation around the CS pattern is δ φ i (t) = φ i (t) − s m (t) for i ∈ C m and m = 0, 1, ..., N. Next, we want to obtain Eq. (C6) in the cluster-based coordinate defined in Eq. (C1). The transverse variations can be written as η N ] . Then, the variational equations transverse to the sync-cluster C 0 reaḋ for κ = 2, ..., N. As shown in Ref. 35, the cluster-based coordinates can block-diagonalize the adjacency matrix B (c) according to the CS pattern so that the block corresponding to the sync-cluster C 0 can be represented by the matrix diag(λ κ for κ = 2, ..., N can be chosen to be the eigenvectors of the adjacency matrix 35,61 . Hence, the variational equations transverse to the sync-manifold are given bẏ where the lower-right block corresponds to the sync-cluster C 0 and we obtain λ (0) κ = −1 for all κ for our global intra-and inter-population topology. Hence, if we consider the summation term in Eq. (C8) as an external forcing field 17 , then it gives approximated values of the (N − 1)-fold degenerate transverse LEs in Eq. (9).
To estimate the Lyapunov exponent along the syncmanifold for the synchronized population Λ (0) perturb , the perturbation should be performed along the sync-manifold. This means we obtain the variational equation when the small perturbation s 0 (t) → s 0 (t) + δ s 0 (t) where |δ s 0 (t)| 1 is applied to Eq. (C4).
For the Kuramoto-Sakaguchi phase oscillators in the nonlocal intra-population network, we use the same ansatz introduced above where we treated the chimera state as a CS pattern dynamics. We again start the analysis with the governing equation that, however, now contain the nonlocal adjacency matrix for i = 1, ..., 2N where F(φ ), K i j , and H(x) are the same as defined in Eq. (C2). The matrix B (n) ∈ R 2N×2N , which defines the global inter-and nonlocal intra-population network, is given by where A ∈ R N×N is defined in Eq. (15) and J N ∈ R N×N is the unit matrix whose elements are all 1. In this ansatz, the quotient adjacency matrix is given bỹ wherein the termsÃ (n) 00 = N − 2 andÃ (n) i j = A mm for m, m = 1, ..., N ensuring that the intra-population topology is not global but nonlocal. FromÃ (n) , we obtain the quotient dynamics according to the CS pattern describing our chimeras with the variables s 0 (t) = φ i (t) (sync.) and s m (t) = φ i+N (t) (incoh.) for i = m = 1, ..., N: for the synchronized population, and for κ = 2, ..., N. Therefore, Eq. (C14) yields two different groups of degenerate Lyapunov exponents transverse to the sync-manifold provided that Z is treated as an external forcing field. Also, there is another LE of the synchronized population, which arises from a perturbation along the sync-manifold. Here, a small perturbation s 0 → s 0 + δ s 0 is imposed on Eq. (C12) where |δ s 0 | 1. This perturbation gives Λ (0) sync = − ν N Z < 0 strongly depending on the motion of the incoherent oscillators.

Stuart-Landau Planar Oscillators
Let us consider the spectra corresponding to the amplitude DOFs in more detail. Using the corresponding ansatz as above, the evolution of the amplitude DOFs can be expressed as for i = 1, ..., 2N, where F (amp) (r) = ε −1 (1 − r 2 )r + µ N rcosα and H (amp) (r) = r. Here, we regard the phase variables as external forcing functions, which means we define the coupling weight in Eq. (C16) as K if i, j belong to the same population and K (amp) i j = ν N cos(φ j − φ i − α) if i, j belong to the different populations. This equation is a Pecora-type equation (cf. Eq. (B1)), and the amplitude Lyapunov exponents can be approximated as follows.
According to the chimera CS pattern dynamics introduced in Sec. C 2, we denote the amplitude degrees of freedom by r i (t) = R 0 (t) = 1 for the synchronized population and r i+N (t) = R m (t) for the incoherent one, and, correspondingly, the phase DOFs by s 0 (t) = φ i (t) (sync.) and s m (t) = φ i+N (t) (incoh.) for i = m = 1, ..., N. Then, the quotient dynamics of the amplitude DOFs for the synchronized population with the quotient adjacency matrix in Eq. (C3) is governed by Considering a small deviation around the CS dynamics, i.e., δ r i (t) = r i (t) − R m (t) for i ∈ C m for m = 0, 1, ..., N, we obtain the coupled variational equations as where the last term is zero and since the adjacency matrix is block-diagonalizd in the cluster-based coordinates κ δ κκ for κ = 2, ..., N. Hence, the N − 1 variational equations transversal to the syncmanifold are given bẏ Here, the λ (0) κ = −1 since they are the same as for the global intra-population network (Eq. (C8)). Thus, with Eq. (C18) we obtain the approximate values of the (N − 1)-fold degenerate transverse Lyapunov exponents in the amplitude DOFs as Λ (amp,0) trans,κ ≈ ε −1 (1 − 3R 2 0 ) < 0 in Eq. (23) for κ = 2, ..., N. Next, to estimate the Lyapunov exponent associated with the perturbation along the sync-manifold in the amplitude DOFs, we perform a small perturbation along the syncmanifold R 0 (t) → R 0 (t) + δ R 0 (t) with |δ R 0 | 1 in Eq. (C17) and obtain Hence, it gives a slightly greater Lyapunov exponent than the transverse ones which shows that Λ (amp,0) trans Λ (amp,0) perturb . As for the other negative exponents, we guess that the other stable Lyapunov exponents of the amplitude DOFs are linked to the incoherent oscillators governed by the quotient dynamics in Eq. (B2) for m = 1, ..., N.
Next, we deal with the phase degrees of freedom of the Stuart-Landau oscillators ensemble. Here, we also exploit the network structure with appropriately defined coupling weights. With this approach, the governing equations for the phase DOFs read for i = 1, ..., 2N, where the uncoupled dynamics is governed by F (ph) (φ i ) = −σ r 2 i − µ N sinα and the coupling function is defined as H(x) = sin(x − α). The coupling weights are defined by K Here, η κi δ φ i (t) for κ = 2, ..., N and the deviation along the CS dynamics is δ φ i (t) = φ i (t) − s 0 (t) for i ∈ C 0 , and the eigenvalues λ (0) κ = −1 for all κ. In addition, the LE in the sync population coming from a perturbation along the sync-manifold has the value of Λ (0) perturb = − ν NZ < 0 and is expected to be found in the synchronized phase DOFs.
Furthermore, we can also rationalize the eigenvalue branches of the synchronized oscillators that were already discussed in the continuum limit in Ref. 26 as follows. We again consider the real-valued coordinate of the SL variables in the vector form as x k (t) = (a k (t), b k (t)) ∈ R 2 where a k and b k are defined in Eq. (22). Then, the SL oscillators evolve according to i j and K i j are defined in Eq. (C2), the uncoupled dynamics is governed by and the coupling function is written as for i = 1, ..., 2N. If we also regard the chimera state as a CS pattern dynamics: x i (t) = s 0 (sync.
which gives Λ 1 ∼ −2ε −1 corresponding to the amplitude DOF branch and Λ 2 0 corresponding to the phase DOF branch for the synchronized oscillators. This result and our previous analysis strongly suggest that the negative branch indeed arises from the amplitude DOFs and the near-zero branch comes from the phase DOFs including slow and stable Lyapunov exponents, and both render the Poisson chimeras attracting.
Appendix D: Concurrent dynamical and topological variations: Stuart-Landau oscillators on nonlocal intra-population topology Here, both the topological and dynamical variations are introduced simultaneously. Thus, we consider Stuart-Landau amplitude oscillators in the nonlocal intra-population network topology, and focus on weak coupling with (ε = 0.01). Starting from PIC, we observe chimera states that are similar to those in Sec. III A. Hence, the Poisson chimeras with the parameters A = 0.2 and A = 0.35 follow the similar incoherent dynamics as in Fig. 7.
The Lyapunov analysis for the nonlocal Stuart-Landau oscillators obviously results in the properties dictated by the κ = −2 for κ = N/2+2, ..., N. Also, the negative LE corresponding to the sync-manifold perturbation is given as Λ (0) perturb = − ν NZ < 0, strongly depending on the collective behavior of the incoherent oscillators. Finally, in the incoherent population, we obtain the same N/2 pairs of the two nearly-degenerate exponents that result from the discrete symmetries of the phase governing equations of the Stuart-Landau oscillators. Therefore, the Poisson chimera trajectories of this system are also attracting more strongly than other cases.
Regarding the amplitude DOFs, the N − 1 transverse Lyapunov exponents also show the two different values of the degenerate exponents approximated as  Fig. 12 (b,d)). Then, we expect to find the sync-manifold perturbation exponent of the amplitude DOFs, Λ (amp,0) perturb ≈ ε −1 (1 − 3R 2 0 ) + µ N (N − 1)cosα which is slightly greater than the transverse exponents. As for the other exponents, we only know that they arise from the incoherent governing equations.
Judging from the above observation, we conclude that the Poisson chimera states are definitely attracting. A comparison with the systems that have only one 'perturbation' compared to the globally coupled phase oscillators, i.e. either the nonlocal coupling topology or the amplitude DOF, suggests that the attraction rate of the phase DOF is mainly determined by the non-local network topology.