Despite their simplicity, networks of coupled phase oscillators can give rise to intriguing collective dynamical phenomena. However, the symmetries of globally and identically coupled identical units do not allow solutions where distinct oscillators are frequency-unlocked—a necessary condition for the emergence of chimeras. Thus, forced symmetry breaking is necessary to observe chimera-type solutions. Here, we consider the bifurcations that arise when full permutational symmetry is broken for the network to consist of coupled populations. We consider the smallest possible network composed of four phase oscillators and elucidate the phase space structure, (partial) integrability for some parameter values, and how the bifurcations away from full symmetry lead to frequency-unlocked weak chimera solutions. Since such solutions wind around a torus they must arise in a global bifurcation scenario. Moreover, periodic weak chimeras undergo a period-doubling cascade leading to chaos. The resulting chaotic dynamics with distinct frequencies do not rely on amplitude variation and arise in the smallest networks that support chaos.

Networks of coupled oscillators occur in a wide range of systems in biology, medicine, and technology. The proper functioning of such systems often relies on the emergence (or absence) of collective modes such as synchronization, where oscillators lock their frequencies and/or phases^{1,2}—a phenomenon that can be studied in the prominent Kuramoto model of phase oscillators.^{3,4} Chimeras are symmetry-breaking solutions where part of the population synchronizes while the other oscillates incoherently,^{5} even if oscillators are identical. While this phenomenon has received much attention in recent years,^{6,7} the bifurcations of chimera solutions in finite networks have remained elusive. Previous studies^{8,9} concerning variants of the Kuramoto model, composed of two populations of identical oscillators with heterogeneous intra- and inter-population coupling/phase lags, found that bifurcations leading to chimera states emerge for phase lags near multiples of $\pi 2$—these are parameters that can be related to resonance points in physical/mechanical oscillator models.^{10,11} Here, we analyze such a variant of the Kuramoto model with two populations with two phase oscillators near these parameter points. We present a detailed analysis of how phase space is organized and elucidate the emergence and bifurcations of chaotic weak chimeras.^{9,12,13}

## I. INTRODUCTION

Coupled phase oscillator networks have been instrumental in understanding collective dynamic phenomena in real-world oscillatory systems.^{1,2} The state of each dynamical node is given by a single phase variable on the circle $T:=R/2\pi Z$. To understand the dynamics, in many cases, it makes sense to assume an idealized system where all nodes are identically connected to all other nodes (such as in the Kuramoto model^{3,4}) and/or all nodes are identical (for example, in the classical master stability approach^{14}). Both of these assumptions together lead to full permutational symmetry, which restricts the possible dynamics: While these systems can still exhibit intriguing collective dynamics,^{15,16} the symmetry restricts the oscillators’ frequencies: For networks of phase oscillators, all nodes have to be frequency synchronized since oscillators’ phases cannot pass each other.^{8}

Dynamics where different oscillators evolve at different (asymptotic) frequencies have attracted significant attention in recent years, a phenomenon commonly known as a *chimera* (see Refs. 7, 17, and 18 for reviews). As such solutions/invariant sets are impossible in fully symmetric phase oscillator networks, breaking system symmetries—also known as forced symmetry breaking^{19}—is necessary for frequency-unlocked solutions to arise. A less symmetric system will have fewer invariant subspaces that restrict the dynamics and create the possibility for frequency-unlocked chimera-type solutions. Note that we focus here on first-order phase oscillators rather than systems with additional degrees of freedom in which chimera phenomena have been observed. These include oscillators with amplitude variation^{20–22} or second-order phase equations;^{23} due to the larger phase space, chimera-type dynamics^{24,25} can occur in these extended systems also for fully symmetric networks.

A classical way to break full permutational symmetry is to consider coupled populations of phase oscillators: While oscillators are still all-to-all coupled, distinct populations form when coupling within populations can differ from the coupling between populations. Such a coupling structure for $L=MN$ phase oscillators of Kuramoto-type, where the coupling is mediated by a single harmonic, arises by introducing distinct coupling strengths $Ks,Kn>0$ and phase lags $\alpha s,\alpha n\u2208T$ within and between populations. Specifically, the phase $\theta \sigma ,k$ of oscillator $k\u2208{1,\u2026,N}$ of population $\sigma \u2208{1,\u2026,M}$ evolves according to

In other words, coupling between oscillators is mediated by the coupling functions $gs(\varphi )=Kssin\u2061(\varphi \u2212\alpha s)$, $gn(\varphi )=Knsin\u2061(\varphi \u2212\alpha n)$ that determine the *self-coupling* within populations and *coupling to other populations* (its “neighbors”), respectively. While system (1) is not symmetric with respect to any permutation of oscillator indices, it still has a wreath product symmetry that permutes oscillators within populations or permutes the populations^{26} (see also Ref. 27). In particular, the populations in (1) are interchangeable in contrast to networks with distinct populations.^{28}

Splitting the oscillators into just $M=2$ populations can lead to frequency-unlocked chimera dynamics. Indeed, the classical work by Abrams *et al.*^{29} demonstrated the emergence of chimera solutions for two coupled populations with identical phase lags $\alpha s=\alpha n$ but disparate coupling strength $Kn\u2260Ks$ in the limit of $N\u2192\u221e$ oscillators. More generic interactions with nonidentical phase lags $\alpha s\u2260\alpha n$ in the limit of large networks may lead to chaotic collective dynamics.^{9,13} Chimera solutions not only arise in the mean-field limit but also for finite networks^{6} and general sinusoidal coupling is sufficient to obtain chaotic chimera dynamics.^{13} While it was already indicated in Ref. 13 that chaotic dynamics are possible even in the smallest networks of $M=2$ populations of $N=2$ oscillators, the bifurcations that lead to chaotic chimera dynamics remained elusive.

The main focus of this paper is to analyze small networks that consist of $M=2$ populations of $N=2$ oscillators with generic sinusoidal coupling of Kuramoto-type (1) [see Fig. 1(a)]. This is the smallest network of this type in which chaotic dynamics can occur. Even as the full permutational symmetry is broken, phase space is organized by special structures for specific parameter: We first analyze degenerate cases where the dynamical system is either a gradient system or Hamiltonian-like with conserved quantities. These special cases allow us to understand how frequency-unlocked dynamics arise in these small networks. Perturbing away from the special parameter values yields global bifurcation scenarios that give rise to frequency-unlocked limit cycle solutions—they correspond to *weak chimeras*. As such, these solutions have nonzero winding number, wrap around the torus, and have nontrivial homology. How (some of) these solutions bifurcate further is shown in the numerical bifurcation diagram in Fig. 1(b) in terms of the phase difference $\psi 1:=\theta 1,1\u2212\theta 1,2$ (note that further branches co-exist): The dynamics eventually undergo a period-doubling cascade that leads to chaotic weak chimeras. Finally, we highlight that for certain parameter values, the system has simultaneously conservative and dissipative dynamics in different regions of phase space.

The paper is organized as follows. In Sec. II, we collect basic information about the governing equation (1), identify invariant subspaces and equilibria, and analyze their linear stability. In Sec. III, we consider the case of pure sine coupling (*odd coupling*) for which the equations of motion are a gradient system. In Sec. IV, we consider the case of pure cosine coupling (*even coupling*) that leads to the emergence of conserved quantities; the phase space structure elucidates the emergence or frequency-unlocked solutions corresponding to weak chimeras. Subsequently, in Sec. V, we show how these weak chimeras bifurcate. Specifically, we detail the bifurcations in Fig. 1 that lead to the emergence of chaotic dynamics. In Sec. VI, we briefly comment on coexisting conservative–dissipative dynamics for combined pure sine/cosine coupling. We conclude with a discussion in Sec. VII; here, we also briefly consider the dynamics of small networks with $M=3$ populations.

## II. FROM GLOBAL COUPLING TO COUPLED POPULATIONS OF KURAMOTO OSCILLATORS

Before we analyze coupled phase oscillator networks, we briefly recall some notions related to dynamical systems that are equivariant with respect to the action of a group of symmetries. Let $F:X\u2192TX$ be a smooth vector field on $X$ where $TX$ denotes the tangent bundle. Suppose that a group $\Gamma $ acts on $X$. The vector field $F$ is $\Gamma $-equivariant if

for all $\gamma \u2208\Gamma $ where $\gamma ^$ is the induced action on the tangent space. A $\Gamma $-equivariant vector field $F$ defines a $\Gamma $-*equivariant dynamical system*,

on $X$.^{30,31} The group $\Gamma x:={\gamma \u2208\Gamma |\gamma x=x}$ denotes the *stabilizer* or *isotropy subgroup* of $x\u2208X$. While the isotropy subgroup describes the symmetries of a single point, the *symmetries of a set $X\u2282X$* are $\Sigma (X):={\gamma \u2208\Gamma |\gamma X\u2282X}$ where $\gamma X={\gamma x|x\u2208X}$. If $H\u2282\Gamma $ is a subgroup, then the set $Fix\u2061(H):={x\u2208X|\gamma x=x\u2200\gamma \u2208H}$ is invariant under the flow induced by (2).

### A. Networks of *L* oscillators

While the later analysis focuses on small networks, we first collect general properties of the general network (1) consisting of $M$ populations of $N$ phase oscillators. We consider the coupling of Kuramoto type with $gs(\varphi )=Kssin\u2061(\varphi \u2212\alpha s)$ and $gn(\varphi )=Knsin\u2061(\varphi \u2212\alpha n)$ throughout the paper but the symmetry analysis in this section is independent of this specific choice. Since the coupling is through phase differences only, the system is $T$-equivariant: There is a continuous $T$ symmetry where $\gamma \u2208T$ acts on the phase space $TL$ by a common phase shift $\theta \sigma ,k\u21a6\theta \sigma ,k+\gamma $ to all oscillators. For Kuramoto-type coupling, note that by rescaling time, we can set $Ks+Kn=1$ without loss of generality and parameterize the coupling strength through the disparity parameter $A:=Ks\u2212Kn$ (see also Ref. 29).

#### 1. Globally and identically coupled oscillators

If the coupling within and between populations is identical, that is, $A=0$ ($Ks=Kn$) and $\alpha s=\alpha n=:\alpha $, then we have a globally coupled network of $L=MN$ identical oscillators. Omitting the population index $\sigma $, write $\theta =(\theta 1,\u2026,\theta L)$ and the evolution of the phase $\theta k$ of oscillator $k\u2208{1,\u2026,L}$ is

System (3) is $SL$-equivariant, where $SL$ denotes the symmetric group on $L$ symbols, which acts on $TL$ by permuting the oscillator indices. This implies that the sets $Pmn:={\theta |\theta m=\theta n}$ are dynamically invariant as fixed-point sets of the transposition that swaps oscillators $m$ and $n$. Thus, the dynamics can be reduced to the canonical invariant region

that is bounded by $Pnm$ (see Refs. 32 and 33 for details). Indeed, since the coupling function only contains a first harmonic, the canonical invariant region is foliated into dynamically invariant sheets on which the dynamics are effectively two-dimensional.^{34}

The symmetries and type of coupling imply that certain phase configurations are dynamically invariant. The phase configuration

where all oscillators are *phase synchronized* is dynamically invariant as the intersection of all $Pmn$. Phase synchrony is typically quantified by the *Kuramoto order parameter*

with $i:=\u22121$: We have $|Z|=1$ if and only if $\theta \u22080$. Define the *antiphase* or *incoherent phase configurations* as

The set $I$ is dynamically invariant for globally coupled networks and is a union of manifolds (see Ref. 33).

The symmetries also constrain the frequencies of the oscillators. For a solution $\theta (t)$ of (3) with initial condition $\theta (0)=\theta 0$, define the *asymptotic average frequency*^{24,35} as

Since the sets $Pmn$ are dynamically invariant, the phase difference between oscillators is bounded. This implies that

for all $m,n\u2208{1,\u2026,L}$.

While no precise mathematical definition of a chimera exists (but attempts to classify them^{36}), a distinguishing feature of a *weak chimera* (as defined in Ref. 33) is the separation of frequencies: A trajectory $\theta (t)$ with initial condition $\theta (0)=\theta 0$ converges to a weak chimera if there are distinct $m,n,n\u2032$ such that $\Omega m(\theta 0)\u2260\Omega n(\theta 0)=\Omega n\u2032(\theta 0)$. Thus, weak chimeras cannot exist in fully symmetric systems.

#### 2. Breaking full symmetry: *M* populations of *N* oscillators

Breaking the symmetry is necessary to observe weak chimeras in a network of phase oscillators: Here, we now consider the dynamics of $M$ populations of $N$ oscillators each. Equation (1) that determines the evolution is $SN\u2240SM=(SN)M\u22caSM\u2282SL$ equivariant, where $\u22ca$ denotes a semidirect product. Each copy of $SN$ acts on $TL$ by permuting the oscillator indices within a given population $\sigma $ and $SM$ acts by permuting the populations (i.e., the population indices $\sigma $); the symmetries within and between populations do not necessarily commute and, hence, the semi direct product.

The system has fewer symmetries than the fully symmetric system. For example—in the notation of Sec. II A 1—the sets $Pmn$ are only dynamically invariant if oscillators $m,n$ belong to the same population. Note that the oscillators in each population are still fully symmetric. This implies that for (1) there are invariant sets $P\sigma ,mn:={\theta |\theta \sigma ,m=\theta \sigma ,n}$ for $m,n\u2208{1,\u2026,N}$ and any $\sigma \u2208{1,\u2026,M}$.

The (global) Kuramoto order parameter $Z$ is defined as in (5); it quantifies coherence of all oscillators in the network. Naturally, we also define phase synchrony and incoherence for each population: The (local) Kuramoto order parameter for population $\sigma $ is

and we have $Z=1M\u2211\sigma =1MZ\sigma $.

### B. Networks of four oscillators

The main focus of our analysis is the smallest multi-population phase oscillator networks that allow for chaotic dynamics. Consider networks (1) that consist of $M=2$ populations of $N=2$ oscillators; the network is sketched in Fig. 1(a). In this case, the governing equations (1) read

with $\omega \u2032=\omega \u2212Ks4sin\u2061\alpha s$. If the coupling is fully symmetric, i.e., $Ks=Kn$ and $\alpha s=\alpha n$, then (10) is $S4$-equivariant. This means that up to the rotational symmetry, the canonical invariant region is a tetrahedron bounded by the dynamically invariant $Pmn$ (as in Sec. II A 1) (see Refs. 33 and 37 for details). The phase configuration $0$ with full phase synchrony corresponds to each of the four corners of the tetrahedron. The incoherent phase configurations $I$ in $C$ are the points with $S2\xd7S2$ isotropy, that is, a line of points parameterized by $(a,b,a+\pi ,b+\pi )\u2208T4$.

For a generic choice of coupling parameters, the dynamical system (10) is $S2\u2240S2\u2261D4$-equivariant, the symmetries of the square [cf. Fig. 1(a)]. Specifically, $D4$ is generated by the rotational symmetry

[a counterclockwise rotation by an angle of $\pi 2$ in Fig. 1(a)] and the mirror symmetry

[a flip about the main diagonal in Fig. 1(a)]. General properties of $D4$-equivariant phase oscillator networks can be found in Ref. 32. The codimension-1 invariant subspaces are $P1,12={\theta 1,1=\theta 1,2}$ and $P2,12={\theta 2,1=\theta 2,2}$. Since we only have two oscillators per population, we will write $P1:=P1,12$, $P2:=P2,12$ for simplicity.

Three values of the disparity parameter stand out where the network structure is qualitatively different. As already discussed, $A=0$ corresponds to $Kn=Ks$, which is the strength of interactions within and between populations is the same. If in addition also $\alpha n=\alpha s$, the oscillators have full permutational symmetry and we have the Kuramoto–Sakaguchi equations for identical oscillators. If $A=1$, we have $Kn=0$ and the populations are uncoupled; this means that the green links are absent in Fig. 1(a). If $A=\u22121$, we have $Ks=0$ and the oscillators within populations are uncoupled. This configuration is equivalent to a ring of oscillators with nearest-neighbor coupling as illustrated in Fig. 1(a): If the purple links are absent, for a given oscillator, the two oscillators of the other population are the direct neighbors.

#### 1. Reduced dynamics and symmetries

The phase-shift symmetry of (10) means that the system dynamics is effectively three-dimensional. We exploit this by rewriting the system in terms of phase differences $\psi 1:=\theta 1,1\u2212\theta 1,2,\psi 2:=\theta 1,2\u2212\theta 2,1,\psi 3:=\theta 2,1\u2212\theta 2,2$. This yields the reduced system

If the coupling is fully symmetric, i.e., $A=0$ ($Ks=Kn$) and $\alpha s=\alpha n$, then the phase space is organized by invariant planes $Pnm$, the canonical invariant region, and its images under the symmetry. In the reduced system (12), these invariant planes are given by $\psi 1=0$, $\psi 3=0$, $\psi 2=0$, $\psi 1+\psi 2=0$, $\psi 2+\psi 3=0$, $\psi 1+\psi 2+\psi 3=0$, which split $T3$ into six dynamically invariant regions, the canonical invariant regions and its images. The geometry is illustrated in Fig. 3(a) where one invariant region is bounded by coordinate planes and a light gray plane $\psi 3=2\pi \u2212\psi 1\u2212\psi 2$. The intersection of all planes corresponds to full phase synchrony. As defined above, the set of incoherent phase configurations are the points where the (global) Kuramoto order vanishes; in phase difference coordinates, these are given by

For arbitrary parameters $\alpha s$, $\alpha n$, $A$, the network (10) has dihedral symmetry $S2\u2240S2\u2261D4$. In phase differences (12), the generators of the symmetry action (11) are

that correspond to the rotation and reflection mentioned above. Solutions transform in a nontrivial way along variable $\psi 2$ that describes the state of the two populations relative to one another under $\gamma r$, $\gamma m$. The system also has parameter symmetries

where (17) is time reversing for $A<0$.

#### 2. Equilibria, invariant subspaces, and linear stability

The four oscillator network (12) supports phase-synchronized (or coherent) solutions. Specifically, there are fixed points

The first fixed point corresponds to the phase configuration where all oscillators are in phase (the point of full isotropy where all phases are equal) [cf. Fig. 2(a)]. The second corresponds to the configuration where both populations are in phase but antiphase with respect to each other [cf. Fig. 2(b)].

Linear stability of the solutions $0$ and $\pi $ is determined by the eigenvalues

respectively; the signs are for $0$ and $\pi $, respectively.

The incoherent phase configurations in $I$ can be parameterized as three lines

such that $I=I0\u222aI1\u222aI2$. The set $I0$ is the phase configuration where each population is incoherent (that is, the phase difference of the two oscillators is $\pi $), and the sets $I1$ and $I2$ correspond to phase configurations where oscillators in distinct populations have a phase difference of $\pi $ [cf. Figs. 2(c) and 2(d)].

The line $I0$ is a continuum of equilibria. Linearizing the equations, we obtain that linear stability is determined by the eigenvalues,

Writing $a=Ks4cos\u2061\alpha s$, $b=Kn4cos\u2061(\alpha n+\phi )$, and $c=Kn4cos\u2061(\alpha n\u2212\phi )$ we have that $\lambda 2,3=2(a\u2213bc)$. This yields the linear stability properties of the equilibria: Phase configurations in $I0$ are linearly stable if $a<0$ and $bc<a2$ or if $a<min{0,\u2212bc}$. How the equilibria bifurcate depends on $\phi $: The eigenvalues $\lambda 2,\lambda 3$ are a complex conjugated pair if $\phi \u2208(|\alpha n\u2212\pi 2|\u2212\pi ,\u2212|\alpha n\u2212\pi 2|)\u222a(|\alpha n\u2212\pi 2|,\pi \u2212|\alpha n\u2212\pi 2|),$ which suggests a (degenerate) Hopf bifurcation for $Kscos\u2061\alpha s=0$. For other $\phi $, the eigenvalues are real. In the limiting case of global coupling ($Ks=Kn:=1$, $\alpha s=\alpha n=\alpha $), the eigenvalues evaluate to

Now consider the incoherent phase configurations $I1$ and $I2$. These are invariant lines of (12) for any parameter values. The dynamics on each of these lines is given by

Thus, $I1$ and $I2$ are continua of nonisolated fixed points if $\alpha s=\alpha n=\xb1\pi 2$ as in the traditional Kuramoto–Sakaguchi equations. Linear stability of these fixed points for $\alpha s=\alpha n=\alpha =\xb1\pi 2$ is determined by the nontrivial eigenvalues

This implies that the set of incoherent phase configurations $I=I0\u222aI1\u222aI2$ are surrounded by sets of periodic orbits when $\alpha =\xb1\pi 2$ and $KnKs>0$ ($|A|<1$). If $\alpha s,\alpha n\u2248\xb1\pi 2$, the invariant lines persist and are surrounded by spiral trajectories.

Finally, we consider the set of *two-cluster synchronized phase configurations* $S$, which are the points $S2\xd7S2$ isotropy in the globally coupled system: These are the points where there are two clusters of two oscillators each. In phase difference coordinates on $T3$ for (12), we define

and have $S=S0\u222aS1\u222aS2$ [see Figs. 2(f) and 2(g)]. The set $S0$ consists of the phase configurations where each of the populations is phase synchronized—this means that $S0=P1\u2229P2={\psi 1=0}\u2229{\psi 3=0}$ (cf. Fig. 3). The sets $S1,S2$ are phase configurations where two oscillators of distinct populations are phase synchronized. The set $S$ is dynamically invariant for any values of $\alpha s$, $\alpha n$, and $A$; $S0$ is the fixed-point set of a subgroup of $D4$. Now, the dynamics on $S0$ are determined by

For $\alpha n\u2260\xb1\pi 2$, there are exactly two fixed points on this line that corresponds to $0$ and $\pi $, respectively. For $\alpha n=\xb1\pi 2$, the set $S0$ is a continuum of fixed points. Linearization yields the eigenvalues

Note that the invariant planes $P1={\psi 1=0}$ and $P2={\psi 3=0}$ bound the phase differences within each population (see Fig. 2 for the corresponding phase configurations). Thus, $\Omega 11,12=0$ and $\Omega 21,22=0$. By contrast, the phase difference $\psi 2$ between the two populations may not be bounded—this corresponds to weak chimeras. In Secs. III–VI, we will elucidate the dynamical mechanisms that lead to such solutions.

## III. GRADIENT DYNAMICS FOR ODD COUPLING

If the coupling function is a pure sine function ($\alpha s=\alpha n=0$), system (10) is a gradient system.^{38} More precisely, with the potential

Eq. (10) can be written as

for $\sigma ,k\u2208{1,2}$. The system has a parameter symmetry given by the action of $\gamma (A)$ as defined in (17). This allows us to restrict the parameter range to $A\u2208(\u22121,1)$ and makes the bifurcation behavior of the system for the special parameter values $A=\xb11$ (corresponding to uncoupled populations and a ring topology, respectively) more transparent. Note that for these parameter values $\gamma (\u2212)$ is an additional $Z2$-symmetry.

### A. Equilibria and their stability

We first analyze the equilibria and their stability for $\alpha s=\alpha n=0$. According to (20), the stability of the coherent equilibrium $0$ is determined by the eigenvalues $\lambda 1=(A\u22121)/2,$ $\lambda 2,3=\u22121/2.$ Thus, $0$ is an attractor for $A<1$ and a saddle equilibrium otherwise. Similarly, the stability of the equilibrium $\pi $ is determined by $\lambda 1=(1\u2212A)/2,$ $\lambda 2,3=\u2212A/2$. Thus, $\pi $ is a source for $A\u22640$, a saddle for $A\u2208(0,1)$, and a sink for $A\u22651$ (with neutral linear stability along $I0$ for $A=1$). Finally, the linearization of the vector field at the incoherent invariant line $I0={(\pi ,\phi ,\pi )|\phi \u2208T}$ has eigenvalues (22), which evaluate to

Therefore, different points $I0$ have different transverse stability for certain values of the parameter $A$. Specifically, for $A<\u22121$, the equilibria $(\pi ,\phi ,\pi )$ are stable if

where $\phi \u2217=arccos(1+A|1\u2212A|).$ Decreasing $A$ from zero, the equilibria on $I0$ are (transversely) stable around $\phi =\xb1\pi 2$; the part of $I0$ with transverse stability increases as $A$ is decreased. Thus, there is multistability of the equilibrium $0$ and segments of the manifold $I0$ if $A<\u22121$. Conversely, for $A\u2208(\u22121,0)$ equilibria on $I0$ are repellers if $\phi \u2208(\phi \u2217\u2212\pi ,\u2212\phi \u2217)\u222a(\phi \u2217,\pi \u2212\phi \u2217)$. For $A>0$, the entire line $I0$ is repelling. In all other cases, the points of $I0$ are of saddle type.

The system has a number of equilibria apart from $0,\pi $, and the continuum $I0$. First, the points $(\pi ,0,0)$, $(0,0,\pi )$, $(\pi ,\pi ,0)$, and $(0,\pi ,\pi )$—two-cluster configurations of one and three oscillators [cf. Figs. 2(j) and 2(k)]—are equilibria for any $A$ when $\alpha s=\alpha n=0$. Linear stability is determined by the eigenvalues

These equilibria are of saddle type since $Kn>0$ implies $\lambda 1\lambda 2<0$, $Kn<0$ implies $\lambda 1\lambda 3<0$, and $Kn=0$ implies $\lambda 2\lambda 3<0$.

Second, for $A>0$, the system has four symmetric equilibria $(\u22122\psi \u2217,\psi \u2217,0)$, $(2\psi \u2217,\u2212\psi \u2217,0)$, $(0,\psi \u2217,\u22122\psi \u2217)$, and $(0,\u2212\psi \u2217,2\psi \u2217)$ with $\psi \u2217=arccos\u2061(A\u22121A+1)$. These equilibria are of saddle type and move along a straight line as the parameter $A$ is varied. Pairs of these saddle points (in each plane $\psi 1=0$, $\psi 3=0$) disappear in pitchfork bifurcations with $\pi $ at $A=0$.

### B. Bifurcations and integrability

The system bifurcates at $A=\u22121$, $A=0$, and $A=1$. In each of these cases, the dynamics have additional symmetries and there can be conserved quantities.

If $A=0$, the system is fully symmetric and the well-known Kuramoto model for identical oscillators. The system has a global attractor $0$ and repellers $I=I0\u222aI1\u222aI2$ consisting entirely of equilibria. Other equilibria with phase difference coordinates $0$ and $\pi $ are saddles in this case. Breaking the full permutational symmetry for $A\u22600$ leads to the disappearance of fixed points from invariant varieties $I1$ and $I2$ but these lines remain invariant.

For $A=\xb11$, the system has conserved quantities. If $A=1$, the oscillator populations are uncoupled and there is a conserved quantity. We have

The proof is a computation that is given in Appendix A. Moreover, for the specific gradient cases, we have an additional constant of motion; thus, trajectories are completely determined.

If $A=\u22121$, there is no coupling between distinct populations and we have a network of identical Kuramoto oscillators on a ring. This system also has a conserved quantity.

Note that the constants of motion $H(1,\u22c5)$ and $H(\u22121,0)$ do not depend on the interpopulation phase difference $\psi 2$. Thus, all solutions belong to cylinder-like surfaces. Figure 4 shows the projections of surfaces determined by $H(1,\u22c5)(\psi 1,\psi 3)=C$ and $H(\u22121,0)(\psi 1,\psi 3)=C$ for fixed $C\u2208[\u22121,1]$ onto the plane $(\psi 1,\psi 3)$ for different values of parameter $C$. In addition, for $A=1$, all solutions of the system belong to the parallel planes $H(1,0)=\beta $, $\beta \u2208T$: Each regular trajectory starts from a point of manifold $I0$ and tends to a point of the set $\psi 1=\psi 3=0$ [see Fig. 4(c)]. The manifold $S$ consists completely of fixed points in this case. For $A=\u22121$, the whole manifold $I$ consists of saddle equilibria (transverse to $I$). In addition, the system has two more invariant lines of degenerate saddle points given by $\psi 1=\pi $, $2\psi 2+\psi 3=0$ and $\psi 1+2\psi 2=0$, $\psi 3=\pi $. Stable and unstable one-dimensional invariant manifolds of each saddle belong are contained in the same cylinder $H(\u22121,0)=C$ for fixed $C$ (again, the saddles are neutral in the third direction). All trajectories that are not equilibria or separatrices connect the source $\pi $ to the attractor $0$ [see Fig. 4(d)]. In particular, these trajectories on each cylinder $H(\u22121,0)=C$ are bounded by one-dimensional invariant manifolds of two saddle points.

To summarize, for the gradient case, we observe the following stability/bifurcation behavior. For $A<0$, the system has multistability of $0$ and two parts of the manifold $I0$ described in Sec. III A. Moreover, the equilibrium $\pi $ is a repeller. For $A\u2208(\u22121,0)$, the system has only one attractor $0$ and repellers $\pi $ and parts of $I0$. For $A>0$ the entire line $I0$ becomes a repeller. For $A\u2208(0,1)$, the equilibrium $0$ is the only attractor. Finally, if $A>1$, the equilibrium $\pi $ is an attractor and $0$ is repelling. One-dimensional manifolds of equilibria only exist at the bifurcation values of $A$; except for $I0$, these do not persist as the degeneracies are broken.

## IV. FROM CONSERVATIVE DYNAMICS FOR EVEN COUPLING TO CHIMERAS

Weak chimeras as solutions where the two populations have distinct frequencies arise for phase-lag parameters close to $\xb1\pi 2$ (cf. Refs. 8, 9, and 12). Here, we elucidate the mechanisms that lead to such solutions. We mainly consider the singular case $\alpha s=\alpha n=\alpha =\xb1\pi 2$: then, the vector field is conservative (Hamiltonian-like). Indeed, if the coupling function is even (this is the case if $\alpha s=\alpha n=\xb1\pi 2$), system (12) is divergence-free, that is, if $G$ denotes the right hand side of (12), then

for any values of $Ks$ and $Kn$ (see also Ref. 33). This means that the system does not have any attractors or repellers and there are one-parameter families of fixed points and two-parametric families of periodic orbits or families of homo/heteroclinic cycles. Moreover, the system has the time-reversing symmetry given by

with fixed point space $Fix\u2061(\gamma (t))={(\psi 1,\psi 2,\psi 3)|\psi 1=\psi 3}$. System (12) with parameters $A$, $\alpha s=\xb1\pi 2$, $\alpha s=\u2213\pi 2$ is equivalent to the same system with parameters $A\u22121$, $\alpha s=\alpha n=\xb1\pi 2$ due to the parameter symmetry

### A. Phase space and integrability

Consider the case that the phase-lag parameters $\alpha n,\alpha s$ are identical, that is, $\alpha :=\alpha n=\alpha s=\xb1\pi 2$. For $A=0$ ($Kn=Ks$), the system is fully symmetric and the phase space is organized into the canonical invariant region and its symmetric copies (cf. Sec. II B). Specifically, the sets $S1$, $S2$ are dynamically invariant as points with isotropy $S2\xd7S2$. For $\alpha =\xb1\pi 2$, these sets remain invariant even for nonidentical coupling strength $A\u22600$ (this is not true in general): They form continua of equilibria whose linear stability is determined by the eigenvalues,

The zero eigenvalue corresponds to the direction along $S1$ or $S2$, respectively, and the equilibria are degenerate saddles or degenerate centers depending on the sign $KsKn$ (or $1\u2212A2$).

By direct calculation, one can verify the existence of a first integral (see also the constants of motion in Refs. 6 and 34).

The existence of the conserved quantity implies that the phase space is organized by invariant sets

parameterized by $C\u2208(0,1)$. These are cylindrical surfaces if lifted to $R3$ and two-dimensional tori in $T3$; in the following, we will simply refer to $L(C)$ as *invariant cylinders*. Note that we can also parameterize the invariant cylinders by their diameter $d=42arccos\u2061(C)$ (at $\psi 3=\psi 1$). One of these cylinders is shown in Fig. 3(c) and typical projections of such cylinders on the $(\psi 1,\psi 3)$-plane for different values of $C$ are shown in Fig. 3(d). In the limiting case of $C=1$, we have $L(C)=I0$. If $C=0$, it corresponds to the “square” cylinder of invariant planes $\psi 1=0$, $\psi 1=\pi $, $\psi 3=0$, and $\psi 3=\pi $. For $C\u2208(0,1)$, the set $L(C)$ contains eight equilibria, which correspond to the intersection of $L(C)$ with the sets $I1$, $I2$ and $S1$, $S2$: With the parameterization of $I,S$ as above, the intersection is at $\phi =2arcsin\u2061(C)$ and we have

The intersections are shown in Fig. 3(c): Each invariant line $I\u2113$, $\u2113=1,2$, intersects the cylinder $L(C)$ twice, at the point $I\u2113$ and the symmetric point $I\u2113\u2032=\u2212I\u2113$ (or $I\u2113\u2032=2\pi \u2212I\u2113$ in $R3$) for given $C\u2208(0,1)$. The invariant line $S\u2113$ intersects the cylinder $L(C)$ at points $S\u2113$ and its symmetric counterpart $S\u2113\u2032=\u2212S\u2113$. In the limiting case $C=0$, the points $I1$, $I1\u2032$, $I2$, $I2\u2032$ correspond to $\pi $ (and its symmetric copies) while $S1$, $S1\u2032$, $S2$, $S2\u2032$ correspond to the fully synchronized phase configuration $0$.

### B. The emergence of chimera-like solutions

The global dynamics are determined by the dynamics on the invariant cylinders. To obtain coordinates on the invariant cylinder (torus), write the variables $\psi 1$, $\psi 3$ in polar coordinates $\rho $, $\varphi $. Now, the two-dimensional dynamics on $L(C)$ can be expressed in coordinates $(\psi 2,\varphi )$; these depend on the value of $C$ and may bifurcate as $C$ is varied. Note that the dynamics for any $C$ is present since $C$ parameterizes the family of invariant tori that foliates phase space rather than being a system parameter.^{39,40} First consider the special case $C=0$; in this case, $L(C)=L(0)$ corresponds to the degenerate cylinder that consists of the invariant planes. The dynamics on $L(0)$ depend on the coupling parameter $A$ as shown in Figs. 5(a)–5(c): For $A\u2208(0,1)$, there are continua of heteroclinic trajectories from equilibria in $S0$ as shown in Fig. 5(a); these are bounded by homoclinic trajectories from $0$ to itself (on the torus) and degenerate to the equilibrium $\pi $. Moreover, there are families of periodic orbits with a nontrivial winding number (shaded areas). For $A=0$, these families of periodic orbits disappear as the stable and unstable manifolds of $0$ form a heteroclinic “web” on $L(0)$ [see Fig. 5(b)]. Finally, for $A\u2208(\u22121,0)$, there are homoclinic/heteroclinic trajectories from (points in) $S0$ to itself. For fixed $A\u2208(0,1)$, the bifurcation scenario is similar with $0$ replaced by $Sj,Sj\u2032$, $j=1,2$ and $\pi $ replaced by $Ij,Ij\u2032$, $j=1,2$ [see Figs. 5(d)–5(f)]. Define the critical cylinder size as

(Here, we used the symmetry $\gamma (A,t)$ to obtain the formula for $A\u22651$.) For $C\u2208(0,C\u2217)$, we have families of periodic orbits where the populations rotate relative to each other. For $C=C\u2217$, we have a heteroclinic web between $Sj,Sj\u2032$. Finally, for $C\u2208(C\u2217,1)$, there are periodic frequency-locked solutions.

### C. No coupling between populations

If $A=1$ ($Kn=0$), the network consists of two uncoupled populations that evolve independently of one another. In this case, $\gamma (A,t)$ acts as a symmetry of the system since it maps $A\u21a6A\u22121$. The system (12) reduces to

and its solutions form a continuum of lines

parameterized by the initial conditions $(\psi 1(0),\psi 2(0),\psi 3(0))=(\psi 10,\psi 20,\psi 30)$. The dynamics are schematically shown in Fig. 3(b). For most initial conditions, the phase difference $\psi 2$ between the populations is increasing or decreasing monotonously. Specifically, we have that each population is frequency synchronized and they rotate relative to one another with $\Omega 2,j\u2212\Omega 1,k=cos\u2061(\psi 3(0))\u2212cos\u2061(\psi 1(0))$, $j,k=1,2$. Thus, $\psi 2$ increases monotonically if $|\psi 3|>|\psi 1|$ for $\psi 1,\psi 3\u2208[\u2212\pi ,\pi ]$ and $\psi 2$ decreases monotonically if $|\psi 3|<|\psi 1|$ for $\psi 1,\psi 3\u2208[\u2212\pi ,\pi ]$. The direction of the flow is shown in Fig. 3(b) in the cube $[0,2\pi ]3$. On the planes $\psi 3=\xb1\psi 1$, the dynamics are trivial (there are nonisolated fixed points). Note that the invariant lines $I0$ and $S0$ are the intersections of these planes.

For coupling parameters $A$ ($Kn$) close to the singular case of uncoupled populations, the coupling between populations is weak. Frequency unlocked solutions persist but now deviate from straight lines. This gives the frequency-unlocked solutions on the invariant cylinders $L(C)$ for $C<C\u2217=C\u2217(A)$ and given $A\u2208(0,1)$ described above. For $|A|>1$, the parameter symmetry $\gamma (A,t)$ implies that system (12) has the same phase portraits but with the shift $\psi 2\u21a6\psi 2+\pi $ and a potential time reversal. The symmetry also means that the saddles $S\u2113$, $S\u2113\u2032$ swap places with centers $I\u2113$, $I\u2113\u2032$. Local bifurcations occur at each fixed point on the invariant manifolds $I0$, $S0$ of the system (two eigenvalues of each point transform from pure imaginary to real symmetric or vice versa). More globally, all invariant manifolds of the fixed points $S\u2113$ swap with those of $I\u2113$, $\u2113=1,2$. Figure 6 schematically shows the bifurcation that occurs in the narrow stripe of phase space $(\psi 2,\varphi )\u2208L(C)$ at the bifurcation point $A=1$: There is a heteroclinic cycle between the saddle $S\u2113$ ($S\u2113\u2032$) and its copy with coordinate $\psi 2+2\pi $, which bounds a family of periodic orbits around the center $I\u2113$ ($I\u2113\u2032$) [see Fig. 6(a)]. At the bifurcation point $A=1$, the phase space is foliated by straight lines as discussed above [see Fig. 6(b)]. For $A>1$, a new heteroclinic loop appears between the saddle $I\u2113$ ($I\u2113\u2032$) that bounds periodic orbits centered at $S\u2113$ ($S\u2113\u2032$). Frequency-unlocked chimera solutions exist for all $A>1$.

We can now describe the fate of the frequency-unlocked chimera solutions in more detail. At the bifurcation point $A=1$, the whole phase space is filled with frequency-unlocked chimera solutions except on the invariant planes $\psi 3=\xb1\psi 1$. Indeed, there are two “chimera tubes” with opposite directions of trajectories. As the parameter $A$ deviates from one, narrow cylinders with phase-locked periodic orbits appear. The heteroclinic trajectories form the boundaries between the region of frequency locking and a region of frequency-unlocked chimera solutions. Tracing these trajectories along the invariant cylinders $L(C)$ that foliate phase space for $C\u2208(0,C\u2217)$, we obtain a tube of phase-unlocked chimera solutions (cf. Fig. 7). Since the surface bounding the tube resembles a snake, we will refer to it as a *serpentine chimera region*. Figure 7(c) shows boundary surface of the tube schematically (“skin of the chimera-snake”) that consists of piecewise smooth surfaces. As $A\u21920$, the tube degenerates to a single heteroclinic connection with a nontrivial winding number. Given the parameter symmetry $\gamma (A,t)$, we can conclude that there are phase-unlocked chimera solutions for any $A>0$.

### D. No coupling within populations

Finally, consider the case $A=\u22121$ when the oscillators within each population are uncoupled ($Ks=0$)—they still remain coupled indirectly through the interaction with the other population. In this case, $\gamma (A,t)$ acts as a time-reversing symmetry since it maps $A\u21a6A\u22121$. A direct calculation shows:

Note that $H(\u22121,\pi 2)=H(1,0)$ in Proposition 2.

Thus, for $A=\u22121$, there are two conserved quantities, $H(\u22c5,\pi 2)$ and $H(\u22121,\pi 2)$, and, thus, solutions lie on the intersection of the invariant cylinders $L(C)$ with the plane $H(\u22121)=C(\u22121)$ for $C\u2208(0,1)$, $C(\u22121)\u2208T$, given by

Similar to the case $A=1$, local and global bifurcations occur when $A=\u22121$. The time reversing symmetry organizes the solutions on the invariant cylinders (28): Invariant lines of saddle fixed points $S\u2113$ mutually change the structure of their fixed points with the lines of center fixed points $I\u2113$, $\u2113=1,2$. Global bifurcation with phase-locked periodic orbits (lines in phase space $(\psi 2,\varphi )\u2208R2$) for $A=\u22121$ is of the same type that the global bifurcation with the appearance of the family of straight lines in the described case $A=1$. In the case $A=\u22121$, global bifurcations occur on the square $K(0,0)\u2208T3$ (with vertices at the points $0$, $\pi $ and two their copies) instead of $S0$ as for $A=1$.

## V. DISSIPATIVE DYNAMICS: BIFURCATIONS AND TRANSITIONS TO CHAOS

In the conservative case, a limit cycle solution with a nontrivial winding number along the phase difference $\psi 2$ (the phase difference between the two populations) arises in phase space. How do these solutions bifurcate as parameters are varied and the additional structures are broken? In the following, we will refer to any limit cycle solution with a nontrivial winding number as a weak chimera.

### A. From conservative to dissipative dynamics

Before we turn to chimeras, we consider the overall organization of phase space by invariant subsets and qualitatively describe what dynamics are possible. The conservative dynamics for $\alpha s=\xb1\alpha n=\xb1\pi 2$ correspond to a bifurcation point. Recall that for these parameter values, the sets $I1$, $I2$ and $S1$, $S2$ are continua of equilibria; in particular, the linearization at each equilibrium has a zero eigenvalue with an eigenvector in the direction of the corresponding set. Away from the bifurcation point, there are only two equilibria for each set: These are ${\pi ,(\pi ,0,\pi )}\u2282I1$, ${\pi ,(\pi ,\pi ,\pi )}\u2282I2$, ${0,(\pi ,0,\pi )}\u2282S1$, and ${0,(\pi ,\pi ,\pi )}\u2282S2$. Before the bifurcation point, one of the two equilibria is repelling along the set while the other one is attracting. After the bifurcation, the situation is reversed. Hence, the bifurcation corresponds to reversing the flow along the four sets in a degenerate bifurcation. The stability transverse to the invariant set determines the dynamics of trajectories nearby: If there are complex conjugate eigenvalues, nearby trajectories spiral around the set as the conservative structure is broken.

Away from the bifurcation point $\alpha s,\alpha n=\xb1\pi 2$, the invariant cylinders [Fig. 3(c)] break up, which allows for more intricate dynamics. For small deviations from the bifurcation parameter $\alpha s=\xb1\pi 2\u2212\epsilon $, there is a slow drift transverse to the cylinders $L(C)$ (which are not invariant anymore). Since the families of heteroclinic orbits that separate frequency-locked and frequency-unlocked solutions also disconnect, it is possible for trajectories to come close to the continuum of equilibria $I0$ as well as the degenerate cylinder $L(0)$ that consists of the invariant planes $\psi 1=0$ and $\psi 3=0$.

Taken together, typical trajectories can exhibit dynamics that explore a large region of phase space beyond the bifurcation point. Trajectories can move from the vicinity of the points $\pi $ ($0$), moving along the invariant line $I1$ or $I2$ ($S1$ or $S2$) in a spiraling fashion to approach $I0$, before then drifting along the level sets $L(C)$ into the region where the phase difference between populations increases and approaching a frequency-unlocked solution on the invariant set $L(0)$. In the following, we will consider the bifurcations of these asymptotic solutions.

### B. Bifurcations on invariant manifolds and flat chimeras

We first consider the weak chimeras on the invariant planes $P1,P2$ defined by $\psi 1=0$ or $\psi 3=0$ (these correspond to the limiting cylinder $L(0)$ in Sec. IV B); we will refer to these solutions as *flat chimeras*. These correspond to in-phase synchronization of one population while the other population is approximately in anti-phase.^{8,12} The relative phase difference between the two populations increases without bound as $\Omega 1,2\u22600$. Without loss of generality, we assume that $\psi 3=0$, that is, the second population is phase synchronized. Dynamics (12) on the invariant subspace restrict to a two-dimensional system; therefore, no chaotic flat chimeras are possible. Recall that families of flat chimeras arise for $\alpha s=\xb1\pi 2$, $\alpha n=\xb1\pi 2$ as described above [cf. Fig. 5(a)]. This family of flat chimeras is usually bounded by homoclinic or heteroclinic cycles [Fig. 5(b)]. These flat chimeras are neurally stable transverse to the invariant plane $\psi 3=0$ due to the (partial) integrability of the system.

Figure 8 elucidates the bifurcation structure for parameter values where the conservative structure is broken. Note that because of the symmetries, the system can have only an even number of fixed points on the plane $\psi 3=0$. The system has 2, 4, or 6 fixed points depending on parameter values. Saddle-node and saddle-connection bifurcations must occur in pairs simultaneously. Specifically, the panels in Figs. 8(a)–8(p) show all possible phase portraits for the dynamics of $(\psi 2,\psi 1)\u2208T2$; these are arranged to show bifurcation transitions from (a) to (b), from (b) to (c), and so on. The system exhibits the following bifurcations on the plane: A pitchfork bifurcation of $0$ or $\pi $ transverse to the invariant line $\psi 1=0$ [Figs. 8(c)–8(d), 8(k)–8(l), 8(l)–8(m), and 8(o)–8(p)]; two simultaneous saddle-node bifurcations on the flat chimera limit cycle [Figs. 8(a)–8(c) and 8(i)–8(k)]; two simultaneous saddle-node bifurcations away from the chimera [Figs. 8(o)–8(n)]; a subcritical pitchfork bifurcation of limit cycles [Figs. 8(h)–8(g)]; a supercritical pitchfork bifurcation of limit cycles [Figs. 8(g)–8(h)]; two simultaneous saddle-connection bifurcations [Figs. 8(f)–8(g) and 8(m)–8(n)].

We highlight the bifurcations that relate to the creation and destruction of stable flat chimeras. First, there are two simultaneous saddle-node bifurcations on an invariant circle that lead to a limit cycle solution with nonzero winding number, a flat chimera [Figs. 8(a)–8(c) and 8(k)–8(m)]. Second, there are pitchfork bifurcations of limit cycles, both sub- [Figs. 8(h) and 8(i)] and supercritical [Figs. 8(p) and 8(q)]; these bifurcations can stabilize flat chimeras (within the invariant subspace). Third, there are simultaneous saddle-connection bifurcations of saddle equilibria that lead to the emergence of two symmetric limit cycles with unbounded phase difference $\psi 2$ between populations; the resulting limit cycles can be stable [Figs. 8(f) and 8(g)] or unstable [Figs. 8(m) and 8(n)].

Note that the stability of flat chimeras in the full system (12) depends on the transverse stability with respect to the third direction.

### C. Chaotic weak chimeras

Two populations of two phase oscillators with Kuramoto–Sakaguchi coupling support chaotic weak chimeras; the bifurcation scenario, obtained numerically,^{41} is indicated in Fig. 1. In the following, we describe the bifurcations that lead to the emergence of chaotic weak chimeras in more detail. Since $|D4|=8$, every point with trivial isotropy has eight images under the symmetry action. Recall that $\Sigma (A)$ are the symmetries that preserve the set $A$. If $A\eta $ is a family of limit sets indexed by a parameter $\eta $ and $\Sigma (A\eta )$ changes as $\eta $ is varied, we have a *symmetry-increasing bifurcation*.^{42}

We describe the bifurcations as a single parameter is varied; for concreteness, we fix $A=0.7$ and $\alpha n=0.44$ and vary $\alpha s$. For these parameters, there exists a *flat chimera* (limit cycle) on $P1\u222aP2$ that consists of points with nontrivial isotropy as the plane is invariant under a reflection: Consider the stable limit cycle shown in Figs. 8(k) and 8(n) within $P\sigma $. A numerical calculation shows that it is also stable transverse to $P\sigma $. This limit cycle loses (transverse) stability in a pitchfork bifurcation of limit cycles which leads to two stable limit cycles with trivial isotropy that are mapped to each other through the reflectional symmetry. [Note that other stable limit cycles on $P\sigma $, such as those shown in Fig. 8(h), which are unstable transversely and, therefore, lead to saddle limit cycles in a pitchfork bifurcation.] Since the bifurcations happen simultaneously on $P1$ and $P2$, there is a total of four such solutions. The projections of four chimera states into the $(\psi 1,\psi 3)$-plane are shown in Fig. 9(a). Note that each of the resulting limit cycles has a setwise reflectional symmetry. This symmetry is broken in a supercritical pitchfork bifurcation of limit cycles at $\alpha s\u22481.616$, which leads to a total of eight nonsymmetric weak chimera limit cycle solutions [cf. Figs. 9(a) and 9(b)].

Each stable limit cycle loses stability in a period-doubling bifurcation [see Figs. 9(b) and 9(c)]. The resulting limit cycle bounds a Möbius strip that wraps around the torus in the $\psi 2$ direction; the original limit cycle is of saddle type after the bifurcation and lies at the center of the strip. Note that the width of the Möbius strip varies along variable $\psi 2$ and for different parameter values. This indicates the heterogeneity of the attraction strength along the periodic solution that will have a further impact on the formation of chaos. The first period-doubling bifurcation is followed by a chain of period-doubling bifurcations: The $n$th period-doubling bifurcation leads to a $2n$-limit cycle as shown in Figs. 9(c)–9(f). At each period-doubling bifurcation, the geometry of the Möbius strips becomes more elaborate, allowing trajectories to follow different directions as the trajectory wraps around the torus in $\psi 2$ direction.

A chaotic attractor with nontrivial winding number (i.e., a chaotic weak chimera) emerges as a result of a period-doubling cascade as shown in Fig. 9(g). We estimate the fractal dimension of the chaotic attractor to be slightly larger than two. Roughly speaking, as trajectories wind around $\psi 2$ on the attractor, they can either take a direct path, closer to the original limit cycle, or make an “excursion” in the $\psi 1$ or $\psi 3$ direction. The Poincaré section shown in Fig. 10(a) shows the finer structure of the attractor. The attractor undergoes a symmetry-increasing bifurcation as two symmetrically related attractors merge [Figs. 9(g) and 9(h)]; see also Ref. 24 and note that a similar effect has also been observed in oscillators with amplitude variations.^{22,25} The enlargement of the chaotic attractor can also be seen in the Poincaré section in Fig. 10(b); note that the section does not necessarily respect the attractor symmetry.

As the parameter $\alpha s$ is increased beyond $\alpha s=1.64166$, there is multistability between the equilibrium $0$ and the four stable symmetric chaotic weak chimeras. Moreover, there are two saddle equilibria $S\u2113$ on each of the invariant planes $P1,P2$, as shown in Fig. 8(h) or Fig. 8(m). The stable manifold of each saddle is two-dimensional, intersects the appropriate invariant plane transversely, and separates the basin of attraction of $0$ and the chaotic attractors. Finally, the chaotic attractor is destroyed through a homoclinic tangency^{43–46} and trajectories eventually converge to $0$. The transients in the vicinity of the former chaotic attractor can be very long before the stable equilibrium is approached.

### D. Minichimerapedia for networks of two populations of two phase oscillators

To summarize, we have given an overview of the solutions of (12) that correspond to weak chimeras for two coupled populations of phase oscillators (10). All of these solutions share the property that the phase difference between the populations ($\psi 2$ in the reduced dynamics) is unbounded as time evolves. As they wrap around the torus, the solutions are nonhomological to zero and must arise in global bifurcations as described above.

Limit cycle solutions on the invariant planes $P1={\psi 1=0}$ or $P2={\psi 3=0}$ [Figs. 8(k)–8(n)], referred to as flat chimeras. The situation corresponds to phase synchronization of one of the populations with local order parameter $|Z1|=1$ (or $|Z2|=1$) and $|Z2|=|Z2(t)|\u2208(0,1)$ [or $|Z1|=|Z1(t)|\u2208(0,1)$].

A one-parameter family of neutral periodic orbits on the invariant plane for $\alpha s=\xb1\pi 2$, $\alpha n=0$, $\alpha n=\xb1\pi 2$, $\alpha n=\pi $ [Fig. 11(a)].

The limit cycles with setwise symmetry [Fig. 9(a)].

The limit cycles without symmetry [Figs. 9(b)–9(f)].

A two-parameter family of neutrally stable periodic orbits (parameters as in item 2) [cf. Figs. 3(c), 5(a), and 5(d)].

The eight nonsymmetric chaotic attractors [Figs. 9(g), 10(a), and 10(c)].

The four symmetric chaotic attractors that appear as the result of the symmetry-increasing bifurcation [Figs. 9(h), 10(b), and 10(d)].

Note that trajectories close to homoclinic/heteroclinic orbits that are nonhomological to zero [see Fig. 8(f)] can also show (transient) frequency-unlocked dynamics between the populations.

## VI. COEXISTENCE OF CONSERVATIVE AND DISSIPATIVE DYNAMICS

Finally, we remark that the system also has an interesting and unusual dynamics when $\alpha s=\xb1\pi 2$, $\alpha n\u2208{0,\pi}$—this corresponds to the function of the coupling function $gs$ that determines the coupling with populations being even and the interaction function $gn$ between populations being odd. In this case, the system has the first integral given in Proposition 3; for simplicity, we use the same notation here.

For $\alpha s=\pi 2$, $\alpha n=0$ and arbitrary $A$ the system has the first integral $H(\u22121,0)$ as defined in (25).

Thus, the dynamics of the system evolve on the level sets $H(\u22121,0)(\psi 1,\psi 3)=C$ for fixed $C\u2208[\u22121,1]$ [see Fig. 4(b)]. In the cases $C=\xb11$, we have a planar system on the invariant planes $P1={\psi 3=0}$ and $P2={\psi 1=0}$. As $|C|$ is decreased, the level sets deform continuously and limit to dynamics on the half-planes $\psi 3=\psi 1$ for $\psi 1\u2208[0,\pi ]$ and $\psi 3=2\pi \u2212\psi 1$ for $\psi 1\u2208[\pi ,2\pi ]$ (or their symmetric counterparts) for $C=0$.

This system has coexistence of conservative (Hamiltonian-like) and dissipative dynamics in phase space for $\alpha s=\pi 2$, $\alpha n=0$, and $A\u2208[A\u2217,1/A\u2217]$ with $A\u2217\u22480.38146$. Conservative–dissipative dynamics are often related to the presence of time-reversing symmetries; such dynamics has been described, for example, in a three-dimensional laser system,^{47} globally coupled superconducting Josephson junction arrays,^{48} chains of locally coupled phase oscillators,^{49} and circulant networks of phase oscillators with skew-symmetric coupling.^{50} Here, the system has the conservative region filled with a two-parameter family of neutral periodic orbits (weak chimeras) and this region is bounded with the surface of one-parameter family of heteroclinic cycles forming a “heteroclinic tube.” Outside the Hamiltonian-like region, the dynamics is dissipative with attractor $0$ and repeller $\pi $ when $A<1$ or attractor $\pi $ and repeller $0$ when $A>1$. The boundary surface of the three-dimensional conservative region has a structure similar to that shown in Fig. 7(a). The system has two one-parameter families of saddle points (which are curved lines compared to the case $\alpha s=\alpha n=\xb1\pi 2$). There are heteroclinic (homoclinic in $T3$) cycles that consist of these saddles and their stable and unstable invariant manifolds on the same level surface (each saddle is neutral in transversal towards the level surface). Note that while for fixed $\alpha s=\pi 2$ the dynamics for $\alpha n=\pi 2$ (see above) and $\alpha n=0$ both take place on neutrally stable cylindrical surfaces, these surfaces have essentially different shapes [compare Figs. 3(d) and 4(b)]: Invariant lines that consist of neutral saddles are no longer straight lines (as $I1$ and $I2$ for $\alpha n=\pi 2$) but they deform when changing the parameter $A$.

The theory of bifurcations without parameters^{39,40} can also be used to study the dynamics of the current case. We can consider two-dimensional dynamics on surfaces $H(\u22121,0)=C$ using $C$ as a parameter. This allows us to analyze conservative–dissipative dynamics on an individual level set and then extrapolate to the entire phase space $T3$ as $C$ parameterizes the level sets and attractors and repellers are located on the common boundary of all level sets $\psi 1=\psi 3=0$. Figure 11 shows schematic phase portraits on surfaces $H(\u22121,0)=C$ for different $A$ and $C$. The conservative region on the cylindric surface is the biggest for $C=\xb11$ for any $A=1$ and it shrinks with decreasing $|C|$ to the heteroclinic cycle [transition from (a) to (c) in Fig. 11]. Figure 11(d) corresponds to the limit case $C=0$ with degenerate (green) line $I0$ of saddles in the entire phase space $T3$. There are two conservative regions close to the straight lines $(0,\psi 2,\pi )$ and $(\pi ,\psi 2,0)$ and also one common dissipative region in $T3$.

For $A=1$, the three-dimensional conservative region corresponds to the whole phase space $T3$. Finally, the conservative region is destroyed in a saddle-connection bifurcation when the connections of stable and unstable invariant manifolds of each saddle are simultaneously broken. This happens either when one of the equalities $\alpha s=\xb1\pi 2$, $\alpha n=0$ is violated or when the parameter $A$ leaves the interval $[A\u2217,1/A\u2217]$: If $A$ reaches the bifurcation value $A\u2217$ or $1/A\u2217$ [Fig. 8(c)], the conservative region collapses onto the one-dimensional heteroclinic cycle (between two saddles that belong to $\psi 1=0$ or $\psi 3=0$).

## VII. DISCUSSION

Here, we considered the properties of phase oscillator networks that consist of $M$ populations of $N$ phase oscillators. While we analyzed networks for $M=N=2$ with sinusoidal coupling in details, some of the observations hold also for larger networks with more general coupling. We, therefore, briefly discuss a few general properties of (1). First, system (1) has dihedral symmetry $DMN$ for any $M$ and $N$. Second, (1) is a gradient system for any odd coupling functions $gs$, $gn$ for all $M$ and $N$. For example, if $gs(\varphi )=gn(\varphi )=g(\varphi )$ satisfies $g(\varphi )=\u2212g(\u2212\varphi )$ then the system has the potential

where $h(\varphi )$ is an even function such that $h\u2032(\varphi )=g(\varphi )$. The gradient system (23) is a special case. Third, in the case of an even coupling function $gs(\varphi )=gn(\varphi )=g(\varphi )$ with $g(\varphi )=g(\u2212\varphi )$, system (1) is divergence-free (this generalizes results in Ref. 33 for $Ks=Kn$). Other properties will be discussed elsewhere.

Our results also shed light on the phase space structure in higher dimensions, for example, $M=3$ populations of $N=2$ oscillators. In this case, the system has a continuous set of neutral chimera solutions for coupling $gs(x)=gn(x)=\xb1cos\u2061(x)$ (or $\alpha s=\alpha n=\xb1\pi 2$). These solutions can be periodic, quasi-periodic or chaotic (similar to ABC flows,^{51,52} chaos that fills an entire torus). Figure 12 shows an example of such chaotic dynamics for $M=3$ populations of $N=2$ oscillators. This suggests that the situation with $(N\u22122)$-parametric neutral chimeras also occurs in the case of arbitrary $M$, $N$ for $\alpha s=\alpha n=\xb1\pi 2$.

The system in (1) with $M=2$ populations and $N>3$ oscillators has been the subject of a number of studies^{17} with small and large finite oscillator systems,^{6,13,53} including the continuum limit^{9,13,29,54} of infinitely many oscillators, some based on low-dimensional descriptions obtained via the Watanabe–Strogatz or Ott–Antonsen reductions.^{34,55,56} For future work, it would be interesting to investigate how our findings extend to the dynamic behavior of such larger system sizes, particularly regarding invariant subspaces and how they are affected by the symmetry-breaking mechanism leading to weak chimera states. Finally, our results are, of course, also relevant for more general oscillator systems that allow for amplitude variation. If oscillators are weakly coupled so that the system possesses an attracting torus such that (1) provides an adequate phase reduction, then one would expect some persistence of the bifurcation scenarios. As nonchaotic weak chimeras have been found in real-world experiments with electrochemical oscillators,^{12} it would be interesting to see the bifurcations outlined here in such systems.

Previous studies found that bifurcation curves leading to chimera states are organized around points in parameter space where $|\alpha |=\pi 2$;^{8,9,12} here, we explained the symmetry-breaking bifurcations generating weak chimeras near these points. It is interesting to note that experiments using mechanical oscillators^{10,11,17} indicate that chimera states emergence in a scenario near resonance—this conjecture is confirmed as one can rigorously show^{11} that resonance in such systems is related to phase lags being $\alpha s,\alpha n=\xb1\pi 2$.

## ACKNOWLEDGMENTS

O.B. acknowledges support of the National Research Foundation of Ukraine (Project No. 2020.02/0089).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Oleksandr Burylko:** Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Erik A. Martens:** Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal). **Christian Bick:** Conceptualization (equal); Investigation (equal); Writing – original draft (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

### APPENDIX A: PROOFS

The proofs of Propositions 1–6 are very similar and rely on direct computation with trigonometric identities. We need to show that the respective constant of motion $H$ does not change over time, that is, $dHdt=0$ along solutions. For completeness, we include the proof of Proposition 1.

#### Proof of Proposition 1

*Proof of Proposition 1*

## REFERENCES

*Chemical Oscillations, Waves, and Turbulence*, Springer Series in Synergetics Vol. 19 (Springer, Berlin, 1984).

*Pattern Formation in Continuous and Coupled Systems. The IMA Volumes in Mathematics and its Applications*, edited by M. Golubitsky, D. Luss, and S. H. Strogatz (Springer, New York, 1999), Chap. 10, pp. 121–135.

*The Symmetry Perspective*, Progress in Mathematics Vol. 200 (Birkhäuser Verlag, Basel, 2002).

*Dynamics and Symmetry*, ICP Advanced Texts in Mathematics Vol. 3 (Imperial College Press, 2007).

*Dynamics Beyond Uniform Hyperbolicity*, Encyclopaedia of Mathematical Sciences (Mathematical Physics) Vol. 102 (Springer Verlag, 2004).

*Sur la topologie des écoulements stationnaires des fluides parfaits*, (Springer, 1965), Chap. Vladimir I. Arnold—Collected Works, Vol. 2, pp. 15–18.