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 approach based on 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–10} and are believed to be important for certain biological manifestations, such as unihemispheric sleep of some animals or so-called 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–15} A detailed analysis of their dynamics is much facilitated with a simple topology, the simplest one consisting of two coupled populations.^{16–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 that the same type of bifurcations exists.^{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 of a Poisson chimera, which gives a natural way to classify the chimera states arising from different initial conditions. The main methods employed are Lyapunov analysis^{29–34} and network symmetry.^{35–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–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 intra-coupling 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}

This simpler structure has been exploited in numerous studies.^{4,16–27,40} In many of them, the continuum limit was considered.^{16,18,19} Furthermore, in order to address the robustness of chimera states, heterogeneities have been introduced,^{21,23,41} or non-complete networks of oscillators were considered with a static^{22} or time varying^{24} network structure. Besides phase oscillators, also, planar oscillators were studied.^{25,26}

Studies with finite-sized populations revealed a strong dependence of the chimera states on initial conditions (ICs).^{19–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 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.

We will exploit that the dynamics of globally coupled phase oscillators with sinusoidal coupling can be reduced to low-dimensional phase space. Two situations should be distinguished; For finite-sized ensembles, Watanabe and Strogatz (WS)^{45} showed that such a system consisting of $N$ oscillators possesses $N\u22123$ constants of motion, the remaining degrees of freedom describing the global behavior in terms of one radial and two angular variables denoted below as $\rho $, $\Psi $, and $\Phi $, respectively. In the continuum limit of infinitely many oscillators, Ott and Antonsen^{43,46} discovered that for a suitable initial condition, the dynamics lives on a two-dimensional manifold and captures the evolution of the macroscopic order parameter. In Refs. 28 and 44, it was demonstrated that when choosing an appropriate initial condition, namely, PICs, the OA manifold corresponds to the uniformly distributed constants of motion in the WS theory, whereas trajectories starting from n-PICs correspond to a nonuniform distribution of constants of motion in the WS theory.

After considering two-population networks and the relation between the microscopic and macroscopic dynamics, we introduce two simple ways that render Poisson chimera states stable in the sense that they attract nearby trajectories that start from n-PIC and evolve toward close vicinity of the Poisson submanifold.^{28,42,43,47,48} 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 nonlocal intra-population network topologies are used.

Our main method to access the properties of the various chimera trajectories is the Lyapunov spectral analysis, which yields the spectra of the Lyapunov exponents (LEs) and the covariant Lyapunov vectors (CLVs).^{29–33} The analysis reveals 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–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.^{49–51} The details of the background theories are compiled in Appendixes A and 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.

## II. POISSON AND NON-POISSON CHIMERAS

### A. Model and observable dynamics

In this section, we consider a set of identical Kuramoto–Sakaguchi (KS) phase oscillators arranged in the two-population network topology with global inter- and intra-population 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 $\varphi i\u2208T=[\u2212\pi ,\pi )$ for $i=1,\u2026,2N$. The governing equations of the oscillators in the first population are

with $i=1,\u2026,N$ and those of the second population are

with $i=1,\u2026,N$. Notice that all the oscillators are identical; i.e., they have the same intrinsic frequency $\omega =0$ and the same Sakaguchi phase-lag parameter $\alpha =\pi /2\u2212\beta $, where $\beta $ is small enough such that chimera states exist.^{13,15}$\nu $ and $\mu $ are the inter- and intra-population coupling strengths (see Fig. 1). We rescale time such that $\mu +\nu =1$ and define $A=\mu \u2212\nu $. Throughout this work, we set $\beta =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 $r1(t)ei\Theta 1(t)=1N\u2211j=1Nei\varphi j(t)$ and $r2(t)ei\Theta 2(t)=1N\u2211j=1Nei\varphi j+N(t)$. Chimera states in a two-population network have one population consisting of perfectly synchronized oscillators with $rsync(t)=1$ and the other one being composed of incoherent oscillators with $0<rincoh(t)<1$.^{17}

Numerical solutions of Eqs. (1) and (2) suggest that for each parameter set $A$ and $\beta $, 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 Sec. II B), 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 Figs. 2(b) and 2(d), the magnitude of the order parameter $rincoh(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 $rincoh(t)$ is virtually indistinguishable from one of the OA solutions in the continuum limit. For small population sizes $N$, as in Figs. 2(a) and 2(c), $rincoh(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, $rincoh(t)$ shows a more complicated motion, strongly depending on the given initial conditions (Fig. 3). This initial condition dependence of $rincoh(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 Sec. II B. Furthermore, we explain 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 the 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,46} To obtain PICs, one first has to solve the two-dimensional Ott–Antonsen reduced equations for the incoherent population, which for the stable stationary chimera state with the parameter set $A=0.2$ and $\beta =0.08$ results in $\rho 0=0.69998$ and $\phi 0=6.11918$, where $\phi 0=\phi 1\u2212\phi 2$ and $\phi i$ for $i=1,2$ is the OA phase variable for each population, respectively.^{16} Then, consider the Poisson kernel

where $a0=\rho 0e\u2212i\phi 0$ and its inverse cumulative distribution function (inverse CDF). For our finite-size chimeras, we want the initial incoherent phase distribution ${\varphi i+N(0)}i=1N$ 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,\u2026,N$. For the synchronized population, the initial phases ${\varphi i(0)}i=1N$ are picked from the delta distribution $f(1)(\varphi )=\delta (\varphi \u2212\varphi 0)$, which manifests that this population consists of the perfectly synchronized oscillators.

Simulations of the governing equations (1) and (2) can also be initiated from an n-PIC. In this work, n-PIC consists of initial phases ${\varphi i(0)}i=12N$ that are randomly and independently from each other picked from the uniform distribution within $[\u2212\pi ,\pi )$.

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$, $rincoh(t)$ is virtually indistinguishable from one of the continuum limits, which is a solution of the OA reduced dynamics.^{16} For small $N$, the motion of $rincoh(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 $rincoh(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) and (2) with PICs reveals that the instantaneous velocity of each incoherent oscillator ${\varphi \u02d9i+N(t)}i=1N$ 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^{52–54} [see Fig. 4(e)]. Numerically, the period of the instantaneous velocity $T$ has a value $T\u224823.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 $\varphi \u02d9i(t\u2212jNT)=\varphi \u02d9i+j(t)$ for an arbitrary $j\u2208{1,\u2026,N}$, which gives $\varphi i(t\u22121NT)=\varphi i+1(t)+W$ for $i=N+1,\u2026,2N$ with $\varphi 2N+1\u2261\varphi N+1$ where $W\u2208R$ is a common constant. Plugging the expression for ${\varphi i+N(t)}i=1N$ in the definition of the order parameter, we obtain

for all $t\u2208R$. Thus, $rincoh(t)$ in Eq. (5) is indeed a periodic function, and its period $\tau =T/N$ is continuously decreasing as $N$ increases. In Fig. 4(b), the simulations 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 $\tau (N)$ of the order parameter oscillations is 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 $rincoh(t)$ also decreases with increasing $N$. To explain this, we here consider the Watanabe–Strogatz reduced dynamics $\rho 2(t)$, $\Phi 2(t)$, and $\Psi 2(t)$ for the incoherent population.^{17,45} These quantities are related to the Kuramoto order parameter according to^{17,27,28}

where

and ${\psi k(2)}k=1N$ 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 $\psi k(2)=2\pi kN$ for $k=1,\u2026,N$.^{17,27} For $\gamma 2$, one can obtain^{28,42}

Simulations suggest that the values of the radial variable $\rho 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 $\rho 2(t)=\rho 0$ in Eq. (7). Then, the Kuramoto order parameter can be rewritten as

where $O(t;\rho 0,\Psi 2)=eiN(2\pi N\u2212\Psi 2(t))1\u2212(\u2212\rho 0)NeiN(2\pi N\u2212\Psi 2(t))$. 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\u2192\u221e$. Thus, the periodic behavior of $rincoh(t)$ gradually disappears with increasing $N$ such that $rincoh(t)\u2192\rho 0$ as $N\u2192\u221e$. Regarding the small-size breathing chimera state, $rincoh(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 $\rho 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 $\rho 2(t)<1$ for $\u2200t\u22650$, which makes $(1\u2212\rho 2(t))(\u2212\rho 2(t))N\u21920$ as $N\u2192\u221e$, and the dynamics of the chimera states approach one of the continuum limits.

Our analysis has revealed that both period and amplitude of the secondary oscillation of $rincoh(t)$ continuously decrease as the system size increases. From approximately $N\u227324$, the secondary oscillation is not discernible anymore. Alternatively, $rincoh(t)$ displays a motion indistinguishable from one of the OA dynamics in the continuum limit. We, therefore, classify chimeras with population sizes $N\u227324$ as large-size chimeras and 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\u2192\u221e$.

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 $\beta =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. Figure 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 $rincoh(t)$ is more complicated in all four cases. In particular, the fluctuations of $rincoh(t)$ do not disappear for the large-size chimeras, and the overall motion of $rincoh(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 ${\varphi \u02d9i+N(t)}i=1N$ do not form a splay state-like behavior but rather their shapes differ from an 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,46} For the finite-sized chimeras initially starting from PIC, we can deduce from the splay form of $\varphi \u02d9i(t\u2212\tau )=\varphi \u02d9i+1(t)$ that at least at $t=n\tau $ for $n\u2208N$, 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 $\varphi i(t\u2212\tau )=\varphi 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)=|\u27e8ei\varphi (t)\u27e92\u2212\u27e8e2i\varphi (t)\u27e9|$, where $\u27e8\u22c5\u27e9$ 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)\u2248O(10\u22125)$] 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 Figs. 4(c) and 4(g)] and after a long enough transient time $E(t)\u2248O(10\u22121)$. 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 Figs. 3(c) and 3(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 ${\varphi i(t)}i=12N$ of a given ensemble of oscillators satisfy the following three dynamical characteristics:

The sync-population is perfectly synchronized and invariant.

The incoherent phase distribution of Poisson chimeras remains in the Poisson kernel or at least in close vicinity of the Poisson submanifold.

Large-size Poisson chimeras are characterized by an incoherent order parameter being close to one of the continuum limits and the small-size Poisson chimeras by an incoherent order parameter whose motion is a superposition of one of large-size Poisson chimeras and a secondary oscillation that continuously disappears through an increasing frequency and vanishing amplitude as $N\u2192\u221e$.

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 $\varphi \u02d9i(t\u2212jNT)=\varphi \u02d9i+j(t)$ for an arbitrary $j\u2208{1,\u2026,N}$ and for $i=N+1,\u2026,2N$ with $\varphi 2N+1(t)\u2261\varphi 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 the $(N\u22123)$-parameter family of invariant subspaces determined by $N\u22123$ constants of motion based on the WS framework (see Fig. 8 in Ref. 45). 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 $MPoisson$ 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 $Mincoh$ and general non-uniform constants of motion. In Ref. 45, 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 Sec. II C.

### 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 the 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 approach 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 symmetry-induced cluster pattern dynamics is compiled in Appendix C.

In Fig. 5, panels (a) and (b) display numerically determined Lyapunov spectra along Poisson chimera trajectories for stationary (a) and breathing (b) chimeras. Details about the numerical method used can be found in Refs. 29–33 and are summarized in Appendix A. The Lyapunov spectrum of stationary chimera states [Fig. 5(a)] is composed of four groups of exponents: (i) $(N\u22121)$-fold degenerate zero exponents denoted by $\Lambda zero(incoh)=0$, (ii) $(N\u22121)$-fold degenerate negative exponents denoted by $\Lambda trans(0)$, and (iii) and (iv) two individual negative LEs, denoted by $\Lambda perturb(0)$ and $\Lambda \rho (incoh)$. The spectrum obtained from a breathing Poisson chimera trajectory [Fig. 5(b)] exhibits a similar partition of the exponents; however, there is just one individual non-degenerate negative exponent, $\Lambda perturb(0)$, 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$.

#### 1. Synchronized population: $\Lambda trans(0)$ and $\Lambda perturb(0)$

In Fig. 5, there are $(N\u22121)$-fold degenerate transverse Lyapunov exponents denoted by $\Lambda trans(0)$. The approximate analytical expressions of them are given as

for $\kappa =2,\u2026,N$ (indicating the indices for the $N\u22121$ transverse directions) where $Z=\u2211m\u2032=1Ncos(sm\u2032\u2212s0\u2212\alpha )$ is treated as an external forcing field and ${sm}m=0N$ 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 $\u2212\mu cos\alpha \u226a\u2212\nu NZ<0$. It also follows from the numerics that the covariant Lyapunov vectors corresponding to the LEs in Eq. (9) have the form

for $\kappa =2,\u2026,N$ where $\varphi ch(t)\u2208T2N$ stands for the given chimera trajectory and $T\varphi ch(t)(T2N)$ is the tangent space at the point along such a chimera trajectory. These numerical CLVs have $\u2211i=1Nv\kappa i(trans)=0$, which ascertains that these LEs correspond indeed to LEs transverse to the sync-manifold of the synchronized population.

We also discover in Figs. 5(a) and 5(b) another negative exponent $\Lambda perturb(0)$ for the synchronized population. The approximated value of it is given as

where $Z$ is again considered 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 ${\varphi i+N(t)=sm(t)|i=m=1,\u2026,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 $\Lambda perturb(0)$ has the form $vperturb(0)=[v,\u2026,v,v1(incoh),\u2026,vN(incoh)]\u22a4\u2208T\varphi ch(t)(T2N)$, where $\u2211j=1Nvj(incoh)\u22600$. 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) and (2). Note that the Lyapunov exponents corresponding to the sync-population obtained here in Eqs. (9) and (11) are consistent with the 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\u22121)$-fold degenerate $\Lambda trans,\kappa (0)$ for $\kappa =2,\u2026,N$ and $\Lambda perturb(0)$ 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.

#### 2. Incoherent population: $\Lambda zero(incoh)$ and $\Lambda \rho (incoh)$

Next, we turn to the $(N\u22121)$-fold degenerate zero Lyapunov exponents $\Lambda zero(incoh)=0$ and the negative exponent $\Lambda \rho (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,45}

where $a=1,2$ denotes the population index and $\psi i(a)$ are the constants of motion determined by the initial condition. This transformation leads to the six-dimensional reduced set of equations,^{17}

for $a=1,2$. The mean-field forcing $Ha$ is given by

where $\gamma a$ is defined by the same way in Eq. (6) for each population. The six-dimensional reduced dynamics in Eq. (13) with the tangent space dynamics along the corresponding chimera reference trajectory [$\rho 1(t)=1$ and $\rho 2(t)<1$] is associated with six Lyapunov exponents, which can be determined numerically. In Figs. 5(c) and 5(d), their values are shown vs the index for the same parameters, which were used in the calculations of the full Lyapunov spectra depicted in Figs. 5(a) and 5(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\u22123$ constants of motion, i.e., by the initial condition^{45} (here, PICs and the uniform distribution of the constants of motion consistent with the Poisson submanifold), which yield $N\u22123$ neutral directions, i.e., $N\u22123$ zero LEs. In addition, there are two further zero exponents associated with the incoherent population that come from the two angular variables ($\Phi 2$, $\Psi 2$) in the reduced dynamics.^{55} Hence, we obtain in total $N\u22121$ zero exponents. Apart from these zero LEs, there exists one negative LE that corresponds to the stable fixed point of the radial variable $\rho 2(t)\u2248\rho 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 time-shift 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}

#### 3. 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 $IPR(i)\u2208[(2N)\u22121,1]$ and $vj(i)$ is the $j$th component of the CLV $v(i)\u2208T\varphi ch(t)(T2N)$ corresponding to a given exponent denoted by $\Lambda i(N)$ defined in Eq. (A2) for $i=1,\u2026,2N$. By definition, $IPR(i)(N)$ is close to $1$ if the given vector is well localized but close to $12N$ if the vector components spread out through all the oscillators. Therefore, a CLV is a collective mode if $IPR(i)(N)\u223c1N$ as $N$ increases, whereas a CLV is localized when $IPR(i)(N)\u223cconst.$ as $N$ increases.^{32,33}

In Figs. 5(e) and 5(f), the numerically obtained IPRs of the CLVs corresponding to $\Lambda perturb(0)$ and $\Lambda \rho (incoh)$ of the stationary Poisson chimera are plotted vs the system size. The proportionality of $IPR(N)\u223c1N$ 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 ($\Lambda perturb(0)$ and $\Lambda \rho incoh$) are collective (non-localized) throughout all the oscillators and not restricted to the incoherent oscillator population.

#### 4. Non-Poisson chimeras

Finally, we turn to the Lyapunov exponents of the non-Poisson chimera trajectories that start from a given n-PIC. Two examples of the temporal evolution of magnitude of the order parameter of non-Poisson chimera trajectories that were obtained from different n-PIC but otherwise identical parameters in the governing equations are depicted in Figs. 6(a) and 6(b) together with the corresponding numerically determined Lyapunov spectra (c) and (d). In line with our discussion 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\u22121)$-fold degenerate $\Lambda trans,\kappa (0)$ for $\kappa =2,\u2026,N$ of the synchronized population given by Eq. (9). Likewise, the numerical CLV analysis confirms that 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 $\Lambda perturb(0)$ 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\u22121)$-fold degenerate $\Lambda zero(incoh)=0$ are also found from the WS reduced dynamics. However, since $\Lambda \rho (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 Figs. 6(c) and 6(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 toward at least a close neighborhood of the OA manifold and the Poisson submanifold for the continuum limit and the 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–43,47,48} 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 and 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

The governing equations of the Kuramoto–Sakaguchi phase oscillators in the nonlocal intra-population topology are for the first population

for $i=1,\u2026,N$ and

for $i=1,\u2026,N$ for the second one. $Aij$ 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 $\u2208MPoisson$, 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 Figs. 4(c) and 4(g). Additionally, the nonlocal Poisson chimera also shows 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 is 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.

### B. Lyapunov analysis of Poisson chimeras in the nonlocal intra-population network

#### 1. Synchronized population: $\Lambda trans(0)$ and $\Lambda perturb(0)$

For the Poisson chimeras on the nonlocal topology, the Lyapunov spectrum of the nonlocal Poisson chimeras is qualitatively different from one of the global Poisson chimeras, as can be seen in Fig. 8. There are $N\u22121$ transverse Lyapunov exponents consisting of two different values. This splitting of the values of $\Lambda trans(0)$ is due to the fact that the transversal variational equations include two different eigenvalues of the adjacency matrix corresponding to the same synchronized cluster according to the nonlocal network symmetry [compare Eqs. (C14) and (C15)]. The analytical approximate expressions of the $N\u22121$ transverse Lyapunov exponents to the sync-manifold are

provided that $Z$ is treated as an external forcing field. The simulation 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\kappa (0)=[v\kappa 1(trans),\u2026,v\kappa N(trans),0,\u2026,0]\u22a4\u2208T\varphi ch(t)(T2N)$, while $\u2211i=1Nv\kappa i(trans)=0$ for $\kappa =2,\u2026,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 $\Lambda perturb(0)=\u2212\nu NZ<0$ that strongly depends on the motion of the incoherent oscillators. Hence, we conclude that in the nonlocal intra-population topology, the synchronized population of Poisson chimera states is also stable in both the directions transverse and parallel to the sync-manifold.

#### 2. Incoherent population: Pairs 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 the 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 $sm(t)$ and $sm+N/2(t)$ are characterized by the same evolution equation. It is also known that such discrete symmetries cause near-degeneracy in the Lyapunov spectrum.^{34} Thus, $N/2$ pairs of two nearly identical exponents occur in the incoherent population [see Figs. 8(b) and 8(d)]. Therefore, unlike the globally coupled Poisson chimeras, which are neutrally stable, the incoherent population of nonlocal Poisson chimeras is stable (see Fig. 9), as suggested by the fact that all pairs of the incoherent Lyapunov exponents are definitely negative except for the two zero exponents that are connected to the continuous symmetries: the phase shift [$vps=(\delta \varphi 0,\u2026,\delta \varphi 0)\u22a4$ where $|\delta \varphi 0|\u226a1$] and the time shift [$vts\u221d\varphi \u02d9ch=f(\varphi ch)$], respectively, which, in fact, do not affect the stability of the trajectory.^{33} For large-size Poisson chimeras in Figs. 8(c) and 8(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 $sm$ and $sm+N/2$ should be included in the uncoupled term outside the summation and Eq. (C13) becomes

with $\omega ~m(t)=\u2212\mu Nsin(sm+N/2\u2212sm\u2212\alpha )\u2248O(N\u22121)$. Thus, we can interpret $\omega ~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,47,48}

### C. Dynamical variation: Stuart–Landau oscillators

As the second way to obtain attracting Poisson chimeras, we consider Stuart–Landau (SL) planar oscillators. This two-population network of SL oscillators has been studied recently in the continuum limit,^{25,26} in which attracting chimera 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 $\varphi i(t)\u2208[\u2212\pi ,\pi )$ and an amplitude $ri(t)\u2208(0,\u221e)$ variable. The governing equations are

for $i=1,\u2026,N$, which depicts the evolution of the amplitude variables of the oscillators in the first oscillator population and

for $i=1,\u2026,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: $\sigma =0.2$ and $\omega =0$. Notice that as $\epsilon \u21920$, the system approaches the evolution equations (1) and (2) of the phase-only oscillators whose amplitude $ri\u21921$ for all $i=1,\u2026,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 $ri(0)=1$ for $i=1,\u2026,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 $ri(t)=1$ for $i=1,\u2026,N$, and the amplitudes of the oscillators in the other population show some distribution with the degree of variation depending on the parameter $\epsilon $.

For the SL oscillators, the coupling strength $\epsilon $ acts as a bifurcation parameter. For weak coupling strength, i.e., sufficiently small $\epsilon $ (here, we use $\epsilon =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; $rincoh(t)$ is stationary for $A=0.2$ and exhibits a breathing motion for $A=0.35$. However, when increasing $\epsilon $ at constant $A=0.2$, the stationary chimera undergoes eventually a Hopf bifurcation, giving rise to breathing chimeras, which are observed, e.g., for $\epsilon =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,\u2026,2N$, where $ak$ and $bk$ are real-valued functions of time. Thus, the perturbation vectors in the tangent space are written in the form $v(i)=(a1,\u2026,aN,aN+1,\u2026,a2N,b1,\u2026,bN,bN+1,\u2026,b2N)\u22a4\u2208Txch(t)(R4N)$. This coordinate transformation is a unitary transformation; hence, it can uphold the information on Lyapunov exponents. Such a coordinate transformation is useful to investigate the covariant Lyapunov vectors. In addition, we also calculated the LEs in the original coordinate system (phase and amplitude) and obtained the same result. This fact and some analytical considerations in Ref. 26 and in Appendix C2 allow the phase and amplitude LEs to be distinguished.

#### 1. Amplitude degrees of freedom

In Fig. 10, the numerically obtained Lyapunov spectra of chimera states for strong [$\epsilon =0.1$ (a) and (b) and $0.15$ (c) and (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 and the latter in the middle one.

First, consider the amplitude DOFs of the synchronized oscillators. They have $(N\u22121)$-fold degenerate strongly negative Lyapunov exponents, which are transverse to the sync-manifold. The approximate values of these Lyapunov exponents are

for $\kappa =2,\u2026,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:

where $\u2211i=1Na\kappa i(amp,0)=\u2211i=1Nb\kappa i(amp,0)=0$ for $\kappa =2,\u2026,N$. In Figs. 10(b) and 10(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) and (C20). For the analytical value of it, one obtains

which is a slightly greater Lyapunov exponent than the transverse ones $\Lambda trans(amp,0)\u2272\Lambda perturb(amp,0)$, in line with the numerical observations in Fig. 10. The numerical CLV analysis reveals that this LE has the form $vperturb(amp,0)=(a,\u2026,a,a1(inc),\u2026,aN(inc),b,\u2026,b,b1(inc),\u2026,bN(inc))\u22a4\u2208Txch(t)(R4N)$, where $a,b$$\u2208R$ are constant and $\u2211j=1Naj(inc)\u22600$ and $\u2211j=1Nbj(inc)\u22600$. 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 ($\epsilon =0.01$), depicted in Fig. 11, show the same behavior.

#### 2. Phase degrees of freedom

In the phase degrees of freedom, the synchronized oscillators also have the $(N\u22121)$-fold degenerate transverse Lyapunov exponents in Figs. 10 and 11, whose analytical approximate expressions are

for $\kappa =2,\u2026,N$, where $Z~=\u2211m\u2032=1NRm\u2032R0cos(sm\u2032\u2212s0\u2212\alpha )$ should be considered as an external forcing field. In addition, the LE in the sync group coming from a perturbation along the sync-manifold has the value of $\Lambda perturb(0)=\u2212\nu NZ~<0$ and is expected to be found in the synchronized phase DOFs. The numerical CLV analysis also confirms that $\Lambda trans,\kappa (0)$ and $\Lambda perturb(0)$ associated with the synchronized population are indeed transverse and parallel to the sync-manifold, respectively.

What makes Poisson chimeras of SL oscillators attractive are the incoherent LEs $\Lambda (incoh)$ in Fig. 10 (see the inset) and Figs. 11(a) and 11(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) and (C5) except for the amplitude variables that can be considered a small self-organized heterogeneity $\Omega ~m(t)$ in the phase governing equations, Eq. (C23).^{25} For the strongly coupled systems with $\epsilon =0.1$ and $0.15$ as in Figs. 10(a) and 10(c), there are clearly negative Lyapunov exponents in the incoherent phase DOFs. For the stationary chimera ($\epsilon =0.1$), we have, in addition, two zero exponents for the breathing chimera ($\epsilon =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–43,47,48} For the weak coupling case $\epsilon =0.01$ in Figs. 11(a) and 11(c), the amplitude fluctuations are not that strong ($Rm\u2248R0=1$ for $m=1,\u2026,N$), and as a result, Eq. (C23) can be approximated by Eqs. (C4) and (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 Figs. 11(a) and 11(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, and 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 ($\epsilon =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 stronger negative Lyapunov exponents than the globally coupled SL oscillators rather similar to the phase-only system in Fig. 8. Therefore, the simultaneous perturbations also cause the Poisson chimeras to evolve toward a close neighborhood of the Poisson submanifold, i.e., the Poisson chimera trajectory (see Fig. 9).

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 $\epsilon =0.15$. As can be seen in Figs. 10(e)–10(g), for at least six Lyapunov modes, the IPR shows the tendency to decrease according to $IPR(i)(N)\u223c1N$ 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 $\epsilon =0.1$ shows the same collective modes, not shown here). In PARTs 2 and 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, $\Lambda 2N+1$ in PART 2 and $\Lambda 4N$ in PART 3, are not localized but affect all the oscillators collectively (not restricted only on the incoherent group, but also 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 two-population 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–46} 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 dynamics 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 neither 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. Alternatively, 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,^{45,55} the stationary chimera states in a two-population network with global intra- and inter-population coupling topology, whether it is a Poisson or non-Poisson chimera, are neutrally stable in $N\u22121$ 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 intra-population 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,47} Considering the WS transformation, it was also shown that a finite-sized system is evolving toward at least a close vicinity of the Poisson submanifold when the system has a suitable heterogeneity, such as nonidentical natural frequencies, noisy oscillators, or experiences a heterogeneous mean-field forcing.^{28,48} 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 toward 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.

## 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 (DFG) (Project No. KR1189/18 “Chimera States and Beyond”).

## DATA AVAILABILITY

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

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### APPENDIX A: LYAPUNOV ANALYSIS

To study the spectral properties of a chimera trajectory in the 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 $x\u02d9(t)=f(x(t))$ with an initial condition $x(0)=x0$, where $x(t)\u2208Rn$ is the dynamical variable, $f$ is the vector field, and $n$ is the dimension of the state space. A reference trajectory $xref(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 $\delta x(t)$ resides; i.e., $\delta x(t)\u2208Txref(t)(Rn)$. Those perturbation vectors are governed by the Jacobian matrix of the vector field, evaluated along the reference trajectory, which can be represented as $\delta x\u02d9(t)=J(t;xref(t))\delta x(t)$ where the Jacobian matrix is defined by $(J)ij=\u2202x\u02d9i\u2202xj|xref(t)$. From this, we consider the fundamental matrix solution such that $O\u02d9(t)=J(t;xref(t))O(t)$ with $O(0)=In$; this solution defines $the tangent linear propagator$, $M(t0,t)=O(t)O\u22121(t0)$, of the perturbation vector from a given point in time point to the future time so that $\delta x(t)=M(t0,t)\delta x(t0)$.^{30}

Oseledets’ theorem^{29,56} tells us that the limits (A1) exist and share the same real positive eigenvalues denoted by $\mu 1>\mu 2>\cdots >\mu n$ (Here, we only consider the nondegenerate case). The forward and backward Oseledets matrices are, respectively, defined by

where $\u22a4$ stands for the transpose of a matrix and $\u2212\u22a4$ 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\xb1(i)(t)$. However, these vectors are not covariant under the dynamics; i.e., it does not bear any information on the local expansion/contraction of the perturbation vectors. Nevertheless, we can construct Oseledets’ splitting that decomposes the tangent space according to the local expansion/contraction behavior along the reference trajectory. We define nested subspaces, which construct Oseledets’ splitting in the following way: $(\Gamma (i)(t))+=\u2a01j=in(U(j)(t))+$ and $(\Gamma (i)(t))\u2212=\u2a01j=1i(U(j)(t))\u2212$, where $(U(j)(t))\xb1$ are the eigenspaces of the forward/backward Oseledets matrices spanned by ${d\xb1(j)(t)}j=1n$.^{29} Therefore, we have the decomposition of the tangent space such that $Txref(t)(Rn)=\u2a01j=1n\Omega (j)(t)$, where $\Omega (i)(t)=(\Gamma (i)(t))+\u2229(\Gamma (i)(t))\u2212$ is called Oseledets’ splitting. This Oseledets’ splitting is covariant under the given dynamics in the sense that $\Omega (i)(t)=M(t0,t)\Omega (i)(t0)$. The spanning vectors ${v(i)(t)}i=1n$ 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

for $U(t0)\u2208(\Gamma (i)(t0))+\u2216(\Gamma (i+1)(t0))+$ where the nested subspaces are $Rn=(\Gamma (1)(t))+\u2283(\Gamma (2)(t))+\u2283\cdots \u2283(\Gamma (n)(t))+$. Hence, the Lyapunov exponents characterize the exponential asymptotic growth rate $||M(t0,t)v(i)(t0)||\u223c||v(i)(t0)||exp(\Lambda it)$, 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 intra-population cases. Furthermore, the discrete network symmetries are represented by the automorphism group of a given network.^{35,36,49,50} 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,51,52,57–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 the 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 $\sigma $ of the set of nodes that preserve the adjacency relation among the nodes in the way that $Aij=A\sigma (i)\sigma (j)$.^{49} Consider the group action under a subgroup $G\u2264Aut(G)$. Then, an orbit partition of a given network $G$ under the subgroup $G$ is a set of orbits defined by $\phi (G,i)={\sigma (i)|\sigma \u2208G}$, which defines a mathematical partition such that $\phi (G,i)=\phi (G,j)$ for all $j\u2208\phi (G,i)$ and $\phi (G,i)\u2229\phi (G,j)=\u2205$ if $j\u2209\phi (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–37,51} 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,51} and the other one the Kuramoto-type equation, which describes diffusively coupled oscillators,^{35,52,60}

for $i=1,\u2026,N$ where $xi(t)\u2208Rn$ 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 a CS pattern, we consider the set of all clusters (orbits) ${\phi (i,G)}i=1N={Cm}m=1M$, 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 ${sm(t)=xi(t)|i\u2208Cm,1\u2264m\u2264M}$ under the quotient adjacency matrix $A~mm\u2032=\u2211j\u2208Cm\u2032Aij$ for an arbitrary node $i\u2208Cm$, which is nothing but the number of links from an arbitrary node in $Cm$ to all the nodes in $Cm\u2032$. Hence, the quotient dynamics of the CS pattern is given by

for $m=1,\u2026,M$.

The set of $N$-dimensional orthonormal vectors ${u\kappa (m)}\kappa =1|Cm|$ for $m=1,\u2026,M$ called the $cluster-based coordinates$ is defined by the following rules:^{35} (i) $u\kappa i(m)=0$ if $i\u2209Cm$, (ii) for $\kappa =1$, all the nonzero elements of $u1(m)$ should be $1/|Cm|$ that defines the cluster sync-manifold, and (iii) the other vectors ${u\kappa (m)}\kappa =2|Cm|$ are mutually orthogonal and also to $u1(m)$. The cluster-based coordinate transformation can block-diagonalize 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,51}

To study the spectral properties of the dynamics on each cluster, for the moment, we only consider the Pecora-type equation in Eqs. (B1) and (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 clusters,

for $i=1,\u2026,N$, where $\delta xi(t)=xi(t)\u2212sm(t)$ for $i\u2208Cm$ 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 $\eta \kappa (m)=\u2211i\u2208Cmu\kappa i(m)\delta xi$ for $m=1,\u2026,M$ and $\kappa =2,\u2026,|Cm|$, where $\eta \kappa (m)$ for $\kappa \u22652$ represents the perturbation of the transverse direction to the cluster $Cm$, we get the variational equation for that cluster, independent of the other clusters, provided that the given cluster $Cm$ is non-intertwined with the others [for non-intertwined clusters, see Ref. 36 or for independently synchronizable cluster sets (ISC sets), see the supplementary material in Ref. 35]. Therefore, the $|Cm|\u22121$ transversal variational equations of the cluster $Cm$ both for the Pecora-type and the Kuramoto-type are given by^{35}

for $\kappa =2,\u2026,|Cm|$, where $\lambda \kappa (m)$ is the eigenvalue of the adjacency matrix corresponding to the cluster $Cm$ with the eigenvector $U\kappa (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

#### 1. Kuramoto–Sakaguchi phase oscillators

As a first step, we identify the cluster-synchronization (CS) pattern corresponding to the chimera state on the two-population 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 $C0$, whereas each incoherent oscillator in the other population is treated as a trivial cluster denoted by $Cm$ with $m=1,\u2026,N$. This gives us the corresponding cluster-based coordinates $U\u22a4=[u1(0),u1(1),\u2026,u1(N),u2(0),\u2026,uN(0)]$ for the chimera pattern.^{35} Here, $u1(0)$ indicates the direction along the synchronized cluster $C0$ of the chimera state so that $u1j(0)=1N$ for $j\u2208C0$ and $u1j(0)=0$ if $j\u2209C0$. For the transverse directions, we obtain $\u2211j\u2208C0u\kappa j(0)=0$ and $u\kappa j(0)=0$ if $j\u2209C0$ for $\kappa =2,\u2026,N$. Finally, for the incoherent trivial clusters, we have $u1j(m)=1$ if $j\u2208Cm$ and $u1j(m)=0$ otherwise, with $m=1,\u2026,N$. Note that all the cluster-based coordinate vectors should be mutually orthonormalized. An example of a possible candidate of the cluster-based coordinates is^{61}

where the first column

The cluster-based coordinates decouple the variational equations according to the given CS pattern, as demonstrated in Appendix B and in the supplementary 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 ${Cm}m=0N$, the Lyapunov exponents corresponding to the synchronized cluster $C0$ can be analytically estimated. According to this approach, the governing equation can be written as

for $i=1,\u2026,2N$ where the uncoupled dynamics is $F(\varphi )=\u2212\mu Nsin\alpha $ (here, just a constant) and the coupling function is $H(x)=sin(x\u2212\alpha )$. This is nothing but the Kuramoto-type equation discussed in Eq. (B1). The adjacency matrix $Bij(c)\u2208R2N\xd72N$ stands for the complete graph with $2N$ nodes, and the coupling weights are defined by $Kij=\mu N$ if $i,j$ belong to the same population and $Kij=\nu N$ if $i,j$ belong to different populations, respectively, for $i,j=1,\u2026,2N$. From the CS pattern ${Cm}m=0N$, the quotient adjacency matrix is given as

where

for the synchronized cluster ($C0$) where the quotient adjacency matrix $A~00=N\u22121$ and $A~0m\u2032=1$ for $m\u2032=1,\u2026,N$, and $H(0)=\u2212sin\alpha $. The quotient governing equations of the $N$ trivial clusters ($C1,\u2026,CN$) for the incoherent population read

for $m=1,\u2026,N$. From the quotient dynamics, we consider the variational equations of the synchronized oscillators around the CS pattern as

for each $i\u2208C0$ where the deviation around the CS pattern is $\delta \varphi i(t)=\varphi i(t)\u2212sm(t)$ for $i\u2208Cm$ and $m=0,1,\u2026,N$. Next, we want to obtain Eq. (C6) in the cluster-based coordinate defined in Eq. (C1). The transverse variations can be written as $\eta \kappa (0)(t)=\u2211i\u2208C0u\kappa i(0)\delta \varphi i(t)$ with $U=[u1(0),u1(1),\u2026,u1(N),u2(0),\u2026,uN(0)]\u22a4$. Then, the variational equations transverse to the sync-cluster $C0$ read

for $\kappa =2,\u2026,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 $C0$ can be represented by the matrix $diag(\lambda 2(0),\lambda 3(0),\u2026,\lambda N(0))\u2208R(N\u22121)\xd7(N\u22121)$ and the off-diagonal blocks are zero. This, in turn, means that the last term in Eq. (C7) should be zero and $\u2211i\u2208C0\u2211k\u2208C0u\kappa i(0)Bik(c)u\kappa \u2032k(0)=\lambda \kappa (0)\delta \kappa \kappa \u2032$ for $\kappa ,\kappa \u2032=2,\u2026,N$ where $\lambda \kappa (0)$ are the eigenvalues of the adjacency matrix since $u\kappa (0)$ for $\kappa =2,\u2026,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 by

for $\kappa =2,\u2026,N$. Notice that for the global intra- and inter- population network, the eigenvalues in Eq. (C8)$\lambda \kappa (0)=\u22121$ for all $\kappa =2,\u2026,N$. Consider as an example the system with $N=4$. Its block-diagonalized adjacency matrix reads^{61}

where the lower-right block corresponds to the sync-cluster

^{17}then it gives approximated values of the $(N\u22121)$-fold degenerate transverse LEs in Eq. (9).

To estimate the Lyapunov exponent along the sync-manifold for the synchronized population $\Lambda perturb(0)$, the perturbation should be performed along the sync-manifold. This means that we obtain the variational equation when the small perturbation $s0(t)\u2192s0(t)+\delta s0(t)$ where $|\delta s0(t)|\u226a1$ is applied to Eq. (C4),

Then, we obtain Eq. (11) provided that $Z=\nu N\u2211m\u2032=1Ncos(sm\u2032\u2212s0\u2212\alpha )$ is regarded as an external forcing function.

For the Kuramoto–Sakaguchi phase oscillators in the nonlocal intra-population network, we use the same approach 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,\u2026,2N$, where $F(\varphi )$, $Kij$, and $H(x)$ are the same as defined in Eq. (C2). The matrix $B(n)\u2208R2N\xd72N$, which defines the global inter- and nonlocal intra-population network, is given by

where

wherein the terms

for the synchronized population and

where $\omega ~m(t)=\u2212\mu Nsin(sm+N/2\u2212sm\u2212\alpha )$ with $A~0m(n)=1$ and $A~m0(n)=N$ for $m=1,\u2026,N$ for the incoherent trivial clusters.

As seen in Sec. III B, there are $N\u22121$ transverse Lyapunov exponents consisting of two different values. This splitting of the values of $\Lambda trans(0)$ is due to the two different eigenvalues of the nonlocal adjacency matrix. As clear from Eq. (B4), one has to consider the eigenvalues of the adjacency matrix associated with the cluster-based vector, which are the eigenvectors of the adjacency matrix $U\kappa (m)$, to obtain the transverse variational equations. For the global topology discussed above, these eigenvalues $\lambda \kappa (0)=\u22121$ are the same for $\kappa =2,\u2026,N$. In contrast, the nonlocal adjacency matrix has two different eigenvalues: $\lambda \kappa (0)=0$ for $\kappa =2,\u2026,N/2+1$ and $\lambda \kappa (0)=\u22122$ for $\kappa =N/2+2,\u2026,N$.^{61} This leads to two different variational equations with the same method in Eqs. (C7) and (C8),

for $\kappa =2,\u2026,N$. Therefore, Eq. (C14) yields two different groups of degenerate Lyapunov exponents transverse to the sync-manifold

provided that $Z$, defined in Eq. (9), 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 $s0\u2192s0+\delta s0$ is imposed on Eq. (C12) where $|\delta s0|\u226a1$. This perturbation gives $\Lambda sync(0)=\u2212\nu NZ<0$ strongly depending on the motion of the incoherent oscillators.

#### 2. Stuart–Landau planar oscillators

Let us consider the spectra corresponding to the amplitude DOFs in more detail. Using the corresponding approach as above, the evolution of the amplitude DOFs can be expressed as

for $i=1,\u2026,2N$, where $F(amp)(r)=\epsilon \u22121(1\u2212r2)r+\mu Nrcos\alpha $ and $H(amp)(r)=r$. Here, we regard the phase variables as external forcing functions, which means that we define the coupling weight in Eq. (C16) as $Kij(amp)=\mu Ncos(\varphi j\u2212\varphi i\u2212\alpha )$ if $i,j$ belong to the same population and $Kij(amp)=\nu Ncos(\varphi j\u2212\varphi i\u2212\alpha )$ 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 1, we denote the amplitude degrees of freedom by $ri(t)=R0(t)=1$ for the synchronized population and $ri+N(t)=Rm(t)$ for the incoherent one and, correspondingly, the phase DOFs by $s0(t)=\varphi i(t)$ (sync.) and $sm(t)=\varphi i+N(t)$ (incoh.) for $i=m=1,\u2026,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., $\delta ri(t)=ri(t)\u2212Rm(t)$ for $i\u2208Cm$ for $m=0,1,\u2026,N$, we obtain the coupled variational equations as

for each $i\u2208C0$ and $Cm\u2032m=cos(sm\u2032\u2212sm\u2212\alpha )$ for $m,m\u2032=0,\u2026,N$. Then, viewing these in the cluster-based coordinates with $\xi \kappa (0)(t)=\u2211i\u2208C0u\kappa i(0)\delta ri(t)$ for $\kappa =2,\u2026,N$, the transversal variational equations in Eq. (B4) read

where the last term is zero and since the adjacency matrix is block-diagonalizd in the cluster-based coordinates $\u2211i\u2208C0\u2211k\u2208C0u\kappa i(0)Bik(c)u\kappa \u2032k(0)=\lambda \kappa (0)\delta \kappa \kappa \u2032$ for $\kappa =2,\u2026,N$. Hence, the $N\u22121$ variational equations transversal to the sync-manifold are given by

Here, $\lambda \kappa (0)=\u22121$ 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\u22121)$-fold degenerate transverse Lyapunov exponents in the amplitude DOFs as $\Lambda trans,\kappa (amp,0)\u2248\epsilon \u22121(1\u22123R02)<0$ in Eq. (23) for $\kappa =2,\u2026,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 sync-manifold $R0(t)\u2192R0(t)+\delta R0(t)$ with $|\delta R0|\u226a1$ in Eq. (C17) and obtain

Hence, it gives a slightly greater Lyapunov exponent than the transverse ones

which shows that $\Lambda trans(amp,0)\u2272\Lambda perturb(amp,0)$.

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,\u2026,N$.

Next, we deal with the phase degrees of freedom of the Stuart–Landau oscillator 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,\u2026,2N$, where the uncoupled dynamics is governed by $F(ph)(\varphi i)=\u2212\sigma ri2\u2212\mu Nsin\alpha $ and the coupling function is defined as $H(x)=sin(x\u2212\alpha )$. The coupling weights are defined by $Kij(ph)=\mu Nrjri$ if $i,j$ belong to the same population and $Kij(ph)=\nu Nrjri$ otherwise, provided that the amplitude variables are treated as external forcing functions. Thus, the resulting equation is the Kuramoto-type equation of Eq. (B2). From the quotient adjacency matrix defined in Eq. (C3), the quotient dynamics of the synchronized and incoherent populations in phase DOFs are obtained as

where $\Omega ~m(t)=\u2212\sigma Rm2(t)$ for $m=1,\u2026,N$. The quotient dynamics reveal that the phase DOFs of the SL oscillator ensemble in the synchronized population also have $(N\u22121)$-fold degenerate transverse Lyapunov exponents with

for $\kappa =2,\u2026,N$, where $Z~=\u2211m\u2032=1NRm\u2032R0cos(sm\u2032\u2212s0\u2212\alpha )$ should be considered an external forcing field, which follows from the transversal variational equations [Eq. (B4)],

Here, $\eta \kappa (0)(t)=\u2211i\u2208C0u\kappa i(0)\delta \varphi i(t)$ for $\kappa =2,\u2026,N$; the deviation along the CS dynamics is $\delta \varphi i(t)=\varphi i(t)\u2212s0(t)$ for $i\u2208C0$ and the eigenvalues $\lambda \kappa (0)=\u22121$ for all $\kappa $. In addition, the LE in the sync-population coming from a perturbation along the sync-manifold has the value of $\Lambda perturb(0)=\u2212\nu 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 $xk(t)=(ak(t),bk(t))\u22a4\u2208R2$ where $ak$ and $bk$ are defined in Eq. (22). Then, the SL oscillators evolve according to

for $i=1,\u2026,2N$, where $Bij(c)$ and $Kij$ are defined in Eq. (C2), and the uncoupled dynamics is governed by

and the coupling function is written as

for $i=1,\u2026,2N$. If we also regard the chimera state as a CS pattern dynamics, $xi(t)=s0$ (sync.) and $xi+N(t)=sm(t)$ (incoh.) for $i=m=1,\u2026,N$, then the variational equations transversal to the synchronized cluster $C0$ in the cluster-based coordinates are given by

for $\kappa =2,\u2026,N$ where the Jacobians of the dynamical functions in Eqs. (C27) and (C28) read

and

Since for the synchronized SL oscillators we have $rkei\varphi k=ei\varphi 0=12(a0+ib0)$ for $k=1,\u2026,N$, we can rewrite the transversal variational equations in the following form:

for all the directions transverse to the sync-manifold. Notice that the eigenvalues of the adjacency matrix $\lambda \kappa (0)=\u22121$ for all $\kappa $. If we consider $\varphi 0$ as an external forcing function, then the eigenvalues of the matrix $Jtrans(0)$ are

which gives $\Lambda 1\u223c\u22122\epsilon \u22121$ corresponding to the amplitude DOF branch and $\Lambda 2\u22720$ 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 (see Sec. III A and Fig. 1) and dynamical (see Sec. III C) variations are introduced simultaneously. Thus, we consider Stuart–Landau amplitude oscillators in the nonlocal intra-population network topology and focus on weak coupling with ($\epsilon =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 given nonlocal topology of the network. From the same method discussed in Appendixes. A– C, we obtain the two different values of the degenerate transverse Lyapunov exponents in the synchronized group of phase DOFs,

where $\lambda \kappa (0)=0$ for $\kappa =2,\u2026,N/2+1$ and $\lambda \kappa (0)=\u22122$ for $\kappa =N/2+2,\u2026,N$. Also, the negative LE corresponding to the sync-manifold perturbation is given as $\Lambda perturb(0)=\u2212\nu 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\u22121$ transverse Lyapunov exponents also show the two different values of the degenerate exponents approximated as

since for the nonlocal network $\lambda \kappa (0)=0$ for $\kappa =2,\u2026,N/2+1$ and $\lambda \kappa (0)=\u22122$ for $\kappa =N/2+2,\u2026,N$ [distinguished by the gray dashed line in Figs. 12(b) and 12(d)]. Then, we expect to find the sync-manifold perturbation exponent of the amplitude DOFs, $\Lambda perturb(amp,0)\u2248\epsilon \u22121(1\u22123R02)+\mu N(N\u22121)cos\alpha $, 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 non-local coupling topology or the amplitude DOF, suggests that the attraction rate of the phase DOF is mainly determined by the non-local network topology.

## REFERENCES

*Kudos 2009, Equitable, PA 2009*(2009).