# Exact finite-dimensional reduction for a population of noisy oscillators
and its link to Ott-Antonsen and Watanabe-Strogatz theories

Rok Cestnik  Department of Physics and Astronomy, University of Potsdam, Karl-
Liebknecht-Strasse 24/25, 14476, Potsdam-Golm, Germany  Arkady Pikovsky
Department of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-
Strasse 24/25, 14476, Potsdam-Golm, Germany

###### Abstract

Populations of globally coupled phase oscillators are described in the
thermodynamic limit by kinetic equations for the distribution densities, or
equivalently, by infinite hierarchies of equations for the order parameters.
Ott and Antonsen [Chaos 18, 037113 (2008)] have found an invariant finite-
dimensional subspace on which the dynamics is described by one complex
variable per population. For oscillators with Cauchy distributed frequencies
or for those driven by Cauchy white noise, this subspace is weakly stable and
thus describes the asymptotic dynamics. Here we report on an exact finite-
dimensional reduction of the dynamics outside of the Ott-Antonsen subspace. We
show, that the evolution from generic initial states can be reduced to that of
three complex variables, plus a constant function. For identical noise-free
oscillators, this reduction corresponds to the Watanabe-Strogatz system of
equations [Phys. Rev. Lett. 70, 2391 (1993)]. We discuss how the reduced
system can be used to explore the transient dynamics of perturbed ensembles.

> Large ensembles of globally coupled oscillators can be described by means of
> kinetic equations for the evolution of the distribution densities. These
> equations take the simplest form if the oscillators are described by their
> phases only. Alternatively, one can write an infinite set of ordinary
> differential equations for the set of order parameters (Fourier modes of the
> distribution). For oscillators driven by Cauchy white noise or for a Cauchy
> distribution of natural frequencies, the system of equations for the order
> parameters (for coupling in the first harmonics of the phase), has a
> remarkable property first discovered by Ott and Antonsen in 2008: it
> possesses an invariant manifold on which the dynamics reduces to just one
> complex equation. In this paper we extend this result by showing, that for
> arbitrary initial conditions the dynamics reduces to that of three complex
> variables. In the noise-free case of identical oscillators, our equations
> are equivalent to the system derived by Watanabe and Strogatz in 1993. The
> finite-dimensional reduction allows exact calculation of transients to the
> attracting Ott-Antonsen states, by solving a simple six-dimensional system
> of equations.

##  I Introduction

Ensembles of coupled oscillators is a popular object in studies of complex
systems, with a wide range of applications; from physical systems (lasers
Nixon  _et al._ (2013), Josephson junctions Wiesenfeld and Swift (1995);
Wiesenfeld, Colet, and Strogatz (1998), chemical reactions Totz  _et al._
(2018)), to engineering (pedestrians on a bridge Eckhardt  _et al._ (2007))
and life sciences (neurons Luke, Barreto, and So (2013); Laing (2014), nephron
cells Holstein-Rathlou  _et al._ (2001), genetic circuits Prindle  _et al._
(2012)). A common theoretical approach includes different levels of reductions
and idealizations. If the units are self-sustained periodic oscillators, and
the coupling is weak, one can perform a phase reduction, neglecting variations
of the oscillators’ amplitudes that appear in the higher orders in coupling
strength Kuramoto (1984). As a result, each oscillator is described by just
one variable on a unit circle – the phase, which enormously simplifies the
analysis. Another idealization, which is appropriate for large ensembles, is
the thermodynamic limit of an infinite number of units. This allows for a
formulation of the evolution in terms of kinetic equations for the
distribution of the phases. An important class of models are those with global
(or mean-field) coupling. Such models appear naturally, e.g., for Josephson
junctions with a common load and for pedestrians on a bridge; in other cases
(e.g., for neural ensembles) they are justified by a huge number of
interconnections between the units.

Among the setups for ensembles of globally coupled phase oscillators, the
paradigmatic Kuramoto model Kuramoto (1975a) and its generalizations Sakaguchi
and Kuramoto (1986); Acebrón  _et al._ (2005) are particularly popular. Here
one assumes a relatively simple coupling, where the dynamics of the
oscillator’s phase depends only on the first harmonics of the phase itself. To
define the coupling, one introduces mean fields which are the circular moments
of the phase distribution. Different setups with identical deterministic
units, as well as ones having different natural frequencies and/or being
driven by noise have been considered in the literature.

One of the striking properties of the Kuramoto-type models is the possibility
to reduce the dynamics to a finite-dimensional one. Watanabe and Strogatz
Watanabe and Strogatz (1993, 1994) (WS) have demonstrated that ensembles of
identical, noise-free units can be exactly reduced to three dynamical
equations (plus constants of motion). Ott and Antonsen Ott and Antonsen (2008)
(OA) found a particular family of phase distributions (wrapped Cauchy
distribution) that is invariant under the dynamical evolution. This holds not
only for identical units, but also for ones with a Cauchy distribution of
natural frequencies, and for ones driven by white Cauchy noise Tanaka (2020);
Tönjes and Pikovsky (2020). In contradistinction to WS theory, the OA
reduction is not valid for arbitrary initial states - they should belong to
the OA invariant manifold. However, because there are arguments that the OA
manifold is attracting (although not in a trivial sense, see discussion in Ott
and Antonsen (2009); Pietras and Daffertshofer (2016); Engelbrecht and Mirollo
(2020)), the OA equations correctly describe the asymptotic in time regimes.

The goal of this paper is to fill, at least partially, the gap between WS and
OA theories. We will develop, in the thermodynamic limit, a low-dimensional
description of the Kuramoto-type phase ensembles with Cauchy noise and/or
Cauchy distribution of natural frequencies, valid for arbitrary initial
conditions. Of course, this reduction contains WS and OA equations as
particular cases.

The paper is organized as follows. In section II we formulate the problem. In
section III we introduce our basic tools (generating functions) and define a
family of finite-dimensional invariant manifolds (these results have been also
presented in a short communication Cestnik and Pikovsky (2022)). Section IV
contains the main result - we show how the evolution of generic states can be
reduced to three complex variables plus a constant function. Here we also
discuss different possibilities of introducing these variables based on
initial conditions. In section V we demonstrate stability of the OA manifold
in the presence of noise. In section VI we consider identical noise-free
oscillators, and demonstrate that the dynamics reduces to the WS equations.
Section VII is devoted to the implications for the spectrum of the Lyapunov
exponents. In section VIII we discuss how our approach allows for finding the
evolution outside of the OA manifold. We conclude and discuss possible further
developments in section IX. Many technical details are shifted from the main
text to appendices.

##  II Problem formulation

In this paper we consider populations of phase oscillators, subject to global
coupling or to a global common force, in the thermodynamic limit of an
infinite number of units. Consequently, the proper description is in terms of
the phase distribution functions. Our theory is valid for a restricted class
of systems: (i) important is that the coupling/forcing is proportional to the
first harmonics of the phase only, and (ii) the oscillators can differ from
each other only in additive terms in their phase dynamics, which are either
Cauchy-distributed white noise terms, or Cauchy-distributed frequency
constants, or a combination of both. In this section we introduce these
models.

###  II.1 Ensemble of phase oscillators with independent Cauchy noise forces

We consider an ensemble of noisy phase oscillators coupled in the first
harmonic:

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

Here $\omega$ is a combination of a natural frequency and a real-valued
additive force, and $h(t)$ is a complex-valued force that includes the first
harmonic of the phase. Both these quantities can potentially depend on the
mean fields of the population, thus readily describing global coupling. There
is no restriction on these forces, e.g., they can include noise which is then
the common noise for all elements of the population (cf. Refs. Braun  _et al._
, 2012; Gong  _et al._ , 2019). The terms $\xi_{j}(t)$ represent independent,
normalized Cauchy white noise forces, with $\gamma$ being the real and
positive noise strength Chechkin  _et al._ (2003); Toenjes, Sokolov, and
Postnikov (2013); Tönjes and Pikovsky (2020); Tanaka (2020).

We consider the thermodynamic limit of infinitely many oscillators. In this
case it is natural to describe the state with the phase density function
$P(\varphi,t)$, and express the original dynamics in terms of the continuity
equation, a partial differential equation (PDE) where the Cauchy noise begets
a term with a fractional derivative on the right-hand side:

$\frac{\partial}{\partial
t}P+\frac{\partial}{\partial\varphi}\big{(}\dot{\varphi}P\big{)}=\gamma\bigg{|}\frac{\partial}{\partial\varphi}\bigg{|}P\;.$
(2)

With $\left|\frac{\partial}{\partial\varphi}\right|^{\alpha}$ one denotes an
operator, which in the Fourier representation reduces to a multiplication with
$|n|^{\alpha}$:
$\left|\frac{\partial}{\partial\varphi}\right|^{\alpha}e^{in\varphi}=-|n|^{\alpha}e^{in\varphi}$,
cf. Ref. Toenjes, Sokolov, and Postnikov, 2014. In this representation, the
usual Gaussian noise corresponds to $\alpha=2$, while the Cauchy noise
corresponds to $\alpha=1$ (other non-integer values of $\alpha$ describe
different $\alpha$-stable distributions).

The phase density is commonly expressed as a Fourier series:

$\begin{split}P(\varphi,t)&=\frac{1}{2\pi}\Big{(}1+\sum\limits_{n=1}^{\infty}Z_{n}(t)e^{-\mathrm{i}n\varphi}+\text{c.c.}\Big{)}\;,\\\
\qquad Z_{n}(t)&=\langle
e^{\mathrm{i}n\varphi}\rangle=\int\limits_{0}^{2\pi}d\varphi\,e^{\mathrm{i}n\varphi}P(\varphi,t)\;.\end{split}$
(3)

Quantities $Z_{n}$ represent complex order parameters, also known as Kuramoto-
Daido order parameters Kuramoto (1975a); Daido (1996). These circular moments
of the phase distribution are in fact the “mean fields” which may govern the
ensemble. In terms of these order parameters, the dynamics is represented as
an infinite set of ordinary differential equations (ODE):

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

(here it suffices to consider positive $n$ only, so we replace $|n|$ in the
noisy term by $n$). These equations have been discussed in Ref. Tönjes and
Pikovsky, 2020 and represent the exact dynamics of system (1) in the
thermodynamic limit, without any approximation or assumption. Normalization of
the phase density implies $Z_{0}\equiv 1$.

###  II.2 Ensemble with a Cauchy distribution of natural frequencies

Equivalent equations can also be derived to represent the case of Cauchy
distributed natural frequencies, the situation widely considered starting from
the initial formulation by Kuramoto Kuramoto (1975b, 1984). In this case we
consider the terms $\xi_{j}$ in Eq. (1) as constants with a normalized Cauchy
distribution $g(\xi)=\pi^{-1}(1+\xi^{2})^{-1}$. The total additive force
$\omega(t)+\gamma\xi_{j}$ can be interpreted as an instantaneous frequency of
oscillator $j$. If $\omega=\omega_{0}$ is a constant, then
$\omega_{0}+\gamma\xi_{j}$ is the natural frequency of oscillator $j$.
Parameter $\gamma$ describes the width of the distribution of natural
frequencies (as we will see below, this parameter plays the same role as the
strength of the Cauchy noise above, thus we use the same letter to describe
it). In our further derivation, we follow the way presented recently in Ref.
Engelbrecht and Mirollo, 2020. One introduces the parameter $\xi$ into the
distribution of phases $P(t,\varphi;\xi)$, and the equation for this
distribution (2) then reads

$\frac{\partial}{\partial
t}P\big{\rvert}_{\xi}+\frac{\partial}{\partial\varphi}\left(\left[\omega+\gamma\xi-\mathrm{i}he^{-\mathrm{i}\varphi}+\mathrm{i}h^{*}e^{\mathrm{i}\varphi}\right]P\big{\rvert}_{\xi}\right)=0\;,$
(5)

where we used compact notation $P\big{\rvert}_{\xi}\equiv P(\varphi,t;\xi)$.
Of interest are the order parameters (circular moments), averaged over the
additions to the frequency $\xi$:

$Z_{n}(t)=\int\limits_{-\infty}^{\infty}d\xi\int\limits_{0}^{2\pi}d\varphi\,e^{\mathrm{i}n\varphi}P\big{\rvert}_{\xi}g(\xi)\;.$
(6)

The main assumption allowing for explicit equations for these order parameters
is analyticity of the distribution $P(\varphi,t;\xi)$ in the upper halfplane
of complex $\xi$. This assumption has been first introduced by Ott and
Antonsen in their seminal paper Ott and Antonsen (2008). The main reason
behind it is the possibility to calculate the integrals via residue
integration. Indeed, employing the residue theorem for a contour closing the
upper halfplane in (6) and taking the only pole at $\xi=\mathrm{i}$, one
reduces (6) to $Z_{n}(t)=\int
d\varphi\,e^{\mathrm{i}n\varphi}P(\varphi,t;\mathrm{i})$. Now, let us multiply
(5) with $e^{\mathrm{i}n\varphi}g(\xi)$ and integrate in $\xi$ and $\varphi$.
The only additional integral to be calculated (again by virtue of the residue
method) is

$\int\limits_{0}^{2\pi}d\varphi\int\limits_{-\infty}^{\infty}d\xi\;\xi
e^{\mathrm{i}n\varphi}P\big{\rvert}_{\xi}g(\xi)=\mathrm{i}\int\limits_{0}^{2\pi}d\varphi\,e^{\mathrm{i}n\varphi}P\big{\rvert}_{\mathrm{i}}=\mathrm{i}Z_{n}\;.$

This yields the system of equations

$\dot{Z}_{n}=\mathrm{i}n(\omega+\mathrm{i}\gamma)Z_{n}+nhZ_{n-1}-nh^{*}Z_{n+1}\;,$

which coincides with (4).

We end this section with two remarks. First, the validity of Eqs. (4) for
Cauchy independent noises is unconditional, while for the Cauchy distributed
constant additions to the frequency, an extra assumption of analyticity has to
be adopted; the validity of this assumption is commonly assumed in the OA
theory and its applications, but one can construct distributions of the phase
which, at least during some time interval, violate this assumption Pikovsky
and Rosenblum (2011a).

The second remark is that if one has both Cauchy-distributed constant and
noisy additions to the frequency, with intensities $\gamma_{1}$ and
$\gamma_{2}$, then one can use Eqs. (4) with the total intensity
$\gamma=\gamma_{1}+\gamma_{2}$.

##  III Generating functions and finite-dimensional reductions of the dynamics

###  III.1 Ordinary and exponential generating functions

In our treatment of the infinite system (4) we will make use of generating
functions, which are formal power series. We will use both the ordinary
generating function (OGF) 111The sum defining ordinary generating functions is
typically considered from index $n=0$, but in our context it is convenient to
start from $n=1$; note that we always have $f_{0}=1$ for normalization
reasons.

$\mathcal{F}(k)=\sum_{n=1}^{\infty}f_{n}k^{n}\;,$

and the exponential generating function (EGF)

$\mathsf{F}(k)=\sum_{n=0}^{\infty}f_{n}\frac{k^{n}}{n!}\;.$

There is no simple relation between these functions for the same sequence
$\\{f_{n}\\}$, and in different situations we will use different generating
functions.

###  III.2 Finite-dimensional reductions of the infinite system for circular
moments

In this section we briefly introduce the finite-dimensional reductions
described in our recent letter Cestnik and Pikovsky (2022). First, we
characterize the state with complex order parameters $Z_{n}$ by introducing
the complex-valued EGF (because $|Z_{n}|\leq 1$, this series converges for all
$k$):

$\mathsf{Z}(k,t)=\sum\limits_{n=0}^{\infty}Z_{n}(t)\frac{k^{n}}{n!}\;.$ (7)

Then the dynamics (4) are recast as a single PDE (see Appendix A for the
derivation), which in contrast to Eq. (2) is generally complex (prime denotes
derivative with respect to $k$):

$\dot{\mathsf{Z}}=(\mathrm{i}\omega-\gamma)k\mathsf{Z}^{\prime}+hk\mathsf{Z}-h^{*}k\mathsf{Z}^{\prime\prime}\;.$
(8)

The normalization condition is $\mathsf{Z}(0,t)=1$. The structure of this
equation allows for a particular solution with the exponential ansatz
$\mathsf{Z}(k,t)=e^{kQ(t)}$, revealing a single ODE for the complex variable
$Q(t)$:

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

This is commonly known as the Ott-Antonsen ansatz Ott and Antonsen (2008),
which reveals a two-dimensional invariant manifold in the infinite system (4).
In this case, higher circular moments are powers of the first one:
$Z_{n}=Q^{n}$. The distribution of the phases is the wrapped Cauchy
distribution (a.k.a. Poisson kernel):

$P(\varphi,t)=\frac{1}{2\pi}\frac{1-|Q|^{2}}{|1-Qe^{-\mathrm{i}\varphi}|^{2}}\;.$
(10)

We recently generalized this solution with an ansatz allowing for an
additional function Cestnik and Pikovsky (2022):

$\mathsf{Z}(k,t)=e^{kQ(t)}\mathsf{B}(k,t)\;,$ (11)

in which case we obtain, in addition to (9), another PDE for the newly
introduced function $\mathsf{B}(k,t)$:

$\dot{\mathsf{B}}=\big{(}\mathrm{i}\omega-\gamma-2h^{*}Q\big{)}k\mathsf{B}^{\prime}-h^{*}k\mathsf{B}^{\prime\prime}\;.$
(12)

Although at first glance this equation is similar to Eq. (8), it does not
contain a term without $k$-derivative of $\mathsf{B}(k,t)$, and thus allows
for a more general dimensionality reduction. Namely, we expand the function
$\mathsf{B}(k)$ as an EGF (we will see below that coefficients $\beta_{n}$
grow not faster than $\sim\text{const}^{n}$, thus this series converges for
all $k$):

$\mathsf{B}(k,t)=\sum\limits_{n=0}^{\infty}\beta_{n}(t)\frac{k^{n}}{n!}\;,$
(13)

($\beta_{0}\equiv 1$ due to normalization), thus introducing new dynamical
variables $\beta_{n}(t)$, that describe the dynamics with an infinite set of
ODEs (plus one ODE (9) for $Q$):

$\displaystyle\dot{Q}$
$\displaystyle=(\mathrm{i}\omega-\gamma)Q+h-h^{*}Q^{2}\;\;,$ (14a)
$\displaystyle\frac{1}{n}\dot{\beta}_{n}$
$\displaystyle=(\mathrm{i}\omega-\gamma-2h^{*}Q)\beta_{n}-h^{*}\beta_{n+1}\;,\quad
n\geq 1\;,$ (14b)

(see Appendix A for the relation of (14b) to (12)). Notice how the right-hand
side of (14b) only contains terms proportional to $\beta_{n}$ and
$\beta_{n+1}$, but no term with $\beta_{n-1}$ is present. This means that if
the system is truncated at a finite number $N$ of variables $\beta_{n}$ (i.e.
assuming that all higher terms vanish identically: $\beta_{n\geq N}=0$), the
dynamics is exactly described by the first $N$ equations of system (14b) for
all times. These truncations represent dynamically invariant finite-
dimensional manifolds. The $\beta_{n}$ variables relate to the Kuramoto-Daido
order parameters $Z_{n}$ via a modified binomial transform (cf. Refs.
Haukkanen, 1993; Prodinger, 1994):

$\begin{split}Z_{n}(t)&=\sum\limits_{m=0}^{n}\binom{n}{m}\beta_{m}(t)\big{[}Q(t)\big{]}^{n-m}\;,\\\
\beta_{n}(t)&=\sum\limits_{m=0}^{n}\binom{n}{m}Z_{m}(t)\big{[}-Q(t)\big{]}^{n-m}\;.\end{split}$
(15)

For example, the first three order parameters are expressed with the newly
introduced variables as

$\displaystyle Z_{1}$ $\displaystyle=Q+\beta_{1}\;,\qquad
Z_{2}=Q^{2}+2Q\beta_{1}+\beta_{2}\;,$ $\displaystyle\quad
Z_{3}=Q^{3}+3Q^{2}\beta_{1}+3Q\beta_{2}+\beta_{3}\;.$

##  IV Reduction of the dynamics to three complex variables

As outlined in the previous section III, there are many finite-dimensional
invariant manifolds (with a finite number of additional variables $\beta_{n}$)
beyond the OA two-dimensional manifold (which corresponds to vanishing
$\beta_{n}$ for $n\geq 1$). However, as already mentioned in Ref. Cestnik and
Pikovsky, 2022, it is not excluded that different $\beta_{n}$ could be
dependent. Below we show that this is indeed the case, and the dynamics of the
whole (even infinite) hierarchy of variables $\beta_{n}(t)$ can be reduced to
two complex equations.

###  IV.1 Six-dimensional reduction

We now introduce two new complex variables $y(t),s(t)$ and new dynamical
equations:

$\displaystyle\dot{Q}$
$\displaystyle=(\mathrm{i}\omega-\gamma)Q+h-h^{*}Q^{2}\;,$ (16a)
$\displaystyle\dot{y}$ $\displaystyle=(\mathrm{i}\omega-\gamma-2h^{*}Q)y\;,$
(16b) $\displaystyle\dot{s}$ $\displaystyle=h^{*}y\;.$ (16c)

Our goal below is to demonstrate, that these equations are equivalent to the
infinite system (14) and therefore to the original system (4). At this point
we would like to mention, that Eqs. (16) look like a skew system: variable $Q$
appears on the r.h.s. of (16b), and variable $y$ appears on the r.h.s. of
(16c), but variables $y,s$ do not appear on the r.h.s. of (16a). However, in
most applications one describes a population with global coupling, where
$\omega,h$ depend on the order parameters $Z_{n}$, and thus on all dynamical
variables $Q,y,s$ (cf. Eq. (23) below).

To show how this system represents dynamics (14), we first introduce
additional auxiliary variables $\alpha_{n}(t)$ by transforming $\beta_{n}(t)$:

$\beta_{n}(t)=y^{n}(t)\alpha_{n}(t)\;.$ (17)

We take the time derivative of this relation and divide both sides by
$n\beta_{n}$

$\frac{\frac{1}{n}\dot{\beta}_{n}}{\beta_{n}}=\frac{\dot{y}}{y}+\frac{\frac{1}{n}\dot{\alpha}_{n}}{\alpha_{n}}\;,$

and then insert the dynamics of $\beta_{n}$ (14b) and $y$ (16b) :

$(\mathrm{i}\omega-\gamma-2h^{*}Q)-h^{*}\frac{\beta_{n+1}}{\beta_{n}}=(\mathrm{i}\omega-\gamma-2h^{*}Q)+\frac{\frac{1}{n}\dot{\alpha}_{n}}{\alpha_{n}}\;.$

Notice how the majority of the terms cancel, including all the effects of
frequency $\omega(t)$ and noise $\gamma$. As a result, the dynamics of the
variables $\alpha_{n}(t)$ simplifies to:

$\frac{1}{n}\dot{\alpha}_{n}=-h^{*}y\;\alpha_{n+1}\;.$ (18)

Now let us introduce the OGF of the variables $\alpha_{n}(t)$:

$\mathcal{A}(k,t)=\sum\limits_{n=1}^{\infty}\alpha_{n}(t)k^{n}\;,$ (19)

and express the dynamics (18) in terms of this OGF (see Appendix A for the
derivation):

$\dot{\mathcal{A}}=-h^{*}y\Big{[}\mathcal{A}^{\prime}-\frac{1}{k}\mathcal{A}\Big{]}\;.$
(20)

Next we introduce yet another set of variables $\mu_{n}$. This time we launch
them not directly, but via an expression of the corresponding OGF
$\mathcal{M}(k,t)=\sum\limits_{n=1}^{\infty}\mu_{n}k^{n}$ in terms of
$\mathcal{A}(k,t)$:

$\frac{\mathcal{M}(k)}{k}=\frac{\mathcal{A}(k+s)}{k+s}$ (21)

(where $s$ is the variable in (16c)). By taking the time derivative of Eq.
(21), we obtain:

$\frac{\dot{\mathcal{M}}(k)}{k}=\frac{\dot{\mathcal{A}}(k+s)}{k+s}+\dot{s}\Big{[}\frac{\mathcal{A}^{\prime}(k+s)}{k+s}-\frac{\mathcal{A}(k+s)}{(k+s)^{2}}\Big{]}\;.$

Now we insert the dynamics of $s$ according to Eq. (16c), as well as the
dynamics of $\mathcal{A}$ according to Eq. (20) and behold, the right-hand
side of this relation is zero. This means that the OGF $\mathcal{M}$ and the
corresponding variables $\mu_{n}$ are constant in time: $\dot{\mu}_{n}=0$,
$\dot{\mathcal{M}}(k)=0$. In other words, the variables $\mu_{n}$ are
integrals of motion. This completes the proof that Equations (16) are
equivalent to Equations (14).

The radius of convergence of the constant function $\mathcal{M}(k)$ is
determined by the asymptotic behavior of the coefficients $\mu_{n}$. Below we
show that these coefficients initially coincide with the order parameters
$Z_{n}(0)$. Since they are bounded $|Z_{n}|\leq 1$, the radius of convergence
of $\mathcal{M}$ is at least one. During dynamical evolution, the relevant
argument of function $\mathcal{M}$ is $s(t)$, and in all our simulations we
never observed it getting larger than one in absolute value: $|s(t)|<1$.

It is instructive to rephrase the relation (21), formulated above in terms of
OGFs, to the level of the variables (where it corresponds to a modified
binomial transform, see Appendix B for the derivation):

$\begin{split}\mu_{n}&=\sum\limits_{m=n}^{\infty}\binom{m-1}{n-1}\alpha_{m}(t)\big{[}s(t)\big{]}^{m-n}\;,\\\
\alpha_{n}(t)&=\sum\limits_{m=n}^{\infty}\binom{m-1}{n-1}\mu_{m}\big{[}-s(t)\big{]}^{m-n}\;.\end{split}$
(22)

Notice that we do not write the time argument of $\mu_{n}$ because these
quantities are constants.

Using this relation, as well as how the variables $\alpha_{n}(t)$ relate to
$\beta_{n}(t)$ (17), and then how $\beta_{n}(t)$ relate to the order
parameters $Z_{n}(t)$ (15), we can express the order parameters in terms of
the constant function $\mathcal{M}(k)$ and the three dynamic variables
$Q(t),y(t),s(t)$ (see Appendix C for the derivation) for all times (here we
omit the time dependence in notation for convenience):

$Z_{n}=Q^{n}-\sum\limits_{m=1}^{n}\binom{n}{m}Q^{n-m}y^{m}\sum\limits_{d=0}^{m-1}\frac{s^{d-m}}{d!}\mathcal{M}^{(d)}(-s)\;,$
(23)

where $\mathcal{M}^{(d)}$ denotes the $d^{\text{th}}$ derivative of
$\mathcal{M}$ with respect to $k$. In particular, the first circular moment
(the Kuramoto order parameter) is expressed as:

$Z_{1}=Q-y\frac{\mathcal{M}(-s)}{s}\;.$ (24)

Notice how at $s=0$ one has to take the limit $\lim\limits_{\varepsilon\to
0}\frac{\mathcal{M}(-\varepsilon)}{\varepsilon}=-\mu_{1}$. For all moment
expressions (23) we have to consider similar limits:
$\lim\limits_{\varepsilon\to
0}\sum\limits_{d=0}^{m-1}\frac{\varepsilon^{d-m}}{d!}\mathcal{M}^{(d)}(-\varepsilon)=-\mu_{m}$.
When performing numerical integration, one thus requires an expansion of the
above expression (23) for small $Q$ and small $s$:

$Z_{n}=y^{n}\left[\mu_{n}+n\mu_{n-1}\frac{Q}{y}-n\mu_{n+1}s\right]+\mathcal{O}(Q^{2},Qs,s^{2})\;,$
(25)

we remind that $\mu_{0}\equiv 1$ for normalization reasons. The need for
expansion (25) can also be avoided by considering a different definition of
the constant function (21), see Appendix E \- it also simplifies the
expression for moments (23).

###  IV.2 Initial conditions

The new set of variables $Q,y,s,\mu_{n}$ is not uniquely determined by the
initial order parameters $Z_{n}$, and in this section we discuss possible
variants of determining them. Different choices for the initial conditions of
$Q,y,s,\mu_{n}$ define the constant function $\mathcal{M}(k)$ differently. We
first illustrate this with the simplest example of the OA manifold.

####  IV.2.1 Different choices of variables for OA initial conditions

As discussed above, on the OA manifold the EGF reads
$\mathsf{Z}(k,t)=\exp[kZ(t)]$ and the order parameter $Z$ obeys
$\dot{Z}=(\mathrm{i}\omega-\gamma)Z+h-h^{*}Z^{2}$. Suppose, having a set of
initial moments $Z_{n}(0)=Z^{n}(0)$, we want to introduce new variables
$Q,\beta_{n}$ according to (11). One can immediately see that a set
$\beta_{n}(0)=\beta^{n}(0)$ is admissible if $Q(0)+\beta(0)=Z(0)$, in this
case $\mathsf{B}(k,0)=\exp[k\beta(0)]$. Furthermore, relation (17) allows for
different choices of $y(0)$ and $\alpha_{n}(0)$. For any choice of $y(0)$, we
obtain $\alpha_{n}(0)=\alpha^{n}$, with $\alpha=\beta(0)/y(0)$. It is easy to
see that independently of the choice of $\beta(0)$ and $y(0)$, the dynamics of
the order parameters is the same. Indeed, in this case $\mathcal{M}(k)=\alpha
k/(1-\alpha k)$ and the main order parameter according to (24) is
$Z_{1}=Q+\alpha y/(1+\alpha s)$. Calculation of the derivative $\dot{Z}_{1}$
from the general dynamical equations (16) yields the correct equation
$\dot{Z}_{1}=(\mathrm{i}\omega-\gamma)Z_{1}+h-h^{*}Z_{1}^{2}$, i.e. the system
remains on the OA manifold. In this specific case of pure OA dynamics, it is
natural to consider $\alpha=0$ such that $\mathcal{M}=0$ and the only relevant
equation is the OA equation (16a) Pikovsky and Rosenblum (2008, 2011b).

####  IV.2.2 Variant A: A simple choice of initial variables

Here we present possibly the simplest choice of initial conditions for
$Q,y,s,\mu_{n}$. We initially set $Q$ and $s$ to zero and set $y$ to 1, so
that all the variable sets $\mu_{n},\alpha_{n},\beta_{n}$ and $Z_{n}$
coincide:

$\begin{split}&Q(0)=0\;,\\\ &y(0)=1\;,\\\ &s(0)=0\;,\\\
&\mu_{n}=\alpha_{n}(0)=\beta_{n}(0)=Z_{n}(0)\;.\end{split}$ (26)

Then the constant function can simply be determined by the initial order
parameters:

$\mathcal{M}(k)=\sum\limits_{n=1}^{\infty}Z_{n}(t=0)\,k^{n}\;.$ (27)

It is instructive to express this function in terms of the initial
distribution of the phases

$\mathcal{M}(k)=\int_{0}^{2\pi}d\varphi\;P(\varphi,0)\frac{ke^{\mathrm{i}\varphi}}{1-ke^{\mathrm{i}\varphi}}\;.$
(28)

As discussed above, this is not the only possible choice of initial conditions
and in some cases might not be optimal. Notice that for this choice, the OA
manifold corresponds to the constants being powers of $Z_{1}(0)$:
$\mu_{n}=Z^{n}_{1}(0)$, so they do not vanish, as is commonly considered
Pikovsky and Rosenblum (2008, 2011b). In terms of dynamics, the two
descriptions are equivalent.

Next we list specific functions $\mathcal{M}(k)$ for some examples of the
initial states in this variant A of the initial conditions (26). We stress
here, that a representation via an initial distribution density is valid for
the interpretation of the system with identical oscillators under Cauchy white
noise (section II.1). For the case of non-identical oscillators with a
distribution of frequencies (section II.2), one should operate with the order
parameters directly.

  * •

A uniform distribution $P(\varphi,0)=\frac{1}{2\pi}$ corresponds to
$\mathcal{M}(k)=0$. Here all the moments vanish, which is a trivial invariant
state of the dynamics (2).

  * •

A delta distribution of the phases $P(\varphi,0)=\delta(\varphi-\varphi_{0})$
corresponds to $Z_{n}=\exp(\mathrm{i}n\varphi_{0})$ and thus
$\mathcal{M}(k)=\frac{e^{\mathrm{i}\varphi_{0}}\,k}{1-e^{\mathrm{i}\varphi_{0}}\,k}$.

  * •

A wrapped Cauchy distribution
$P(\varphi,0)=\frac{1}{2\pi}\frac{1-|\mu|^{2}}{|1-\mu
e^{-\mathrm{i}\varphi}|^{2}}$ with a complex parameter $\mu\in\mathbb{C}$,
corresponds to $\mathcal{M}(k)=\frac{\mu k}{1-\mu k}$. The moments are powers
of the parameter $\mu$: $\mu_{n}=\mu^{n}$, which means this state is on the OA
manifold Ott and Antonsen (2008).

  * •

A Kato-Jones distribution Kato and Jones (2015)
$P(\varphi,0)=\frac{1}{2\pi}(1+c\frac{\mu e^{-\mathrm{i}\varphi}}{1-\mu
e^{-\mathrm{i}\varphi}}+\text{c.c.})$ with complex parameters
$c,\mu\in\mathbb{C}$, corresponds to $\mathcal{M}(k)=c\frac{\mu k}{1-\mu k}$.
Its moments are described by a power series multiplied with a complex
constant: $\mu_{n}=c\,\mu^{n}$, making it a skewed/asymmetric generalization
of the wrapped Cauchy distribution.

  * •

Distributions with a finite number of moments
$P(\varphi,0)=\frac{1}{2\pi}(1+\sum\limits_{n=1}^{N}\mu_{n}e^{-\mathrm{i}n\varphi}+\text{c.c.})$
correspond to polynomials $\mathcal{M}(k)=\sum\limits_{n=1}^{N}\mu_{n}k^{n}$.

  * •

Distributions with binomial moments $\mu_{n}=\binom{n}{m}\mu^{n}$ for $m\geq
1$ correspond to rational functions $\mathcal{M}(k)=\frac{(\mu k)^{m}}{(1-\mu
k)^{m+1}}$. Their phase density reads:
$P(\varphi,0)=\frac{1}{2\pi}(1+\frac{\mu^{m}e^{-\mathrm{i}m\varphi}}{(1-\mu
e^{-\mathrm{i}\varphi})^{m+1}}+\text{c.c.})$.

  * •

A half-uniform distribution

$P(\varphi,0)=\begin{cases}\frac{1}{\pi}&\text{if}\
\varphi\in(\varphi_{0},\varphi_{0}+\pi)\;,\\\ 0&\text{else}\;,\end{cases}$
(29)

corresponds to
$\mathcal{M}(k)=\frac{2\mathrm{i}}{\pi}\operatorname{arctanh}(e^{\mathrm{i}\varphi_{0}}k)$.
Here odd moments are fractions:
$\mu_{2n-1}=\frac{2\mathrm{i}}{\pi}\frac{\exp(\mathrm{i}(2n-1)\varphi_{0})}{2n-1}$,
and even ones are equal to zero: $\mu_{2n}=0$, $n\geq 1$.

  * •

A sawtooth distribution
$P(\varphi,0)=\frac{1}{\pi}(1-\frac{\varphi-\varphi_{0}}{\pi}+\lfloor\frac{\varphi-\varphi_{0}}{\pi}\rfloor)$
corresponds to
$\mathcal{M}(k)=-\frac{\mathrm{i}}{\pi}\log(1-e^{\mathrm{i}2\varphi_{0}}k^{2})$.
Here even moments are fractions:
$\mu_{2n}=\frac{2\mathrm{i}}{\pi}\frac{\exp(\mathrm{i}2n\varphi_{0})}{2n}$,
and odd ones are equal to zero: $\mu_{2n-1}=0$, $n\geq 1$.

In practice, the phase density can only be observed empirically, and therefore
it may not be clear how to choose $\mathcal{M}(k)$. In such situations it can
always be approximated with a finite truncation of its Taylor series:
$\mathcal{M}(k)\approx\sum\limits_{n=1}^{N}Z_{n}(0)k^{n}$, one only needs to
estimate the first few initial moments $Z_{n}(0)$. This just corresponds to
approximating the initial state with a finite Fourier representation. In
Appendix F we present a numerical example where we test the convergence of
such approximations.

We end this subsection with the following remark: if the initial distribution
of the phases is a weighted sum of “elementary” distributions
$P(\varphi,0)=\sum_{m}c_{m}P_{m}(\varphi,0)$ with real weights
$c_{m}\in\mathbb{R}$ that add up to 1: $\sum_{m}c_{m}=1$ (additionally one
should ensure that $P(\varphi,0)\geq 0$), then the constant generating
function is the weighted sum of the corresponding “elementary” generating
functions $\mathcal{M}(k)=\sum_{m}c_{m}\mathcal{M}_{m}(k)$. In particular, in
Refs. Engelbrecht and Mirollo, 2020; Ichiki and Okumura, 2020 a superposition
of several wrapped Cauchy distributions has been considered as an initial
state; in terms of the approach above this corresponds to
$\mathcal{M}(k)=\sum_{m}c_{m}\frac{\varkappa_{m}k}{1-\varkappa_{m}k}$, where
complex parameters $\varkappa_{m}$ characterize partial distributions.

####  IV.2.3 Variant B: Initial conditions based on the OA manifold

Often initial states that are close to the OA manifold are of interest.
Suppose that the order parameters are well described as powers of a complex
constant, with minor perturbations:

$Z_{n}(0)=R^{n}+\varepsilon_{n}\;.$ (30)

In this case a different initial condition appears natural:

$\begin{split}&Q(0)=R\;,\\\ &y(0)=1\;,\\\ &s(0)=0\;,\\\
&\mu_{n}=\alpha_{n}(0)=\beta_{n}(0)=\sum\limits_{m=1}^{n}\binom{n}{m}\varepsilon_{m}(-R)^{n-m}\;,\end{split}$
(31)

and the constant function $\mathcal{M}(k)$ is expressed as:

$\mathcal{M}(k)=\sum\limits_{n=1}^{\infty}k^{n}\sum\limits_{m=1}^{n}\binom{n}{m}\varepsilon_{m}(-R)^{n-m}\;.$
(32)

This function is small if values of $\varepsilon_{n}$ are small. Notice
however, that this “perturbation” approach is actually global, because
smallness of $\varepsilon_{n}$ is not needed, and $\mathcal{M}(k)$ need not be
small for this description to be valid. Note that the definition of function
$\mathcal{M}(k)$ depends on the choice of initial conditions (cf. (27) and
(32)), so the list of specific $\mathcal{M}(k)$ functions in Section IV.2.2
does not apply here. The moments are still described with Eq. (23) but for
numerical integration one needs to expand them beyond (25) for only small $s$
(such expansions can be avoided by considering an alternative constant
function to $\mathcal{M}$, see Appendix E).

A simple specific example where the perturbation to the OA manifold in (30),
is one where only the first harmonic term is perturbed: $\varepsilon_{n}=0$
for $n>1$. In this case $\mathcal{M}(k)=\varepsilon_{1}\frac{k}{(1+Rk)^{2}}$
and so the first moment is given by
$Z_{1}=Q+\varepsilon_{1}\frac{y}{(1-Rs)^{2}}$. The dynamics follow (16).

We mention here that in some cases, an extension of the set of variables might
be appropriate. As an example, we show in Appendix D a possibility to describe
an initial state similar to (30) with a system of four complex variables.

###  IV.3 Finite ensemble numerical comparison

Figure 1: Comparison of the finite ensemble dynamics of $10\,000$ Josephson
phases (33) with the low-dimensional dynamics in the thermodynamic limit (16).
In panel $(a)$ the first two circular moments are compared in the complex
plane, and in panel $(b)$ their absolute values are compared as functions of
time. Solid lines: simulation of a finite ensemble; dotted lines (which
practically overlap with solid lines): solution of exact equations in the
thermodynamic limit.

Here we numerically compare the derived low-dimensional dynamics (16), which
is exact in the thermodynamic limit, with a simulation of a finite ensemble.
We take the example already explored in Ref. Cestnik and Pikovsky, 2022: an
array of overdamped noisy Josephson junctions coupled via a resistive load
Watanabe and Strogatz (1994). The equations for the Josephson phases read

$\dot{\varphi}_{j}=1+a\sin(\varphi_{j})+\frac{\varepsilon}{N}\sum_{n=1}^{N}\sin(\varphi_{n})+\gamma\xi_{j}(t)\;.$
(33)

In terms of the basic model (1), this corresponds to the choice of $\omega,h$:
$\omega=1+\varepsilon\,\text{Im}[Z_{1}]$ and $h=-\frac{a}{2}$. We consider
parameters: $a=-0.7$, $\varepsilon=1.5$, and noise strength $\gamma=0.02$. To
highlight the advantage of the new derivation, we consider the initial
distribution of phases to be half-uniform (29), thus starting far from the OA
manifold, and also far from initial states that are easy to represent in terms
of the $Q,\beta_{n}$ hierarchy Cestnik and Pikovsky (2022). For the finite
ensemble we consider $10\,000$ phases, which are initially randomly sampled in
the interval $[0,\pi)$. For our exact reduction in the thermodynamic limit, we
take initial conditions from variant A (26): $Q(0)=s(0)=1-y(0)=0$, and
therefore the constant function
$\mathcal{M}(k)=\frac{2\mathrm{i}}{\pi}\operatorname{arctanh}(k)$. The
comparison of the trajectories in Fig. 1 shows a very close match, and we
expect it to be even closer for larger ensembles. The discrepancies appear to
be only due to finite size effects, we stress that in the limit of an
infinitely large ensemble our reduction (16) is exact.

##  V Global stability of the OA manifold

Stability of the OA manifold has been discussed in Refs. Ott and Antonsen,
2009; Pietras and Daffertshofer, 2016; Engelbrecht and Mirollo, 2020. Here we
demonstrate how these results are reproduced in our approach. To show the
attractiveness of the OA manifold it is enough to demonstrate that the
variable $y$ tends to zero $y\to 0$. Indeed, for $y=0$ we have from (17)
$\beta_{n}=0$, $n\geq 1$, and from (11) it follows that the solution is on the
OA manifold.

Let us introduce two new variables $Y,S$ according to relations

$y=(1-|Q|^{2})Y\;,\qquad s=Q^{*}S\;.$ (34)

The equations for these variables read

$\displaystyle\dot{Y}$
$\displaystyle=\left(\mathrm{i}\omega+hQ^{*}-h^{*}Q-\gamma\frac{1+|Q|^{2}}{1-|Q|^{2}}\right)Y\;,$
(35a) $\displaystyle\dot{S}$
$\displaystyle=\Big{(}\mathrm{i}\omega+hQ^{*}-h^{*}Q+\gamma\Big{)}S+\left(h^{*}Q-\frac{h^{*}}{Q^{*}}\right)(S-Y)\;.$
(35b)

Here we focus on the equation for $Y$, and will use Eq. (35b) in Section VI
below. From (35a) we obtain the following evolution of $|Y|$:

$\frac{d|Y|}{dt}\frac{1}{|Y|}=-\gamma\frac{1+|Q|^{2}}{1-|Q|^{2}}\leq-\gamma\;.$
(36)

The latter inequality on the r.h.s. follows from the property
$0\leq|Q|^{2}<1$. Indeed, the equation for the evolution of $|Q|^{2}$ reads

$\frac{d}{dt}|Q|^{2}=-2\gamma|Q|^{2}+(hQ^{*}+h^{*}Q)(1-|Q|^{2})$

and at $|Q|=1$ the derivative $\frac{d}{dt}|Q|^{2}=-2\gamma$ is negative for
$\gamma>0$, thus this boundary is not reachable.

Integrating inequality (36) yields

$|Y(t)|\leq|Y(0)|\,e^{-\gamma t}\;,$ (37)

and consequently, since $|y|\leq|Y|$:

$|y(t)|\leq\frac{|y(0)|}{1-|Q(0)|^{2}}\,e^{-\gamma t}\;,$ (38)

which means that $|Y|$ and $|y|$ vanish exponentially fast in time. This
proves the attractiveness of the OA manifold for system (4) for $\gamma>0$.

Figure 2: Time series of the coupled Josephson junctions array (33) for three
different noise strengths; $\gamma=0$ depicted with red, $\gamma=10^{-4}$
depicted with blue and $\gamma=2\times 10^{-4}$ depicted with green. Panel
$(a)$: the norm of the first order parameter $Z_{1}$; panel $(b)$: the norm of
variable $y$. The initial conditions for all three cases are the same:
$Q(0)=s(0)=1-y(0)=0$ (variant A) and constant function $\mathcal{M}(k)=0.4k$,
which means the initial distribution has only one harmonic. In the case of no
noise, chaotic dynamics is observed, while in the noisy case there is an
initial chaotic stage which is followed by a clear exponential decay of the
$y$ variable. In panel $(b)$ the upper bound (38) for all three cases is shown
with a dashed black line.

###  V.1 Example

Here we illustrate stability of the OA manifold numerically. We take the
already explored example of Josephson junctions (33), with the same parameters
$a=-0.7$, $\varepsilon=1.5$. In Ref. Cestnik and Pikovsky, 2022 it was
demonstrated that for $\gamma=0$, the dynamics outside of the OA manifold is
chaotic. In Fig. 2 we again demonstrate chaotic behavior in this system for
$\gamma=0$, and a transition to regular dynamics for $\gamma=10^{-4}$ and
$\gamma=2\times 10^{-4}$. The exponential decay of $|y|$, which is bounded by
(38) is evident at large times.

##  VI Noise-free case and a relation to the Watanabe-Strogatz theory

Watanabe and Strogatz Watanabe and Strogatz (1993, 1994) demonstrated that a
population of identical noiseless oscillators can be reduced to three real
dynamical variables plus constants of motion. To see that this case is
included in our theory, let us consider identical oscillators and no noise,
thus taking $\gamma=0$.

It is instructive to start with Eqs. (35). It is easy to see, that for
$\gamma=0$, the manifold $S=Y$ is invariant:

$\frac{d}{dt}(S-Y)=\left[\mathrm{i}\omega+hQ^{*}-\frac{h^{*}}{Q^{*}}\right](S-Y)\;.$

Using the definition of variables $Y,S$ (34), the manifold is described in
terms of $y$ and $s$ as:

$Y=\frac{y}{1-|Q|^{2}}=\frac{s}{Q^{*}}=S\;.$ (39)

Since we initially set $Q(0)=s(0)=1-y(0)=0$ for arbitrary states, one can
always set initial conditions on this manifold such that $Y(0)=S(0)=1$. Thus,
two equations (35) reduce to one. Moreover, because of (36), for $\gamma=0$
the variable $Y$ remains on the unit circle for all times, and we can
introduce an angle variable $\theta(t)=\text{arg}(Y(t))$. Its evolution
follows from (35):

$\dot{\theta}=\omega-\mathrm{i}(hQ^{*}-h^{*}Q)\;.$ (40)

This angle variable $\theta$ just corresponds to the WS angle variable, while
$Q$ is the WS order parameter Watanabe and Strogatz (1993, 1994). The two
equations (16a) (with $\gamma=0$) and (40) then represent the exact evolution.
We explicitly derive the equivalence with the WS approach in Appendix G.

In the WS theory, the relation between the original phases $\varphi_{j}(t)$
and constant phases $\psi_{j}$ is given by the Möbius transform Watanabe and
Strogatz (1994); Marvel, Mirollo, and Strogatz (2009):

$e^{\mathrm{i}\varphi_{j}}=\frac{e^{\mathrm{i}(\psi_{j}+\theta)}+Q}{1+Q^{*}e^{\mathrm{i}(\psi_{j}+\theta)}}\;,\qquad
e^{\mathrm{i}(\psi_{j}+\theta)}=\frac{e^{\mathrm{i}\varphi_{j}}-Q}{1-Q^{*}e^{\mathrm{i}\varphi_{j}}}\;.$
(41)

Our choice of the initial conditions $Q(0)=\theta(0)=0$ (it corresponds to the
“identity conversion” in terms of WS, cf. Eq. (5.10) in Ref. Watanabe and
Strogatz, 1994) means that $\psi_{j}=\varphi_{j}(0)$. Thus, because the
integrals $\mu_{n}$ are defined as $\mu_{n}=Z_{n}(0)$, these quantities are
the circular moments of the transformed constant phase variables in the WS
approach $\mu_{n}=\langle e^{\mathrm{i}n\psi}\rangle$.

Watanabe and Strogatz have shown that these transformations are also valid for
a finite number of oscillators, but this case is not covered by our approach.
We mention here, that in the WS formalism there is also a freedom in choosing
the order parameter $Q$ and the phase variable $\theta$; this freedom is
similar to the one discussed in Section IV.2 above.

##  VII Lyapunov spectrum

Our theory describes the dynamics outside of the OA manifold, and is thus
suitable for consideration of small perturbations transversal to this
manifold. Such perturbations define the Lyapunov spectrum of the dynamics,
together with the perturbations tangential to this manifold. The system of
equations (14) is most suitable for this analysis. The OA manifold corresponds
to vanishing $\beta_{n}$, therefore Eqs. (14b) define the transversal
perturbations. Since these equations are a skew system, each $\beta_{n}$
defines two Lyapunov exponents (because $\beta_{n}$ are complex). One can
straightforwardly derive from (14b), omitting the skew term $\sim\beta_{n+1}$
on the r.h.s., the averaged evolution for the magnitude of a perturbation:

$\frac{1}{2n}\left\langle\frac{d}{dt}\ln|\beta_{n}|^{2}\right\rangle=\left\langle-\gamma-h^{*}Q-hQ^{*}\right\rangle=\Lambda\;.$

Thus, the Lyapunov spectrum consists of the exponents within the OA manifold
(which are calculated using linearised Eq. (9)), and of doubly degenerated
values $n\Lambda$, $n=1,2,3,\ldots$.

##  VIII Response of the Ott-Antonsen regime to a resetting

As has been already discussed in the literature Ott and Antonsen (2009);
Pietras and Daffertshofer (2016); Engelbrecht and Mirollo (2020) and in
Section V, in the system of equations (16) the OA manifold is attracting if
$\gamma>0$ (at least in the weak sense, but because we follow only the moments
of the phase distribution, such an attraction is enough). In terms of
variables $Q,y,s$ with a nontrivial constant function $\mathcal{M}(k)$, this
corresponds to $y\to 0$ as $t\to\infty$ (see Section V). For the conservative
case $\gamma=0$, see Section VI above.

The approach above allows for calculating the evolution from an arbitrary
state to the OA manifold via solutions of (16). One can reformulate such a
problem as a resetting one: One starts with the dynamics on the OA manifold;
then an instant “resetting” to a state outside of this manifold is performed.
The evolution of (16) then shows what will be the final state after re-
attraction to the OA manifold. A particular question of interest depends on
the type of the attractors on the OA manifold. If there is only one global
attractor, then the trajectory returns to it. If this attractor is periodic or
quasiperiodic, the returning trajectory will be phase shifted with respect to
the unperturbed one (in the quasiperiodic case one expects phase shifts in
every direction of independent oscillations). Here one speaks about a phase
resetting or a phase response curve (PRC) Canavier (2006); Smeal, Ermentrout,
and White (2010). For a chaotic global attractor, generally one does not
expect a resetting to have a drastic effect (although for strange attractors
with well-defined phase variables a phase resetting similar to the periodic
case can be defined Schwabedal  _et al._ (2012); it can lead to phase
synchronization of chaos if periodically repeated Pikovsky  _et al._ (1997)).
In the case of multistability, the most drastic effect of resetting would be a
jump to another basin of attraction, so that the final state will be another
attractor on the OA manifold (in case of multistable periodic attractors one
can additionally follow the phase response Grines, Osipov, and Pikovsky
(2018)). Below we consider several examples, for small and large resettings.

###  VIII.1 Perturbation theory in terms of an (infinitesimal) PRC

Suppose we have a state on the OA manifold with a complex order parameter $R$,
so that $\langle e^{\mathrm{i}n\varphi}\rangle=R^{n}$. Let us apply to all the
phases a transformation

$\varphi\to\varphi+\varepsilon f(\varphi)\;,$ (42)

where $f(\varphi)=\sum_{m}f_{m}e^{\mathrm{i}m\varphi}$ is a PRC function
(given by its Fourier representation) and $\varepsilon\ll 1$ is assumed to be
small. Let us calculate the circular moments just after a resetting, in order
$\varepsilon$:

$\displaystyle Z_{n}=\left\langle e^{\mathrm{i}n(\varphi+\varepsilon
f(\varphi))}\right\rangle\approx\left\langle
e^{\mathrm{i}n\varphi}(1+\mathrm{i}n\varepsilon f(\varphi))\right\rangle=$
$\displaystyle=R^{n}+\mathrm{i}n\varepsilon\sum\limits_{m=-\infty}^{\infty}f_{m}\left\langle
e^{\mathrm{i}(n+m)\varphi}\right\rangle\;.$

Since $m$ can be negative, calculation of the latter average is not a simple
expression, because

$\left\langle
e^{\mathrm{i}(n+m)\varphi}\right\rangle=\begin{cases}R^{n+m}&n+m\geq 0\;,\\\
(R^{*})^{|n+m|}&n+m<0\;,\end{cases}$

therefore we restrict ourselves to two simplest cases.

####  VIII.1.1 First harmonics resetting

In this case
$f(\varphi)=f_{1}e^{\mathrm{i}\varphi}+f_{1}^{*}e^{-\mathrm{i}\varphi}$. For
$n\geq 1$ we have $n+m\geq 0$ and therefore for both $m=\pm 1$ we can write
$\left\langle e^{\mathrm{i}(n+m)\varphi}\right\rangle=R^{n+m}$. Thus

$Z_{n}=R^{n}(1+\mathrm{i}n\varepsilon(f_{1}R+f_{1}^{*}R^{-1}))\;,\qquad n\geq
0\;.$

Calculation of the EGF yields

$\mathsf{Z}(k,0)=e^{kR}(1+\mathrm{i}\varepsilon k(f_{1}R^{2}+f_{1}^{*}))\;,$

where we used $\sum_{n}n\frac{x^{n}}{n!}=xe^{x}$.

Let us now transform to variables $Q,y,s$ and take $Q(0)=R$ (like in variant
B, Section IV.2.3). This means that the EGF $\mathsf{B}$ is

$\mathsf{B}(k,0)=1+\mathrm{i}\varepsilon k(f_{1}R^{2}+f_{1}^{*})\;.$

We come to the conclusion, that only one variable $\beta_{1}$ is non-zero, and
the system can be directly and exactly solved with variables $Q,\beta_{1}$ by
virtue of Eqs. (14); there is no need to go to the full system (16).
Alternatively, one can consider only the first two equations of system (16)
and function $\mathcal{M}(k)=\beta_{1}(0)k$ and then the third variable $s$
does not matter. If we rewrite Eqs. (14) in terms of variables
$(Z_{1},\beta_{1})$, we obtain

$\displaystyle\dot{Z}_{1}$
$\displaystyle=(\mathrm{i}\omega-\gamma)Z_{1}+h-h^{*}Z_{1}^{2}+h^{*}\beta_{1}^{2}\;,$
$\displaystyle\dot{\beta}_{1}$
$\displaystyle=(\mathrm{i}\omega-\gamma-2h^{*}(Z_{1}-\beta_{1}))\beta_{1}\;.$

One can see that the correction to the standard OA equation is
$h^{*}\beta_{1}^{2}\sim\varepsilon^{2}$. Thus, in the first order in
$\varepsilon$, inclusion of the additional variable $\beta_{1}$ is irrelevant
and the resetting is well described within the OA equation.

####  VIII.1.2 Second harmonics resetting

In this case
$f(\varphi)=f_{2}e^{\mathrm{i}2\varphi}+f_{2}^{*}e^{-\mathrm{i}2\varphi}$, and
we have $\sum_{m}f_{m}\left\langle
e^{\mathrm{i}(n+m)\varphi}\right\rangle=f_{2}\left\langle
e^{\mathrm{i}(n+2)\varphi}\right\rangle+f_{2}^{*}\left\langle
e^{\mathrm{i}(n-2)\varphi}\right\rangle$. Thus, the term with $n=1$ reads
$f_{2}R^{3}+f_{2}^{*}R^{*}$, while all the higher-order terms $n\geq 2$ can be
written in a unified way $f_{2}R^{n+2}+f_{2}^{*}R^{n-2}$. Rewriting the term
with $n=1$ as $f_{2}R^{3}+f_{2}^{*}R^{-1}+[f_{2}^{*}R^{*}-f_{2}^{*}R^{-1}]$,
we obtain

$Z_{n}=R^{n}+\mathrm{i}\varepsilon
nR^{n}[f_{2}R^{2}+f_{2}^{*}R^{-2}]+\mathrm{i}\varepsilon\delta_{n,1}[f_{2}^{*}R^{*}-f_{2}^{*}R^{-1}]\;,$

where $\delta_{n,1}$ is the Kronecker delta. This yields the following EGF

$\mathsf{Z}(k,0)=e^{kR}(1+\mathrm{i}\varepsilon
k(f_{2}R^{3}+f_{2}^{*}R^{-1}))+\mathrm{i}\varepsilon
k[f_{2}^{*}R^{*}-f_{2}^{*}R^{-1}]\;.$

Now the EGF $\mathsf{B}(k,0)$ is nontrivial, choosing $Q(0)=R$:

$\mathsf{B}(k,0)=1+\mathrm{i}\varepsilon
k(f_{2}R^{3}+f_{2}^{*}R^{-1})+\mathrm{i}\varepsilon
ke^{-kR}[f_{2}^{*}R^{*}-f_{2}^{*}R^{-1}]\;.$

This allows for obtaining a closed expression of the constant function (for
choice $y(0)=1,\;s(0)=0$) as

$\mathcal{M}(k)=\mathrm{i}\varepsilon
k\left[f_{2}R^{3}+f_{2}^{*}\frac{R^{*}+2k+Rk^{2}}{(1+Rk)^{2}}\right]\;.$

After this, system (16) is to be solved.

###  VIII.2 Large resettings

Unfortunately, a transformation of the type (42) is hardly tractable for large
$\varepsilon$. Here we discuss another way of resetting, which leads to closed
expressions even for large changes of the phases. This approach is applicable
to identical oscillators subject to Cauchy white noise, but not for the
distribution of natural frequencies.

Suppose, in the OA state with order parameters $Z_{n}=R^{n}$, we randomly
choose a portion $\varepsilon$ of all oscillators and reset them completely
(this means that they “forget” their old states), cf. Ref. Sarkar and Gupta,
2022. We consider two variants below.

####  VIII.2.1 Random resetting

Here we assume that the new phases in the affected set become uniformly
distributed in the interval $[0,2\pi)$. These oscillators do not contribute to
new order parameters which thus take the values $Z_{n}=(1-\varepsilon)R^{n}$.
This corresponds to an offset Cauchy distribution, or specifically, a weighted
superposition of the Cauchy distribution and the uniform distribution. If one
uses variant A of initial conditions, then evolution starts from the initial
values $Q(0)=s(0)=0$, $y(0)=1$ and $\mathcal{M}(k)=(1-\varepsilon)Rk/(1-Rk)$
is determined via Eq. (27). Alternatively, adopting variant B, one can start
from the initial conditions $Q(0)=R$, $y(0)=1$, $s(0)=0$ and then
$\mathcal{M}(k)=-\varepsilon Rk/(1+Rk)$ is determined via Eq. (32). Then the
system evolves according to Eqs. (16).

However, due to the simplicity of this example, there is an even easier way of
treating this situation. Notice how moments can be viewed as a superposition
of two OA contributions, referred to as Poisson kernels by Ref Engelbrecht and
Mirollo, 2020:

$Z_{n}=(1-\varepsilon)Q_{1}+\varepsilon Q_{2}\;,$ (43)

where initially $Q_{1}(0)=R$ and $Q_{2}(0)=0$. In this case therefore, one can
evolve the system by considering two OA equations (9), which only interact
through the forcing $h$, and the solution maintains the form (43) for all
times.

####  VIII.2.2 Coherent resetting

Consider now that reset phases are not distributed uniformly but rather take
on another distribution $P^{\text{(res)}}(\varphi)$. If this distribution is a
wrapped Cauchy (which includes the uniform and the delta distribution) the
setting can again be treated simply as a superposition of the OA modes (a.k.a.
Poisson kernels) 222Resetting with a Kato-Jones Kato and Jones (2015)
distribution (asymmetric generalization of the wrapped Cauchy) can be
accomplished with two additional OA modes Cestnik and Pikovsky (2022)..
However, the reset distribution can generally have a different form. Below we
consider two cases.

Case (i): Partially coherent resetting. Here the reset phases are distributed
according to a single harmonic density:
$P^{\text{(res)}}(\varphi)=\frac{1}{2\pi}\left[1+2c\cos(\varphi-\varphi_{0})\right]$,
$c,\varphi_{0}\in\mathbb{R}$. The reset distribution therefore has only one
non-zero moment: $Z_{1}^{\text{(res)}}=ce^{\mathrm{i}\varphi_{0}}$, and the
full distribution after a portion $\varepsilon$ of phases are reset is
described by: $Z_{n}=(1-\varepsilon)R^{n}+\varepsilon
ce^{\mathrm{i}\varphi_{0}}\delta_{n,1}$. If using variant A initial
conditions, the variables initialize as $Q(0)=s(0)=0$, $y(0)=1$ and
$\mathcal{M}(k)=(1-\varepsilon)\frac{Rk}{1-Rk}+\varepsilon
ce^{\mathrm{i}\varphi_{0}}k$ according to Eq. (27). Alternatively if choosing
variant B, then $Q(0)=R$, $y(0)=1$, $s(0)=0$ and
$\mathcal{M}(k)=-\varepsilon\left[\frac{Rk}{1+Rk}+\frac{ce^{\mathrm{i}\varphi_{0}}k}{(1+Rk)^{2}}\right]$
according to Eq. (32). Then the system evolves following Eqs. (16). One could
also treat this setting as one OA mode and one general contribution described
by full Eqs. (16), we provide such a description in Appendix D.

Case (ii): Fully coherent resetting. Here the reset phases take the same value
$\varphi_{0}$. The new order parameters are thus
$Z_{n}=(1-\varepsilon)R^{n}+\varepsilon e^{\mathrm{i}n\varphi_{0}}$. Again
both variants of initializing the variables after the reset are possible. In
variant A one sets $Q(0)=s(0)=0$, $y(0)=1$ and
$\mathcal{M}(k)=(1-\varepsilon)\frac{Rk}{1-Rk}+\varepsilon\frac{e^{\mathrm{i}\varphi_{0}}k}{1-e^{\mathrm{i}\varphi_{0}}k}$
is determined via Eq. (27), while following variant B, we can set $Q(0)=R$,
$y(0)=1$, $s(0)=0$ and
$\mathcal{M}(k)=\varepsilon\frac{(e^{\mathrm{i}\varphi_{0}}-R)k}{1-(e^{\mathrm{i}\varphi_{0}}-R)k}$
is determined via Eq. (32). As mentioned before, since the delta distribution
is a special case of the wrapped Cauchy, this case could also be described by
just two OA modes (43).

The expressions above can be readily extended to a setup where several
randomly chosen subpopulations of oscillators are reset with different
distributions. Such an approach has been discussed in the context of
application to de-synchronization of neurons for Parkinson patients Tass
(1999).

Figure 3: Switching domains for large resetting in the bistable system (44).
Both the synchronous and asynchronous regimes are stable. Starting from the
synchronous regime on the OA manifold, we reset an $\varepsilon$ portion of
phases according to three example distributions: uniform (red domain), case
(i) – single harmonic density with amplitude $c=0.5$ (blue domain), and case
(ii) – delta distribution (green domain). The shaded regions mark areas that
induce a switch to asynchrony.

####  VIII.2.3 Numerical example for large resettings

As an example we consider a simple Kuramoto-type system with a synchrony-
asynchrony bistability Pikovsky and Rosenblum (2009). In this setup
$\omega=\text{const}$ (and one can without loss of generality set this
parameter to zero), and force is

$h=Z_{1}\exp[\mathrm{i}\theta_{0}+\mathrm{i}\theta_{1}|Z_{1}|^{2}]\;.$ (44)

We set $\gamma=0.1$, $\theta_{0}=0.8\pi$ and $\theta_{1}=4$. For these
parameters the states with $Z_{1}=0$ and $|Z_{1}|\approx 0.948$ are both
stable. We start with the latter state of a nearly synchronized ensemble, and
apply the three types of resetting as described above. By solving the reduced
six-dimensional equations (16), we obtain the domain of parameters
$\varepsilon,\varphi_{0}$ for which the resettings lead to a transition to the
asynchronous state $Z_{1}=0$.

For the random resetting, there is no dependence on $\varphi_{0}$, and the
corresponding domain is $\varepsilon>0.196$ (above red line in Fig. 3). For
the coherent resetting, we consider two cases discussed above, (i) and (ii),
and their corresponding basins are depicted in Fig. 3 with blue and green
domains, respectively. One can see that in all three cases, a finite
perturbation is needed to suppress synchrony. For the coherent resettings,
there is an optimal combination of $\varepsilon$ and $\varphi_{0}$; the
coherent subpopulation should be phase shifted around $\pi$ relative to the
phase of the mean field of the non-reset units. For case (ii) we also see that
if $\varepsilon$ is too large, the reset units form a new cluster and the
synchrony remains.

##  IX Conclusion

First, we summarize the approach and findings of this paper. Our starting
point is an infinite system of equations for the circular moments (order
parameters). These equations contain damping due to either Cauchy white noise,
or a Cauchy distribution of natural frequencies. By virtue of several
transformations, which are formulated in terms of generating functions, we
reduce this system to three complex equations. Additionally, a complex-valued
function of one variable is defined, which remains constant during the
evolution. The order parameters at each moment of time are represented through
this function and the three complex dynamical variables.

The original set of equations for the order parameters have the same form in
two situations: if the phase oscillators are subject to a Cauchy white noise,
and if the natural frequencies are Cauchy distributed (but time-independent).
Only in the former case there is a simple unique correspondence between the
order parameters and the distribution of the phases. In the latter situation,
one can calculate the order parameters from the distribution of the phases
(under the assumption of analyticity of the density in the upper complex plane
of frequencies), but it appears impossible to reconstruct this distribution
from the order parameters without further assumptions. Therefore, the results
of the paper are fully applicable to noisy ensembles, but some approaches
(e.g., phase resetting) are not suitable for the oscillators with distributed
frequencies.

The theory includes both the WS description (noise-free identical oscillators)
and the OA manifold (on which the dynamical variable $y$ vanishes). In the
framework of our approach, one can simply demonstrate that the dynamical
variable $y(t)$ tends to zero, which corresponds to the weak stability of the
OA manifold discussed in the literature. Therefore, our approach is an
essential improvement compared to OA theory, if a transient evolution from an
initial state outside of the OA manifold is important. In particular, it
allows for a calculation of the full basins of different attractors lying on
the OA manifold.

In this paper we operated with the phase equations. In some cases it is
convenient to transform the phase equations to other variables (e.g., theta-
neurons, equations which belong to class (1) Luke, Barreto, and So (2013);
Laing (2014), can be transformed to so-called quadratic integrate-and-fire
neurons Laing (2015); Montbrió, Pazó, and Roxin (2015); Bick  _et al._
(2020)). An extension of the theory to quadratic integrate-and-fire neurons
will be presented elsewhere Pietras, Cestnik, and Pikovsky (2022).

###### Acknowledgements.

We thank L. Smirnov and R. Toenjes for useful discussions. The work was
supported by DFG (grant No. PI 220/21-1).

##  Appendix A From the dynamics of the moments to the PDEs for the generating
functions

Consider a sequence of variables $\\{f_{n}\\}$ (the only condition is that the
generating functions below do exist; because we apply the theory to bounded
circular moments, this appears to always be the case). Also, in our case
$f_{0}\equiv 1$.

We distinguish between two types of generating functions: exponential
generating function (EGF) defined as

$\mathsf{F}(k)=\sum\limits_{n=0}^{\infty}f_{n}\frac{k^{n}}{n!}\;,$ (45)

and ordinary generating function (OGF) defined as

$\mathcal{F}(k)=\sum\limits_{n=1}^{\infty}f_{n}k^{n}\;.$ (46)

Here we show how the dynamics of the variables $f_{n}$, given by an infinite
set of ODEs, can be translated to the dynamics in terms of generating
functions $\mathsf{F}$ and $\mathcal{F}$, given by a single partial
differential equation (PDE). Suppose the dynamics of $f_{n}$ is as follows:

$\dot{f}_{n}=n\left(a\,f_{n-1}+b\,f_{n}+c\,f_{n+1}\right)\;,\quad n\geq 1\;,$
(47)

where $a,b,c$ are arbitrary complex quantities. Then the PDE for the EGF
$\mathsf{F}$ reads:

$\dot{\mathsf{F}}=a\,k\mathsf{F}+b\,k\mathsf{F}^{\prime}+c\,k\mathsf{F}^{\prime\prime}\;.$
(48)

In fact, in the case of an EGF $\mathsf{F}(k,t)$, any term $f_{n+m}$ for
$m\geq-1$ on the right-hand side of (47) just corresponds to a term of the
form $k\mathsf{F}^{(m+1)}$ in equation (48).

To prove this, it is sufficient to express derivatives of $\mathsf{F}$ as
formal series:

$\displaystyle\dot{\mathsf{F}}$
$\displaystyle=\sum\limits_{n=0}^{\infty}\dot{f}_{n}\frac{k^{n}}{n!}\;,$ (49)
$\displaystyle\mathsf{F}$
$\displaystyle=\sum\limits_{n=0}^{\infty}f_{n}\frac{k^{n}}{n!}=\frac{1}{k}\sum\limits_{n=0}^{\infty}(n+1)\,f_{n}\frac{k^{n+1}}{(n+1)!}=$
$\displaystyle=\frac{1}{k}\sum\limits_{n=-1}^{\infty}(n+1)\,f_{n}\frac{k^{n+1}}{(n+1)!}=\frac{1}{k}\sum\limits_{n=0}^{\infty}n\,f_{n-1}\frac{k^{n}}{n!}\;,$
$\displaystyle\mathsf{F}^{\prime}$
$\displaystyle=\sum\limits_{n=1}^{\infty}f_{n}\frac{k^{n-1}}{(n-1)!}=\frac{1}{k}\sum\limits_{n=1}^{\infty}n\,f_{n}\frac{k^{n}}{n!}=\frac{1}{k}\sum\limits_{n=0}^{\infty}n\,f_{n}\frac{k^{n}}{n!}\;,$
$\displaystyle\mathsf{F}^{\prime\prime}$
$\displaystyle=\sum\limits_{n=2}^{\infty}f_{n}\frac{k^{n-2}}{(n-2)!}=\frac{1}{k}\sum\limits_{n=2}^{\infty}(n-1)\,f_{n}\frac{k^{n-1}}{(n-1)!}=$
$\displaystyle=\frac{1}{k}\sum\limits_{n=1}^{\infty}n\,f_{n+1}\frac{k^{n}}{n!}=\frac{1}{k}\sum\limits_{n=0}^{\infty}n\,f_{n+1}\frac{k^{n}}{n!}\;,$
$\displaystyle\mathsf{F}^{(m)}$
$\displaystyle=\sum\limits_{n=m}^{\infty}f_{n}\frac{k^{n-m}}{(n-m)!}=$
$\displaystyle=\frac{1}{k}\sum\limits_{n=m}^{\infty}(n-m+1)\,f_{n}\frac{k^{n-m+1}}{(n-m+1)!}=$
$\displaystyle=\frac{1}{k}\sum\limits_{n=1}^{\infty}n\,f_{n+m-1}\frac{k^{n}}{n!}=\frac{1}{k}\sum\limits_{n=0}^{\infty}n\,f_{n+m-1}\frac{k^{n}}{n!}\;.$

Notice that the sum index changes in the derivation, e.g. if we have the
summed terms proportional to $n$ and the sum starts with index “1”:
$\sum\limits_{n=1}n\,(...)$, we can just add the zero term (since it vanishes)
and write the sum from index “0”. The way we have derived the above
expressions, the right-most forms of (49) all correspond to a sum starting at
$n=0$ and have a factor of $\frac{k^{n}}{n!}$, therefore we can just compare
the terms under the sum to see that the dynamics (47) corresponds to the PDE
(48).

More generally, a term in the dynamics of the variables $f_{n}$ having the
form

$\dot{f}_{n}=...+c\,n\,f_{n+m}+...\;,\quad m\geq-1\;,$ (50)

in terms of the EGF $\mathsf{F}(k,t)$ corresponds to the partial derivative in
the equation for $\mathsf{F}$:

$\dot{\mathsf{F}}=...+c\,k\mathsf{F}^{(m+1)}+...$ (51)

Now let us consider the OGF $\mathcal{F}(k,t)$. This generating function is
applicable to systems with only the homogeneous term $f_{n}$ and one higher
term $f_{n+1}$ on the r.h.s., i.e. to the case $a=0$ in (47). Let us write the
relevant sums:

$\displaystyle\dot{\mathcal{F}}$
$\displaystyle=\sum\limits_{n=1}^{\infty}\dot{f}_{n}k^{n}\;,$ (52)
$\displaystyle k\mathcal{F}^{\prime}$
$\displaystyle=\sum\limits_{n=1}^{\infty}n\,f_{n}k^{n}\;,$
$\displaystyle\mathcal{F}^{\prime}-\frac{1}{k}\mathcal{F}$
$\displaystyle=\sum\limits_{n=1}^{\infty}n\,f_{n}k^{n-1}-f_{n}k^{n-1}=$
$\displaystyle=\sum\limits_{n=1}^{\infty}(n-1)\,f_{n}k^{n-1}=\sum\limits_{n=1}^{\infty}n\,f_{n+1}k^{n}\;.$

Again, by comparing the right-most expressions of (52), we can see that the
dynamics

$\dot{f}_{n}=b\,n\,f_{n}+c\,n\,f_{n+1}\;,$ (53)

in terms of the OGF $\mathcal{F}$, corresponds to

$\dot{\mathcal{F}}=b\,k\mathcal{F}^{\prime}+c\,\Big{[}\mathcal{F}^{\prime}-\frac{1}{k}\mathcal{F}\Big{]}\;.$
(54)

##  Appendix B Transformation (21) in terms of variables

We rewrite here relation (21):

$\frac{\mathcal{M}(k)}{k}=\frac{\mathcal{A}(k+s)}{k+s}\;.$ (55)

Expressed as an OGF series (with a different index for the right-hand side),
it reads:

$\sum\limits_{n=1}^{\infty}\mu_{n}k^{n-1}=\sum\limits_{r=1}^{\infty}\alpha_{r}(k+s)^{r-1}\;.$

Next, we expand the binomial term on the r.h.s.:

$\sum\limits_{n=1}^{\infty}\mu_{n}k^{n-1}=\sum\limits_{r=1}^{\infty}\alpha_{r}\sum\limits_{d=0}^{r-1}\binom{r-1}{d}k^{r-d-1}s^{d}\;,$

gather together terms for which $r-d-1=n-1$, and use the same index $n$ as on
the l.h.s.:

$\sum\limits_{n=1}^{\infty}\mu_{n}k^{n-1}=\sum\limits_{n=1}^{\infty}k^{n-1}\sum\limits_{m=n}^{\infty}\binom{m-1}{m-n}\alpha_{m}s^{m-n}\;.$

Now it is clear that this transformation corresponds to (22):

$\mu_{n}=\sum\limits_{m=n}^{\infty}\binom{m-1}{m-n}\alpha_{m}s^{m-n}\;,$ (56)

where the binomial coefficient $\binom{m-1}{m-n}$ can also be written as
$\binom{m-1}{n-1}$.

##  Appendix C Expressing moments $Z_{n}$ with the constant function
$\mathcal{M}(k)$

The $d^{\text{th}}$ derivative in $k$ of $\mathcal{M}(k)$ is expressed as:

$\mathcal{M}^{(d)}(k)=\frac{d!}{k^{d}}\sum\limits_{m=1}^{\infty}\binom{m}{d}\mu_{m}k^{m}\;.$
(57)

Next, using the binomial relation:
$\binom{m-1}{n-1}=(-1)^{n-1}\sum\limits_{d=0}^{n-1}(-1)^{d}\binom{m}{d}$ and
the relation between $\alpha_{n}$ and $\mu_{n}$ variables (22), we can write:

$\alpha_{n}=-\sum\limits_{d=0}^{n-1}\frac{s^{d-n}}{d!}\mathcal{M}^{(d)}(-s)\;.$
(58)

Then we can use relation (17) to express $\beta_{n}$ as:

$\beta_{n}=-y^{n}\sum\limits_{d=0}^{n-1}\frac{s^{d-n}}{d!}\mathcal{M}^{(d)}(-s)\;,\quad
n\geq 1\;,$ (59)

and we can use (15) to express the moments:

$Z_{n}=Q^{n}-\sum\limits_{m=1}^{n}\binom{n}{m}Q^{n-m}y^{m}\sum\limits_{d=0}^{m-1}\frac{s^{d-m}}{d!}\mathcal{M}^{(d)}(-s)\;.$
(60)

For initial condition variant A: $Q(0)=s(0)=0$, $y(0)=1$, to see that
expression (60) is correct initially at $t=0$, we have to evaluate the
following limits:

$\lim\limits_{\varepsilon\to
0}\sum\limits_{d=0}^{m-1}\frac{\varepsilon^{d-m}}{d!}\mathcal{M}^{(d)}(-\varepsilon)=-\mu_{m}\;,$
(61)

which confirm that $Z_{n}(0)=\mu_{n}$.

##  Appendix D Alternative perturbation around the OA manifold

There are different ways one can write states close to the OA manifold. An
alternative to the main text (30) is:

$Z_{n}=(1-\varepsilon)Z^{n}+\varepsilon\,p_{n}\;,\quad|\varepsilon|\ll 1\;,$
(62)

Notice that due to the linearity of Eq. (8), we can view this simply as
splitting the state into two parts: one on and one off the OA manifold:
$Z_{n}^{(\text{OA})}=Z^{n}$ and $Z_{n}^{(\text{pert})}=p_{n}$. Each of the two
parts can then be treated separately, the only thing connecting them is the
force $h$. We know that the part on the OA manifold requires only one complex
equation (16a) by considering $\mathcal{M}^{(\text{OA})}=0$, while the rest
can be treated as before with the full set of three equations (16) and the
simplest general initial conditions (26). The complete dynamics is then
written as:

$\begin{split}\dot{Z}&=(\mathrm{i}\omega-\gamma)Z+h-h^{*}Z^{2}\;,\\\
\dot{Q}&=(\mathrm{i}\omega-\gamma)Q+h-h^{*}Q^{2}\;,\\\
\dot{y}&=(\mathrm{i}\omega-\gamma-2h^{*}Q)y\;,\\\
\dot{s}&=h^{*}y\;,\end{split}$ (63)

where the perturbations $Z_{n}^{(\text{pert})}=p_{n}$ express with $Q,y,s$ via
Eq. (23). Global order parameters are simply the weighted sum of both
contributions, e.g. the first order parameter is expressed as:
$Z_{1}=(1-\varepsilon)Z+\varepsilon\Big{(}Q-y\frac{\mathcal{M}(-s)}{s}\Big{)}$,
where $\mathcal{M}=\sum\limits_{n=1}^{\infty}p_{n}(t=0)k^{n}$. Notice that
this perturbation description is actually global, $\varepsilon$ does not need
to be small.

##  Appendix E Alternative constant function
$\mathcal{N}(k)=\frac{\mathcal{M}(-k)}{-k}$

Instead of transformation (21) we can consider the following:

$\mathcal{N}(k)=\frac{\mathcal{A}(s-k)}{s-k}=\frac{\mathcal{M}(-k)}{-k}\;.$
(64)

The OGF $\mathcal{N}(k)$ is also a constant function in time, its time
derivative yields:

$\dot{\mathcal{N}}(k)=\frac{\dot{\mathcal{A}}(s-k)}{s-k}+\dot{s}\Big{[}\frac{\mathcal{A}^{\prime}(s-k)}{s-k}-\frac{\mathcal{A}(s-k)}{(s-k)^{2}}\Big{]}\
=0\;.$

The benefit of this function is that it simplifies the expression for
Kuramoto-Daido order parameters (23):

$Z_{n}=Q^{n}\Big{[}1-\sum\limits_{m=1}^{n}\binom{n}{m}\Big{(}-\frac{y}{Q}\Big{)}^{m}\frac{\mathcal{N}^{(m-1)}(s)}{(m-1)!}\Big{]}\;,$
(65)

the first three moments explicitly:

$\displaystyle Z_{1}$ $\displaystyle=Q+y\mathcal{N}(s)\;,$ $\displaystyle
Z_{2}$ $\displaystyle=Q^{2}+2Qy\mathcal{N}(s)-y^{2}\mathcal{N}^{\prime}(s)\;,$
$\displaystyle Z_{3}$
$\displaystyle=Q^{3}+3Q^{2}y\mathcal{N}(s)-3Qy^{2}\mathcal{N}^{\prime}(s)+y^{3}\frac{\mathcal{N}^{\prime\prime}(s)}{2}\;,$

and no expansions for small $s$ like (25) are needed.

If we consider this function $\mathcal{N}(k)$ as the OGF of some variables
$\nu_{n}$ we now have to also consider the zero${}^{\text{th}}$ term:

$\mathcal{N}(k)=\sum_{n=0}^{\infty}\nu_{n}k^{n}\;.$ (66)

(notice $n$ starting from 0). Constant quantities $\nu_{n}$ express with
$\mu_{n}$ as:

$\nu_{n}=(-1)^{n}\mu_{n+1}\;.$ (67)

Considering variant A initial conditions (26): $Q(0)=s(0)=0,\ y(0)=1$, these
quantities can be determined as:

$\nu_{n}=(-1)^{n}Z_{n+1}(0)\;,$ (68)

and so the constant function $\mathcal{N}(k)$ can be expressed as:

$\mathcal{N}(k)=\sum_{n=0}^{\infty}(-1)^{n}Z_{n+1}(0)k^{n}\;.$ (69)

##  Appendix F Approximating $\mathcal{M}(k)$ with a finite series

Figure 4: The switching domains for resetting the bistable system (44) with a
$\delta$-distribution, for different approximations the constant function
$\mathcal{M}(k)$ (the black cross-shaded region corresponds to the green
region in Fig. 3). The considered approximations are truncations of its Taylor
series at different numbers of terms:
$\mathcal{M}_{N}=\sum\limits_{n=1}^{N}Z_{n}(0)k^{n}$.

For empirically observed phase densities it might not always be clear how to
determine the constant function $\mathcal{M}(k)$, but one can always estimate
the first few moments $Z_{n}(0)$ numerically and then approximate the function
with a series: $\mathcal{M}(k)\approx\sum\limits_{n=1}^{N}Z_{n}(0)k^{n}$. Here
in Fig. 4 we show an example of that. We consider the bistable system (44) and
just like in Fig. 3 compute the switching domain for resetting with a
$\delta$-distribution. Several different truncations of the series are
considered: $N=1,2,3,4$, and the corresponding domains are depicted in Fig. 4.
We see a fast convergence of the approximations.

##  Appendix G Relation to the Watanabe-Strogatz theory

In this appendix we demonstrate explicitly that the derived equations coincide
with the Watanabe-Strogatz equations for identical oscillators. Our basic
equation for the oscillator dynamics (1) in the noise-free case reads

$\dot{\varphi}=\omega+\text{Im}[2he^{-\mathrm{i}\varphi}]=\omega+2\text{Im}[h]\cos\varphi-2\text{Re}[h]\sin\varphi\;.$

Watanabe and Strogatz (Eq. (5.3) in Ref. Watanabe and Strogatz, 1994) write
this equation in slightly different notations

$\dot{\varphi}=F+G\cos\varphi+H\sin\varphi\;,$

so that $F=\omega$, $G=2\text{Im}[h]$, $H=-2\text{Re}[h]$. After a
transformation (Eq. (5.4) in Ref. Watanabe and Strogatz, 1994)

$\tan\left(\frac{\varphi-\Phi}{2}\right)=\sqrt{\frac{1+\eta}{1-\eta}}\tan\left(\frac{\psi-\Psi}{2}\right)\;,$

they obtain a set of equations for the variables $\eta,\Phi,\Psi$ (Eq. (5.9)
in Ref. Watanabe and Strogatz, 1994)

$\displaystyle\dot{\eta}$
$\displaystyle=-(1-\eta^{2})(G\sin\Phi-H\cos\Phi)\;,$ (70)
$\displaystyle\eta\dot{\Phi}$ $\displaystyle=\eta F-G\cos\Phi-H\sin\Phi\;,$
$\displaystyle\eta\dot{\Psi}$
$\displaystyle=-\sqrt{1-\eta^{2}}(G\cos\Phi+H\sin\Phi)\;.$

The variables $\psi$ are constants of motion.

Let us introduce a new variable $\rho$ according to

$\frac{1-\rho}{1+\rho}=\sqrt{\frac{1+\eta}{1-\eta}}\;.$

Then $\eta=-\frac{2\rho}{1+\rho^{2}}$,
$1-\eta^{2}=\frac{(1-\rho^{2})^{2}}{(1+\rho^{2})^{2}}$,
$\sqrt{1-\eta^{2}}=\frac{1-\rho^{2}}{1+\rho^{2}}$, and
$\dot{\eta}=-2\dot{\rho}\frac{1-\rho^{2}}{(1+\rho^{2})^{2}}$. We substitute
these relations in (70) and obtain

$\displaystyle\dot{\rho}$
$\displaystyle=\frac{1-\rho^{2}}{2}(G\sin\Phi-H\cos\Phi)\;,$ (71)
$\displaystyle\dot{\Phi}$
$\displaystyle=F+\frac{1+\rho^{2}}{2\rho}(G\cos\Phi+H\sin\Phi)\;,$
$\displaystyle\dot{\Psi}$
$\displaystyle=\frac{1-\rho^{2}}{2\rho}(G\cos\Phi+H\sin\Phi)\;.$

Let us now introduce variables $Q=\rho e^{\mathrm{i}\Phi}$ and
$\theta=\Phi-\Psi$ and assume that these variables fulfill Eqs. (16a) (with
$\gamma=0$), (40) above:

$\displaystyle\dot{Q}$ $\displaystyle=\mathrm{i}\omega Q+h-h^{*}Q^{2}\;,$ (72)
$\displaystyle\dot{\theta}$
$\displaystyle=\omega-\mathrm{i}(hQ^{*}-h^{*}Q)\;.$

The equations for $\rho,\Phi,\Psi$ following from (72) read:

$\displaystyle\dot{\rho}$
$\displaystyle=(1-\rho^{2})(\text{Re}[h]\cos\Phi+\text{Im}[h]\sin\Phi)\;,$
(73) $\displaystyle\dot{\Phi}$
$\displaystyle=\omega+\frac{1+\rho^{2}}{\rho}(\text{Im}[h]\cos\Phi-\text{Re}[h]\sin\Phi)\;,$
$\displaystyle\dot{\Psi}$
$\displaystyle=\frac{1-\rho^{2}}{\rho}(\text{Im}[h]\cos\Phi-\text{Re}[h]\sin\Phi)\;.$

Substituting here $\omega=F$, $\text{Re}[h]=-\frac{H}{2}$ and
$\text{Im}[h]=\frac{G}{2}$, we obtain exactly system (71), which proves the
equivalence of our equations (16a) (with $\gamma=0$), (40) and the WS
equations (71).

## References

  * Nixon  _et al._ (2013) M. Nixon, E. Ronen, A. A. Friesem, and N. Davidson, “Observing geometric frustration with thousands of coupled lasers,” Phys. Rev. Lett. 110, 184102 (2013).
  * Wiesenfeld and Swift (1995) K. Wiesenfeld and J. W. Swift, “Averaged equations for Josephson junction series arrays,” Phys. Rev. E 51, 1020–1025 (1995).
  * Wiesenfeld, Colet, and Strogatz (1998) K. Wiesenfeld, P. Colet, and S. Strogatz, “Frequency locking in Josephson arrays: Connection with the Kuramoto model,” Phys. Rev. E 57, 1563–1569 (1998).
  * Totz  _et al._ (2018) J. F. Totz, J. Rode, M. R. Tinsley, K. Showalter, and H. Engel, “Spiral wave chimera states in large populations of coupled chemical oscillators,” Nat. Phys. 14, 282 (2018).
  * Eckhardt  _et al._ (2007) B. Eckhardt, E. Ott, S. H. Strogatz, D. M. Abrams, and A. McRobie, “Modeling walker synchronization on the Millennium Bridge,” Phys. Rev. E 75, 021110 (2007).
  * Luke, Barreto, and So (2013) T. B. Luke, E. Barreto, and P. So, “Complete classification of the macroscopic behavior of a heterogeneous network of theta neurons,” Neural Comput. 25, 3207–3234 (2013).
  * Laing (2014) C. R. Laing, “Derivation of a neural field model from a network of theta neurons,” Phys. Rev. E 90, 010901(R) (2014).
  * Holstein-Rathlou  _et al._ (2001) N. H. Holstein-Rathlou, K. P. Yip, O. V. Sosnovtseva, and E. Mosekilde, “Synchronization phenomena in nephron-nephron interaction,” Chaos 11, 417–426 (2001).
  * Prindle  _et al._ (2012) A. Prindle, P. Samayoa, I. Razinkov, T. Danino, L. S. Tsimring, and J. Hasty, “A sensing array of radically coupled genetic “biopixels”,” Nature 481, 39–44 (2012).
  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves and Turbulence_ (Springer, Berlin, 1984).
  * Kuramoto (1975a) Y. Kuramoto,  _in International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki, Vol. 39 (Springer, New York, 1975) p. 420.
  * Sakaguchi and Kuramoto (1986) H. Sakaguchi and Y. Kuramoto, “A soluble active rotator model showing phase transition via mutual entrainment,” Prog. Theor. Phys. 76, 576–581 (1986).
  * Acebrón  _et al._ (2005) J. A. Acebrón, L. L. Bonilla, C. J. P. Vicente, F. Ritort, and R. Spigler, “The Kuramoto model: A simple paradigm for synchronization phenomena,” Rev. Mod. Phys. 77, 137–175 (2005).
  * Watanabe and Strogatz (1993) S. Watanabe and S. H. Strogatz, “Integrability of a globally coupled oscillator array,” Phys. Rev. Lett. 70, 2391–2394 (1993).
  * Watanabe and Strogatz (1994) S. Watanabe and S. H. Strogatz, “Constants of motion for superconducting Josephson arrays,” Physica D 74, 197–253 (1994).
  * Ott and Antonsen (2008) E. Ott and T. M. Antonsen, “Low dimensional behavior of large systems of globally coupled oscillators,” Chaos 18, 037113 (2008).
  * Tanaka (2020) T. Tanaka, “Low-dimensional dynamics of phase oscillators driven by cauchy noise,” Phys. Rev. E 102, 042220 (2020).
  * Tönjes and Pikovsky (2020) R. Tönjes and A. Pikovsky, “Low-dimensional description for ensembles of identical phase oscillators subject to Cauchy noise,” Phys. Rev. E 102, 052315 (2020).
  * Ott and Antonsen (2009) E. Ott and T. M. Antonsen, “Long time evolution of phase oscillator systems,” Chaos 19, 023117 (2009).
  * Pietras and Daffertshofer (2016) B. Pietras and A. Daffertshofer, “Ott-Antonsen attractiveness for parameter-dependent oscillatory systems,” Chaos 26, 103101 (2016).
  * Engelbrecht and Mirollo (2020) J. R. Engelbrecht and R. Mirollo, “Is the Ott-Antonsen manifold attracting?” Phys. Rev. Research 2, 023057 (2020).
  * Cestnik and Pikovsky (2022) R. Cestnik and A. Pikovsky, “Hierarchy of exact low-dimensional reductions for populations of coupled oscillators,” Phys. Rev. Lett. 128, 054101 (2022).
  * Braun  _et al._ (2012) W. Braun, A. Pikovsky, M. A. Matias, and P. Colet, “Global dynamics of oscillator populations under common noise,” EPL 99, 20006 (2012).
  * Gong  _et al._ (2019) C. C. Gong, C. Zheng, R. Toenjes, and A. Pikovsky, “Repulsively coupled Kuramoto-Sakaguchi phase oscillators ensemble subject to common noise,” Chaos 29, 033127 (2019).
  * Chechkin  _et al._ (2003) A. V. Chechkin, J. Klafter, V. Y. Gonchar, R. Metzler, and L. V. Tanatarov, “Bifurcation, bimodality, and finite variance in confined Lévy flights,” Phys. Rev. E 67, 010102 (2003).
  * Toenjes, Sokolov, and Postnikov (2013) R. Toenjes, I. M. Sokolov, and E. B. Postnikov, “Nonspectral relaxation in one dimensional Ornstein-Uhlenbeck processes,” Phys. Rev. Lett. 110, 150602 (2013).
  * Toenjes, Sokolov, and Postnikov (2014) R. Toenjes, I. M. Sokolov, and E. B. Postnikov, “Spectral properties of the fractional fokker-planck operator for the Lévy flight in a harmonic potential,” The European Physical Journal B 87, 1–11 (2014).
  * Daido (1996) H. Daido, “Onset of cooperative entrainment in limit-cycle oscillators with uniform all-to-all interactions: bifurcation of the order function,” Physica D 91, 24–66 (1996).
  * Kuramoto (1975b) Y. Kuramoto, “Self-entrainment of a population of coupled nonlinear oscillators,” in  _International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki (Springer Lecture Notes Phys., v. 39, New York, 1975) p. 420.
  * Pikovsky and Rosenblum (2011a) A. Pikovsky and M. Rosenblum, “Dynamics of heterogeneous oscillator ensembles in terms of collective variables,” Physica D 240, 872–881 (2011a).
  * Note (1) The sum defining ordinary generating functions is typically considered from index $n=0$, but in our context it is convenient to start from $n=1$; note that we always have $f_{0}=1$ for normalization reasons.
  * Haukkanen (1993) P. Haukkanen, “Formal power series for binomial sums of sequences of numbers,” Fibonacci Quart. 31, 28–31 (1993).
  * Prodinger (1994) H. Prodinger, “Some information about the binomial transform,” Fibonacci Quart. 32, 412–415 (1994).
  * Pikovsky and Rosenblum (2008) A. Pikovsky and M. Rosenblum, “Partially integrable dynamics of hierarchical populations of coupled oscillators,” Phys. Rev. Lett. 101, 264103 (2008).
  * Pikovsky and Rosenblum (2011b) A. Pikovsky and M. Rosenblum, “Dynamics of heterogeneous oscillator ensembles in terms of collective variables,” Physica D 240, 872–881 (2011b).
  * Kato and Jones (2015) S. Kato and M. C. Jones, “A tractable and interpretable four-parameter family of unimodal distributions on the circle,” Biometrika 102, 181–190 (2015).
  * Ichiki and Okumura (2020) A. Ichiki and K. Okumura, “Diversity of dynamical behaviors due to initial conditions: Extension of the Ott-Antonsen ansatz for identical Kuramoto-Sakaguchi phase oscillators,” Phys. Rev. E 101, 022211 (2020).
  * Marvel, Mirollo, and Strogatz (2009) S. A. Marvel, R. E. Mirollo, and S. H. Strogatz, “Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action,” Chaos 19, 043104 (2009).
  * Canavier (2006) C. C. Canavier, “Phase response curve,” Scholarpedia 1, 1332 (2006).
  * Smeal, Ermentrout, and White (2010) R. M. Smeal, G. B. Ermentrout, and J. A. White, “Phase-response curves and synchronized neural networks,” Philosophical Transactions of the Royal Society B: Biological Sciences 365, 2407–2422 (2010).
  * Schwabedal  _et al._ (2012) J. Schwabedal, A. Pikovsky, B. Kralemann, and M. Rosenblum, “Optimal phase description of chaotic oscillators,” Phys. Rev. E 85, 026216 (2012).
  * Pikovsky  _et al._ (1997) A. Pikovsky, M. Rosenblum, G. Osipov, and J. Kurths, “Phase synchronization of chaotic oscillators by external driving,” Physica D 104, 219–238 (1997).
  * Grines, Osipov, and Pikovsky (2018) E. Grines, G. Osipov, and A. Pikovsky, “Describing dynamics of driven multistable oscillators with phase transfer curves,” Chaos 28, 106323 (2018).
  * Sarkar and Gupta (2022) M. Sarkar and S. Gupta, “Synchronization in the Kuramoto model in presence of stochastic resetting,” arXiv , 2203.00339 (2022).
  * Note (2) Resetting with a Kato-Jones Kato and Jones (2015) distribution (asymmetric generalization of the wrapped Cauchy) can be accomplished with two additional OA modes Cestnik and Pikovsky (2022).
  * Tass (1999) P. A. Tass,  _Phase Resetting in Medicine and Biology. Stochastic Modelling and Data Analysis._ (Springer-Verlag, Berlin, 1999).
  * Pikovsky and Rosenblum (2009) A. Pikovsky and M. Rosenblum, “Self-organized partially synchronous dynamics in populations of nonlinearly coupled oscillators,” Physica D 238(1), 27–37 (2009).
  * Laing (2015) C. R. Laing, “Exact neural fields incorporating gap junctions,” SIAM J. Appl. Dyn. Syst. 14, 1899–1929 (2015).
  * Montbrió, Pazó, and Roxin (2015) E. Montbrió, D. Pazó, and A. Roxin, “Macroscopic description for networks of spiking neurons,” Phys. Rev. X 5, 021028 (2015).
  * Bick  _et al._ (2020) C. Bick, M. Goodfellow, C. R. Laing, and E. A. Martens, “Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review,” J. Math. Neurosci. 10, 9 (2020).
  * Pietras, Cestnik, and Pikovsky (2022) B. Pietras, R. Cestnik, and A. Pikovsky, “Exact finite-dimensional description for networks of globally coupled spiking neurons,” arXiv 2209.00922 (2022).