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 systems1,2 where regions of synchrony and asynchrony form spontaneously.3 They were observed in diverse experiments4–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 analysis29–34 and network symmetry.35–37 

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 static22 or time varying24 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.

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 the opposite one.

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 the opposite one.

Close modal

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 N3 constants of motion, the remaining degrees of freedom describing the global behavior in terms of one radial and two angular variables denoted below as ρ, Ψ, and Φ, respectively. In the continuum limit of infinitely many oscillators, Ott and Antonsen43,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.

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 ϕiT=[π,π) for i=1,,2N. The governing equations of the oscillators in the first population are

(1)

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

(2)

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 r1(t)eiΘ1(t)=1Nj=1Neiϕj(t) and r2(t)eiΘ2(t)=1Nj=1Neiϕ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 β, 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.

FIG. 2.

The magnitudes of Kuramoto order parameters r(t) of the coherent and incoherent populations of chimera states in the two-population network starting from PICs after transients have died out (t105). For each figure, the gray solid line indicates the order parameter for the perfectly synchronized population [r(t)=1] and the black solid line the incoherent population [r(t)<1]: (a) and (b) Stationary chimera states with A=0.2 and (c) and (d) breathing chimera states with A=0.35 for the system sizes N=6 (left) and N=60 (right), respectively.

FIG. 2.

The magnitudes of Kuramoto order parameters r(t) of the coherent and incoherent populations of chimera states in the two-population network starting from PICs after transients have died out (t105). For each figure, the gray solid line indicates the order parameter for the perfectly synchronized population [r(t)=1] and the black solid line the incoherent population [r(t)<1]: (a) and (b) Stationary chimera states with A=0.2 and (c) and (d) breathing chimera states with A=0.35 for the system sizes N=6 (left) and N=60 (right), respectively.

Close modal
FIG. 3.

The magnitudes of Kuramoto order parameters r(t) of the coherent and incoherent populations of chimera states in the two-population network starting from n-PICs after transients have died out (t105). For each figure, the gray solid line indicates the order parameter for the perfectly synchronized population [r(t)=1] and the black solid line the incoherent population [r(t)<1]: (a) and (b) A=0.2 (for which with PICs, stationary chimeras are obtained) and (c) and (d) A=0.35 (for which with PICs, breathing chimeras are obtained) for the system sizes N=6 (left) and N=60 (right), respectively.

FIG. 3.

The magnitudes of Kuramoto order parameters r(t) of the coherent and incoherent populations of chimera states in the two-population network starting from n-PICs after transients have died out (t105). For each figure, the gray solid line indicates the order parameter for the perfectly synchronized population [r(t)=1] and the black solid line the incoherent population [r(t)<1]: (a) and (b) A=0.2 (for which with PICs, stationary chimeras are obtained) and (c) and (d) A=0.35 (for which with PICs, breathing chimeras are obtained) for the system sizes N=6 (left) and N=60 (right), respectively.

Close modal

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 β=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

(3)

where a0=ρ0eiφ0 and its inverse cumulative distribution function (inverse CDF). For our finite-size chimeras, we want the initial incoherent phase distribution {ϕ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

(4)

for i=1,,N. For the synchronized population, the initial phases {ϕi(0)}i=1N are picked from the delta distribution f(1)(ϕ)=δ(ϕϕ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 {ϕi(0)}i=12N that are randomly and independently from each other picked from the uniform distribution within [π,π).

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 fluctuation17 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 {ϕ˙i+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 state52–54 [see Fig. 4(e)]. Numerically, the period of the instantaneous velocity T has a value T23.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 ϕ˙i(tjNT)=ϕ˙i+j(t) for an arbitrary j{1,,N}, which gives ϕi(t1NT)=ϕi+1(t)+W for i=N+1,,2N with ϕ2N+1ϕN+1 where WR is a common constant. Plugging the expression for {ϕi+N(t)}i=1N in the definition of the order parameter, we obtain

(5)

for all tR. Thus, rincoh(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 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 τ(N) of the order parameter oscillations is indeed decreasing with N according to T/N.

FIG. 4.

(a) Oscillations of the magnitude of the order parameter for A=0.2 and different system sizes: N=4 (red), N=8 (blue), N=16 (green), and N=32 (black). (b) Period of rincoh(t) as determined numerically (red) and predicted from Eq. (5) as a function of the system size N. (c) Snapshot of the sorted incoherent phases in the numerical order with N=100 as a function of the rescaled index after a time t106 for a Poisson chimera (black dots), a non-Poisson chimera (gray diamonds), and the theoretical curve of the inverse CDF of the Poisson kernel (red solid curve). (d) Magnitude of the secondary oscillation as a function of system size. (e) and (f) Instantaneous frequencies of the incoherent oscillators of the system N=8 for a Poisson chimera (e) and a non-Poisson chimera (f). (g) Snapshot of the incoherent phase distribution for a Poisson chimera (red) and the non-Poisson chimera (blue) for t106 with N=100 oscillators. Each solid line indicates the theoretical Poisson kernel curve corresponding to ρ0 within an appropriate rotating frame.

FIG. 4.

(a) Oscillations of the magnitude of the order parameter for A=0.2 and different system sizes: N=4 (red), N=8 (blue), N=16 (green), and N=32 (black). (b) Period of rincoh(t) as determined numerically (red) and predicted from Eq. (5) as a function of the system size N. (c) Snapshot of the sorted incoherent phases in the numerical order with N=100 as a function of the rescaled index after a time t106 for a Poisson chimera (black dots), a non-Poisson chimera (gray diamonds), and the theoretical curve of the inverse CDF of the Poisson kernel (red solid curve). (d) Magnitude of the secondary oscillation as a function of system size. (e) and (f) Instantaneous frequencies of the incoherent oscillators of the system N=8 for a Poisson chimera (e) and a non-Poisson chimera (f). (g) Snapshot of the incoherent phase distribution for a Poisson chimera (red) and the non-Poisson chimera (blue) for t106 with N=100 oscillators. Each solid line indicates the theoretical Poisson kernel curve corresponding to ρ0 within an appropriate rotating frame.

Close modal

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 ρ2(t), Φ2(t), and Ψ2(t) for the incoherent population.17,45 These quantities are related to the Kuramoto order parameter according to17,27,28

where

(6)

and {ψ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 ψk(2)=2πkN for k=1,,N.17,27 For γ2, one can obtain28,42

(7)

Simulations 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

(8)

where O(t;ρ0,Ψ2)=eiN(2πNΨ2(t))1(ρ0)NeiN(2πNΨ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. Thus, the periodic behavior of rincoh(t) gradually disappears with increasing N such that rincoh(t)ρ0 as N. 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 ρ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 ρ2(t)<1 for t0, which makes (1ρ2(t))(ρ2(t))N0 as N, 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 N24, 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 N24 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.

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. 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 {ϕ˙i+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 ϕ˙i(tτ)=ϕ˙i+1(t) that at least at t=nτ for nN, 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)=|eiϕ(t)2e2iϕ(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(105)] 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)O(101). 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 {ϕi(t)}i=12N of a given ensemble of oscillators satisfy the following three dynamical characteristics:

Condition 1

The sync-population is perfectly synchronized and invariant.

Condition 2

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

Condition 3

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.

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(tjNT)=ϕ˙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 the (N3)-parameter family of invariant subspaces determined by N3 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.

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) (N1)-fold degenerate zero exponents denoted by Λzero(incoh)=0, (ii) (N1)-fold degenerate negative exponents denoted by Λtrans(0), and (iii) and (iv) two individual negative LEs, denoted by Λperturb(0) and Λρ(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, Λ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.

FIG. 5.

(a) and (b) Lyapunov spectra of the full dynamics of Poisson chimera states with N=12 for A=0.2 (stationary chimeras) (a) and A=0.35 (breathing chimeras) (b). For the meaning of the Λs, see the text. (c) and (d) Lyapunov spectra for the six-dimensional Watanabe–Strogatz reduced dynamics of the chimera states in (a) and (b), respectively. The exponents marked by the black dashed lines indicate the LE corresponding to the radial WS variable. (e) and (f) Inverse participation ratio of the covariant Lyapunov vectors of the stationary chimera in (a) corresponding to Λρ(incoh) and Λperturb(0), respectively, as a function of system size N. The gray dashed lines indicate IPR(i)(N)1N.

FIG. 5.

(a) and (b) Lyapunov spectra of the full dynamics of Poisson chimera states with N=12 for A=0.2 (stationary chimeras) (a) and A=0.35 (breathing chimeras) (b). For the meaning of the Λs, see the text. (c) and (d) Lyapunov spectra for the six-dimensional Watanabe–Strogatz reduced dynamics of the chimera states in (a) and (b), respectively. The exponents marked by the black dashed lines indicate the LE corresponding to the radial WS variable. (e) and (f) Inverse participation ratio of the covariant Lyapunov vectors of the stationary chimera in (a) corresponding to Λρ(incoh) and Λperturb(0), respectively, as a function of system size N. The gray dashed lines indicate IPR(i)(N)1N.

Close modal

1. Synchronized population: Λtrans(0) and Λperturb(0)

In Fig. 5, there are (N1)-fold degenerate transverse Lyapunov exponents denoted by Λtrans(0). The approximate analytical expressions of them are given as

(9)

for κ=2,,N (indicating the indices for the N1 transverse directions) where Z=m=1Ncos(sms0α) 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 μcosανNZ<0. It also follows from the numerics that the covariant Lyapunov vectors corresponding to the LEs in Eq. (9) have the form

(10)

for κ=2,,N where ϕch(t)T2N stands for the given chimera trajectory and Tϕch(t)(T2N) is the tangent space at the point along such a chimera trajectory. These numerical CLVs have i=1Nvκ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 Λperturb(0) for the synchronized population. The approximated value of it is given as

(11)

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 {ϕi+N(t)=sm(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 Λperturb(0) has the form vperturb(0)=[v,,v,v1(incoh),,vN(incoh)]Tϕch(t)(T2N), where j=1Nvj(incoh)0. 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 (N1)-fold degenerate Λtrans,κ(0) for κ=2,,N and Λ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: Λzero(incoh) and Λρ(incoh)

Next, we turn to the (N1)-fold degenerate zero Lyapunov exponents Λzero(incoh)=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,45

(12)

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

(13)

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

where γ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 [ρ1(t)=1 and ρ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 N3 constants of motion, i.e., by the initial condition45 (here, PICs and the uniform distribution of the constants of motion consistent with the Poisson submanifold), which yield N3 neutral directions, i.e., N3 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 N1 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 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 to32,33

(14)

where q=2 and IPR(i)[(2N)1,1] and vj(i) is the jth component of the CLV v(i)Tϕch(t)(T2N) corresponding to a given exponent denoted by Λi(N) defined in Eq. (A2) for i=1,,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)1N as N increases, whereas a CLV is localized when IPR(i)(N)const. as N increases.32,33

In Figs. 5(e) and 5(f), the numerically obtained IPRs of the CLVs corresponding to Λperturb(0) and Λρ(incoh) of the stationary Poisson chimera are plotted vs the system size. The proportionality of IPR(N)1N 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 (Λperturb(0) and Λρ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 (N1)-fold degenerate Λtrans,κ(0) for κ=2,,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 Λperturb(0) strongly depends on the motion of the incoherent oscillators through Z in Eq. (11), which is determined by the initial condition.

FIG. 6.

(a) and (b) Temporal evolution of the magnitude of the Kuramoto order parameter obtained from non-Poisson chimera times series starting from different n-PICs after a time t106 for N=12 and A=0.2. (c) and (d) Lyapunov spectra corresponding to the dynamics of (a) and (b).

FIG. 6.

(a) and (b) Temporal evolution of the magnitude of the Kuramoto order parameter obtained from non-Poisson chimera times series starting from different n-PICs after a time t106 for N=12 and A=0.2. (c) and (d) Lyapunov spectra corresponding to the dynamics of (a) and (b).

Close modal

Concerning the LEs in the incoherent population, (N1)-fold degenerate Λzero(incoh)=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 Figs. 6(c) and 6(d).

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.

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

(15)

where the

i
th oscillator is disconnected to the (i+N2)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

(16)

for i=1,,N and

(17)

for i=1,,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 MPoisson, 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.

FIG. 7.

The Kuramoto order parameters of the phase oscillators governed by the nonlocal intra-group coupling. (a) and (b) Chimera states with the parameter A=0.2 and β=0.08 corresponding to stationary chimeras for the system sizes N=6 and N=60, respectively. (c) and (d) Chimera states with A=0.35 corresponding to the breathing chimera states. Gray line: synchronized group [r(t)=1]. Black line: incoherent order parameter [r(t)<1].

FIG. 7.

The Kuramoto order parameters of the phase oscillators governed by the nonlocal intra-group coupling. (a) and (b) Chimera states with the parameter A=0.2 and β=0.08 corresponding to stationary chimeras for the system sizes N=6 and N=60, respectively. (c) and (d) Chimera states with A=0.35 corresponding to the breathing chimera states. Gray line: synchronized group [r(t)=1]. Black line: incoherent order parameter [r(t)<1].

Close modal

1. Synchronized population: Λtrans(0) and Λ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 N1 transverse Lyapunov exponents consisting of two different values. This splitting of the values of Λ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 N1 transverse Lyapunov exponents to the sync-manifold are

(18)

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κ(0)=[vκ1(trans),,vκN(trans),0,,0]Tϕch(t)(T2N), while i=1Nvκi(trans)=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.

FIG. 8.

(a) Schematic drawing of the two-population oscillator network with nonlocal coupling for N=6. The same color in the incoherent group indicates that the oscillators marked by the same color are characterized by the same evolution dynamics. (b) and (c) The Lyapunov spectra for A=0.2 with N=6 and N=60, respectively. (d) and (e) The Lyapunov spectra for A=0.35 with N=6 and N=60, respectively.

FIG. 8.

(a) Schematic drawing of the two-population oscillator network with nonlocal coupling for N=6. The same color in the incoherent group indicates that the oscillators marked by the same color are characterized by the same evolution dynamics. (b) and (c) The Lyapunov spectra for A=0.2 with N=6 and N=60, respectively. (d) and (e) The Lyapunov spectra for A=0.35 with N=6 and N=60, respectively.

Close modal

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 Λperturb(0)=ν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=(δϕ0,,δϕ0) where |δϕ0|1] and the time shift [vtsϕ˙ch=f(ϕ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).

FIG. 9.

Schematic representation of chimera trajectories in the invariant manifold Mincoh (sphere) and the Poisson submanifold MPoissonMincoh (red curve). Each line schematically represents a trajectory of a chimera state from a given initial condition. The arrow indicates the time flow in the incoherent phase space. (a) For the Kuramoto–Sakaguchi phase oscillators on the global intra-group coupling; the incoherent trajectories of the chimeras dwell in the neutrally stable manifold Mincoh. The Poisson chimera trajectories reside in the invariant and also a neutrally stable Poisson submanifold only if the trajectory starts from PIC; the non-Poisson chimera from n-PIC dwells in the manifold outside the Poisson submanifold. Thus, the non-Poisson chimeras exhibit various incoherent motions according to the given n-PIC. (b) Attracting Poisson chimeras for the nonlocal intra-group topology or Stuart–Landau oscillators. The trajectories starting even from n-PIC eventually settle down on or close to the Poisson trajectory.

FIG. 9.

Schematic representation of chimera trajectories in the invariant manifold Mincoh (sphere) and the Poisson submanifold MPoissonMincoh (red curve). Each line schematically represents a trajectory of a chimera state from a given initial condition. The arrow indicates the time flow in the incoherent phase space. (a) For the Kuramoto–Sakaguchi phase oscillators on the global intra-group coupling; the incoherent trajectories of the chimeras dwell in the neutrally stable manifold Mincoh. The Poisson chimera trajectories reside in the invariant and also a neutrally stable Poisson submanifold only if the trajectory starts from PIC; the non-Poisson chimera from n-PIC dwells in the manifold outside the Poisson submanifold. Thus, the non-Poisson chimeras exhibit various incoherent motions according to the given n-PIC. (b) Attracting Poisson chimeras for the nonlocal intra-group topology or Stuart–Landau oscillators. The trajectories starting even from n-PIC eventually settle down on or close to the Poisson trajectory.

Close modal

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

(19)

with ω~m(t)=μNsin(sm+N/2smα)O(N1). 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,47,48

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 ϕi(t)[π,π) and an amplitude ri(t)(0,) variable. The governing equations are

(20)

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

(21)

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) and (2) of the phase-only oscillators whose amplitude ri1 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 ri(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 ri(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; rincoh(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.

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

(22)

for k=1,,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,,aN,aN+1,,a2N,b1,,bN,bN+1,,b2N)Txch(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 [ε=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.

FIG. 10.

Lyapunov exponents vs index of the strongly coupled SL oscillators with global intra-group coupling for A=0.2, N=20 and (a) and (b) ε=0.1 (stationary) and (c) and (d) ε=0.1 (breathing). (a) and (c) Lyapunov exponents of phase DOFs. The insets show a magnification of the Lyapunov exponents corresponding to the incoherent phase DOFs. (b) and (d) Lyapunov exponents corresponding to amplitude DOFs. (e)–(g) IPR vs system size N for ε=0.15 and Lyapunov modes corresponding to the exponents in PART 1 (e), PART 2 (f), and PART 3 (g). The black dashed guidelines indicate 1/N.

FIG. 10.

Lyapunov exponents vs index of the strongly coupled SL oscillators with global intra-group coupling for A=0.2, N=20 and (a) and (b) ε=0.1 (stationary) and (c) and (d) ε=0.1 (breathing). (a) and (c) Lyapunov exponents of phase DOFs. The insets show a magnification of the Lyapunov exponents corresponding to the incoherent phase DOFs. (b) and (d) Lyapunov exponents corresponding to amplitude DOFs. (e)–(g) IPR vs system size N for ε=0.15 and Lyapunov modes corresponding to the exponents in PART 1 (e), PART 2 (f), and PART 3 (g). The black dashed guidelines indicate 1/N.

Close modal

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

(23)

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:

where i=1Naκi(amp,0)=i=1Nbκi(amp,0)=0 for κ=2,,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

(24)

which is a slightly greater Lyapunov exponent than the transverse ones Λtrans(amp,0)Λ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,,a,a1(inc),,aN(inc),b,,b,b1(inc),,bN(inc))Txch(t)(R4N), where a,bR are constant and j=1Naj(inc)0 and j=1Nbj(inc)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.

FIG. 11.

Full Lyapunov spectra of N=20 Stuart–Landau oscillators for weak coupling ε=0.01 and (a) and (b) stationary Poisson chimeras with A=0.2 and (c) and (d) breathing Poisson chimeras with A=0.35. Left column: phase DOFs. Right column: amplitude DOFs.

FIG. 11.

Full Lyapunov spectra of N=20 Stuart–Landau oscillators for weak coupling ε=0.01 and (a) and (b) stationary Poisson chimeras with A=0.2 and (c) and (d) breathing Poisson chimeras with A=0.35. Left column: phase DOFs. Right column: amplitude DOFs.

Close modal

2. Phase degrees of freedom

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

(25)

for κ=2,,N, where Z~=m=1NRmR0cos(sms0α) 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 Λperturb(0)=νNZ~<0 and is expected to be found in the synchronized phase DOFs. The numerical CLV analysis also confirms that Λtrans,κ(0) and Λ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 Λ(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 Ω~m(t) in the phase governing equations, Eq. (C23).25 For the strongly coupled systems with ε=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 (ε=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–43,47,48 For the weak coupling case ε=0.01 in Figs. 11(a) and 11(c), the amplitude fluctuations are not that strong (RmR0=1 for m=1,,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 (ε=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 ε=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)1N 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 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, Λ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 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.

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 N1 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.

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”).

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

The authors have no conflicts to disclose.

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˙(t)=f(x(t)) with an initial condition x(0)=x0, where x(t)Rn 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 δx(t) resides; i.e., δx(t)Txref(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 δx˙(t)=J(t;xref(t))δx(t) where the Jacobian matrix is defined by (J)ij=x˙ixj|xref(t). From this, we consider the fundamental matrix solution such that O˙(t)=J(t;xref(t))O(t) with O(0)=In; this solution defines the tangent linear propagator, M(t0,t)=O(t)O1(t0), of the perturbation vector from a given point in time point to the future time so that δx(t)=M(t0,t)δx(t0).30 

Oseledets’ theorem29,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

(A1)

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 vectorsd±(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: (Γ(i)(t))+=j=in(U(j)(t))+ and (Γ(i)(t))=j=1i(U(j)(t)), where (U(j)(t))± are the eigenspaces of the forward/backward Oseledets matrices spanned by {d±(j)(t)}j=1n.29 Therefore, we have the decomposition of the tangent space such that Txref(t)(Rn)=j=1nΩ(j)(t), where Ω(i)(t)=(Γ(i)(t))+(Γ(i)(t)) is called Oseledets’ splitting. This Oseledets’ splitting is covariant under the given dynamics in the sense that Ω(i)(t)=M(t0,t)Ω(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

(A2)

for U(t0)(Γ(i)(t0))+(Γ(i+1)(t0))+ where the nested subspaces are Rn=(Γ(1)(t))+(Γ(2)(t))+(Γ(n)(t))+. Hence, the Lyapunov exponents characterize the exponential asymptotic growth rate ||M(t0,t)v(i)(t0)||||v(i)(t0)||exp(Λit), and the covariant Lyapunov vectors indicate the stable/unstable directions of the perturbation vectors in the state space.34 

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 σ of the set of nodes that preserve the adjacency relation among the nodes in the way that Aij=Aσ(i)σ(j).49 Consider the group action under a subgroup GAut(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)= 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 network35–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 equation36,37,51 and the other one the Kuramoto-type equation, which describes diffusively coupled oscillators,35,52,60

(B1)

for i=1,,N where xi(t)Rn 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) {φ(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)|iCm,1mM} under the quotient adjacency matrix A~mm=jCmAij for an arbitrary node iCm, which is nothing but the number of links from an arbitrary node in Cm to all the nodes in Cm. Hence, the quotient dynamics of the CS pattern is given by

(B2)

for m=1,,M.

The set of N-dimensional orthonormal vectors {uκ(m)}κ=1|Cm| for m=1,,M called the cluster-based coordinates is defined by the following rules:35 (i) uκi(m)=0 if iCm, (ii) for κ=1, all the nonzero elements of u1(m) should be 1/|Cm| that defines the cluster sync-manifold, and (iii) the other vectors {uκ(m)}κ=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,

(B3)

for i=1,,N, where δxi(t)=xi(t)sm(t) for iCm 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 ηκ(m)=iCmuκi(m)δxi for m=1,,M and κ=2,,|Cm|, where ηκ(m) for κ2 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|1 transversal variational equations of the cluster Cm both for the Pecora-type and the Kuramoto-type are given by35 

(B4)

for κ=2,,|Cm|, where λκ(m) is the eigenvalue of the adjacency matrix corresponding to the cluster Cm 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.

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,,N. This gives us the corresponding cluster-based coordinates U=[u1(0),u1(1),,u1(N),u2(0),,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 jC0 and u1j(0)=0 if jC0. For the transverse directions, we obtain jC0uκj(0)=0 and uκj(0)=0 if jC0 for κ=2,,N. Finally, for the incoherent trivial clusters, we have u1j(m)=1 if jCm and u1j(m)=0 otherwise, with m=1,,N. Note that all the cluster-based coordinate vectors should be mutually orthonormalized. An example of a possible candidate of the cluster-based coordinates is61 

(C1)

where the first column

u1(0)=[1N,,1N,0,,0]
indicates the sync-manifold direction, D=diag(1,,1)RN×N indicating the incoherent trivial clusters, each O is a zero-matrix, and PRN×N1 representing the directions transverse to C0, and it 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 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

(C2)

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

(C3)

where

A(c)RN×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: A~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 s0(t)=ϕi(t) (synchronized, C0) and sm(t)=ϕi+N(t) (incoherent, Cm) for m=i=1,,N,

(C4)

for the synchronized cluster (C0) where the quotient adjacency matrix A~00=N1 and A~0m=1 for m=1,,N, and H(0)=sinα. The quotient governing equations of the N trivial clusters (C1,,CN) for the incoherent population read

(C5)

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

(C6)

for each iC0 where the deviation around the CS pattern is δϕi(t)=ϕi(t)sm(t) for iCm 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 ηκ(0)(t)=iC0uκi(0)δϕi(t) with U=[u1(0),u1(1),,u1(N),u2(0),,uN(0)]. Then, the variational equations transverse to the sync-cluster C0 read

(C7)

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 C0 can be represented by the matrix diag(λ2(0),λ3(0),,λN(0))R(N1)×(N1) and the off-diagonal blocks are zero. This, in turn, means that the last term in Eq. (C7) should be zero and iC0kC0uκi(0)Bik(c)uκk(0)=λκ(0)δκκ for κ,κ=2,,N where λκ(0) are the eigenvalues of the adjacency matrix since uκ(0) 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 by

(C8)

for κ=2,,N. Notice that for the global intra- and inter- population network, the eigenvalues in Eq. (C8)λκ(0)=1 for all κ=2,,N. Consider as an example the system with N=4. Its block-diagonalized adjacency matrix reads61 

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

C0
⁠, 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 (N1)-fold degenerate transverse LEs in Eq. (9).

To estimate the Lyapunov exponent along the sync-manifold for the synchronized population Λ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)s0(t)+δs0(t) where |δs0(t)|1 is applied to Eq. (C4),

(C9)

Then, we obtain Eq. (11) provided that Z=νNm=1Ncos(sms0α) 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

(C10)

for i=1,,2N, where F(ϕ), Kij, and H(x) are the same as defined in Eq. (C2). The matrix B(n)R2N×2N, which defines the global inter- and nonlocal intra-population network, is given by

where

ARN×N
is defined in Eq. (15) and JNRN×N is the unit matrix whose elements are all 1. In this approach, the quotient adjacency matrix is given by

(C11)

wherein the terms

A~00{(n)}=N2
and A~ij{(n})=Amm for m,m=1,,N ensuring that the intra-population topology is not global but nonlocal. From A~(n), we obtain the quotient dynamics according to the CS pattern describing our chimeras with the variables s0(t)=ϕi(t) (sync.) and sm(t)=ϕi+N(t) (incoh.) for i=m=1,,N,

(C12)

for the synchronized population and

(C13)

where ω~m(t)=μNsin(sm+N/2smα) with A~0m(n)=1 and A~m0(n)=N for m=1,,N for the incoherent trivial clusters.

As seen in Sec. III B, there are N1 transverse Lyapunov exponents consisting of two different values. This splitting of the values of Λ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κ(m), to obtain the transverse variational equations. For the global topology discussed above, these eigenvalues λκ(0)=1 are the same for κ=2,,N. In contrast, the nonlocal adjacency matrix has two different eigenvalues: λκ(0)=0 for κ=2,,N/2+1 and λκ(0)=2 for κ=N/2+2,,N.61 This leads to two different variational equations with the same method in Eqs. (C7) and (C8),

(C14)

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

(C15)

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 s0s0+δs0 is imposed on Eq. (C12) where |δs0|1. This perturbation gives Λsync(0)=ν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

(C16)

for i=1,,2N, where F(amp)(r)=ε1(1r2)r+μNrcosα 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)=μNcos(ϕjϕiα) if i,j belong to the same population and Kij(amp)=νNcos(ϕ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. C1, 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)=ϕi(t) (sync.) and sm(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

(C17)

Considering a small deviation around the CS dynamics, i.e., δri(t)=ri(t)Rm(t) for iCm for m=0,1,,N, we obtain the coupled variational equations as

for each iC0 and Cmm=cos(smsmα) for m,m=0,,N. Then, viewing these in the cluster-based coordinates with ξκ(0)(t)=iC0uκi(0)δri(t) for κ=2,,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 iC0kC0uκi(0)Bik(c)uκk(0)=λκ(0)δκκ for κ=2,,N. Hence, the N1 variational equations transversal to the sync-manifold are given by

(C18)

Here, λκ(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 (N1)-fold degenerate transverse Lyapunov exponents in the amplitude DOFs as Λtrans,κ(amp,0)ε1(13R02)<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 sync-manifold R0(t)R0(t)+δR0(t) with |δR0|1 in Eq. (C17) and obtain

(C19)

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

(C20)

which shows that Λtrans(amp,0)Λ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),

(C21)

for m=1,,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

(C22)

for i=1,,2N, where the uncoupled dynamics is governed by F(ph)(ϕi)=σri2μNsinα and the coupling function is defined as H(x)=sin(xα). The coupling weights are defined by Kij(ph)=μNrjri if i,j belong to the same population and Kij(ph)=ν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

(C23)

where Ω~m(t)=σRm2(t) for m=1,,N. The quotient dynamics reveal that the phase DOFs of the SL oscillator ensemble in the synchronized population also have (N1)-fold degenerate transverse Lyapunov exponents with

(C24)

for κ=2,,N, where Z~=m=1NRmR0cos(sms0α) should be considered an external forcing field, which follows from the transversal variational equations [Eq. (B4)],

(C25)

Here, ηκ(0)(t)=iC0uκi(0)δϕi(t) for κ=2,,N; the deviation along the CS dynamics is δϕi(t)=ϕi(t)s0(t) for iC0 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 Λperturb(0)=ν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))R2 where ak and bk are defined in Eq. (22). Then, the SL oscillators evolve according to

(C26)

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

(C27)

and the coupling function is written as

(C28)

for i=1,,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,,N, then the variational equations transversal to the synchronized cluster C0 in the cluster-based coordinates are given by

(C29)

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

(C30)

and

(C31)

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

(C32)

for all the directions transverse to the sync-manifold. Notice that the eigenvalues of the adjacency matrix λκ(0)=1 for all κ. If we consider ϕ0 as an external forcing function, then the eigenvalues of the matrix Jtrans(0) are

(C33)

which gives Λ12ε1 corresponding to the amplitude DOF branch and Λ20 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.

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 (ε=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 λκ(0)=0 for κ=2,,N/2+1 and λκ(0)=2 for κ=N/2+2,,N. Also, the negative LE corresponding to the sync-manifold perturbation is given as Λperturb(0)=ν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 N1 transverse Lyapunov exponents also show the two different values of the degenerate exponents approximated as

(D1)

since for the nonlocal network λκ(0)=0 for κ=2,,N/2+1 and λκ(0)=2 for κ=N/2+2,,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, Λperturb(amp,0)ε1(13R02)+μN(N1)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.

FIG. 12.

Full Lyapunov spectra of the Stuart–Landau oscillators with nonlocal intra-population topology for (a) the phase degrees of freedom and (b) the amplitude degrees of freedom. The parameter set used here is N=6, A=0.2, and ε=0.01. (c) and (d) The breathing chimera states with A=0.35. Note that the Lyapunov exponents in (a) follow the same behavior as the stationary chimera state of the phase-only system (compare Fig. 8).

FIG. 12.

Full Lyapunov spectra of the Stuart–Landau oscillators with nonlocal intra-population topology for (a) the phase degrees of freedom and (b) the amplitude degrees of freedom. The parameter set used here is N=6, A=0.2, and ε=0.01. (c) and (d) The breathing chimera states with A=0.35. Note that the Lyapunov exponents in (a) follow the same behavior as the stationary chimera state of the phase-only system (compare Fig. 8).

Close modal

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.

1.
A.
Pikovsky
,
M.
Rosenblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
Cambridge
,
2001
).
2.
S. H.
Strogatz
,
Sync
(
Hyperion
,
New York
,
2003
).
3.
Y.
Kuramoto
and
D.
Battogtokh
, “
Coexistence of coherence and incoherence in nonlocally coupled phase oscillators
,”
Nonlinear Phenom. Complex Syst.
5
,
380
(
2002
).
4.
E. A.
Martens
,
A. F. S.
Thutupalli
, and
O.
Hallatschek
, “
Chimera states in mechanical oscillator networks
,”
Proc. Natl. Acad. Sci. U.S.A.
110
,
10563
(
2013
).
5.
L.
Schmidt
,
K.
Schönleber
, and
K.
Krischer
, “
Coexistence of synchrony and incoherence in oscillatory media under nonlinear global coupling
,”
Chaos
24
,
013102
(
2014
).
6.
J. D.
Hart
,
K.
Bansal
,
T. E.
Murphy
, and
R.
Roy
, “
Experimental observation of chimera and cluster states in a minimal globally coupled network
,”
Chaos
26
,
094801
(
2016
).
7.
M.
Wickramasinghe
and
I. Z.
Kiss
, “
Spatially organized dynamical states in chemical oscillator networks: Synchronization, dynamical differentiation, and chimera patterns
,”
PLoS One
8
,
e80586
(
2013
).
8.
J. F.
Totz
,
J.
Rode
,
M. R.
Tinsley
,
K.
Showalter
, and
H.
Engel
, “
Spiral wave chimera states in large populations of coupled chemical oscillators
,”
Nat. Phys.
14
,
282
(
2018
).
9.
A. M.
Hagerstrom
,
T. E.
Murphy
,
R.
Roy
,
P.
Hövel
,
I.
Omelchenko
, and
E.
Schöll
, “
Experimental observation of chimeras in coupled-map lattices
,”
Nat. Phys.
8
,
658
(
2012
).
10.
M. R.
Tinsley
,
S.
Nkomo
, and
K.
Showalter
, “
Chimera and phase-cluster states in populations of coupled chemical oscillators
,”
Nat. Phys.
8
,
662
(
2012
).
11.
C. R.
Laing
and
C. C.
Chow
, “
Stationary bumps in networks of spiking neurons
,”
Neural Comput.
13
,
1473
(
2001
).
12.
N.
Rattenborg
,
C.
Amlaner
, and
S.
Lima
, “
Behavioral, neurophysiological and evolutionary perspectives on unihemispheric sleep
,”
Neurosci. Biobehav. Rev.
24
,
817
(
2000
).
13.
O. E.
Omel’chenko
, “
The mathematics behind chimera states
,”
Nonlinearity
31
,
R121
(
2018
).
14.
O. E.
Omel’chenko
, “
Coherence–incoherence patterns in a ring of non-locally coupled phase oscillators
,”
Nonlinearity
26
,
2469
(
2013
).
15.
M. J.
Panaggio
and
D. M.
Abrams
, “
Chimera states: Coexistence of coherence and incoherence in networks of coupled oscillators
,”
Nonlinearity
28
,
R67
(
2015
).
16.
D. M.
Abrams
,
R.
Mirollo
,
S. H.
Strogatz
, and
D. A.
Wiley
, “
Solvable model for chimera states of coupled oscillators
,”
Phys. Rev. Lett.
101
,
084103
(
2008
).
17.
M. J.
Panaggio
,
D. M.
Abrams
,
P.
Ashwin
, and
C. R.
Laing
, “
Chimera states in networks of phase oscillators: The case of two small populations
,”
Phys. Rev. E
93
,
012218
(
2016
).
18.
D.
Pazó
and
E.
Montbrió
, “
Low-dimensional dynamics of populations of pulse-coupled oscillators
,”
Phys. Rev. X
4
,
011099
(
2014
).
19.
E.
Montbrió
,
J.
Kurths
, and
B.
Blasius
, “
Synchronization of two interacting populations of oscillators
,”
Phys. Rev. E
70
,
056125
(
2004
).
20.
E. A.
Martens
,
M. J.
Panaggio
, and
D. M.
Abrams
, “
Basin of attraction for chimera states
,”
New J. Phys.
18
,
022002
(
2016
).
21.
C. R.
Laing
, “
Chimera states in heterogeneous networks
,”
Chaos
19
,
013113
(
2009
).
22.
C. R.
Laing
,
K.
Rajendran
, and
I. G.
Kevrekidis
, “
Chimeras in random non-complete networks of phase oscillators
,”
Chaos
22
,
013132
(
2012
).
23.
E. A.
Martens
,
C.
Bick
, and
M. J.
Panaggio
, “
Chimera states in two populations with heterogeneous phase-lag
,”
Chaos
26
,
094819
(
2016
).
24.
A.
Buscarino
,
M.
Frasca
,
L. V.
Gambuzza
, and
P.
Hövel
, “
Chimera states in time-varying complex networks
,”
Phys. Rev. E
91
,
022817
(
2015
).
25.
C. R.
Laing
, “
Chimeras in networks of planar oscillators
,”
Phys. Rev. E
81
,
066221
(
2010
).
26.
C. R.
Laing
, “
Dynamics and stability of chimera states in two coupled populations of oscillators
,”
Phys. Rev. E
100
,
042211
(
2019
).
27.
A.
Pikovsky
and
M.
Rosenblum
, “
Partially integrable dynamics of hierarchical populations of coupled oscillators
,”
Phys. Rev. Lett.
101
,
264103
(
2008
).
28.
A.
Pikovsky
and
M.
Rosenblum
, “
Dynamics of heterogeneous oscillator ensembles in terms of collective variables
,”
Physica D
240
,
872
(
2011
).
29.
F.
Ginelli
,
H.
Chaté
,
R.
Livi
, and
A.
Politi
, “
Covariant Lyapunov vectors
,”
J. Phys. A: Math. Theor.
46
,
254005
(
2013
).
30.
P. V.
Kuptsov
and
U.
Parlitz
, “
Theory and computation of covariant Lyapunov vectors
,”
J. Nonlinear Sci.
22
,
727
(
2012
).
31.
J. P.
Eckmann
, “
Ergodic theory of chaos and strange attractors
,”
Rev. Mod. Phys.
57
,
3
(
1985
).
32.
K. A.
Takeuchi
and
H.
Chaté
, “
Collective Lyapunov modes
,”
J. Phys. A: Math. Theor.
46
,
254007
(
2013
).
33.
K.
Höhlein
,
F. P.
Kemeth
, and
K.
Krischer
, “
Lyapunov spectra and collective modes of chimera states in globally coupled Stuart-Landau oscillators
,”
Phys. Rev. E
100
,
022217
(
2019
).
34.
A.
Pikovsky
and
A.
Politi
,
Lyapunov Exponents: A Tool to Explore Complex Dynamics
(
Cambridge University Press
,
Cambridge
,
2016
).
35.
Y. S.
Cho
,
T.
Nishikawa
, and
A. E.
Motter
, “
Stable chimeras and independently synchronizable clusters
,”
Phys. Rev. Lett.
119
,
084101
(
2017
).
36.
L. M.
Pecora
,
F.
Sorrentino
,
A. M.
Hagerstrom
,
T. E.
Murphy
, and
R.
Roy
, “
Cluster synchronization and isolated desynchronization in complex networks with symmetries
,”
Nat. Commun.
5
,
4079
(
2014
).
37.
F.
Sorrentino
,
L.
Pecora
,
A. M.
Hagerstrom
,
T. E.
Murphy
, and
R.
Roy
, “
Complete characterization of the stability of cluster synchronization in complex dynamical networks
,”
Sci. Adv.
2
,
e1501737
(
2016
).
38.
M.
Wolfrum
,
O. E.
Omel’chenko
,
S.
Yanchuk
, and
Y. L.
Maistrenko
, “
Spectral properties of chimera states
,”
Chaos
21
,
013112
(
2011
).
39.
M.
Wolfrum
and
O. E.
Omel’chenko
, “
Chimera states are chaotic transients
,”
Phys. Rev. E
84
,
015201
(
2011
).
40.
T.
Bountis
,
V. G.
Kanas
,
J.
Hizanidis
, and
A.
Bezerianos
, “
Chimera states in a two-population network of coupled pendulum-like elements
,”
Eur. Phys. J. Spec. Top.
223
,
721
(
2014
).
41.
C. R.
Laing
, “
The dynamics of chimera states in heterogeneous Kuramoto networks
,”
Physica D
238
,
1569
(
2009
).
42.
I. V.
Tyulkina
,
D. S.
Goldobin
,
L. S.
Klimenko
, and
A.
Pikovsky
, “
Dynamics of noisy oscillator populations beyond the Ott-Antonsen ansatz
,”
Phys. Rev. Lett.
120
,
264101
(
2018
).
43.
E.
Ott
and
T. M.
Antonsen
, “
Long time evolution of phase oscillator systems
,”
Chaos
19
,
023117
(
2009
).
44.
S. A.
Marvel
,
R. E.
Mirollo
, and
S. H.
Strogatz
, “
Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action
,”
Chaos
19
,
043104
(
2009
).
45.
S.
Watanabe
and
S. H.
Strogatz
, “
Constants of motion for superconducting Josephson arrays
,”
Physica D
74
,
197
(
1994
).
46.
E.
Ott
and
T. M.
Antonsen
, “
Low dimensional behavior of large systems of globally coupled oscillators
,”
Chaos
18
,
037113
(
2008
).
47.
B.
Pietras
and
A.
Daffertshofer
, “
Ott-Antonsen attractiveness for parameter-dependent oscillatory systems
,”
Chaos
26
,
103101
(
2016
).
48.
V.
Vlasov
,
M.
Rosenblum
, and
A.
Pikovsky
, “
Dynamics of weakly inhomogenous oscillator populations: Perturbation theory on top of Watanabe-Strogatz integrability
,”
J. Phys. A: Math. Theor.
49
,
31LT02
(
2016
).
49.
S.
Kudose
, “Equitable partitions and orbit partitions,” in Kudos 2009, Equitable, PA 2009 (2009).
50.
B. D.
MacArthur
,
R. J.
Sánchez-García
, and
J. W.
Anderson
, “
Symmetry in complex networks
,”
Discret. Appl. Math.
156
,
3525
(
2008
).
51.
Y. S.
Cho
, “
Concurrent formation of nearly synchronous clusters in each intertwined cluster set with parameter mismatches
,”
Phys. Rev. E
99
,
052215
(
2019
).
52.
S.
Lee
,
Y. S.
Cho
, and
H.
Hong
, “
Twisted states in low-dimensional hypercubic lattices
,”
Phys. Rev. E
98
,
062221
(
2018
).
53.
S.
Nichols
and
K.
Wiesenfeld
, “
Ubiquitous neutral stability of splay-phase states
,”
Phys. Rev. A
45
,
8430
(
1992
).
54.
D. G.
Aronson
,
M.
Golubitsky
, and
J.
Mallet-Paret
, “
Ponies on a merry-go-round in large arrays of Josephon junctions
,”
Nonlinearity
4
,
903
(
1991
).
55.
A.
Pikovsky
and
M.
Rosenblum
, “
Self-organized partially synchronous dynamics in populations of nonlinearly coupled oscillators
,”
Physica D
238
,
27
(
2009
).
56.
V.
Oseledets
, “
A multiplicative ergodic theorem. Characteristic Liapunov, exponents of dynamical systems
,”
Trans. Mosc. Math. Soc.
19
,
197
(
1968
).
57.
V.
Nicosia
,
M.
Valencia
,
M.
Chavez
, and
A.
Díaz-Guilera
, “
Remote synchronization reveals network symmetries and functional modules
,”
Phys. Rev. Lett.
110
,
174102
(
2013
).
58.
M. T.
Schaub
,
N.
O’Clery
,
Y. N.
Billeh
,
J.
Delvenne
,
R.
Lambiotte
, and
M.
Barahona
, “
Graph partitions and cluster synchronization in networks of oscillators
,”
Chaos
26
,
094821
(
2016
).
59.
F. P.
Kemeth
,
S. W.
Haugland
, and
K.
Krischer
, “
Symmetries of chimera states
,”
Phys. Rev. Lett.
120
,
214101
(
2018
).
60.
S.
Lee
and
Y. S.
Cho
, “
Stable chimeras of non-locally coupled Kuramoto–Sakaguchi oscillators in a finite array
,”
J. Korean Phys. Soc.
78
,
476
(
2020
).
61.
See https://github.com/tnishi0/grouping-clusters/ for information about the cluster synchronization pattern and independently synchronizable clusters and cluster sets.