# Chimera Dynamics of Generalized Kuramoto-Sakaguchi Oscillators in Two-
population Networks

Seungjae Lee & Katharina Krischer  Physik-Department, Technische Universität
München, James-Franck-Straße 1, 85748 Garching, Germany<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Chimera dynamics, an intriguing phenomenon of coupled oscillators, is
characterized by the coexistence of coherence and incoherence, arising from a
symmetry-breaking mechanism. Extensive research has been performed in various
systems, focusing on a system of Kuramoto-Sakaguchi (KS) phase oscillators. In
recent developments, the system has been extended to the so-called generalized
Kuramoto model, wherein an oscillator is situated on the surface of an
$M$-dimensional unit sphere, rather than being confined to a unit circle. In
this paper, we exploit the model introduced in [New. J. Phys. 16, 023016
(2014)] where the macroscopic dynamics of the system was studied using the
extended Watanabe-Strogatz transformation both for real and complex spaces.
Considering two-population networks of the generalized KS oscillators in 2D
complex spaces, we demonstrate the existence of chimera states and elucidate
different motions of the order parameter vectors depending on the strength of
intra-population coupling. Similar to the KS model on the unit circle,
stationary and breathing chimeras are observed for comparatively strong intra-
population coupling. Here, the breathing chimera changes their motion upon
decreasing intra-population coupling strength via a global bifurcation
involving the completely incoherent state. Beyond that, the system exhibits
periodic alternation of the two order parameters with weaker coupling
strength. Moreover, we observe that the chimera state transitions into a
componentwise aperiodic dynamics when the coupling strength weakens even
further. The aperiodic chimera dynamics emerges due to the breaking of
conserved quantities that are preserved in the stationary, breathing and
alternating chimera states. We provide a detailed explanation of this scenario
in both the thermodynamic limit and for finite-sized ensembles. Furthermore,
we note that an ensemble in 4D real spaces demonstrates similar behavior.

  *

##  1 Introduction

Within the framework of the concept ‘More is different’ [1] or ‘More than the
sum’ [2], collective behavior of systems of nonlinearly coupled objects has
been extensively explored in diverse interdisciplinary fields of science and
mathematics, with a particular focus on systems of coupled oscillators [3, 4].
One notable phenomenon in systems of numerous oscillators is the emergence of
chimera states [5, 6, 7], which are characterized by the simultaneous presence
of coherence and incoherence through the breaking of symmetry. Chimera states
were originally observed in a system of coupled oscillators on a ring geometry
as a spatiotemporal pattern where the system exhibits both a coherently
oscillating part and an incoherently oscillating part along the ring [8].
Later this coexistence pattern was dubbed a chimera state to highlight its
peculiar characteristics [9].

For the exploration of the essential properties of chimera states, Abrams et
al. exploited a simplified model where identical Kuramoto-Sakaguchi (KS) phase
oscillators are stationed in two-population networks [10]. Note that before
their work, Montbrió, Kurths, and Blasius observed chimera states in two-
population networks, however, with heterogeneous KS oscillators [11]. In these
studies, the intra- and inter-population connections are considered to be all-
to-all, with different coupling strengths where the former is stronger than
the latter. This configuration effectively emulates a nonlocal coupling
scenario on a ring geometry, and since their introduction extensive research
has been conducted on the collective dynamics in two-population networks of
coupled oscillators. Thereby, the two-population network was considered with
identical [10, 12] and heterogeneous [13, 14, 15] KS oscillators. Also,
numerous other variations were taken into account, including phase oscillators
under higher-order interaction and planar oscillators [16, 17, 18, 19, 20, 21,
22, 23]. In addition, chimera states have been ceaselessly studied by adopting
three- and multi-population networks as underlying system topology [24, 25,
26, 27, 28, 29, 30, 31].

Investigations of the collective dynamics of sinusoidally coupled oscillators,
such as Kuramoto oscillators on the unit circle, have relied heavily on
dimension reduction methods. Notably, the Ott-Antonsen (OA) ansatz [32, 33]
and the Watanabe-Strogatz (WS) transformation [34, 35, 36] were extensively
employed, including the exploration of chimera states [37, 38, 39, 40]. The
introduction of the WS or OA ansatz allows for the description of the system’s
macroscopic dynamics, providing an alternative to the Kuramoto order parameter
which is defined from the microscopic phase information of the system. For
details, see Sec. 2.

Meanwhile, Kuramoto oscillators, originally defined on the unit circle in
$\mathbb{C}^{1}$ or $\mathbb{R}^{2}$, have been extended to so-called
generalized Kuramoto oscillators [41, 42, 43, 44, 45, 46, 47]. In these
models, each oscillator is represented as a unit vector confined to the
surface of a higher-dimensional unit sphere, and their collective dynamics,
such as synchronization of the oscillators, has been studied. Furthermore,
dimension reduction methods like the OA and the WS ansatz were also developed
and employed for the study of the generalized Kuramoto oscillators [48, 49,
50, 51, 52, 53, 54, 55]. In this paper, we exploit, in particular, the model
suggested by T. Tanaka in Ref. [55] where the author proposed the generalized
Watanabe-Strogatz transformation as a vector form of a linear fractional
transformation. This projection map allows for the investigation of the
macroscopic dynamics of the system in the thermodynamic limit, considering
both real ($\mathbb{R}^{M}$) and complex ($\mathbb{C}^{M}$) spaces.

The rest of this paper is organized as follows. In Sec. 2, we revisit the
generalized Kuramoto model and the generalized WS transformation introduced in
Ref. [55], necessary for our study. In Sec. 3, we introduce a suitable
coupling matrix that determines the Benjamin-Feir instability point and
corresponds to the phase-lag parameter in the standard Kuramoto-Sakaguchi
model. Then, in Sec. 4, we discuss observable chimera states of the
generalized KS oscillators in two-population networks with a particular focus
on the 2D complex spaces ($\mathbb{C}^{2}$). Therein, the chimera states are
investigated in the thermodynamic limit, using the WS macroscopic variables,
and in the finite-sized ensembles, using the Kuramoto order parameter,
respectively. Furthermore, in Sec. 5, we also provide the scenario of
emergence of chimera states in the 4D real spaces ($\mathbb{R}^{4}$) in the
thermodynamic limit as a comparison to the 2D complex spaces. Finally, we
summarize our results in Sec. 6. Some useful information will be given through
A-C.

##  2 Dimension Reduction Methods

A system of identical Kuramoto-Sakaguchi phase oscillators is governed by

$\displaystyle\dot{\phi}_{j}(t)$
$\displaystyle=\omega+2~{}\text{Im}\bigg{[}g(t)e^{-i\phi_{j}}\bigg{]}$
$\displaystyle=\omega+\frac{2}{N}\sum_{k=1}^{N}\sin(\phi_{k}-\phi_{j}-\alpha)$
(1)

for $j=1,...,N$. Each oscillator $\phi_{j}(t)\in\mathbb{T}:=[0,2\pi]$ is
defined on a unit circle of $\mathbb{C}^{1}$, i.e., by the argument of a
complex number $e^{i\phi_{j}(t)}$ with unit modulus. The oscillator is
influenced sinusoidally by the mean-field forcing
$g(t):=e^{-i\alpha}m(t)\in\mathbb{C}^{1}$. Here, the Kuramoto order parameter
$m(t)$ is defined by

$\displaystyle
m(t):=\frac{1}{N}\sum_{k=1}^{N}e^{i\phi_{k}(t)}\in\mathbb{C}^{1}$ (2)

which serves as the centroid of the phases on the unit circle[56]. The phase-
lag parameter $\alpha\in[0,2\pi]$ determines the Benjamin-Feir instability
point [8, 57], i.e., $\alpha_{\text{BF}}=\frac{\pi}{2}$: for
$\alpha<\alpha_{\text{BF}}$, the synchronized state
$\\{\phi_{j}(t)=\phi_{0}+\Omega t\\}_{j=1}^{N}$ is stable whereas it becomes
unstable for $\alpha>\alpha_{\text{BF}}$, where $\Omega$ is the common locked
frequency. On the other hand, desynchronized states, e.g., a splay state
$\\{\phi_{j}(t)=\frac{2\pi}{N}j+\Omega t\\}_{j=1}^{N}$, is stable for
$\alpha>\alpha_{\text{BF}}$ whereas it becomes unstable for
$\alpha<\alpha_{\text{BF}}$ [35]. All the oscillators have the same intrinsic
frequency $\omega\in\mathbb{R}$ that can be set to zero due to rotational
symmetry.

To investigate the collective behavior of the system, one can exploit the
Watanabe-Strogatz (WS) transformation [35, 37, 36], which is a linear
fractional transformation [58] that reads

$\displaystyle
e^{i\phi_{j}(t)}=M_{t}(e^{i\theta_{j}}):=\frac{b(t)+e^{i\theta_{j}}e^{i\varphi(t)}}{1+\overline{b(t)}e^{i\theta_{j}}e^{i\varphi(t)}}$
(3)

where $\varphi(t)\in\mathbb{R}$ and $b(t)\in\mathbb{C}$ are called the WS
variables characterizing the macroscopic dynamics of the system. The symbol
with an overbar represents the complex conjugate. In particular, the WS
variable $b(t)$ quantifies the degree of coherence of the system in a similar
manner (although not precisely the same) as the Kuramoto order parameter, and
fulfills $|b(t)|<1$ for all $t$ [38, 36]. Here, $\\{\theta_{j}\\}_{j=1}^{N}$
is a set of constants of motion obeying three constraints:
$\sum_{k=1}^{N}e^{i\theta_{k}}=0$ in order for $b(t)=0$ to be consistent with
the incoherent state $m(t)=0$, and the third one can be written as
$\sum_{k=1}^{N}\theta_{k}=0$. They are determined by a given initial condition
$\\{\phi_{j}(0)\\}_{j=1}^{N}$ (see Ref. [38]). The WS macroscopic variables
are governed by $\dot{b}=-\overline{g}b^{2}+g$ and
$\dot{\varphi}=i\overline{g}b-ig\overline{b}$ [36].

In the thermodynamic limit ($N\rightarrow\infty$) with uniform constants of
motion, one can approach the Ott-Antonsen (OA) manifold [32, 33] where the WS
variable $b(t)$ exactly coincides with the Kuramoto order parameter $m(t)$
[37, 38]. To see this, we first consider the uniform distribution of constants
of motion $d\mu(\theta)=\frac{1}{2\pi}d\theta$ together with the Watanabe-
Strogatz transformation $\phi=T(\theta):=-i\log M_{t}(e^{i\theta})$ for a
fixed $t$. This transformation pushes-forward a measure $\mu$ to $T_{*}\mu$
and gives $d(T_{*}\mu)(\phi)=f(\phi)d\phi$ where the phase density function

$\displaystyle f(\phi)$
$\displaystyle:=\frac{1}{2\pi}\partial_{\phi}T^{-1}(\phi)$
$\displaystyle=\frac{1}{2\pi}\frac{1-|b|^{2}}{1-2|b|\cos(\phi-\arg
b)+|b|^{2}}$ (4)

is the normalized Poisson kernel [36]. This fact reveals that the oscillators’
phases are distributed according to the normalized Poisson kernel in the OA
manifold [39]. The Poisson kernel has the Kuramoto order parameter given by

$\displaystyle m(t)$
$\displaystyle=\int_{\mathbb{T}}M_{t}(e^{i\theta})d\mu(\theta)=\int_{\mathbb{T}}e^{i\phi}d(T_{*}\mu)(\phi)$
$\displaystyle=\int_{\mathbb{T}}e^{i\phi}\frac{1}{2\pi}\frac{1-|b|^{2}}{1-2|b|\cos(\phi-\arg
b)+|b|^{2}}d\phi$ $\displaystyle=b(t).$ (5)

As a result, it becomes possible to examine the macroscopic behavior (i.e, the
dynamics of the Kuramoto order parameter) of the system by focusing solely on
the dynamics of the WS variable $b(t)$ within the OA manifold [38, 36, 39].

The system of Kuramoto-Sakaguchi oscillators (1) defined on a unit circle in
$\mathbb{C}^{1}$ or $\mathbb{R}^{2}$ has a sinusoidal form and therefore can
be extended to a system of generalized Kuramoto oscillators that are defined
on the surface of a unit sphere in $M$-dimensional real and complex spaces
[55]. To get a brief glimpse of this, we can write Eq. (1) as

$\displaystyle\frac{d}{dt}e^{i\phi_{j}}$
$\displaystyle=i\dot{\phi}_{j}e^{i\phi_{j}}=i\omega
e^{i\phi_{j}}+2ie^{i\phi_{j}}\text{Im}\bigg{[}g(t)e^{-i\phi_{j}}\bigg{]}$
$\displaystyle=i\omega
e^{i\phi_{j}}+e^{i\phi_{j}}\bigg{(}g(t)e^{-i\phi_{j}}-\overline{g(t)}e^{i\phi_{j}}\bigg{)}$
$\displaystyle=-e^{i\phi_{j}}\overline{g(t)}e^{i\phi_{j}}+i\omega
e^{i\phi_{j}}+g(t)$ (6)

for $j=1,...,N$. Let us take
$e^{i\phi_{j}(t)}\in\mathbb{C}^{1}\mapsto\bm{x}_{j}(t)\in\mathbb{K}^{M}$ with
a constraint
$\bm{x}_{j}(t)^{{\dagger}}\bm{x}_{j}(t)=\left\langle\bm{x}_{j}(t)\middle|\bm{x}_{j}(t)\right\rangle=1$
for all $t$ and $j=1,...,N$ where $\mathbb{K}=\mathbb{R}$ or $\mathbb{C}$ is a
ground field, and ${\dagger}$ denotes a Hermitian adjoint [59]. Then, one can
treat $\bm{x}_{j}(t)$ as an oscillator defined on the surface of the unit
sphere denoted by
$\mathbb{S}^{M}:=\\{\bm{x}\in\mathbb{K}^{M}|\bm{x}^{\dagger}\bm{x}=1\\}$
either for $\mathbb{K}=\mathbb{R}$ or $\mathbb{C}$. In this consideration, the
oscillators can be assumed to follow

$\displaystyle\dot{\bm{x}}_{j}=-\bm{x}_{j}\bm{g}^{{\dagger}}\bm{x}_{j}+\Omega\bm{x}_{j}+\bm{g}$
(7)

for $j=1,...,N$. Here, the mean-field forcing $\bm{g}(t)\in\mathbb{K}^{M}$ can
be any arbitrary vector that characterizes an external forcing function or a
global coupling of the system. In this paper, we consider the mean-field
forcing $\bm{g}(t)$ which depends only on the Kuramoto order parameter (see
Sec. 3 for details). The Kuramoto order parameter serves as the center of mass
of the oscillators on $\mathbb{S}^{M}$, which is defined by

$\displaystyle\bm{m}(t):=\frac{1}{N}\sum_{k=1}^{N}\bm{x}_{k}(t)\in\mathbb{K}^{M}.$
(8)

Furthermore, $\Omega\in\mathbb{K}^{M\times M}$ with $\Omega^{\dagger}=-\Omega$
is an anti-hermitian natural frequency matrix that corresponds to $i\omega$ in
the above equation. Throughout this paper, we will set $\bm{\Omega}$ to zero
since we only consider identical oscillators. This model and other variants of
vector-formed oscillators, i.e., generalized Kuramoto models have been so far
investigated for a variety of collective dynamics [49, 50, 51, 42, 43, 44, 45,
46]. In particular, the macroscopic dynamics of the generalized Kuramoto
oscillator system can be investigated based on the generalized WS
transformation introduced in Ref. [55].

Below, we briefly revisit the main focus of Ref. [55] that is necessary for
our study. Therein, the higher-dimensional WS transformation was introduced as
a vector form of a linear fractional transformation

$\displaystyle\bm{x}_{j}(t)=M_{t}(\bm{x}_{0,j}):=\frac{\bm{A}\bm{x}_{0,j}+\bm{b}}{\bm{b}^{{\dagger}}\bm{A}\bm{x}_{0,j}+1}$
(9)

where $\bm{A}(t)\in\mathbb{K}^{M\times M}$ and $\bm{b}(t)\in\mathbb{K}^{M}$
are the WS variables that determine the macroscopic dynamics of the system. In
this context, the initial conditions $\bm{x}_{j}(0)=:\bm{x}_{0,j}$ act as the
constants of motion. To describe the incoherent state consistently, i.e.,
$\bm{m}(t)=\bm{b}(t)=0$, the initial conditions should satisfy
$\frac{1}{N}\sum_{k=1}^{N}\bm{x}_{0,k}=0$. This holds, e.g., if they are
uniformly distributed on $\mathbb{S}^{M}$. Under the constraint
$\left\langle\bm{x}_{j}\middle|\bm{x}_{j}\right\rangle=1$ and with uniform
initial conditions, the WS variables are governed by

$\displaystyle\dot{\bm{A}}$
$\displaystyle=\bm{g}\bm{b}^{\dagger}\bm{A}-\bm{g}^{\dagger}\bm{b}\bm{A}$
$\displaystyle\dot{\bm{b}}$
$\displaystyle=-\bm{b}\bm{g}^{\dagger}\bm{b}+\bm{g}.$ (10)

Furthermore, one can show that $\bm{A}=\bm{H}^{1/2}\bm{U}$ (polar
decomposition) where $\bm{U}\in\mathbb{K}^{M\times M}$ is a unitary matrix and
$\bm{H}^{1/2}=\bm{V}\Sigma^{1/2}\bm{V}^{\dagger}$ with
$\Sigma^{1/2}:=(\sqrt{\Sigma_{ij}})$. Here, $\bm{V}\Sigma\bm{V}^{\dagger}$ is
a singular value decomposition of
$\bm{H}:=(1-|\bm{b}|^{2})I_{M}+\bm{b}\bm{b}^{\dagger}$ where
$I_{M}\in\mathbb{R}^{M\times M}$ is the identity matrix. Thus it follows that
$\bm{H}^{1/2}\bm{b}=\bm{b}$.

In the thermodynamic limit, without loss of generality, it is possible to set
$\bm{U}=I_{M}$ as long as $\bm{x}_{0,j}$ are uniformly distributed on
$\mathbb{S}^{M}$. Then, as T. Tanaka showed [55], the Kuramoto order parameter
can also be expressed in terms of the WS (or OA) variables as

$\displaystyle\bm{m}(t)$
$\displaystyle=\int_{|\bm{x}_{0}|=1}M_{t}(\bm{x}_{0})d\mu(\bm{x}_{0})=\frac{1}{S_{M}}\int_{|\bm{x}_{0}|=1}\bm{x}(t)d\bm{x}_{0}$
$\displaystyle=\frac{M-1}{M}\,_{2}F_{1}(\frac{1}{2},1;\frac{M+2}{2};|\bm{b}|^{2})\bm{b}(t)$
$\displaystyle=h(|\bm{b}|^{2},M)\bm{b}(t)$ (11)

for $\mathbb{K}=\mathbb{R}$. Here, $\,{}_{2}F_{1}$ is the ordinary
hypergeometric function [60], $S_{M}:=\frac{2\pi^{M/2}}{\Gamma(M/2)}$ is the
surface area of the $(M-1)$-sphere, and
$h(|\bm{b}|^{2},M):=\frac{M-1}{M}\,_{2}F_{1}(\frac{1}{2},1;\frac{M+2}{2};|\bm{b}|^{2})$.
For $\mathbb{K}=\mathbb{C}$, the Kuramoto order parameter is given by

$\displaystyle\bm{m}(t)$
$\displaystyle=\frac{1}{S_{2M}}\int_{|\bm{x}_{0}|=1}\bm{x}(t)d\bm{x}_{0}=\bm{b}(t)$
(12)

which means the Kuramoto order parameter exactly coincides with the WS
variable $\bm{b}(t)$. For the details, see Ref. [55] and A.

##  3 Coupling Matrix: Benjamin-Feir Instability Point

To explore chimera states in a system of generalized Kuramoto-Sakaguchi
oscillators in two-population networks, we introduce a suitable coupling
matrix in this section, which corresponds to the phase-lag parameter $\alpha$
in Eq. (1). Moreover, the Benjamin-Feir (BF) instability point is obtained for
this coupling matrix. In the remainder of this paper, we consider the system
below the BF point, in order to follow the previous observations of chimeras
in two-population networks [10, 12, 17, 18, 16].

In line with Eq. (1), we define the forcing function as
$\bm{g}(t)=\bm{K}\bm{m}(t)$ where $\bm{K}\in\mathbb{K}^{M\times M}$ is a
coupling matrix and $\bm{m}(t)$ is the Kuramoto order parameter in Eq. (8).
Then, the microscopic dynamics of the oscillators are given by [55]

$\displaystyle\dot{\bm{x}}_{j}$
$\displaystyle=-\bm{x}_{j}\bm{g}^{\dagger}\bm{x}_{j}+\bm{g}$
$\displaystyle=-\bm{x}_{j}\bm{m}^{\dagger}\bm{K}^{\dagger}\bm{x}_{j}+\bm{K}\bm{m}=\frac{1}{N}\sum_{k=1}^{N}\bigg{(}\bm{K}\bm{x}_{k}-(\bm{x}_{k}^{\dagger}\bm{K}^{\dagger}\bm{x}_{j})\bm{x}_{j}\bigg{)}$
(13)

for $j=1,...,N$. In Eq. (1), the phase-lag parameter induces phase rotations
of $-\alpha$ to each oscillator’s phase on the unit circle. Similarly, we
introduce a rotation matrix as the coupling matrix [46, 54]. First, we
consider the real space $\mathbb{K}=\mathbb{R}$. For a rotation in the real
space, we need to distinguish between even and odd dimensional cases [50, 44].
For odd $M$, we set the $M$-axis (see A) as the rotational axis; then
$M^{\perp}$-plane is isoclinically rotating with the same angle. An example is

$\displaystyle\bm{K}=\begin{pmatrix}\cos\alpha&-\sin\alpha&0&0&0\\\
\sin\alpha&\cos\alpha&0&0&0\\\ 0&0&\cos\alpha&-\sin\alpha&0\\\
0&0&\sin\alpha&\cos\alpha&0\\\ 0&0&0&0&1\end{pmatrix}$ (14)

for $\mathbb{R}^{5}$. For even $M$, there is no rotational axis and hence we
exploit planes of rotation, in particular, an isoclinic rotation with the same
rotational angle for each plane. An example is

$\displaystyle\bm{K}=\begin{pmatrix}\cos\alpha&-\sin\alpha&0&0\\\
\sin\alpha&\cos\alpha&0&0\\\ 0&0&\cos\alpha&-\sin\alpha\\\
0&0&\sin\alpha&\cos\alpha\end{pmatrix}$ (15)

for $\mathbb{R}^{4}$.

To determine the Benjamin-Feir point, we consider the macroscopic WS dynamics
in the thermodynamic limit Eq. (10), but use the notation $\bm{\psi}$ instead
of $\bm{b}$ to denote the WS variable for future reference. The WS dynamics is
then written as

$\displaystyle\dot{\bm{\psi}}$
$\displaystyle=-\bm{\psi}\bm{g}^{\dagger}\bm{\psi}+\bm{g}$
$\displaystyle=h(|\bm{\psi}|^{2},M)\bigg{(}-\bm{\psi}\bm{\psi}^{\dagger}\bm{K}^{\dagger}\bm{\psi}+\bm{K}\bm{\psi}\bigg{)}$
(16)

where the order parameter $\bm{m}=h(|\bm{\psi}|^{2},M)\bm{\psi}$ as in Eq.
(11). Using Eq. (16), the dynamics of the magnitude of the WS variable is
given by

$\displaystyle\partial_{t}|\bm{\psi}|^{2}$
$\displaystyle=\partial_{t}(\bm{\psi}^{\dagger}\bm{\psi})$
$\displaystyle=h(|\bm{\psi}|^{2},M)\bigg{(}(-\bm{\psi}^{\dagger}\bm{K}\bm{\psi}\bm{\psi}^{\dagger}+\bm{\psi}^{\dagger}\bm{K}^{\dagger})\bm{\psi}+\bm{\psi}^{\dagger}(-\bm{\psi}\bm{\psi}^{\dagger}\bm{K}^{\dagger}\bm{\psi}+\bm{K}\bm{\psi})\bigg{)}$
$\displaystyle=2h(|\bm{\psi}|^{2},M)(1-|\bm{\psi}|^{2})\left\langle\bm{\psi}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}\right\rangle$
(17)

where $\partial_{t}:=\frac{d}{dt}$. For even $M$, the last term becomes
$\left\langle\bm{\psi}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}\right\rangle=\cos\alpha|\bm{\psi}|^{2}$.
Thus, Equation (17) is simply written as

$\displaystyle\partial_{t}|\bm{\psi}|^{2}=2h(|\bm{\psi}|^{2},M)(1-|\bm{\psi}|^{2})|\bm{\psi}|^{2}\cos\alpha.$
(18)

There exist two fixed points: $\rho_{*}=|\bm{\psi}|=0$ representing the
completely incoherent state, and $\rho^{*}=|\bm{\psi}|=1$, indicating the
synchronized state. The stability of these fixed points depends on the value
of $\alpha$. Specifically, the synchronized state $\rho^{*}$ is stable while
the incoherent state $\rho_{*}$ is unstable for $\alpha<\frac{\pi}{2}$.
Opposite to this, the incoherent state $\rho_{*}$ is stable whereas the
synchronized state $\rho^{*}$ becomes unstable for $\alpha>\frac{\pi}{2}$.
Therefore, the Benjamin-Feir instability occurs at
$\alpha_{\text{BF}}=\frac{\pi}{2}$ for even $M$.

On the other hand, for odd $M$, the last term reads as
$\left\langle\bm{\psi}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}\right\rangle=|\bm{\psi}|^{2}\cos\alpha+(1-\cos\alpha)x_{M}^{2}$
where $x_{M}$ is the coordinate of $\bm{\psi}$ along the rotational axis,
i.e., the $M$-axis. In this case, the synchronized state exists in two
distinct spaces: either on the $M$-axis or on the $M^{\perp}$-plane. In the
former case, the synchronized state is always stable, regardless of the value
of $\alpha$ since it is governed by
$\partial_{t}|\bm{\psi}|^{2}=2h(|\bm{\psi}|^{2},M)(1-|\bm{\psi}|^{2})|\bm{\psi}|^{2}$
on the $M$-axis. Note $\bm{\psi}=(0,...,0,x_{M})^{\top}$ and
$|\bm{\psi}|=x_{M}$ along the $M$-axis. For the latter case, the stability of
the synchronized state depends on the parameter $\alpha$: Equation (17) on the
$M^{\perp}$-plane can be written as
$\partial_{t}|\bm{\psi}|^{2}=2h(|\bm{\psi}|^{2},M)(1-|\bm{\psi}|^{2})|\bm{\psi}|^{2}\cos\alpha$
(here, $x_{M}=0$ on the $M^{\perp}$-plane), which reveals the synchronized
state is stable for $\alpha<\frac{\pi}{2}$ and unstable for
$\alpha>\frac{\pi}{2}$. Consequently, the Benjamin-Feir instability also
occurs at $\alpha_{\text{BF}}=\frac{\pi}{2}$. Note that the synchronized state
on the $M^{\perp}$-plane is always unstable with respect to a perturbation
along the $M$-axis. See B for the 3D real space. A similar result was reported
in Refs. [46, 54].

Moving forward, let us explore the higher-dimensional complex space
$\mathbb{C}^{M}$. In this setting, we choose the coupling matrix as
$\bm{K}=e^{-i\alpha}I_{M}$ for all $M$. Then, the governing equation reads

$\displaystyle\dot{\bm{\psi}}$
$\displaystyle=-\bm{\psi}\bm{g}^{\dagger}\bm{\psi}+\bm{g}$
$\displaystyle=-\bm{\psi}\bm{\psi}^{\dagger}\bm{K}^{\dagger}\bm{\psi}+\bm{K}\bm{\psi}$
(19)

since the Kuramoto order parameter exactly coincides with the WS variable,
i.e., $\bm{m}(t)=\bm{\psi}(t)$ as in Eq. (12). Consequently, the magnitude of
the WS variable is determined by

$\displaystyle\partial_{t}|\bm{\psi}|^{2}$
$\displaystyle=2(1-|\bm{\psi}|^{2})\left\langle\bm{\psi}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}\right\rangle$
$\displaystyle=2(1-|\bm{\psi}|^{2})|\bm{\psi}|^{2}\cos\alpha$ (20)

Therefore, the synchronized state $|\bm{\psi}|=1$ is stable for
$\alpha<\frac{\pi}{2}$ whereas it is unstable for $\alpha>\frac{\pi}{2}$, and
thus $\alpha_{\text{BF}}=\frac{\pi}{2}$.

##  4 Chimera Dynamics in Two-population Networks for $\mathbb{C}^{2}$

Now, we are prepared to investigate chimera states in two-population networks
of identical generalized Kuramoto-Sakaguchi oscillators. The microscopic
dynamics of the system is given by

$\displaystyle\dot{\bm{x}}^{(1)}_{j}$
$\displaystyle=-\bm{x}^{(1)}_{j}\bm{g}_{1}^{\dagger}\bm{x}^{(1)}_{j}+\bm{g}_{1}$
$\displaystyle\dot{\bm{x}}^{(2)}_{j}$
$\displaystyle=-\bm{x}^{(2)}_{j}\bm{g}_{2}^{\dagger}\bm{x}^{(2)}_{j}+\bm{g}_{2}$
(21)

where the mean-field forcing for each population reads

$\displaystyle\bm{g}_{a}$
$\displaystyle=\mu\bm{K}\bm{m}^{(a)}+\nu\bm{K}\bm{m}^{(b)}=\bm{K}\bigg{(}\frac{\mu}{N}\sum_{k=1}^{N}\bm{x}^{(a)}_{k}(t)+\frac{\nu}{N}\sum_{k=1}^{N}\bm{x}^{(b)}_{k}(t)\bigg{)}$
(22)

and the Kuramoto order parameter of each population is defined by

$\bm{m}^{(a)}(t)=\frac{1}{N}\sum_{k=1}^{N}\bm{x}^{(a)}_{k}(t)$ (23)

for $(a,b)=(1,2)$ or $(2,1)$. Throughout this work, the inter-population
coupling strength is fixed as $\mu=1$ and the intra-population coupling
strength $\nu=1-A$ where $A$ is the control parameter. The rotation matrices
for $\mathbb{K}=\mathbb{R}$ in Sec. 3 is employed as the coupling matrix while
it is $\bm{K}=e^{-i\alpha}I_{M}$ for $\mathbb{K}=\mathbb{C}$. Note that in
this paper, we fix $\alpha=\frac{\pi}{2}-0.005<\alpha_{\text{BF}}$.

###  4.1 Stationary and Breathing Chimeras

In the following, we study the chimera states with a particular focus on the
2D complex space $\mathbb{C}^{2}$. Recall that for the complex space the
Kuramoto order parameter $\bm{m}^{(a)}(t)$ is exactly equal to
$\bm{\psi}_{a}(t)$ in the thermodynamic limit. Then, the forcing field is
written as $\bm{g}_{a}=\mu\bm{K}\bm{\psi}_{a}+\nu\bm{K}\bm{\psi}_{b}$ where
$\bm{K}=e^{-i\alpha}I_{M}$. Therefore, the WS variable $\bm{\psi}_{a}$ is
governed by

$\displaystyle\dot{\bm{\psi}}_{a}$
$\displaystyle=-\bm{\psi}_{a}\bm{g}_{a}^{\dagger}\bm{\psi}_{a}+\bm{g}_{a}$
$\displaystyle=-\bm{\psi}_{a}\bigg{(}\mu\bm{\psi}_{a}^{\dagger}\bm{K}^{\dagger}+\nu\bm{\psi}_{b}^{\dagger}\bm{K}^{\dagger}\bigg{)}\bm{\psi}_{a}+\mu\bm{K}\bm{\psi}_{a}+\nu\bm{K}\bm{\psi}_{b}$
(24)

for $(a,b)=(1,2)$ and $(2,1)$ in the thermodynamic limit.

As an initial approach, we examine the dynamics of the magnitude of the order
parameter vectors:
$|\bm{\psi}_{a}|=\sqrt{\left\langle\bm{\psi}_{a}\middle|\bm{\psi}_{a}\right\rangle}$
for $a=1,2$. From Eq. (24), the dynamics of the order parameter magnitude
reads

$\displaystyle\partial_{t}|\bm{\psi}_{1}|^{2}$
$\displaystyle=2(1-|\bm{\psi}_{1}|^{2})\bigg{(}\mu\cos\alpha|\bm{\psi}_{1}|^{2}+\nu\text{Re}\bigg{[}e^{-i\alpha}\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle\bigg{]}\bigg{)}$
$\displaystyle\partial_{t}|\bm{\psi}_{2}|^{2}$
$\displaystyle=2(1-|\bm{\psi}_{2}|^{2})\bigg{(}\mu\cos\alpha|\bm{\psi}_{2}|^{2}+\nu\text{Re}\bigg{[}e^{-i\alpha}\left\langle\bm{\psi}_{2}\middle|\bm{\psi}_{1}\right\rangle\bigg{]}\bigg{)}$

and the cross term is governed by

$\displaystyle\partial_{t}\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=\mu\big{(}(1-|\bm{\psi}_{1}|^{2})e^{-i\alpha}+(1-|\bm{\psi}_{2}|^{2})e^{i\alpha}\big{)}\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle~{}~{}~{}~{}+\nu
e^{-i\alpha}(|\bm{\psi}_{1}|^{2}-\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle^{2})+\nu
e^{i\alpha}(|\bm{\psi}_{2}|^{2}-\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle^{2}).$
(26)

From numerical integration [61] of Eq. (24), we find that the cross term can
be represented by
$\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle=|\bm{\psi}_{1}||\bm{\psi}_{2}|e^{i\Theta}$
where $\Theta\in\mathbb{R}$ as explained later in this section. Then, setting
$|\bm{\psi}_{a}|=:\rho_{a}$ for $a=1,2$ and $\rho_{2}=1$, we obtain

$\displaystyle\partial_{t}\rho_{1}$
$\displaystyle=\frac{1-\rho_{1}^{2}}{\rho_{1}}\big{(}\mu\cos\alpha\rho_{1}^{2}+\nu\rho_{1}\cos(\Theta-\alpha)\big{)}$
$\displaystyle\partial_{t}\Theta$
$\displaystyle=\mu\sin\alpha(\rho_{1}^{2}-1)+\nu\Big{(}-2\rho_{1}\sin(\Theta+\alpha)+\sin(\Theta-\alpha)\big{(}-\frac{1}{\rho_{1}}-\rho_{1}\big{)}\Big{)}.$
(27)

Note that the above equations are equivalent to Eq. (10) in Ref. [10], i.e.,
the Ott-Antonsen equation of a system of identical KS oscillators (defined on
the unit circle of $\mathbb{C}^{1}$) in two-population networks. From Eq.
(27), we can obtain an overview of some of the observable chimera states.

Figure 1: (a) A bifurcation diagram of stable (solid curve) and unstable
(dashed curve) chimera states from Eq. (27): stationary chimeras (red) and
breathing chimera states (green). LP: limiting point bifurcation, HB:
supercritical Hopf bifurcation, and Hom: homoclinic bifurcation. Lower panels:
Time evolution of the magnitude of the order parameter vectors (the
synchronized population $|\bm{\psi}_{2}|$ in orange and the incoherent
population $|\bm{\psi}_{1}|$ in light blue) in the left column and time
evolution of the modulus of the components of the order parameter vector for
the incoherent population $\bm{\psi}_{1}=(z_{1},z_{2})^{\top}$ (the first
component $|z_{1}|$ in red and the second component $|z_{2}|$ in blue) in the
right column. (b) Stationary chimera states for $A=0.35$. (c) Breathing
chimera states for $A=0.46$. The presented results are based on data obtained
after disregarding the initial transient behavior for $t>10^{5}$.

In Fig. 1 (a), a bifurcation diagram of chimera states is depicted as the
parameter $A$ varies. Stable (red solid) and unstable (red dashed) stationary
chimera states are born/annihilated at a limiting point (LP) bifurcation for a
rather strong intra-population coupling strength. The stable stationary
chimera state undergoes a supercritical Hopf bifurcation (HB) at
$A=A_{\text{HB}}$, and a breathing chimera state (green solid) emerges as a
stable limit-cycle solution. Integrating Eq. (24) numerically, we find
stationary and breathing chimera states for the parameter values according to
the bifurcation diagram depicted in Fig. 1 (a). For simplicity, we assume that
the first population is incoherent and the second population is synchronized,
i.e., $|\bm{\psi}_{1}(t)|<1$ and $|\bm{\psi}_{2}(t)|=1$. Also, we use the
following notations to denote components of each order parameter vector:
$\bm{\psi}_{1}(t)=(z_{1}(t),z_{2}(t))^{\top}$ and
$\bm{\psi}_{2}(t)=(w_{1}(t),w_{2}(t))^{\top}$ where
$z_{i}(t),w_{i}(t)\in\mathbb{C}$ and $\theta_{i}(t)=\arg
z_{i}(t)\in\mathbb{R}$ and $\phi_{i}(t)=\arg w_{i}(t)\in\mathbb{R}$ for
$i=1,2$.

In Fig. 1 (b), the time evolution of the magnitude of the order parameter
vectors for a stationary chimera state is depicted together with the time
evolution of $|z_{1}(t)|$ and $|z_{2}(t)|$. Both $|\bm{\psi}_{a}|$ and the
components $|z_{i}(t)|$ appear as a fixed point solution. Note that for the
synchronized population, $|w_{i}(t)|$ also shows stationary motion.
Furthermore, Equation (24) is invariant under a unitary transformation defined
by

$\displaystyle\bm{u}=\begin{pmatrix}e^{-i\theta_{1}}&0\\\
0&e^{-i\theta_{2}}\end{pmatrix}$ (28)

which means $F_{a}(\bm{u}\bm{\psi}_{a})=\bm{u}F_{a}(\bm{\psi}_{a})$ where
$F_{a}(\bm{\psi}_{a}):=-\bm{\psi}_{a}\bm{g}_{a}^{\dagger}\bm{\psi}_{a}+\bm{g}_{a}$
for $a=1,2$. This unitary transformation is a continuous symmetry
corresponding to the phase shift invariance of the KS model on the unit
circle. Therefore, the phase difference
$\Theta_{i}(t):=\phi_{i}(t)-\theta_{i}(t)$ for $i=1,2$ also shows a stationary
motion as a function of time (not shown here). For larger $A$, Figure 1 (c)
displays the time evolution of $|\bm{\psi}_{a}|$ and $|z_{i}|$ for a given
$A\in(A_{\text{HB}},A_{\text{Hom}})$. Each variable exhibits a periodic motion
as a function of time, and the stable breathing chimera appears indeed as a
limit-cycle solution after a supercritical Hopf bifurcation. Further
increasing $A$, the breathing chimera state disappears via a homoclinic
bifurcation (Hom) as the period of it approaches infinity.

###  4.2 Alternating and Aperiodic Chimeras

Thus far, the scenario of the emergence of chimera states for $\mathbb{C}^{2}$
bears a strong resemblance to that observed in the system of identical KS
oscillators on the unit circle of $\mathbb{C}^{1}$ (e.g., see Refs. [10, 12]).
Nevertheless, in the two-dimensional space $\mathbb{C}^{2}$, the order
parameter consists of two complex components, providing additional
complexities beyond this. Let us start with the parameter regime after the
homoclinic bifurcation. In Fig. 2 (a), a chimera state is shown for a given
$A>A_{\text{Hom}}$. This chimera state can be characterized by a periodic
alternation between the two order parameter vectors. The magnitude of the
order parameter vectors (Fig. 2 (a-1)) seemingly satisfy
$|\bm{\psi}_{1}(t)|=|\bm{\psi}_{2}(t-\frac{T}{2})|$ within the limits of our
numerical capabilities where $T$ denotes their period. Likewise, each
component of the order parameter vectors (Fig. 2 (a-2)) displays
$|z_{i}(t)|=|w_{i}(t-\frac{T}{2})|$ for $i=1,2$. Such an alternating chimera
motion was reported previously for a system of heterogeneous Kuramoto-
Sakaguchi oscillators in two-population networks after the homoclinic
bifurcation in which the breathing chimera state disappeared [15, 14]. A
similar motion was also reported in a system of identical KS oscillators in
three-population networks, which as well was born near a homoclinic
bifurcation [26]. Together the three observations suggest that the periodic
alternating chimera dynamics is somehow linked to the breathing chimera states
through the homoclinic bifurcation (see Fig. 4 (c)).

Figure 2:  Left column: Time evolution of the magnitude of the order
parameters: $|\bm{\psi}_{1}(t)|$ (blue light) and $|\bm{\psi}_{2}(t)|$
(orange). Right column: Time evolution of the first component of the order
parameter vectors: $|z_{1}(t)|$ (red) and $|w_{1}(t)|$ (blue). (a) Periodic
alternating chimera at $A=0.75$. (b) Aperiodic alternating chimera at
$A=0.86$. The presented results are based on data obtained after disregarding
the initial transient behavior for $t>10^{5}$.

Following $A$ further, we observe a seemingly similar motion to the
alternating chimera dynamics in terms of the magnitude of the order
parameters. In Fig. 2 (b), an example trajectory is shown. The magnitude of
the order parameters (Fig. 2 (b-1)) seems similar to that of the periodic
alternating chimera dynamics except for $\text{min}_{t>0}|\bm{\psi}_{a}(t)|$.
However, the componentwise dynamics shows entirely different characteristics.
In Fig. 2 (b-2), we plot an aperiodic time evolution of the first components
of the order parameter for both populations. To investigate the aperiodic
chimera dynamics in more detail, in Fig. 3 (a), the Poincaré map of the
trajectory is plotted in the section defined by $\text{Re}[z_{2}]\equiv 0$.
The aperiodic dynamics depicted in Fig. 2 (b-2) displays scattered points on
the Poincaré section, as anticipated for the aperiodic motion. This conjecture
is further supported by the Lyapunov exponents along the reference trajectory
in the phase space, following the chimera dynamics in Fig. 3 (a) [62, 63, 64].
In Fig. 3 (b) it can be seen that there are two positive Lyapunov exponents
that indicate a sensitive dependence of the trajectory on initial conditions.
Also, there are two zero Lyapunov exponents arising from the two continuous
symmetries: time shift (an autonomous equation) and the unitary transformation
in Eq. (28) corresponding to a phase shift invariance. All the other Lyapunov
exponents are negative.

Figure 3:  (a) Poincaré map of the first components of the order parameters:
$|z_{1}|$ vs. $|w_{1}|$ measured on the section defined by
$\text{Re}[z_{2}]\equiv 0$. (b) Lyapunov exponents measured along the
reference trajectory in (a): Positive exponents (red), zero exponents (blue)
and negative exponents (orange). (c-d) Time evolution of the broken conserved
quantities of the trajectory in (a) (see main text). The parameter is $A=0.86$
and the transient behavior was discarded ($t>10^{5}$).

Now, we elucidate how this aperiodic chimera state emerges from the periodic
alternating chimera state as $A$ is increased. In Fig. 4 (a), the bifurcation
diagram is replotted with parameter values indicating where the periodic
alternating ($A_{\text{Hom}}$) and aperiodic ($A_{c}$) chimera states emerge,
respectively. For $A<A_{c}$, we numerically obtain two conserved quantities
along the chimera trajectory in phase space. As explained before, the system
is invariant under the unitary transformation in Eq. (28). Using this
transformation, we can define the phase difference of each component between
the two order parameter vectors: $\Theta_{i}(t):=\phi_{i}(t)-\theta_{i}(t)$
for $i=1,2$. The first conserved quantity is given by

$\displaystyle
C_{1}(\Theta_{1},\Theta_{2};t):=\sin(\Theta_{1}(t)-\Theta_{2}(t)).$ (29)

For this quantity, $C_{1}(t)=0$ holds for all $t$ along the chimera trajectory
for $A<A_{c}$, i.e., the stationary, breathing and periodic alternating
chimera states all preserve $C_{1}(t)=0$ along each trajectory. Furthermore,
the numerical integration of Eq. (24) for $A<A_{c}$ confirms the relation
between the cross term and the magnitude of the order parameter vectors:
$\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle=|\bm{\psi}_{1}||\bm{\psi}_{2}|e^{i\Theta}$
where $\Theta:=\Theta_{1}=\Theta_{2}$ for $C_{1}=0$. This relation leads to
the second conserved quantity

$\displaystyle
C_{2}(z_{1},z_{2},w_{1},w_{2};t):=\big{(}|z_{1}(t)||w_{2}(t)|-|z_{2}(t)||w_{1}(t)|\big{)}^{2}$
(30)

whereby $C_{2}(t)=0$ for all $t$ along the chimera trajectory as long as
$A<A_{c}$.

Figure 4:  (a) A re-plot of the bifurcation diagram in Fig. 1 (a) with the
periodic alternating chimeras and the chaotic chimeras. (b) Conserved
quantities as a function of the parameter $A$: $t_{0}=20,000$ and
$t_{1}=10^{5}$. (c) Maximum (blue) and minimum (red) values of the magnitude
of the order parameter vectors as a function of the parameter $A$. (d) Time
evolution of $\cos\Theta_{2}(t)$ before (upper, $A=0.52$) and after (lower,
$A=0.53$) the global bifurcation. $A_{c}$ denotes a parameter point from which
on the chaotic chimera emerges.

However, we find that for $A>A_{c}$, $C_{1}$ and $C_{2}$ are not conserved any
more along the chimera trajectory. In Fig. 3 (c-d), the time evolution of
$C_{1}(t)$ and $C_{2}(t)$ is depicted for the chimera trajectory in Fig. 2
(b). Both quantities display irregular deviation from zero as a function of
time. Once again, note that for $A<A_{c}$, $C_{1}(t)=0$ and $C_{2}(t)=0$ for
all $t$ along the chimera orbit. To see this clearly, we numerically measure
$\sup_{t\in[t_{0},t_{1}]}|C_{1}(t)|$ and $\sup_{t\in[t_{0},t_{1}]}|C_{2}(t)|$
with respect to the parameter $A$ where $t_{0}=2\times 10^{4}$ and
$t_{1}=10^{5}$. Figure 4 (b) confirms that for $A<A_{c}$ where the stationary,
breathing and periodic alternating chimeras are observed, the two quantities
remain zero. However, from $A_{c}$ on, the conserved quantities are broken and
begin to show irregular time-evolution along the given chimera trajectory.

Another notable observation is that the breathing chimera state undergoes a
global bifurcation involving the completely incoherent state, i.e.,
$|\bm{\psi}_{1}|=0$ [14]. In Fig. 4 (c), the maximum and minimum values of the
magnitude of the order parameter vector are depicted for the incoherent
population. For small $A$, $\text{min}(|\bm{\psi}_{\text{incoh}}|)$ is
decreasing continuously as $A$ increases up to $A=A_{g}$. Then, it touches
zero value where the completely incoherent state is located. From that point
on, the minimum value of the order parameter magnitude exhibits a continuous
increase, eventually connecting seamlessly to that of the periodic alternating
chimera state for $A>A_{\text{Hom}}$. Further increasing $A$, we reach $A_{c}$
where the aperiodic chimera emerges and
$\text{min}(|\bm{\psi}_{\text{incoh}}|)$ undergoes a discontinuous change to
higher values. At the global bifurcation $A_{g}$, the motion of the breathing
chimera changes. Figure. 4 (d) shows the time evolution of $\cos\Theta(t)$ for
the breathing chimera state before (d-1) and after (d-2) the global
bifurcation. For $A<A_{g}$, the phase variable $\Theta(t)$ of the breathing
chimera states evolves within some interval smaller than
$\mathbb{T}:=[0,2\pi]$. After touching the incoherent state for $A>A_{g}$, the
phase $\Theta(t)$ of the breathing chimera states monotonically increases as a
function of time thereby sweeping all $\mathbb{T}$. Thus, at the global
bifurcation the motion of the order parameter changes from libration to
rotation.

###  4.3 Chimera States in the Microscopic Dynamics

In this section, we investigate the microscopic dynamics of the system of
generalized Kuramoto-Sakaguchi oscillators in two-population networks. Here,
we directly perform numerical integrations of Eq. (21) for $\mathbb{C}^{2}$.

Figure 5: (a-d) Time evolution of the magnitude of the order parameter vectors
defined in Eq. (23) from solving Eq. (21) with $N=30$: the first population
$|\bm{m}^{(1)}(t)|$ (light blue) and the second population $|\bm{m}^{(2)}(t)|$
(orange). (a) The stationary chimera state with $A=0.35$. (b) The breathing
chimera state with $A=0.55$. (c) The periodic alternating chimera state with
$A=0.75$. (d) The chimera state of the broken conserved quantities with
$A=0.86$. Insets of (c-d): Time evolution of $\bm{m}^{(1)}_{1}(t)$ (blue) and
$\bm{m}^{(2)}_{1}(t)$ (red). (e) Snapshot of the oscillators with $N=100$.
Blue dot: synchronized oscillators $\bm{x}^{(2)}_{j}$. Red dots: incoherent
oscillators $\bm{x}^{(1)}_{j}$. Blue curve: $\bm{m}^{(2)}(t)$. Red curve:
$\bm{m}^{(1)}(t)$. (f) The histogram of the distribution of the phases
$\big{\\{}\arg\big{(}\frac{1}{|\bm{\psi}_{1}|}\bm{\psi}_{1}^{\dagger}\bm{x}_{j}^{(1)}\big{)}\big{\\}}_{j=1}^{N}$
for $N=500$ and $A=0.35$. Red curve indicates the normalized Poisson kernel in
Eq. (35). The presented results are based on data obtained after disregarding
the initial transient behavior for $t>10^{4}$.

In Fig. 5 (a-d), the time evolution is depicted for the magnitude of the
Kuramoto order parameter $|\bm{m}^{(a)}(t)|$ defined in Eq. (23) for 30
oscillators in each population $a=1,2$. All the results are obtained from
random initial conditions of $\bm{x}^{(a)}_{j}$ satisfying
$\left\langle\bm{x}^{(a)}_{j}(0)\middle|\bm{x}^{(a)}_{j}(0)\right\rangle=1$
for $j=1,...,N$ and $a=1,2$. For a given $A\in(A_{\text{LP}},A_{\text{HB}})$,
the microscopic dynamics shows a motion of the order parameter vector
corresponding to the stationary chimera state (Fig. 5 (a)). Due to the finite-
size effect and random initial conditions, we observe a fluctuation in the
magnitude of the Kuramoto order parameter. Likewise, in Fig. 5 (b), for
$A\in(A_{\text{HB}},A_{\text{Hom}})$, the order parameter obtained from the
breathing chimera state exhibits small fluctuation superimposed to a simple
periodic variation. Furthermore, for $A\in(A_{\text{Hom}},A_{c})$, the
microscopic dynamics as well forms alternating chimeras, manifesting
themselves in alternating motion of the magnitude and also of each component
of the order parameter vectors between two populations (see Inset of Fig. 5
(c)). In the microscopic dynamics, we numerically found that the chimera
states of these three types also possess conserved quantities along the
trajectory. More precisely, $C_{1}(t)=0$ and $C_{2}(t)=0$ with

$\displaystyle C_{1}(t)$ $\displaystyle:=\sin(\Theta_{1}(t)-\Theta_{2}(t))$
$\displaystyle C_{2}(t)$
$\displaystyle:=\big{(}|m^{(1)}_{1}(t)||m^{(2)}_{2}(t)|-|m^{(2)}_{2}(t)||m^{(1)}_{1}(t)|\big{)}^{2}$
(31)

Here, $\bm{m}^{(1)}=(m^{(1)}_{1},m^{(1)}_{2})^{\top}\in\mathbb{C}^{2}$ and
$\bm{m}^{(2)}=(m^{(2)}_{1},m^{(2)}_{2})^{\top}\in\mathbb{C}^{2}$ are the
Kuramto order parameters in Eq. (23), and $\Theta_{i}=\arg m^{(1)}_{i}-\arg
m^{(2)}_{i}$ for $i=1,2$. As in the thermodynamic limit, increasing the
parameter to be $A>A_{c}$, $C_{1}$ and $C_{2}$ become time-dependent. In Fig.
5 (d), the time evolution of the order parameters is depicted for an example
trajectory for $A>A_{c}$ where the conserved quantities are broken. Also, the
components of the order parameter vectors do not show the periodic alternation
between two populations but rather they exhibit aperiodic dynamics (see Inset
of Fig. 5 (d)). Consequently, our observation of chimera states in the
thermodynamic limit can be also verified in the ensembles of a finite number
of oscillators. In Fig. 5 (e), a snapshot of a stationary chimera state in a
system of 100 oscillators is shown, obtained from Eq. (21). The synchronized
oscillators (blue dots) behave alike altogether along the trajectory of the
order parameter $\bm{m}^{(2)}(t)$ (blue curve). The incoherent oscillators
(red dots) are distributed in the phase space and the order parameter
$\bm{m}^{(1)}(t)$ (red curve) has lower magnitude than $\bm{m}^{(2)}(t)$.

Concerning the distribution of the oscillators in the incoherent population,
one can at least study the distribution of the incoherent oscillators along
the direction of $\bm{\psi}_{1}$ as follows. Here, we assume the stationary
chimera state for a given $A\in(A_{\text{LP}},A_{\text{HB}})$ for which
$\bm{\psi}_{1}$ is saturated at a stationary value after a transient behavior.
First, define an angular variables of the incoherent oscillators along the
order parameter direction:

$e^{i\varphi}:=\frac{1}{|\bm{\psi}_{1}|}\bm{\psi}_{1}^{\dagger}\bm{x}^{(1)}~{}~{}~{}~{}\text{and}~{}~{}~{}~{}e^{i\varphi_{0}}:=\frac{1}{|\bm{\psi}_{1}|}\bm{\psi}_{1}^{\dagger}\bm{x}_{0}^{(1)}$
(32)

in the thermodynamic limit. Using the generalized Watanabe-Strogatz
transformation in Eq. (9), we obtain

$\displaystyle e^{i\varphi}$
$\displaystyle:=\frac{1}{|\bm{\psi}_{1}|}\bm{\psi}_{1}^{\dagger}\bm{x}^{(1)}=\frac{1}{|\bm{\psi}_{1}|}\frac{\bm{\psi}_{1}^{\dagger}\bm{A}\bm{x}^{(1)}_{0}+|\bm{\psi}_{1}|^{2}}{|\bm{\psi}_{1}|e^{i\varphi_{0}}+1}$
$\displaystyle=\frac{1}{|\bm{\psi}_{1}|}\frac{\bm{\psi}_{1}^{\dagger}\bm{x}^{(1)}_{0}+|\bm{\psi}_{1}|^{2}}{|\bm{\psi}_{1}|e^{i\varphi_{0}}+1}=\frac{e^{i\varphi_{0}}+|\bm{\psi}_{1}|}{|\bm{\psi}_{1}|e^{i\varphi_{0}}+1}$
(33)

with $\bm{A}=\bm{H}^{1/2}\bm{U}$ and $\bm{H}^{1/2}\bm{\psi}_{1}=\bm{\psi}_{1}$
where we can set $\bm{U}=I_{M}$ as long as $\bm{x}^{(1)}_{0}$ are uniformly
distributed on $\mathbb{S}^{M}$. Algebraically rearranging Eq. (33) and
following the same argument in Sec. 2, we obtain

$e^{i\varphi_{0}}=\frac{e^{i\varphi}-|\bm{\psi}_{1}|}{1-|\bm{\psi}_{1}|e^{i\varphi}}.$
(34)

Denoting $\varphi=T(\varphi_{0})$ and $\varphi_{0}=T^{-1}(\varphi)$ (inverse
transformation), one obtains $d(T_{*}\mu)(\varphi)=f(\varphi)d\varphi$ where
$f(\varphi):=\frac{1}{2\pi}\partial_{\varphi}T^{-1}(\varphi)$ is the phase
distribution function, and $d\mu(\varphi_{0})=\frac{1}{2\pi}d\varphi_{0}$
since $\bm{x}_{0}^{(1)}$ are uniformly distributed on $\mathbb{S}^{M}$. Then,
we obtain

$\displaystyle f(\varphi)$
$\displaystyle=\frac{1}{2\pi}\partial_{\varphi}T^{-1}(\varphi)=\frac{1}{2\pi
i}\bigg{(}\frac{ie^{i\varphi}}{e^{i\varphi}-|\bm{\psi}_{1}|}+\frac{i|\bm{\psi}_{1}|e^{i\varphi}}{1-|\bm{\psi}_{1}|e^{i\varphi}}\bigg{)}$
$\displaystyle=\frac{1}{2\pi}\frac{1-|\bm{\psi}_{1}|^{2}}{1-2|\bm{\psi}_{1}|\cos\varphi+|\bm{\psi}_{1}|^{2}}$
(35)

which is the normalized Poisson kernel distribution. This reminds us of the
Ott-Antonsen manifold for the KS oscillators on the unit circle of
$\mathbb{C}^{1}$ where oscillators’ phases are distributed according to the
normalized Poisson kernel in Eq. (4). For the higher dimensional case
$\mathbb{C}^{2}$, we also find that the phase distribution along the
generalized WS variable satisfies the normalized Poisson kernel. In Fig. 5
(f), the histogram of the distribution (blue bars) of
$\\{\varphi_{j}:=\frac{1}{|\bm{\psi}_{1}|}\bm{\psi}_{1}^{\dagger}\bm{x}_{j}^{(1)}\\}_{j=1}^{N}$
is shown for $N=500$ and $A=0.35$. The numerical results for the finite-sized
ensemble fits well to the analytical prediction (red curve) in Eq. (35).

From the above results, it is anticipated that the oscillators on
$\mathbb{S}^{M}$ are distributed according to the higher-dimensional Poisson
kernel in $\mathbb{C}^{M}$. Assume that the oscillators on $\mathbb{S}^{M}$
are distributed by

$\displaystyle
f_{\text{complex}}(\bm{x},\bm{\psi};t)=\frac{1}{S_{2M}}\frac{1-|\bm{\psi}|^{2}}{|\bm{\psi}-\bm{x}|^{2M}}$
(36)

for $\bm{\psi}\in\mathbb{C}^{M}$ and for $\bm{x}\in\mathbb{S}^{M}$. Then, we
obtain for the Kuramoto order parameter

$\displaystyle\bm{m}(t)$
$\displaystyle=\int_{|\bm{x}|=1}\bm{x}f(\bm{x},\bm{\psi})d\bm{x}$
$\displaystyle=\frac{1-|\bm{\psi}|^{2}}{S_{2M}}\int_{|\bm{x}|=1}\frac{\bm{x}}{(1+|\bm{\psi}|^{2}-2\text{Re}\left\langle\bm{\psi}\middle|\bm{x}\right\rangle)^{M}}d\bm{x}$
$\displaystyle=\frac{1-|\bm{\psi}|^{2}}{S_{2M}}\int_{|\eta|\leq
1}\int_{|\bm{n}|=\sqrt{1-|\eta|^{2}}}\frac{\eta\hat{\bm{\psi}}+\bm{n}}{(1+|\bm{\psi}|^{2}-2\text{Re}\left\langle\bm{\psi}\middle|\bm{x}\right\rangle)^{M}}d\bm{n}\frac{1}{\sqrt{1-|\eta|^{2}}}d\eta$
(37)

where $\hat{\bm{\psi}}:=\frac{\bm{\psi}}{|\bm{\psi}|}$ is a unit vector and
$\left\langle\bm{n}\middle|\bm{\psi}\right\rangle=0$. Here, we decompose the
unit vector on the sphere into $\bm{x}=\eta\hat{\bm{\psi}}+\bm{n}$ [55]. Then,
the integral can be written as

$\displaystyle\bm{m}(t)$
$\displaystyle=\frac{1-|\bm{\psi}|^{2}}{S_{2M}}\hat{\bm{\psi}}\int_{|\eta|\leq
1}\int_{|\bm{n}|=\sqrt{1-|\eta|^{2}}}\frac{\eta}{(1+|\bm{\psi}|^{2}-2|\bm{\psi}|\text{Re}\eta)^{M}}d\bm{n}\frac{1}{\sqrt{1-|\eta|^{2}}}d\eta$
$\displaystyle=\frac{S_{2M-2}}{S_{2M}}(1-|\bm{\psi}|^{2})\hat{\bm{\psi}}\int_{|\eta|\leq
1}\frac{\eta}{(1+|\bm{\psi}|^{2}-2|\bm{\psi}|\text{Re}\eta)^{M}}d\bm{n}\frac{(1-|\eta|^{2})^{\frac{2M-3}{2}}}{\sqrt{1-|\eta|^{2}}}d\eta$
$\displaystyle=\frac{S_{2M-2}}{S_{2M}}(1-|\bm{\psi}|^{2})\hat{\bm{\psi}}\int_{0}^{1}\int_{0}^{2\pi}\frac{re^{i\theta}}{(1+|\bm{\psi}|^{2}-2|\bm{\psi}|r\cos\theta)^{M}}(1-r^{2})^{M-2}rd\theta
dr$ $\displaystyle=\bm{\psi}(t)$ (38)

which is consistent with the result above, i.e., $\bm{m}(t)=\bm{\psi}(t)$ for
the complex space. Thus, we can expect that the oscillators for the complex
spaces are distributed according to Eq. (36), i.e., the higher-dimensional
normalized Poisson kernel. However, for the real spaces, $\mathbb{R}^{M}$, the
normalized Poisson kernel is given by

$f_{\text{real}}(\bm{x},\bm{\psi};t)=\frac{1}{S_{M}}\frac{1-|\bm{\psi}|^{2}}{|\bm{\psi}-\bm{x}|^{M}}$
(39)

such that the Kuramoto order parameter reads

$\displaystyle\bm{m}(t)$
$\displaystyle=\frac{1-|\bm{\psi}|^{2}}{S_{M}}\int_{|\bm{x}|=1}\frac{\bm{x}}{(1+|\bm{\psi}|^{2}-2\left\langle\bm{\psi}\middle|\bm{x}\right\rangle)^{M/2}}d\bm{x}$
$\displaystyle=\frac{1-|\bm{\psi}|^{2}}{S_{M}}\int_{-1}^{1}\int_{|\bm{n}|=\sqrt{1-\eta^{2}}}\frac{\eta\hat{\bm{\psi}}+\bm{n}}{(1+|\bm{\psi}|^{2}-2\left\langle\bm{\psi}\middle|\bm{x}\right\rangle)^{M/2}}d\bm{n}\frac{d\eta}{\sqrt{1-\eta^{2}}}$
$\displaystyle=\frac{S_{M-1}}{S_{M}}(1-|\bm{\psi}|^{2})\hat{\bm{\psi}}\int_{-1}^{1}\frac{\eta}{(1+|\bm{\psi}|^{2}-2|\bm{\psi}|\eta)^{M/2}}\frac{(1-\eta^{2})^{\frac{M-2}{2}}}{\sqrt{1-\eta^{2}}}d\eta$
$\displaystyle=\bm{\psi}(t)$ (40)

which is inconsistent with Eq. (11). This implies that the distribution of the
oscillators on $\mathbb{S}^{M}$ for $\mathbb{R}^{M}$ is not given by Eq. (39)
with the given $\bm{\psi}$, i.e., the higher-dimensional normalized Poisson
kernel. It was reported that the real oscillators distributed according to the
higher-dimensional Poisson kernel with a given WS variable $\bm{\psi}$ in the
thermodynamic limit satisfy the OA equations introduced in Refs. [49, 50].
Therein, the spherical harmonics expansion was exploited for the oscillator
distribution function, as Ott and Antonsen used Fourier expansion for 2D real
space [32, 33]. For example, in Ref. [49], the distribution of 3D real
oscillators was assumed to be

$\displaystyle f_{\text{real}}(\theta,\phi;t)$
$\displaystyle=\sum_{l=0}^{\infty}\sum_{m=-l}^{l}f_{lm}Y_{lm}(\theta,\phi)$
$\displaystyle=\frac{1}{S_{3}}\bigg{(}1+4\pi\sum_{l=0}^{\infty}\sum_{m=-l}^{l}\rho^{l}\overline{Y_{lm}}(\Theta,\Phi)Y_{lm}(\theta,\phi)\bigg{)}$
(41)

where $f_{lm}=\rho^{l}\overline{Y_{lm}}(\Theta,\Phi)$ is taken as a generalied
OA ansatz and $Y_{lm}(\theta,\phi)$ are spherical harmonics [60]. To obtain
the phase distribution function, one can use

$\displaystyle\frac{4\pi}{2l+1}\sum_{m=-l}^{l}Y_{lm}(\hat{y})\overline{Y_{lm}}(\hat{x})=P_{l}(\hat{x}\cdot\hat{y})~{}~{}~{}~{}\text{and}~{}~{}~{}~{}\sum_{l=0}^{\infty}y^{l}P_{l}(x)=\frac{1}{\sqrt{1+y^{2}-2xy}}$
(42)

where $P_{l}(x)$ are Legendre polynomials. Using the above relations, we can
reach

$\displaystyle f_{\text{real}}(\theta,\phi;t)$
$\displaystyle=\frac{1}{S_{3}}\bigg{(}1+4\pi\sum_{l=1}^{\infty}\frac{2l+1}{4\pi}\rho^{l}P_{l}(\hat{r}\cdot\hat{\rho})\bigg{)}=\frac{1}{S_{3}}\bigg{(}2\sum_{l=0}^{\infty}l\rho^{l}P_{l}(\hat{r}\cdot\hat{\rho})+\sum_{l=0}^{\infty}\rho^{l}P_{l}(\hat{r}\cdot\hat{\rho})\bigg{)}$
$\displaystyle=\frac{1}{S_{3}}\bigg{(}2\rho\frac{\partial}{\partial\rho}\sum_{l=0}^{\infty}\rho^{l}P_{l}(\hat{r}\cdot\hat{\rho})+\sum_{l=0}^{\infty}\rho^{l}P_{l}(\hat{r}\cdot\hat{\rho})\bigg{)}$
$\displaystyle=\frac{1}{S_{3}}\bigg{(}2\rho\frac{\partial}{\partial\rho}\frac{1}{\sqrt{1+\rho^{2}-2\hat{r}\cdot\hat{\rho}}}+\frac{1}{\sqrt{1+\rho^{2}-2\hat{r}\cdot\hat{\rho}}}\bigg{)}$
(43)

where
$\bm{\rho}=\rho(\sin\Theta\cos\Phi,\sin\Theta\sin\Phi,\cos\Theta)^{\top}\in\mathbb{R}^{3}$
is the OA variable (in our notation, $\bm{\psi}(t)=\bm{\rho}(t)$) and the
microscopic oscillator is represented as
$\bm{r}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)^{\top}\in\mathbb{S}^{3}$
in the thermodynamic limit. Finally, the oscillator distribution is given as
[49]

$\displaystyle
f_{\text{real}}(\theta,\phi,t)=\frac{1}{4\pi}\frac{1-\rho^{2}}{(1+\rho^{2}-2\rho\hat{\rho}\cdot\hat{r})^{3/2}}=\frac{1}{S_{3}}\frac{1-|\bm{\rho}|^{2}}{|\bm{\rho}-\bm{r}|^{3}}$
(44)

as in Eq. (39). However, it was reported in Ref. [49, 50] that the OA variable
is then governed by

$\displaystyle\dot{\bm{\rho}}=\bm{\Omega}\bm{\rho}+\frac{1}{2}(1+|\bm{\rho}|^{2})(\bm{K}\bm{\Gamma})-\big{[}\bm{\rho}^{\top}(\bm{K}\bm{\Gamma})\big{]}\bm{\rho}$
(45)

which looks different from Eq. (10). For the higher dimensional real space,
see the OA equation in Ref. [50]. In conclusion, both for the higher-
dimensional complex and real spaces, the oscillators are distributed according
to the higher-dimensional Poisson kernel in the Ott-Antonsen manifold. For the
complex oscillators, the model and the generalized WS transformation in Ref.
[55] give the correct way, and the Kuramoto order parameter exactly coincides
with the OA variable. However, to obtain the Poisson kernel distribution for
the real spaces in the OA manifold, we need to consider the spherical
harmonics expansion and the governing equations reported in Refs. [49, 50] for
the real spaces, and in that manifold, the Kuramoto order parameter exactly
equals to the OA variable. Nevertheless, we will exploit the model in Ref.
[55] for the higher dimensional real spaces below.

##  5 Chimera Dynamics in Two-population Networks for $\mathbb{R}^{4}$

In this section, as a comparison, we investigate the system of identical
generalized Kuramoto-Sakaguchi oscillators in two-population networks for
$\mathbb{K}=\mathbb{R}$ in the thermodynamic limit. The WS variables
$\bm{\psi}_{a}(t)\in\mathbb{R}^{M}$ are governed by

$\displaystyle\dot{\bm{\psi}}_{a}$
$\displaystyle=-\bm{\psi}_{a}\bm{g}_{a}^{\dagger}\bm{\psi}_{a}+\bm{g}_{a}$
$\displaystyle=-\bm{\psi}_{a}\bm{g}_{a}^{\top}\bm{\psi}_{a}+\bm{g}_{a}$ (46)

for $a=1,2$ where ${\dagger}=\top$ (i.e., hermitian adjoint=transpose) for
$\mathbb{K}=\mathbb{R}$. In this system, the mean-field forcing should be
different since the Kuramoto order parameter not exactly coincides with the WS
variable as in Eq. (11), i.e.,
$\bm{m}_{a}(t)=h(|\bm{\psi}_{a}|^{2},M)\bm{\psi}_{a}(t)$. Thus, we have to
consider

$\bm{g}_{a}:=\mu\bm{K}h(|\bm{\psi}_{a}|^{2},M)\bm{\psi}_{a}(t)+\nu\bm{K}h(|\bm{\psi}_{b}|^{2},M)\bm{\psi}_{b}(t)$
(47)

for $(a,b)=(1,2)$ or $(2,1)$. Here, the coupling matrix $\bm{K}$ is a suitable
rotational matrix introduced in Sec. 3. In this section, we use following
notations: $\bm{\psi}_{1}=(x_{1},...,x_{M})^{\top}$ and
$\bm{\psi}_{2}=(y_{1},...,y_{M})^{\top}$ where $x_{i},y_{i}\in\mathbb{R}$ for
$i=1,...,M$.

To get an overview of the observable chimera states, we first investigate the
reduced dynamics, i.e., the dynamics of the magnitude of $\bm{\psi}_{a}$.
Considering Eq. (46), we obtain

$\displaystyle\dot{\rho}_{1}$
$\displaystyle=\frac{1-\rho_{1}^{2}}{\rho_{1}}\big{(}\mu
h(\rho_{1}^{2},M)\rho_{1}^{2}\cos\alpha+\nu h(\rho_{2}^{2},M)\xi\big{)}$
$\displaystyle\dot{\rho}_{2}$
$\displaystyle=\frac{1-\rho_{2}^{2}}{\rho_{2}}\bigg{(}\mu
h(\rho_{2}^{2},M)\rho_{2}^{2}\cos\alpha+\nu h(\rho_{1}^{2},M)\big{(}\xi\cos
2\alpha+\sin 2\alpha\sqrt{\rho_{1}^{2}\rho_{2}^{2}-\xi^{2}}\big{)}\bigg{)}$
$\displaystyle\dot{\xi}$ $\displaystyle=\mu
h(\rho_{1}^{2},M)\bigg{(}\xi(1-\rho_{1}^{2})\cos\alpha+\sin\alpha\sqrt{\rho_{1}^{2}\rho_{2}^{2}-\xi^{2}}\bigg{)}$
$\displaystyle~{}~{}~{}~{}+\mu
h(\rho_{2}^{2},M)\bigg{(}\xi(1-\rho_{2}^{2})\cos\alpha-\sin\alpha\sqrt{\rho_{1}^{2}\rho_{2}^{2}-\xi^{2}}\bigg{)}$
$\displaystyle~{}~{}~{}~{}+\nu
h(\rho_{1}^{2},M)\bigg{(}(\rho_{1}^{2}-\xi^{2})\cos 2\alpha-\sin
2\alpha\sqrt{\rho_{1}^{2}\rho_{2}^{2}-\xi^{2}}\bigg{)}+\nu
h(\rho_{2}^{2},M)(\rho_{2}^{2}-\xi^{2})$ (48)

where $\rho_{a}:=|\bm{\psi}_{a}|$ for $a=1,2$ and
$\xi:=\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle$
indicates the cross term. For the derivation of Eq. (48), see C. Note that the
dynamics of the magnitude of $\bm{\psi}_{a}$ here depends on the dimension $M$
whereas it is independent of $M$ for the complex space as in Eq. (27).

Figure 6: (a) Bifurcation diagram of stable (solid curve) and unstable (dashed
curve) chimera states in two-population networks. Red: stationary chimeras and
Green: breathing chimera states. LP: limiting point bifurcation, HB:
supercritical Hopf bifurcation, Hom: homoclinic bifurcation, and $A_{c}$ from
which on a chaotic chimera appears. (b) Conserved quantities as a function of
the parameter $A$. $t_{0}=20,000$ and $t_{1}=10^{5}$. (c) Maximum and minimum
values of the magnitude of the WS variable $\bm{\psi}_{1}$ as a function of
the parameter $A$. All the results are obtained for $\mathbb{R}^{4}$. (d)
Bifurcation points $A_{\text{LP}}$ and $A_{\text{HB}}$ as a function of the
dimension $M$.

Below, we focus on the dynamics for $\mathbb{R}^{4}$ unless stated otherwise.
In Fig. 6 (a), a bifurcation diagram from Eq. (48) is depicted as the
parameter $A$ varies by assuming again $\rho_{2}=1$. As similar to Fig. 4 (a),
a stable (red, solid) and unstable (red, dashed) stationary chimeras emerge in
a limiting point (LP) bifurcation, however, at rather larger value of
$A_{\text{LP}}$ compared to that of $\mathbb{C}^{2}$. Then, the stable
stationary chimera state undergoes a supercritical Hopf bifurcation, producing
a stable limit-cycle solution, i.e., a breathing chimera state (green, solid).
The breathing chimera undergoes a homoclinic bifurcation as its period
increases indefinitely. The dynamics of the magnitude of the WS variable in
Eq. (48) only shows this scenario. However, we also observe the periodic
alternating chimera state in the full component-dynamics from Eq. (46). For a
given $A\in(A_{\text{Hom}},A_{c})$, a periodic alternating chimera state can
be observed, apparently following
$\bm{\psi}_{1}(t)=\bm{\psi}_{2}(t-\frac{T}{2})$ and
$x_{i}(t)=y_{i}(t-\frac{T}{2})$ for $i=1,...,4$ where $T$ is the period of
$\bm{\psi}_{a}$ for $a=1,2$. Likewise, we numerically find that there are two
conserved quantities along the trajectory of the stationary, breathing and
periodically alternating chimera states. The conserved quantities are similar
to those in the complex space. First, we define angular variables
$\theta_{1}:=\tan^{-1}\frac{x_{2}}{x_{1}}$,
$\theta_{2}:=\tan^{-1}\frac{x_{4}}{x_{3}}$,
$\phi_{1}:=\tan^{-1}\frac{y_{2}}{y_{1}}$,
$\phi_{2}:=\tan^{-1}\frac{x_{4}}{x_{3}}$. Then, Equation (46) is invariant
under the rotational transformation such as

$Q:=\begin{pmatrix}R(-\theta_{1})&\bm{0}\\\
\bm{0}&R(-\theta_{2})\end{pmatrix}~{}~{}~{}~{}\text{where}~{}~{}~{}~{}R(\theta)=\begin{pmatrix}\cos\theta&\sin\theta\\\
-\sin\theta&\cos\theta\end{pmatrix}$

where $Q\in\mathbb{R}^{4\times 4}$, $R(\theta)\in\mathbb{R}^{2\times 2}$, and
$\bm{0}\in\mathbb{R}^{2\times 2}$ is the zero-matrix. Therefore, we can define
the phase difference $\Theta_{1}:=\phi_{1}-\theta_{1}$ and
$\Theta_{2}:=\phi_{2}-\theta_{2}$. The first conserved quantity is

$\displaystyle C_{1}:=\sin\big{(}\Theta_{1}-\Theta_{2}\big{)}$ (49)

and one can obtain $C_{1}(t)=0$ for all $t$ for $A<A_{c}$ along a chimera
trajectory. We also find numerically another conserved quantity that reads

$\displaystyle C_{2}$
$\displaystyle:=\sum_{k=1}^{4}x_{k}^{2}\sum_{k=1}^{4}y_{k}^{2}-\bigg{(}\sqrt{x_{1}^{2}+x_{2}^{2}}\sqrt{y_{1}^{2}+y_{2}^{2}}+\sqrt{x_{3}^{2}+x_{4}^{2}}\sqrt{y_{3}^{2}+y_{4}^{2}}\bigg{)}^{2}$
(50)

which arises from the relation

$\displaystyle\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=|\bm{\psi}_{1}||\bm{\psi}_{2}|\cos\Theta$
$\displaystyle=\sqrt{x_{1}^{2}+x_{2}^{2}}\sqrt{y_{1}^{2}+y_{2}^{2}}\cos\Theta+\sqrt{x_{3}^{2}+x_{4}^{2}}\sqrt{y_{3}^{2}+y_{4}^{2}}\cos\Theta$
(51)

where $\Theta:=\Theta_{1}=\Theta_{2}$ since $C_{1}=0$. In Fig. 6 (b),
$\sup_{t\in[t_{0},t_{1}]}C_{1}(t)$ and $\sup_{t\in[t_{0},t_{1}]}C_{1}(t)$ are
shown as the parameter $A$ varies for $t_{0}=2\times 10^{4}$ and
$t_{1}=10^{5}$. This numerically verifies that the two quantities are indeed
conserved along a chimera trajectory for $A<A_{c}$. On the other hand, we
observe the chimera state that breaks the two conserved quantities for
$A>A_{c}$ (Fig. 6 (b)). This chimera state, similar as for $\mathbb{C}^{2}$
(cf. Fig. 3), also shows aperiodic motion of the components of the WS
variables $\bm{\psi}_{a}$ for $a=1,2$ (not shown here). Hence, for
$\mathbb{R}^{4}$, we observe as well that aperiodic chimera dynamics appears
via breaking the conserved quantities in Eqs. (49-50).

Furthermore, also in the real space, a global bifurcation appears involving
the completely incoherent state. However, this global bifurcation is observed
for the periodic alternating chimera state for $A>A_{\text{Hom}}$ rather than
the breathing chimeras, which is somewhat different from the case of the
complex space (Fig. 4). In Fig. 6 (c), the maximum and the minimum values of
$|\bm{\psi}_{\text{incoh}}|$ are depicted with the parameter $A$ changed. As
$A$ increases, for breathing chimera state,
$\text{min}|\bm{\psi}_{\text{incoh}}|$ continuously decreases, and then
seamlessly connects to that of the periodic alternating chimera state for
$A>A_{\text{Hom}}$. The phase variables $\Theta(t)$ of these chimeras evolve
within some interval smaller than $\mathbb{T}$, similar to Fig. 4 (d-1).
Further increasing $A$, the alternating chimera state touches the completely
incoherent state of $|\bm{\psi}_{1}|=0$ and its phase variable monotonically
increases as a function of time as in Fig. 4 (d-2).

Finally, we note that the dynamics of the magnitude of the WS variables
depends on the dimension $M$. For small $A$, the scenarios of the emergence of
the chimera states in low-dimensional real spaces reveals that the stationary
chimera state is born/annihilated at a limiting point bifurcation
($A_{\text{LP}}$) and it undergoes a supercritical Hopf bifurcation at
$A=A_{\text{HB}}$. However, Figure 6 (d) displays the parameter interval
between the limiting point and the Hopf bifurcation, i.e.,
$(A_{\text{LP}},A_{\text{HB}})$. It decreases as $M$ increases. This indicates
that for the higher-dimensional real systems, it becomes harder to obtain
stationary chimera states. Not only for the stationary chimeras, but other
chimera states also are hardly observed within our numerical ability for
higher dimensional spaces.

##  6 Summary and Outlook

In this paper, we explored a system of identical generalized Kuramoto-
Sakaguchi oscillators defined on the surface of the unit sphere in two-
population networks. We exploited the model proposed in Ref. [55] to take
advantage of the extended Watanabe-Strogatz transformation that describes the
macroscopic dynamics of the system. First, we introduced a suitable coupling
matrix both for the real and complex spaces, in line with the standard KS
model defined on the unit circle and its phase-lag parameter, which becomes
relevant to determine the Benjamin-Feir instability point. For the 2D complex
space $\mathbb{C}^{2}$, particularly in the thermodynamic limit, stationary
chimeras are created/destroyed at a limiting point bifurcation, and the stable
chimera undergoes a supercritical Hopf bifurcation that produces a stable
breathing chimera state. In this system, the breathing chimera states undergo
a global bifurcation involving the completely incoherent state, which changes
their dynamics in terms of the phase variables. Beyond this global
bifurcation, the periodic alternating chimera dynamics appears seamlessly
connecting to the breathing chimeras in terms of the minimum value of the
order parameter magnitude. These three types of chimera states possess two
conserved quantities along their trajectory in phase space. However, when the
coupling strength is further weakened, the chimera trajectory shows
componentwise aperiodic dynamics via the breaking of the conserved quantities.

To back up this further, we also explore a finite-sized ensemble of identical
generalized KS oscillators in two-populations for $\mathbb{C}^{2}$. The
finite-sized system exhibits states corresponding to each chimera state
mentioned above in terms of the Kuramoto order parameter. Furthermore, the
microscopic dynamics conserves the quantities defined by the components of the
Kuramoto order parameters, which then break for larger parameter values of
$A$. Finally, we compared this result of $\mathbb{C}^{2}$ to the dynamics for
$\mathbb{R}^{4}$ in the thermodynamic limit. The chimera states in the real
space again follow a similar scenario, including stationary and breathing
chimeras, conserved quantities, and also the breaking of the latter to induce
componentwise aperiodic chimera dynamics. However, in this case, the global
bifurcation occurs for the alternating chimera states rather than for the
breathing chimera states.

A system of coupled oscillators in two-population networks has been considered
by many researchers for the study of chimera states, as well as many
variations of this topology. Accordingly, one might explore chimera states of
generalized Kuramoto-Sakaguchi oscillators in three- and multi-population
networks. To consider more realistic situations, heterogeneities on the system
could be imposed on the oscillator ensemble. It would be also an interesting
application to study chimeras of heterogeneous generalized KS oscillators. For
the heterogeneous natural frequency distribution, one should consider the
macroscopic dynamics in Refs. [55, 50] for complex and real spaces,
respectively. Furthermore, it is also possible to arrange the generalized
Kuramoto oscillators along a ring geometry with nonlocal coupling, and then
investigate chimera states [65]. Apart from these, there could be many other
applications of the generalized Kuramoto-Sakaguchi model that provides us with
further understandings of chimera dynamics.

##  Appendix A Derivation of Eqs. (11-12)

In Ref. [55], the author obtained Eqs. (11-12) by directly calculating the
definition of the order parameter in the thermodynamic limit. In this
Appendix, we adopt a different approach to obtain the same result, following a
similar method as introduced by Pikovsky and Rosenblum. [37, 38]. Consider the
order parameter defined in Eq. (8) with the WS transformation in Eq. (9):

$\displaystyle\bm{m}(t)$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\bm{x}_{k}(t)=\frac{1}{N}\sum_{k=1}^{N}\frac{\bm{A}\bm{x}_{0,k}+\bm{b}}{\bm{b}^{\dagger}\bm{A}\bm{x}_{0,k}+1}.$

Substituting $\bm{A}=\bm{H}^{1/2}\bm{U}$ and $\bm{H}^{1/2}\bm{b}=\bm{b}$ into
the above equation, we obtain

$\displaystyle\bm{b}^{\dagger}\bm{m}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}(\bm{b}^{\dagger}\bm{H}^{1/2}\bm{U}\bm{x}_{0,k}+|\bm{b}|^{2})\frac{1}{1+\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k}}$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}(\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k}+|\bm{b}|^{2})\sum_{\ell=0}^{\infty}(-1)^{\ell}(\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k})^{\ell}$
$\displaystyle=-\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell=0}^{\infty}(-1)^{\ell+1}(\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k})^{\ell+1}+|\bm{b}|^{2}\frac{1}{N}\sum_{k=1}^{N}\sum_{\ell=0}^{\infty}(-1)^{\ell}(\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k})^{\ell}$
$\displaystyle=|\bm{b}|^{2}\bigg{(}1+(1-|\bm{b}|^{-2})\sum_{\ell=2}^{\infty}(-1)^{\ell}C_{\ell}\bigg{)}$
(52)

where
$C_{\ell}:=\frac{1}{N}\sum_{k=1}^{N}(\bm{b}^{\dagger}\bm{U}\bm{x}_{0,k})^{\ell}$
and $C_{1}=0$. Taking the thermodynamic limit with uniformly distributed
constants of motion, we can set $\bm{U}=I_{M}$ and write $C_{\ell}$ as

$\displaystyle
C_{\ell}=\begin{dcases}\frac{1}{S_{M}}\int_{|\bm{x}_{0}|=1}(\bm{b}^{\dagger}\bm{x}_{0})^{\ell}d\bm{x}_{0},&\text{for}~{}~{}\mathbb{K}=\mathbb{R}\\\
\\\
\frac{1}{S_{2M}}\int_{|\bm{x}_{0}|=1}(\bm{b}^{\dagger}\bm{x}_{0})^{\ell}d\bm{x}_{0},&\text{for}~{}~{}\mathbb{K}=\mathbb{C}\end{dcases}$
(53)

Let us first consider the real space, i.e., $\mathbb{K}=\mathbb{R}$. Using the
spherical coordinate systems, one can write a position of an oscillator on the
surface of the unit ball in $\mathbb{R}^{M}$ as
$x_{1}=\sin\theta_{1}\sin\theta_{2}\cdots\sin\theta_{M-2}\cos\phi$,
$x_{2}=\sin\theta_{1}\sin\theta_{2}\cdots\sin\theta_{M-2}\sin\phi$, …,
$x_{M-2}=\sin\theta_{1}\sin\theta_{2}\cos\theta_{3}$,
$x_{M-1}=\sin\theta_{1}\cos\theta_{2}$, and $x_{M}=\cos\theta_{1}$ where
$\theta_{1},...,\theta_{M-2}\in[0,\pi]$ and $\phi\in[0,2\pi]$. Here, $x_{M}$
corresponds to the $z$-axis, for example, in the 3D space. We call this axis
the $M$-axis and we refer to the plane perpendicular to the $M$-axis as the
$M^{\perp}$-(hyper)plane throughout this paper. Without loss of generality, we
can assume that $\bm{b}$ is aligned along the $M$-axis, i.e.,
$\bm{b}=|\bm{b}|\hat{x}_{M}$. Consequently, we obtain

$\displaystyle C_{\ell}$
$\displaystyle=\frac{1}{S_{M}}|\bm{b}|^{\ell}\int_{|\bm{x}_{0}|=1}\cos^{\ell}\theta_{1}d^{M}V=\frac{|\bm{b}|^{\ell}}{S_{M}}S_{M-1}\int\cos^{\ell}\theta_{1}\sin^{M-2}\theta_{1}d\theta_{1}$
$\displaystyle=|\bm{b}|^{\ell}\frac{1+(-1)^{\ell}}{2\sqrt{\pi}}\frac{\Gamma\big{(}\frac{\ell+1}{2}\big{)}\Gamma\big{(}\frac{M}{2}\big{)}}{\Gamma\big{(}\frac{\ell+M}{2}\big{)}}$

which leads to

$\displaystyle\sum_{\ell=2}^{\infty}(-1)^{\ell}C_{\ell}$
$\displaystyle=|\bm{b}|^{2}\frac{\Gamma\big{(}\frac{M}{2}\big{)}}{2\Gamma\big{(}1+\frac{M}{2}\big{)}}\,_{2}F_{1}(1,\frac{3}{2};\frac{M+2}{2};|\bm{b}|^{2}).$

Finally, Equation (52) can be written as

$\displaystyle\bm{b}^{\dagger}\bm{m}$
$\displaystyle=|\bm{b}|^{2}\bigg{(}1+\frac{|\bm{b}|^{2}-1}{M}\,_{2}F_{1}(1,\frac{3}{2};\frac{M+2}{2};|\bm{b}|^{2})\bigg{)}$
$\displaystyle=\bm{b}^{\dagger}\bm{b}\frac{M-1}{M}\,_{2}F_{1}(\frac{1}{2},1;\frac{M+2}{2};|\bm{b}|^{2}).$

Hence, we can assume that

$\displaystyle\bm{m}(t)=\frac{M-1}{M}\,_{2}F_{1}(\frac{1}{2},1;\frac{M+2}{2};|\bm{b}|^{2})\bm{b}(t).$
(54)

For $\mathbb{K}=\mathbb{C}$, obtaining $\bm{m}(t)=\bm{b}(t)$ is
straightforward as the coefficients $C_{\ell}=0$ for $\ell\neq 0$.

##  Appendix B  Benjamin-Feir Instability Point for $\mathbb{R}^{3}$

Figure 7: (a) A time-parametric plot shows a trajectory initiated close to the
origin on the $M^{\perp}$-plane (blue). This trajectory gradually converges
towards a limit-cycle solution on the unit circle within the
$M^{\perp}$-plane, indicating a synchronized solution. If a slight
perturbation is introduced along the $M$-axis, the trajectory is abruptly
redirected towards the north pole (red). (b) The plot illustrates the time
evolution of the components of the WS variable $\bm{\psi}^{*}_{\perp}(t)$. The
rotation occurs at a tangential speed of $\sin\alpha$. Blue:
$(\bm{\psi}_{\perp}^{*})_{1}$, Orange: $(\bm{\psi}_{\perp}^{*})_{2}$ and
Green: $(\bm{\psi}_{\perp}^{*})_{3}=0$.

In this section, we provide an example of stability analysis for the
synchronized state in $\mathbb{R}^{3}$. First, consider the synchronized state
on the $M$-axis $\bm{\psi}^{*}=(0,...,0,\pm 1)^{\top}$, i.e., the north and
the south poles, respectively. In fact, the dynamics asymptotically approaches
this solution as $t\rightarrow\infty$ since
$h(|\bm{\psi}|^{2},M)=\frac{|\bm{\psi}|-(1-|\bm{\psi}|^{2})\tanh^{-1}|\bm{\psi}|}{|\bm{\psi}|^{3}}$
for $M=3$ and $\tanh^{-1}(1)=\infty$ while $h(|\bm{\psi}|^{2},M)\rightarrow 1$
as $|\bm{\psi}|\rightarrow 1^{-}$. To study the stability of this solution,
let $\bm{\psi}^{*}=(0,0,1-\delta)^{\top}$. Then, the trajectory approaches
this fixed point solution on the north pole: $\lim_{\delta\rightarrow
0^{+}}F(\bm{\psi}^{*})=(0,0,0)^{\top}$ where
$F(\bm{\psi}):=-\bm{\psi}\bm{g}^{\dagger}\bm{\psi}+\bm{g}$. The Jacobian
matrix evaluated at this solution is given by

$\displaystyle\lim_{\delta\rightarrow
0^{+}}J\bigg{|}_{\bm{\psi}^{*}}=\begin{pmatrix}-1+\cos\alpha&-\sin\alpha&0\\\
\sin\alpha&-1+\cos\alpha&0\\\ 0&0&-2\end{pmatrix}$ (55)

and its eigenvalues are $\lambda_{1}=-2,\lambda_{\pm}=-1+e^{\pm i\alpha}$.
Hence, the synchronized solution on the north pole is always a stable solution
regardless of $\alpha$. If the trajectory starts slightly outside the
$M^{\perp}$-plane, the solution asymptotically approaches the north pole (or
the south pole, depending on whether it initiates above or below the plane.)
as time goes on (Fig. 7 (a)).

Next, we find a trajectory starting on the $M^{\perp}$-plane for
$\alpha<\frac{\pi}{2}$ asymptotically approaching a limit-cycle solution
rotating around the great circle on $x_{1}x_{2}$-plane (Fig. 7 (a-b)). The
limit-cycle solution can also be interpreted as the synchronized state, as the
magnitude of the WS variable is equal to unity. First, letting
$\bm{\psi}^{*}_{\perp}(t)=\varepsilon(\cos\theta(t),\sin\theta(t),0)^{\top}$
gives $\lim_{\varepsilon\rightarrow
1^{-}}F(\bm{\psi}^{*}_{\perp})=\sin\alpha(-\sin\theta(t),\cos\theta(t),0)^{\top}$.
This demonstrates that the trajectory asymptotically converges to the
synchronized limit-cycle trajectory, which rotates counterclockwise along the
unit circle with a tangential speed of $\sin\alpha$. Due to its rotational
symmetry, we can set the synchronized solution as a fixed point on the unit
circle like $\bm{\psi}^{*}_{\perp}=\lim_{\varepsilon\rightarrow
1^{-}}\varepsilon(1,0,0)^{\top}$. Subsequently, the eigenvalue of the Jacobian
matrix, evaluated at this solution, with the corresponding eigendirection
along the $M$-axis, is given by $\lambda_{3}=1-\cos\alpha$. Thus, the
synchronized limit-cycle solution is always unstable along the $M$-axis, while
it remains stable on the $M^{\perp}$-plane for $\alpha<\alpha_{\text{BF}}$.

##  Appendix C Derivation of Eq. (48)

To obtain Eq. (48), we consider

$\displaystyle\partial_{t}|\bm{\psi}_{1}|^{2}$
$\displaystyle=\partial_{t}(\bm{\psi}_{1}^{\dagger}\bm{\psi}_{1})$
$\displaystyle=-\mu
h_{1}\big{(}\bm{\psi}_{1}^{\dagger}\bm{K}\bm{\psi}_{1}\bm{\psi}_{1}^{\dagger}-\bm{\psi}_{1}^{\dagger}\bm{K}^{\dagger}\big{)}\bm{\psi}_{1}-\nu
h_{2}\big{(}\bm{\psi}_{1}^{\dagger}\bm{K}\bm{\psi}_{2}\bm{\psi}_{1}^{\dagger}-\bm{\psi}_{2}^{\dagger}\bm{K}^{\dagger}\big{)}\bm{\psi}_{1}$
$\displaystyle~{}~{}~{}~{}-\mu
h_{1}\bm{\psi}_{1}^{\dagger}\big{(}\bm{\psi}_{1}\bm{\psi}_{1}^{\dagger}\bm{K}^{\dagger}\bm{\psi}_{1}-\bm{K}\bm{\psi}_{1}\big{)}-\nu
h_{2}\bm{\psi}_{1}^{\dagger}\big{(}\bm{\psi}_{1}\bm{\psi}_{2}^{\dagger}\bm{K}^{\dagger}\bm{\psi}_{1}-\bm{K}\bm{\psi}_{2}\big{)}$
$\displaystyle=2(1-|\bm{\psi}_{1}|^{2})\bigg{(}\mu
h_{1}\left\langle\bm{\psi}_{1}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}_{1}\right\rangle+\nu
h_{2}\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle\bigg{)}$
$\displaystyle\partial_{t}|\bm{\psi}_{2}|^{2}$
$\displaystyle=\partial_{t}(\bm{\psi}_{2}^{\dagger}\bm{\psi}_{2})$
$\displaystyle=2(1-|\bm{\psi}_{2}|^{2})\bigg{(}\mu
h_{2}\left\langle\bm{\psi}_{2}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}_{2}\right\rangle+\nu
h_{1}\left\langle\bm{\psi}_{1}\middle|\bm{K}^{\dagger}\middle|\bm{\psi}_{2}\right\rangle\bigg{)}$
$\displaystyle\partial_{t}\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=-\mu
h_{2}\big{(}|\bm{\psi}_{1}|^{2}\cos\alpha\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle-\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle\big{)}-\nu
h_{2}\big{(}\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle^{2}-|\bm{\psi}_{2}|^{2}\big{)}$
$\displaystyle~{}~{}~{}~{}-\mu
h_{2}\big{(}|\bm{\psi}_{2}|^{2}\cos\alpha\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle-\left\langle\bm{\psi}_{1}\middle|\bm{K}^{2}\middle|\bm{\psi}_{2}\right\rangle\big{)}$
$\displaystyle~{}~{}~{}~{}-\nu
h_{1}\big{(}\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle\left\langle\bm{\psi}_{1}\middle|\bm{K}^{\dagger}\middle|\bm{\psi}_{2}\right\rangle-|\bm{\psi}_{1}|^{2}\cos
2\alpha\big{)}$ (56)

where $h_{a}:=h(|\bm{\psi}_{a}|^{2},M)$ for $a=1,2$. Then, we consider the
cross term with the coupling matrix as a variable, i.e.,
$\xi:=\left\langle\bm{\psi}_{1}\middle|\bm{K}\middle|\bm{\psi}_{2}\right\rangle=|\bm{\psi}_{1}||\bm{\psi}_{2}|\cos\theta$
leads to

$\displaystyle\left\langle\bm{\psi}_{1}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=|\bm{\psi}_{1}||\bm{\psi}_{2}|\cos(\theta-\alpha)=\xi\cos\alpha+\sin\alpha\sqrt{|\bm{\psi}_{1}|^{2}|\bm{\psi}_{2}|^{2}-\xi^{2}}$
$\displaystyle\left\langle\bm{\psi}_{1}\middle|\bm{K}^{\dagger}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=|\bm{\psi}_{1}||\bm{\psi}_{2}|\cos(\theta-2\alpha)=\xi\cos
2\alpha+\sin 2\alpha\sqrt{|\bm{\psi}_{1}|^{2}|\bm{\psi}_{2}|^{2}-\xi^{2}}$
$\displaystyle\left\langle\bm{\psi}_{1}\middle|\bm{K}^{2}\middle|\bm{\psi}_{2}\right\rangle$
$\displaystyle=|\bm{\psi}_{1}||\bm{\psi}_{2}|\cos(\theta+\alpha)=\xi\cos\alpha-\sin\alpha\sqrt{|\bm{\psi}_{1}|^{2}|\bm{\psi}_{2}|^{2}-\xi^{2}}.$
(57)

Note that for even $M$, we can easily find
$\left\langle\bm{\psi}_{a}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}_{a}\right\rangle=\cos\alpha|\bm{\psi}_{a}|^{2}$.
However, for odd $M$, we find that chimera states only live on
$M^{\perp}$-plane and hence we assume
$\left\langle\bm{\psi}_{a}\middle|\frac{\bm{K}+\bm{K}^{\dagger}}{2}\middle|\bm{\psi}_{a}\right\rangle=\cos\alpha|\bm{\psi}_{a}|^{2}$
for $a=1,2$. Finally, plugging Eq. (57) into Eq. (56), we obtain Eq. (48) in
Sec. 5.

## Acknowledgments

This work has been supported by the Deutsche Forschungsgemeinschaft (DFG
project KR1189/18-2). The authors would like to thank Julius Fischbach for a
fruitful discussion.

## References

## References

  * [1] Philip W Anderson.  More is different: broken symmetry and the nature of the hierarchical structure of science.  Science, 177(4047):393–396, 1972.
  * [2] Steven Strogatz, Sara Walker, Julia M. Yeomans, Corina Tarnita, Elsa Arcaute, Manlio De Domenico, Oriol Artime, and Kwang-Il Goh.  Fifty years of ‘more is different’.  Nature Reviews Physics, 4(8):508–510, Aug 2022.
  * [3] A. Pikovsky, M. Rosenblum, and J. Kurths.  Synchronization: A Universal Concept in Nonlinear Sciences.  Cambridge University Press, Cambridge, 2001.
  * [4] S. H. Strogatz.  Sync.  Hyperion, New York, 2003.
  * [5] O E Omel’chenko.  Coherence-incoherence patterns in a ring of non-locally coupled phase oscillators.  Nonlinearity, 26(9):2469–2498, jul 2013.
  * [6] O E Omel’chenko.  The mathematics behind chimera states.  Nonlinearity, 31(5):R121–R164, apr 2018.
  * [7] Mark J Panaggio and Daniel M Abrams.  Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators.  Nonlinearity, 28(3):R67–R87, feb 2015.
  * [8] Y. Kuramoto and D. Battogtokh.  Coexistence of coherence and incoherence in nonlocally coupled phase oscillators.  Nonlinear Phenom. Complex Syst., 5:380, 2002.
  * [9] Daniel M. Abrams and Steven H. Strogatz.  Chimera States for Coupled Oscillators.  Phys. Rev. Lett., 93:174102, Oct 2004.
  * [10] Daniel M. Abrams, Rennie Mirollo, Steven H. Strogatz, and Daniel A. Wiley.  Solvable model for chimera states of coupled oscillators.  Phys. Rev. Lett., 101:084103, Aug 2008.
  * [11] Ernest Montbrió, Jürgen Kurths, and Bernd Blasius.  Synchronization of two interacting populations of oscillators.  Phys. Rev. E, 70:056125, Nov 2004.
  * [12] Mark J. Panaggio, Daniel M. Abrams, Peter Ashwin, and Carlo R. Laing.  Chimera states in networks of phase oscillators: The case of two small populations.  Phys. Rev. E, 93:012218, Jan 2016.
  * [13] Erik A. Martens, Christian Bick, and Mark J. Panaggio.  Chimera states in two populations with heterogeneous phase-lag.  Chaos, 26(9):094819, 2016.
  * [14] Shuangjian Guo, Mingxue Yang, Wenchen Han, and Junzhong Yang.  Dynamics in two interacting subpopulations of nonidentical phase oscillators.  Phys. Rev. E, 103:052208, May 2021.
  * [15] Carlo R. Laing.  Disorder-induced dynamics in a pair of coupled heterogeneous phase oscillator networks.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 22(4):043104, 2012.
  * [16] Seungjae Lee and Katharina Krischer.  Attracting Poisson chimeras in two-population networks.  Chaos, 31(11):113101, 2021.
  * [17] Carlo R. Laing.  Chimeras in networks of planar oscillators.  Phys. Rev. E, 81:066221, Jun 2010.
  * [18] Carlo R. Laing.  Dynamics and stability of chimera states in two coupled populations of oscillators.  Phys. Rev. E, 100:042211, Oct 2019.
  * [19] Diego Pazó and Ernest Montbrió.  Low-dimensional dynamics of populations of pulse-coupled oscillators.  Phys. Rev. X, 4:011009, Jan 2014.
  * [20] Oleksandr Burylko, Erik A. Martens, and Christian Bick.  Symmetry breaking yields chimeras in two small populations of kuramoto-type oscillators.  Chaos, 32(9):093109, 2022.
  * [21] Simona Olmi.  Chimera states in coupled kuramoto oscillators with inertia.  Chaos, 25(12):123125, 2015.
  * [22] Simona Olmi, Erik A. Martens, Shashi Thutupalli, and Alessandro Torcini.  Intermittent chaotic chimeras for coupled rotators.  Phys. Rev. E, 92:030901, Sep 2015.
  * [23] Christian Bick and Peter Ashwin.  Chaotic weak chimeras and their persistence in coupled populations of phase oscillators.  Nonlinearity, 29(5):1468, mar 2016.
  * [24] Erik A. Martens.  Bistable chimera attractors on a triangular network of oscillator populations.  Phys. Rev. E, 82:016216, Jul 2010.
  * [25] Erik A. Martens.  Chimeras in a network of three oscillator populations with varying network topology.  Chaos, 20(4):043122, 2010.
  * [26] Seungjae Lee and Katharina Krischer.  Chaotic chimera attractors in a triangular network of identical oscillators.  Phys. Rev. E, 107:054205, May 2023.
  * [27] Carlo R. Laing.  Chimeras on a ring of oscillator populations.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(1):013121, 2023.
  * [28] Christian Bick.  Heteroclinic switching between chimeras.  Phys. Rev. E, 97:050201, May 2018.
  * [29] Christian Bick.  Heteroclinic dynamics of localized frequency synchrony: Heteroclinic cycles for small populations.  Journal of Nonlinear Science, 29(6):2547–2570, Dec 2019.
  * [30] Christian Bick, Tobias Böhle, and Christian Kuehn.  Multi-population phase oscillator networks with higher-order interactions.  Nonlinear Differential Equations and Applications NoDEA, 29(6):64, 2022.
  * [31] Seungjae Lee and Katharina Krischer.  Heteroclinic switching between chimeras in a ring of six oscillator populations.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(6), 06 2023.  063120\.
  * [32] Edward Ott and Thomas M. Antonsen.  Low dimensional behavior of large systems of globally coupled oscillators.  Chaos, 18(3):037113, 2008.
  * [33] Edward Ott and Thomas M. Antonsen.  Long time evolution of phase oscillator systems.  Chaos, 19(2):023117, 2009.
  * [34] Shinya Watanabe and Steven H. Strogatz.  Integrability of a globally coupled oscillator array.  Phys. Rev. Lett., 70:2391–2394, Apr 1993.
  * [35] Shinya Watanabe and Steven H. Strogatz.  Constants of motion for superconducting Josephson arrays.  Physica D: Nonlinear Phenomena, 74(3):197–253, 1994.
  * [36] Seth A. Marvel, Renato E. Mirollo, and Steven H. Strogatz.  Identical phase oscillators with global sinusoidal coupling evolve by möbius group action.  Chaos, 19(4):043104, 2009.
  * [37] Arkady Pikovsky and Michael Rosenblum.  Partially integrable dynamics of hierarchical populations of coupled oscillators.  Phys. Rev. Lett., 101:264103, Dec 2008.
  * [38] Arkady Pikovsky and Michael Rosenblum.  Dynamics of heterogeneous oscillator ensembles in terms of collective variables.  Physica D: Nonlinear Phenomena, 240(9):872–881, 2011.
  * [39] Carlo R. Laing.  The dynamics of chimera states in heterogeneous kuramoto networks.  Physica D: Nonlinear Phenomena, 238(16):1569–1588, 2009.
  * [40] Christian Bick, Marc Goodfellow, Carlo R. Laing, and Erik A. Martens.  Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review.  The Journal of Mathematical Neuroscience, 10(1):9, May 2020.
  * [41] Vladimir Jaćimović and Aladin Crnkić.  Low-dimensional dynamics in non-Abelian Kuramoto model on the 3-sphere.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(8), 08 2018.  083105\.
  * [42] Hyun Keun Lee, Hyunsuk Hong, and Joonhyun Yeo.  Improved numerical scheme for the generalized Kuramoto model.  Journal of Statistical Mechanics: Theory and Experiment, 2023(4):043403, 2023.
  * [43] Wei Zou, Sujuan He, D. V. Senthilkumar, and Jürgen Kurths.  Solvable Dynamics of Coupled High-Dimensional Generalized Limit-Cycle Oscillators.  Phys. Rev. Lett., 130:107202, Mar 2023.
  * [44] Sarthak Chandra, Michelle Girvan, and Edward Ott.  Continuous versus Discontinuous Transitions in the $D$-Dimensional Generalized Kuramoto Model: Odd $D$ is Different.  Phys. Rev. X, 9:011002, Jan 2019.
  * [45] Aladin Crnkić, Vladimir Jaćimović, and Marijan Marković.  On synchronization in Kuramoto models on spheres.  Analysis and Mathematical Physics, 11(3):129, 2021.
  * [46] Marcus A. M. de Aguiar.  Generalized frustration in the multidimensional Kuramoto model.  Phys. Rev. E, 107:044205, Apr 2023.
  * [47] Marcus A. M. de Aguiar.  On the numerical integration of the multidimensional kuramoto model, 2023\.
  * [48] M. A. Lohe.  Systems of matrix Riccati equations, linear fractional transformations, partial integrability and synchronization.  Journal of Mathematical Physics, 60(7), 07 2019.  072701\.
  * [49] Ana Elisa D Barioni and Marcus AM de Aguiar.  Complexity reduction in the 3D Kuramoto model.  Chaos, Solitons & Fractals, 149:111090, 2021.
  * [50] Ana Elisa D Barioni and Marcus AM de Aguiar.  Ott-Antonsen ansatz for the D-dimensional Kuramoto model: A constructive approach.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(11):113141, 2021.
  * [51] Max Lipton, Renato Mirollo, and Steven H. Strogatz.  The Kuramoto model on a sphere: Explaining its low-dimensional dynamics with group theory and hyperbolic geometry.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(9), 09 2021.  093113\.
  * [52] MA Lohe.  Higher-dimensional generalizations of the Watanabe-Strogatz transform for vector models of synchronization.  Journal of Physics A: Mathematical and Theoretical, 51(22):225101, 2018.
  * [53] Sarthak Chandra, Michelle Girvan, and Edward Ott.  Complexity reduction ansatz for systems of interacting orientable agents: Beyond the Kuramoto model.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 29(5), 05 2019.  053107\.
  * [54] Guilhermo L. Buzanello, Ana Elisa D. Barioni, and Marcus A. M. de Aguiar.  Matrix coupling and generalized frustration in Kuramoto oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 32(9), 09 2022.  093130\.
  * [55] Takuma Tanaka.  Solvable model of the collective motion of heterogeneous particles interacting on a sphere.  New Journal of Physics, 16(2):023016, feb 2014.
  * [56] Steven H. Strogatz.  From kuramoto to crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1):1–20, 2000.
  * [57] Vladimir García-Morales and Katharina Krischer.  The complex ginzburg–landau equation: an introduction.  Contemporary Physics, 53(2):79–95, 2012.
  * [58] E.M. Stein and R. Shakarchi.  Complex Analysis.  Princeton lectures in analysis. Princeton University Press, 2010.
  * [59] J.J. Sakurai.  Advanced Quantum Mechanics.  A-W series in advanced physics. Addison-Wesley Publishing Company, 1967\.
  * [60] S. Hassani.  Mathematical Physics: A Modern Introduction to Its Foundations.  Springer International Publishing, 2013.
  * [61] Wolfram Research, Inc.  Mathematica, Version 12.0.  Champaign, IL, 2022: Numerical integration was performed in NDSolve with IDA.
  * [62] A. Pikovsky and A. Politi.  Lyapunov Exponents: A Tool to Explore Complex Dynamics.  Cambridge University Press, Cambridge, 2016.
  * [63] F. Ginelli, H. Chaté, R. Livi, and A. Politi.  Covariant Lyapunov vectors.  J. Phys. A: Math. Theor., 46:254005, 2013.
  * [64] P. V. Kuptsov and U. Parlitz.  Theory and computation of covariant Lyapunov vectors.  J. Nonlinear Sci., 22:727, 2012.
  * [65] Ling-Wei Kong and Ying-Cheng Lai.  Short-lived chimera states.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 33(6):063127, 06 2023.