# Hierarchy of Exact Low-Dimensional Reductions for Populations of Coupled
Oscillators

Rok Cestnik  Department of Physics and Astronomy, University of Potsdam, Karl-
Liebknecht-Strasse 24/25, 14476, Potsdam-Golm, Germany  rokcestn@uni-
potsdam.de Arkady Pikovsky  Department of Physics and Astronomy, University of
Potsdam, Karl-Liebknecht-Strasse 24/25, 14476, Potsdam-Golm, Germany

###### Abstract

We consider an ensemble of phase oscillators in the thermodynamic limit, where
it is described by a kinetic equation for the phase distribution density. We
propose an Ansatz for the circular moments of the distribution (Kuramoto-Daido
order parameters) that allows for an exact truncation at an arbitrary number
of modes. In the simplest case of one mode, the Ansatz coincides with that of
Ott and Antonsen [Chaos 18, 037113 (2008)]. Dynamics on the extended manifolds
facilitate higher-dimensional behavior such as chaos, which we demonstrate
with a simulation of a Josephson junction array. The findings are generalized
for oscillators with a Cauchy-Lorentzian distribution of natural frequencies.

A description of complex systems in terms of a few relevant variables (order
parameters) is an indispensable tool in the theoretical analysis of
equilibrium and nonequilibrium dynamics. In many cases, such a reduction is
possible close to a bifurcation point, where a separation of timescales can be
employed to derive closed equations for a few order parameters. In this
context, a seminal breakthrough has been achieved by Ott and Antonsen (OA) Ott
and Antonsen (2008) for populations of coupled oscillators. They found an
exact low-dimensional reduction of the underlying kinetic equations in terms
of the leading Kuramoto order parameter; this reduction is valid globally and
is not restricted to a vicinity of a synchronization transition. One refers to
the configuration described by OA as the OA manifold. Since its discovery, the
OA Ansatz has been adopted in numerous studies of complex systems out of
equilibrium, such as Josephson junctions Marvel and Strogatz (2009); Marvel
_et al._ (2009); Vlasov and Pikovsky (2013), theta neurons and QIF neurons
Luke  _et al._ (2013); Laing (2014); Montbrió  _et al._ (2015); Laing (2015);
Bick  _et al._ (2020); Goldobin  _et al._ (2021), optomechanical arrays Pelka
_et al._ (2020), Kuramoto-Battogtokh chimeras Kuramoto and Battogtokh (2002);
Laing (2009); Bordyugov  _et al._ (2010); Panaggio and Abrams (2015);
Omel’chenko (2018), etc.

This Letter extends the OA approach by constructing a hierarchy of global
exact finite-dimensional reductions for the same setup. It is suitable for all
the above-mentioned applications of OA. This hierarchy includes the OA
manifold as the zero-order case. Although our variables are not the order
parameters of the ensemble, the latter can be easily calculated. Taking an
array of Josephson junctions as an example, we demonstrate that in our
approach the chaotic collective dynamics, which are impossible within the OA
manifold, are straightforwardly recovered.

We consider an ensemble of identical phase oscillators $\varphi_{n}(t)$, with
common complex-valued forcing $h(t)$ and common real-valued frequency
$\omega(t)$:

$\dot{\varphi}_{n}=\omega(t)+\text{Im}\big{[}2h(t)e^{-\mathrm{i}\varphi_{n}}\big{]}\;.$
(1)

In the thermodynamic limit, where the number of oscillators goes to infinity,
it is natural to describe the system with the phase density distribution
$P(\varphi,t)$, which obeys the kinetic equation

$\frac{\partial}{\partial
t}P+\frac{\partial}{\partial\varphi}\Big{(}(\omega-\mathrm{i}he^{-\mathrm{i}\varphi}+\mathrm{i}h^{*}e^{\mathrm{i}\varphi})P\Big{)}=0\;.$
(2)

It is convenient to express the density in terms of Fourier modes
$Z_{n}=\langle e^{\mathrm{i}n\varphi}\rangle$, which are the circular moments
of the distribution:
$P(\varphi,t)=\frac{1}{2\pi}(-1+\sum_{n=0}^{\infty}Z_{n}(t)e^{-\mathrm{i}n\varphi}+\text{c.c.})$.
These moments have a clear physical meaning – they are the Kuramoto-Daido
order parameters Kuramoto (1975); Daido (1996), characterizing order
(synchronization level) in the population. A global coupling in the population
is described by these order parameters (mean fields) via a dependence of the
driving terms $h,\omega$ on $Z_{n}$. For example, in the seminal Kuramoto-
Sakaguchi model Sakaguchi and Kuramoto (1986), frequency $\omega$ is constant
and $h\propto Z_{1}$; for coupled Josephson junctions, parameter $h$ is fixed
but $\omega$ depends on $Z_{1}$ Marvel and Strogatz (2009); Marvel  _et al._
(2009); Vlasov and Pikovsky (2013).

The equations for $Z_{n}$ (with $Z_{0}\equiv 1$ by definition), resulting from
the kinetic equation (2) read

$\frac{1}{n}\dot{Z}_{n}=\mathrm{i}\omega Z_{n}+hZ_{n-1}-h^{*}Z_{n+1}\;,\quad
n\geq 1\;.$ (3)

For the analysis below we introduce the moment exponential generating function
(EGF) $F(k,t)$, defined by the series:

$F(k,t)=\sum_{n=0}^{\infty}Z_{n}(t)\frac{k^{n}}{n!}\;,\quad F(0,t)\equiv 1\;.$
(4)

The infinite system (3) can then be recast as a partial differential equation
(PDE) for $F(k,t)$

$\frac{\partial}{\partial t}F=\mathrm{i}\omega k\frac{\partial}{\partial
k}F+hkF-h^{*}k\frac{\partial^{2}}{\partial k^{2}}F\;.$ (5)

Below we derive a hierarchy of finite-dimensional representations of the
dynamics. We will first represent all the approaches on the basis of the EGF
(4), exploring possible solutions of Eq. (5). After that, we will show what
these representations mean in terms of the moments $Z_{n}$.

The seminal OA Ansatz Ott and Antonsen (2008) corresponds to an exponential
EGF

$F(k,t)=e^{kQ(t)}\;,$ (6)

this yields the OA equation for the complex variable $Q(t)$:

$\dot{Q}=\mathrm{i}\omega Q+h-h^{*}Q^{2}\;.$ (7)

We now consider two new Ansätze that each generalize the exponential solution
(6) to an arbitrarily high-dimensional manifold.

Ansatz 1:  Since the PDE (5) is linear, we can extend the OA Ansatz (6) to an
arbitrarily large sum of exponentials:

$F(k,t)=\sum_{m=1}^{M}\beta_{0,m}e^{kQ_{m}(t)}\;,$ (8)

where $\beta_{0,m}$ are complex constants, satisfying the normalization
condition

$\sum_{m=1}^{M}\beta_{0,m}=1\;.$ (9)

A substitution of Eq. (8) in (5) yields a system of $M$ identical equations of
form (7) for the complex variables $Q_{m}(t)$. This is the first family of
finite-dimensional invariant manifolds.

Ansatz 2:  First, we represent the EGF (4) as a product of an exponential (6)
and a leftover function which we denote by $G(k,t)$:

$F(k,t)=e^{kQ(t)}G(k,t)\;,\quad G(0,t)\equiv 1\;.$ (10)

Next, we suppose that the variable $Q$ satisfies the OA equation (7).
Substituting Eqs. (10),(7) in Eq. (5) then reveals a PDE for $G(k,t)$:

$\frac{\partial}{\partial t}G=\big{[}\mathrm{i}\omega
k-2h^{*}Qk\big{]}\frac{\partial}{\partial
k}G-h^{*}k\frac{\partial^{2}}{\partial k^{2}}G\;.$ (11)

Equation (11) has trivial solutions $G=1$ and $G=e^{k\beta}$, which do not
provide new solutions of the original problem. To find nontrivial solutions,
let us represent the function $G$ as an infinite series
$G(k,t)=\sum_{r=0}^{\infty}\beta_{r}(t)\frac{k^{r}}{r!}$, where $\beta_{r}$
are complex variables ($\beta_{0}\equiv 1$). Substituting this into Eq. (11)
begets an infinite set of equations for $\beta_{r}$:

$\frac{1}{r}\dot{\beta}_{r}=\mathrm{i}\omega\beta_{r}-2h^{*}Q\beta_{r}-h^{*}\beta_{r+1}\;,\quad
r\geq 1\;.$ (12)

Because an equation for $\beta_{r}$ does not contain $\beta_{r-1}$ on the
r.h.s. (only terms $\sim\beta_{r},\beta_{r+1}$ are present), any truncation
$\beta_{r>R}=0\quad\forall R\in\mathbb{N}$ is invariant. The interesting
finite-dimensional solutions are therefore those where $G(k,t)$ is a
polynomial and the EGF has the following form:

$F(k,t)=e^{kQ(t)}\sum\limits_{r=0}^{R}\beta_{r}(t)\frac{k^{r}}{r!}\;.$ (13)

The $Q$ variable evolves according to Eq. (7) and the $\beta_{r}$ variables
according to $R$ equations (12), where $\beta_{R+1}=0$. This is the second
family of finite-dimensional invariant manifolds.

Combining Ansätze 1 and 2:  The two Ansätze (8) and (13) can be combined due
to the linearity of Eq. (5), yielding

$F(k,t)=\sum\limits_{m=1}^{M}e^{kQ_{m}(t)}\sum\limits_{r=0}^{R_{m}}\beta_{r,m}(t)\frac{k^{r}}{r!}\;,$
(14)

with normalization (9). The finite-dimensional set of equations for
$\sum_{m=1}^{M}(1+R_{m})$ complex variables $Q_{m},\beta_{r,m}$ reads:

$\displaystyle\dot{Q}_{m}$ $\displaystyle=\mathrm{i}\omega
Q_{m}+h-h^{*}Q_{m}^{2}\;,$ (15a) $\displaystyle\frac{1}{r}\dot{\beta}_{r,m}$
$\displaystyle=\mathrm{i}\omega\beta_{r,m}-2h^{*}Q_{m}\beta_{r,m}-h^{*}\beta_{r+1,m}\;.$
(15b)

This combines both approaches and constitutes the extended family of finite-
dimensional reductions, which is the main result of this Letter.

We now establish relations between the new variables $Q_{m},\beta_{r,m}$ and
the moments $Z_{n}$. For the OA case, Eqs. (6) and (4) yield $Z_{n}=Q^{n}$.
Similarly, for Ansatz 1 (8), $Z_{n}=\sum_{m=1}^{M}\beta_{0,m}Q_{m}^{n}$. Less
trivial is the expression for moments corresponding to the Ansatz 2 (13). A
factorization of an EGF with an exponential term (10) has already been
explored in mathematical literature Haukkanen (1993); Prodinger (1994), where
this factorization has been shown to correspond to a so-called modified
binomial transform

$Z_{n}=\sum\limits_{r=0}^{R}\binom{n}{r}\beta_{r}Q^{n-r}\;.$ (16)

Combining both expressions, we obtain a general representation for the moments
in terms of new variables:

$Z_{n}=\sum\limits_{m=1}^{M}\sum\limits_{r=0}^{R_{m}}\binom{n}{r}\beta_{r,m}Q_{m}^{n-r}\;.$
(17)

We stress here, that the derivation of the hierarchy above is not based on any
restrictions or assumptions. Indeed, one can prove by direct calculations SM
that expressions (17) and (15) solve the original system (3).

In the Supplemental Material SM  we additionally relate our approach to
circular cumulants, which have been recently suggested for describing the
vicinity of the OA manifold Tyulkina  _et al._ (2018, 2019). In
contradistinction to the variables $Q_{m},\beta_{r,m}$ above, the hierarchy of
cumulants cannot be truncated beyond the first one Goldobin and Dolmatova
(2019).

It is instructive to look at the phase distributions corresponding to
different finite-dimensional invariant subspaces. Using moments (17) to
calculate the density, we can express $P(\varphi)$ as a sum of elemental real-
valued contributions corresponding to individual variables $\beta_{r,m}$:

$P(\varphi)=-\frac{1}{2\pi}+\sum\limits_{m=1}^{M}\sum\limits_{r=0}^{R_{m}}P_{r,m}(\varphi)\;,$
(18)

where

$P_{r,m}(\varphi)=\frac{1}{2\pi}\Big{(}\frac{\beta_{r,m}e^{-\mathrm{i}r\varphi}}{(1-Q_{m}e^{-\mathrm{i}\varphi})^{r+1}}+\text{c.c.}\Big{)}\,.$
(19)

Observe that for $Q_{m}=0$ these simply reduce to Fourier modes. Because of
the constant $-\frac{1}{2\pi}$ term in Eq. (18), modes $P_{r,m}(\varphi)$ can
be interpreted as densities only up to a constant offset, even if the mean
value is nonzero.

Only terms with $r=0$ contribute to normalization (9). They correspond to
constants $\beta_{0,m}$ and have a mean value of
$\frac{\text{Re}[\beta_{0,m}]}{\pi}$, which can be positive or negative. All
other contributions have zero mean. The case with $r=0$ and real-valued
$\beta_{0,m}\in\mathbb{R}$, corresponds to an offset wrapped Cauchy-Lorentz
distribution (WCLD). A single WCLD represents the OA formulation, and recently
superpositions of several WCLDs were also considered in Refs. Tyulkina  _et
al._ (2019); Engelbrecht and Mirollo (2020). The case with $r=0$ and complex
$\beta_{0,m}\in\mathbb{C}$, is an offset Kato-Jones distribution (KJD) Kato
and Jones (2015), which is an asymmetric generalization of the WCLD. Notice
that a single complex $\beta_{0,m}$ is not permitted by condition (9), which
means that a single KJD is not invariant under evolution (for a single moment
in time it can be represented with two modes: $\beta_{0,1}+\beta_{0,2}=1$,
$Q_{1}=0$). Remarkably, a superposition of KJDs has recently been derived for
a seemingly unrelated case of multiharmonic coupling Tönjes and Pikovsky
(2020). Modes (19) with $r>0$ typically have $r$ local maxima, and to the best
of our knowledge, have not been studied before. In Fig. 1 we depict some
examples of modes (19); for more see Supplemental Material SM  where we also
express them in terms of real quantities only (absolute values and arguments).

The overall distribution (18) has to be non-negative and normalized, thus
imposing restrictions on $Q_{m},\beta_{r,m}$ variables. While normalization is
ensured by condition (9), non-negativity presents conditions that are not
easily expressed in general (see special cases in Supplemental Material SM ,
including the case of a single KJD Kato and Jones (2015)).

Figure 1:  Examples of elemental density contributions $P_{r,m}(\varphi)$
(19), for $Q_{m}=0.2$, depicted in polar coordinates. The first row
corresponds to real-valued $\beta_{r,m}$ and its first entry is an offset WCLD
which [for full density (18)] corresponds to the OA manifold. The second row
corresponds to purely imaginary $\beta_{r,m}$ and its first entry is an offset
KJD with mean zero since $\beta_{r,m}$ has no real component. All other
entries have zero mean by definition. In each subplot, different curves
represent different absolute values of the corresponding variable
$|\beta_{r,m}|$. Bold black circle represents zero and the gray disc inside
holds negative values (origin is -0.25). Thin outer circle marks the
$\frac{1}{2\pi}$ value, representing a uniform normalized distribution for
reference.

On the found finite-dimensional manifolds, one can start with proper initial
conditions in terms of variables $Q_{m},\beta_{r,m}$ and follow the finite-
dimensional dynamics according to Eqs. (15). The inverse problem – how to
determine whether a given initial distribution (e.g., given by its order
parameters $Z_{n}$) lies on a low-dimensional manifold – is more difficult,
because one needs to invert the infinite systems of Eqs. (17). In the case of
just Ansatz 2, binomial transform (16) can be inverted:

$\beta_{r}=\sum_{n=0}^{r}\binom{r}{n}Z_{n}(-Q)^{r-n}\;,$ (20)

but the choice of $Q$ is generally not obvious. However, for an initial
distribution with a finite number of nonzero order parameters:
$Z_{n>N}(t=0)=0$ we can set $Q(t=0)=0$ in Eq. (20) and obtain a finite number
of modes $\beta_{r}{\footnotesize(t=0)}=Z_{r}{\footnotesize(t=0)}$. Thus, the
evolution of an $N$-moment distribution can be described with $N$ equations of
system (12) plus one equation for $Q$ (7). However, finite-$N$ distributions
as a class are not invariant under evolution, since the number of nonvanishing
moments $Z_{n}$ becomes infinite immediately when $Q\neq 0$.

Because Fourier modes form a base, we conclude that for any initial
distribution one can find an approximative finite-dimensional description
based on the Fourier approximation of the initial condition. One might
formulate the problem of efficiently approximating a given distribution with a
finite number of general modes (19); however, this problem lies beyond the
scope of this Letter.

Up to now we considered a population of identical oscillators (1). In typical
situations however, one or several parameters of the oscillators and of the
coupling are distributed. Most popular are models with a distribution of
natural frequencies $\omega$, although other cases have also been considered
in the literature Montbrió and Pazó (2011); Pazó and Montbrió (2011); Vlasov
_et al._ (2014); Iatsenko  _et al._ (2014). For definiteness let us consider a
distribution of frequencies $g(\omega)$. All the expressions above are valid
for a distribution $P(\varphi,t\,|\,\omega)$, conditioned on a particular
value of $\omega$: one has $Z_{n}(t\,|\,\omega)$, $Q_{m}(t\,|\,\omega)$,
$\beta_{r,m}(t\,|\,\omega)$, etc. Global order parameters are obtained by
averaging the latter over $\omega$. For example, in the Kuramoto model with
coupling via the first global order parameter, the latter is defined as
$\mathcal{Z}_{1}(t)=\int Z_{1}(t\,|\,\omega)g(\omega)d\omega$. This
expression, together with a set of equations (15) conditioned under set values
of $\omega$, constitute a set of integro-differential equations describing the
problem.

A particular simplification of this system is possible for a Cauchy-Lorentzian
distribution of frequencies
$g(\omega)=\frac{1}{\pi}\frac{\gamma}{(\omega-\omega_{0})^{2}+\gamma^{2}}$. In
this case, one can integrate over $g(\omega)$ assuming analyticity of
variables $Q_{m},\beta_{r,m}$ in the upper complex $\omega$ half-plane, and
the integration (using residue theorem) results in a value at the pole
$\omega=\omega_{0}+\mathrm{i}\gamma$. Here the argumentation is essentially
the same as used by OA Ott and Antonsen (2008). As a result, one obtains an
exact finite-dimensional description in terms of global variables
$\mathcal{Q}_{m}(t)=Q_{m}(t\,|\,\omega_{0}+\mathrm{i}\gamma)$,
$\mathcal{B}_{r,m}(t)=\beta_{r,m}(t\,|\,\omega_{0}+\mathrm{i}\gamma)$ (which
via (17) relate to global moments $\mathcal{Z}_{n}$):

$\displaystyle\dot{\mathcal{Q}}_{m}$
$\displaystyle=(\mathrm{i}\omega_{0}-\gamma)\mathcal{Q}_{m}\;+h-h^{*}\mathcal{Q}_{m}^{2}\;,$
(21) $\displaystyle\frac{1}{r}\dot{\mathcal{B}}_{r,m}$
$\displaystyle=(\mathrm{i}\omega_{0}-\gamma)\mathcal{B}_{r,m}-2h^{*}\mathcal{Q}_{m}\mathcal{B}_{r,m}-h^{*}\mathcal{B}_{r+1,m}\;.$

One can see that in finite representations $\mathcal{B}_{r>R_{m},m}=0$, the
last term $\mathcal{B}_{R_{m},m}$ does not contain driving, so one expects
that it eventually vanishes because of dissipation $\sim\gamma$. Then, there
is no driving for the next to last $\mathcal{B}_{R_{m}-1,m}$, and this term
will vanish as well, and so on for all $\mathcal{B}_{r>0,m}$
($\mathcal{B}_{0,m}$ are constants). What is left are the $\mathcal{Q}_{m}$
variables, which are driven by the same force $h$, and will therefore with
dissipation eventually converge
$\lim_{t\to\infty}(\mathcal{Q}_{m}-\mathcal{Q}_{n})=0,\ \;\forall\;m,n$. Thus,
one expects the OA manifold
$\mathcal{B}_{r>0,m}=0,\;\mathcal{Q}_{m}=\mathcal{Q},\;\mathcal{Z}_{n}=\mathcal{Q}^{n}$,
to be attractive in this situation, as has been already argued in Refs. Ott
and Antonsen (2009); Engelbrecht and Mirollo (2020); Pietras and Daffertshofer
(2016). As a final note we mention that equations (21) are valid not only for
a Cauchy-Lorentzian distribution of frequencies, but also for oscillators
driven by Cauchy noise Tanaka (2020); Tönjes and Pikovsky (2020), where noise
strength takes the role of $\gamma$.

Now we perform a simple numerical simulation to showcase the dynamics on an
extended invariant manifold and compare them to the dynamics on the OA
manifold. We consider an array of overdamped Josephson junctions coupled via a
resistive load Watanabe and Strogatz (1994). It is described by the equations
for the Josephson phases

$\dot{\varphi}_{n}=1+a\sin(\varphi_{n})+\frac{\epsilon}{N}\sum_{m=1}^{N}\sin(\varphi_{m})\;.$
(22)

The OA description of Eq. (22) yields only regular solutions (because the
dimension of the subspace is 2 Martens  _et al._ (2009)). Below we show that
already within the lowest-order finite-dimensional reduction involving 2
complex variables, chaotic regimes are possible. System (22) belongs to the
class (1) with $\omega=1+\epsilon\,\text{Im}\langle Z_{1}\rangle$ and
$h=-\frac{a}{2}$. On the two-variable manifold of Ansatz 2 (13), the equations
in terms of $Q,\beta_{1}\equiv\beta$ are

$\displaystyle\dot{Q}\ =$ $\displaystyle\
\mathrm{i}Q(1+\epsilon\,\text{Im}[Q+\beta])-\frac{a}{2}(1-Q^{2})\;,$ (23)
$\displaystyle\dot{\beta}\ =$ $\displaystyle\
\mathrm{i}\beta(1+\epsilon\,\text{Im}[Q+\beta])+a\,Q\beta\;.$

In Fig. 2 we demonstrate chaotic behavior in this system for $a=1.5$,
$\epsilon=-0.7$ [Runge-Kutta 4th-order method with time step $2\times 10^{-3}$
was used; see Supplemental Material SM  for a corresponding simulation of a
large finite ensemble (22)]. Of the four Lyapunov exponents, two are nonzero
with values $\approx\pm 0.01$. For $\beta=0$ (i.e., on the OA manifold) all
the solutions are periodic.

Figure 2: Poincaré map (orange dots) for dynamics (23) of the Josephson
junctions array for $a=1.5$, $\epsilon=-0.7$ and initial condition
$Z_{1}{\footnotesize(t=0)}=0.4$, $Z_{n>1}{\footnotesize(t=0)}=0$. The section
is made according to the condition $\arg(\beta)=\pi$,
$\frac{d}{dt}\arg(\beta)>0$. Green line: the trajectory on the OA manifold
with the same initial $Z_{1}$.

There are other examples where the extended finite-dimensional dynamics go
beyond the ones possible in the OA formulation. In Supplemental Material SM
we present simulations for a two-population chimera Abrams  _et al._ (2008).
There, on the OA manifold (which is known for such populations of identical
oscillators to be only neutrally stable Engelbrecht and Mirollo (2020)) the
dynamics of the asynchronous population are periodic, while already the
second-order dynamics yield quasiperiodicity (cf. Refs. Pikovsky and Rosenblum
(2008); Tyulkina  _et al._ (2018, 2019)).

Finally, we shortly discuss whether the exact invariant manifolds survive the
addition of independent Gaussian white noise terms to the dynamics of the
phases (as mentioned above, for Cauchy white noise the invariant manifolds
remain valid). In this case, additional terms appear in equations for moments
(3) and in the PDE for the EGF (5). An attempt of Ansatz 2 (10) results in an
infinite system for $\beta_{r}$ which, in contradistinction to system (12),
cannot be truncated, see Supplemental Material SM  for details. Nevertheless,
truncation might yield an approximative finite-dimensional description for
ensembles with noise; this is a subject of a forthcoming research (for an
approximation in terms of circular cumulants see Ref. Tyulkina  _et al._
(2019)).

Summarizing, we have generalized the finite-dimensional description for
populations of coupled oscillators due to Ott and Antonsen Ott and Antonsen
(2008) by constructing a family of exact finite-dimensional invariant
manifolds for the dynamics. Our approach is as general and as restricted as
that of OA Ott and Antonsen (2008): it is applicable to phase dynamics with
harmonic forcing and in the thermodynamic limit only, but can be generalized
to ensembles with distributed frequencies (or other parameters), to
oscillators driven by Cauchy noise, etc. It is fully applicable to systems
previously analyzed in terms of the OA reduction, such as Josephson junctions
Marvel and Strogatz (2009); Marvel  _et al._ (2009); Vlasov and Pikovsky
(2013), theta neurons and QIF neurons Luke  _et al._ (2013); Laing (2014);
Montbrió  _et al._ (2015); Laing (2015); Bick  _et al._ (2020); Goldobin  _et
al._ (2021), and to Kuramoto-Battogtokh chimera Kuramoto and Battogtokh
(2002); Laing (2009); Bordyugov  _et al._ (2010); Panaggio and Abrams (2015);
Omel’chenko (2018). In situations where the OA manifold is stable, the
extended reductions describe relaxation to this manifold in an exact manner.
This description is applicable if the population is driven away from the OA
manifold, e.g., by a phase resetting. In situations where the OA manifold is
only neutrally stable, the extended reductions provide a family of exact
equations beyond the OA Ansatz (there is no condition of a small vicinity,
like in the cumulant description Tyulkina  _et al._ (2018)). Although the
analysis of general initial distributions and determining whether or not they
lie on one of the finite-dimensional manifolds remains a challenging problem,
we have found an important and simple class of initial states that do allow
for an exact finite-dimensional description: these are the states with a
finite number of nonvanishing moments (Kuramoto-Daido order parameters).

There exists an alternative approach to a finite-dimensional description of
oscillator populations, the Watanabe-Strogatz theory (WS) Watanabe and
Strogatz (1993). While its relation to the OA theory has been clarified in the
literature Pikovsky and Rosenblum (2008); Marvel  _et al._ (2009), a
connection to the present approach remains a subject for future studies. As
the WS theory generally predicts partial integrability for populations of
identical oscillators, we can conclude that any high-dimensional dynamics
within the found invariant manifolds have a large number of integrals of
motion.

###### Acknowledgements.

The authors thank Denis Goldobin, Lev Smirnov, Ralf Toenjes, Michael
Rosenblum, Oleh Omel’chenko, and Erik Mau for useful discussions. The work was
supported by DFG (Grant No. PI 220/21-1).

## References

  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, Chaos 18, 037113 (2008).
  * Marvel and Strogatz (2009) S. A. Marvel and S. H. Strogatz, Chaos 19, 013132 (2009).
  * Marvel  _et al._ (2009) S. A. Marvel, R. E. Mirollo, and S. H. Strogatz, Chaos 19, 043104 (2009).
  * Vlasov and Pikovsky (2013) V. Vlasov and A. Pikovsky, Phys. Rev. E 88, 022908 (2013).
  * Luke  _et al._ (2013) T. B. Luke, E. Barreto, and P. So, Neural Comput. 25, 3207 (2013).
  * Laing (2014) C. R. Laing, Phys. Rev. E 90, 010901(R) (2014).
  * Montbrió  _et al._ (2015) E. Montbrió, D. Pazó, and A. Roxin, Phys. Rev. X 5, 021028 (2015).
  * Laing (2015) C. R. Laing, SIAM J. Appl. Dyn. Syst. 14, 1899 (2015).
  * Bick  _et al._ (2020) C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, J. Math. Neurosci. 10, 9 (2020).
  * Goldobin  _et al._ (2021) D. S. Goldobin, M. di Volo, and A. Torcini, Phys. Rev. Lett. 127, 038301 (2021).
  * Pelka  _et al._ (2020) K. Pelka, V. Peano, and A. Xuereb, Phys. Rev. Research 2, 013201 (2020).
  * Kuramoto and Battogtokh (2002) Y. Kuramoto and D. Battogtokh, Nonlinear Phenom. Complex Syst. 5, 380 (2002).
  * Laing (2009) C. R. Laing, Physica D (Amsterdam) 238, 1569 (2009).
  * Bordyugov  _et al._ (2010) G. Bordyugov, A. Pikovsky, and M. Rosenblum, Phys. Rev. E 82, 035205(R) (2010).
  * Panaggio and Abrams (2015) M. J. Panaggio and D. M. Abrams, Nonlinearity 28, R67 (2015).
  * Omel’chenko (2018) O. E. Omel’chenko, Nonlinearity 31, R121 (2018).
  * Kuramoto (1975) Y. Kuramoto,  _in International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki, Vol. 39 (Springer, New York, 1975) p. 420.
  * Daido (1996) H. Daido, Physica D (Amsterdam) 91, 24 (1996).
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, Prog. Theor. Phys. 76, 576 (1986).
  * Haukkanen (1993) P. Haukkanen, Fibonacci Quart. 31, 28 (1993).
  * Prodinger (1994) H. Prodinger, Fibonacci Quart. 32, 412 (1994).
  * (22) See Supplemental Material at https://arxiv.org/src/2108.10143v4/anc/supplemental.pdf for derivation of the kinetic equations, relation to the circular cumulants, additional expression and depiction of density modes, numerical simulation of a chimera state and kinetic equations for the addition of Gaussian noise.
  * Tyulkina  _et al._ (2018) I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, and A. Pikovsky, Phys. Rev. Lett. 120, 264101 (2018).
  * Tyulkina  _et al._ (2019) I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, and A. S. Pikovsky, Radiophys. Quantum Electron. 61, 640 (2019).
  * Goldobin and Dolmatova (2019) D. S. Goldobin and A. V. Dolmatova, Phys. Rev. Research 1, 033139 (2019).
  * Engelbrecht and Mirollo (2020) J. R. Engelbrecht and R. Mirollo, Phys. Rev. Research 2, 023057 (2020).
  * Kato and Jones (2015) S. Kato and M. C. Jones, Biometrika 102, 181 (2015).
  * Tönjes and Pikovsky (2020) R. Tönjes and A. Pikovsky, Phys. Rev. E 102, 052315 (2020).
  * Montbrió and Pazó (2011) E. Montbrió and D. Pazó, Phys. Rev. E 84, 046206 (2011).
  * Pazó and Montbrió (2011) D. Pazó and E. Montbrió, Europhys. Lett. 95, 60007 (2011).
  * Vlasov  _et al._ (2014) V. Vlasov, E. E. N. Macau, and A. Pikovsky, Chaos 24, 023120 (2014).
  * Iatsenko  _et al._ (2014) D. Iatsenko, P. McClintock, and A. Stefanovska, Nat. Commun. 5, 4118 (2014).
  * Ott and Antonsen (2009) E. Ott and T. M. Antonsen, Chaos 19, 023117 (2009).
  * Pietras and Daffertshofer (2016) B. Pietras and A. Daffertshofer, Chaos 26, 103101 (2016).
  * Tanaka (2020) T. Tanaka, Phys. Rev. E 102, 042220 (2020).
  * Watanabe and Strogatz (1994) S. Watanabe and S. H. Strogatz, Physica D (Amsterdam) 74, 197 (1994).
  * Martens  _et al._ (2009) E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P. So, and T. M. Antonsen, Phys. Rev. E 79, 026204 (2009).
  * Abrams  _et al._ (2008) D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A. Wiley, Phys. Rev. Lett. 101, 084103 (2008).
  * Pikovsky and Rosenblum (2008) A. Pikovsky and M. Rosenblum, Phys. Rev. Lett. 101, 264103 (2008).
  * Watanabe and Strogatz (1993) S. Watanabe and S. H. Strogatz, Phys. Rev. Lett. 70, 2391 (1993).