# A unified quantification of synchrony in globally coupled populations with
the Wiener order parameter

Arkady Pikovsky  Institute of Physics and Astronomy, University of Potsdam,
Karl-Liebknecht-Str. 24/25, 14476 Potsdam-Golm, Germany  Michael Rosenblum
Institute of Physics and Astronomy, University of Potsdam, Karl-Liebknecht-
Str. 24/25, 14476 Potsdam-Golm, Germany

###### Abstract

We tackle the quantification of synchrony in globally coupled populations.
Furthermore, we treat the problem of incomplete observations when the
population mean field is unavailable, but only a small subset of units is
observed. We introduce a new order parameter and demonstrate its efficiency
for quantifying synchrony via monitoring general observables, regardless of
whether the oscillations can be characterized in terms of the phases. Under
condition of a significant irregularity in the dynamics of the coupled units,
this order parameter provides a unified description of synchrony in
populations of units of various complexity. The main examples include noise-
induced oscillations, coupled strongly chaotic systems, and noisy periodic
oscillations. Furthermore, we explore how this parameter works for the
standard Kuramoto model of coupled regular phase oscillators. The most
significant advantage of our approach is its ability to infer and quantify
synchrony from the observation of a small percentage of the units and even
from a single unit, provided the observations are sufficiently long.

oscillatory ensembles, synchrony, order parameter, autocorrelation, data
analysis

> Coupling many oscillating systems may result in a synchronization
> transition, at which macroscopic collective motion appears. Examples include
> pedestrians on a bridge, synchronously blinking fireflies, and applause in
> theater halls. In many cases, it is not easy to distinguish synchrony and
> asynchrony. We suggest a method for quantifying synchrony by measuring the
> regularity of the observed fields. Remarkably, not the whole population
> should be observed; in many cases, just observing one unit allows for a
> reliable identification of the synchronization transition in the whole
> ensemble.

##  I Introduction

The emergence of collective activity in large networks of coupled active units
has been known for decades Winfree (1967); Kuramoto (1975); Strogatz and
Stewart (1993); Acebrón  _et al._ (1998); Pikovsky, Rosenblum, and Kurths
(2001), and yet remains a topic of intensive studies in the physics of complex
systems. The relevant applications range from bridge engineering, where
pedestrian synchrony under certain circumstances induces bridge vibration
Eckhardt  _et al._ (2007), to neuroscience, where coordination of neuronal
activity forms macroscopic brain rhythms, vital or pathological Glass and
Mackey (1988); Buzsáki (2006); Lehnertz (2008); Daffertshofer and Pietras
(2020).

This paper discusses the quantification of the degree of unit activity
coordination in globally coupled ensembles, a popular model for highly
interconnected networks. The addressed problem is pertinent both for
theoretical and experimental studies. A well-known solution applies when all
units are accessible, so that a mean field can be obtained. Macroscopic
oscillations of this mean field indicate for collective synchrony, while in
the asynchronous case, the mean field just fluctuates. This suggest using the
variance of the mean field for characterization of synchrony Ginzburg and
Sompolinsky (1994); Hansel and Sompolinsky (1996); Pikovsky, Rosenblum, and
Kurths (1996); Golomb, Hansel, and Mato (2001). A better insight is possible
if the units possess well-defined phases, or the phases of all units can be
inferred from data. In that case, the Kuramoto order parameter (KOP) and
generalized Kuramoto-Daido order parameters are typically used to characterize
the degree of synchrony.

However, the quantification is not that easy in the two outlined cases.

  1. 1.

For some systems, the phases are not defined or are difficult to extract; the
examples detailed below include strongly chaotic systems and noise-induced
oscillations, for which the phase description is hardly possible. Further
relevant cases are, e.g., bursting neurons. Though the mean field observation
quantifies the collective activity in these complicated cases, this measure
can be substituted by a more sensitive one, as shown below.

  2. 2.

The more challenging case is that of an incomplete observation. Suppose we
monitor only a small subpopulation of units - can we quantify the level of
synchrony in the whole, unobserved ensemble? This paper provides a positive
answer for the most frequent case of a regular collective mode.

We tackle the synchrony quantification problem by introducing a novel order
parameter based on Wiener’s lemma and measuring the intensity of the regular
(periodic or quasiperiodic) component in a time series. We demonstrate that
the Wiener order parameter (WOP) (i) successfully differentiates synchronous
and asynchronous states for coupled units; (ii) in the cases where the phases
are accessible and the Kuramoto-type order parameters can be calculated, the
WOP still provides a better discrimination between these states than KOP and
the variance-based measure for a finite-size ensemble; (iii) yields a proper
quantification from an incomplete observation, even from an observation of a
single unit. We illustrate our approach with examples covering noise-induced,
chaotic, noisy, and regular dynamics.

The paper is organized as follows. In Section II we provide theoretical
considerations, assumptions and limitations, define the new order parameter,
and discuss the computational aspects. Sections III-V exemplify the approach
by analysis of ensembles of irregular identical and non-identical units.
Section VI extends the consideration to the case of regular non-identical
units. Finally, Section VII concludes and discusses our findings.

##  II Regularity of mean fields and its characterization

###  II.1 General mean-field coupling

We consider a population of units described by deterministic or noisy
equations. We assume at the beginning that the units are identical, and the
noises (terms $\vec{\xi}_{k}(t)$ in (1)), if present, are independent.
Furthermore, we assume that the systems are coupled via mean fields
$\vec{Y},\vec{Z}$ which are either just global averages of some observables
(global variables $\vec{Y}$) or obey dynamical equations where only the mean
fields are entering (global variables $\vec{Z}$). Then, the equations read

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\vec{x}_{k}$
$\displaystyle=\vec{f}(\vec{x}_{k},\vec{\xi}_{k}(t),\vec{Y},\vec{Z}),\quad
k=1,\ldots,N\;,$ (1) $\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\vec{Z}$
$\displaystyle=\vec{F}(\vec{Z},\vec{Y})\;,$ (2) $\displaystyle\vec{Y}(t)$
$\displaystyle=\frac{1}{N}\sum_{k=1}^{N}\vec{y}(\vec{x}_{k}(t))\;.$ (3)

Notice that dimensions of the vectors $\vec{x},\vec{\xi},\vec{Y},\vec{Z}$
generally differ. The global coupling in the population is organized via a
summation of an observable $\vec{y}$ defined as a function of local variables
over all the units in the ensemble. Quite often, one assumes an “algebraic”
coupling, where only mean fields $\vec{Y}$ are present, but in many
situations, like in coupled clocks or metronomes on a beam Kapitaniak  _et
al._ (2012); Czolczyński  _et al._ (2013), pedestrians on a bridge Eckhardt
_et al._ (2007), and electronic oscillators with a common load Wiesenfeld,
Colet, and Strogatz (1998); Temirbayev  _et al._ (2012) also the global
dynamical equations (2) are present.

We start by assuming that individual units are irregular, i.e., their dynamics
either possess noisy forces $\vec{\xi}_{k}$, or the units operate in a chaotic
regime. Furthermore, we will assume the absence of multistability in
individual units. Below, we will also discuss the potential applicability of
the approach to regular systems.

In the thermodynamic limit $N\to\infty$, it is appropriate to introduce the
probability density of the unit states $w(\vec{x},t)$, which evolves according
to some linear operator that follows from Eq. (1) and is the Liouville
operator in the deterministic case, a Fokker-Plank operator in the case of
Gaussian noises, or some generalized operator if the noise is not white and
Gaussian. Then, instead of system (1-3) we can write

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}w$
$\displaystyle=\hat{L}(\vec{Y},\vec{Z})w\;,$ (4)
$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\vec{Z}$
$\displaystyle=\vec{F}(\vec{Z},\vec{Y})\;,$ (5) $\displaystyle\vec{Y}(t)$
$\displaystyle=\int\mathrm{d}\vec{x}\;w(\vec{x},t)\;\vec{y}(\vec{x})\;.$ (6)

Dependence of the operator $\hat{L}$ on the mean fields, which themselves
depend on the density $w$ via Eqs. (5,6), makes the whole system (4-6)
nonlinear.

In this description, we attribute a steady (time-independent) solution of Eqs.
(4-6) to asynchrony and a time-dependent (periodic or more complex) solution
to synchrony. Indeed, if the mean fields are constants, then the individual
units can be treated as independent identical irregular oscillators. Thus, the
state of each oscillator admits a statistical description with a stationary
invariant distribution, which is a stable stationary solution of Eq. (4). On
the contrary, a synchronous state with oscillating mean fields can be
interpreted as follows. Oscillatory mean fields impose a common forcing on all
units, and due to this, the dynamics of each unit contains two components: an
irregular one (and these components in different units are independent) and a
forcing-induced one (and this component in different units is the same). At
the calculation of the mean field, the irregular components due to the law of
large numbers sum to a constant, while forcing-induced components sum to a
macroscopic time-varying field.

Below, we concentrate on the most common case of regular (periodic or
quasiperiodic) mean fields; a short discussion of irregular mean fields will
be given in Section II.5.

###  II.2 Complete and partial observations

We will suppose that a scalar observable of a specific system $u(\vec{x}_{k})$
is available, possibly not for all units (such a situation often occurs in
neuroimaging Wenzel and Hamm (2021)). Let us define the partial mean fields as

$U_{M}(t)=M^{-1}\sum_{k\in\mathcal{S}_{M}}u(\vec{x}_{k}(t))\;.$ (7)

Here, $1\leq M\leq N$ is the size of the set $\mathcal{S}_{M}$ over which the
averaging is performed. The set of the “observed” units $\mathcal{S}_{M}$ can
be chosen at random; because we assume identical elements, the statistical
properties of $U_{M}$ should not depend on this choice. The case $M=1$ means
the observation of just one oscillator from the population; the case $M=N$
yields the usual mean field to which all the units contribute. We denote this
full observable $U(t)$, dropping the index. Its definition is similar to that
of mean fields $\vec{Y}$ (Eq. (3)), and in fact one can (but this is not
necessary) choose $U$ as one of the components of $\vec{Y}$.

In the thermodynamic limit,

$U(t)=\int\mathrm{d}\vec{x}\;w(\vec{x},t)\;u(\vec{x})\;.$ (8)

In this limit, $U=\text{const}$ in an asynchronous state; while time-dependent
$U(t)$ (typically periodic or quasiperiodic) manifests synchrony in the
population. A straightforward way to distinguish these cases is to calculate
the variance of $U(t)$:

$\text{Var}(U)=\langle(U(t)-\bar{U})^{2}\rangle,\qquad\bar{U}=\langle
U(t)\rangle\;,$

where the brackets denote statistical averaging 111In this paper we assume
ergodicity and stationarity, so practically the averaging is performed over
time.. This quantity vanishes in the asynchronous regime and is finite in the
synchronous one, thus it is widely used as an empirical order parameter
characterizing transition to synchrony Ginzburg and Sompolinsky (1994); Hansel
and Sompolinsky (1996); Pikovsky, Rosenblum, and Kurths (1996); Golomb,
Hansel, and Mato (2001).

###  II.3 Finite-sample fluctuations and their elimination

Even in the thermodynamic limit, for partial observations of sample size $M$,
one cannot expect complete regularity of the observations $U_{M}$. Indeed, a
scalar observable of a specific unit $u(\vec{x}_{k}(t))$ contains a regular
part $u_{reg}(t)$ due to a regular forcing of this unit by the regular mean
fields $\vec{Y}(t),\vec{Z}(t)$, and an irregular part $u_{irr,k}(t)$ from the
irregular specific dynamics. Noteworthy is that the regular part is the same
for all oscillators (here, the assumption of global coupling is important!),
while irregular parts can be considered independent. In the asynchronous case,
$u_{reg}(t)=\text{const}$. Now, the finite-sample average can be represented
as

$U_{M}(t)=\frac{1}{M}\sum_{k\in\mathcal{S}_{M}}u_{irr,k}(t)+u_{reg}(t)\;.$ (9)

In finite-sample observations $U_{M}(t)$, finite-size fluctuations due to the
first term on the r.h.s. of Eq. (9) may be significant for small samples.
These fluctuations lead to a non-vanishing variance $\text{Var}(U_{M})$ in the
whole range of parameters, also for parameters where asynchrony occurs. Here,
one generally expects a law of large numbers $\text{Var}(U_{M})\sim M^{-1}$ to
be valid (cf. Refs. Ginzburg and Sompolinsky, 1994; Hansel and Sompolinsky,
1996; Golomb, Hansel, and Mato, 2001), although there might be more complex
dependencies in a vicinity of the synchronization transition.

We suggest eliminating the fluctuating term in Eq. (9) by virtue of the
autocorrelation analysis. We will rely on the irregularity of the local
dynamics, which ensures that the autocorrelations of $u_{irr,k}(t)$ decay to
zero at large time lags. For the observables $U_{M}(t)$ we define the standard
autocorrelation functions (ACF),

$\Gamma[U_{M}](\tau)=\langle(U_{M}(t)-\bar{U}_{M})(U_{M}(t+\tau)-\bar{U}_{M})\rangle\;,$
(10)

where $\bar{U}_{M}$ is the time average of the observable $U_{M}$. Note that
$\Gamma[U_{M}](0)=\text{Var}(U_{M})$. Because the irregular and regular
components of $U_{M}(t)$ can be considered as independent, the ACF (10)
decomposes in two parts 222In the calculation of the ACF it is usual to assume
that the regular component possesses a random phase shift:

$\Gamma[U_{M}](\tau)=\frac{1}{M}\Gamma_{irr}(\tau)+\Gamma_{reg}(\tau)\;,$ (11)

where

$\Gamma_{irr}(\tau)=\langle(u_{irr}(t)-\bar{u}_{irr})(u_{irr}(t+\tau)-\bar{u}_{irr}),$
$\Gamma_{reg}(\tau)=\langle(u_{reg}(t)-\bar{u}_{reg})(u_{reg}(t+\tau)-\bar{u}_{reg})\;.$

Note that in the asynchronous case $\Gamma_{reg}(\tau)=0$.

Now we utilize different time lag dependence of two parts of
$\Gamma[U_{M}](\tau)$: while the value of $\Gamma_{irr}(\tau)$ decays at
$\tau\to\infty$, the regular part $\Gamma_{reg}(t)$ is a periodic or
quasiperiodic function of $\tau$. In the power spectrum, the regular component
corresponds to delta-peaks, while the irregular component corresponds to a
continuous spectral density. A proper way to quantify the intensity of the
regular component in the autocorrelation function $\Gamma(\tau)$ is delivered
by the Wiener formula Wiener (1930)

$W(\Gamma)=\lim_{\mathcal{T}\to\infty}\frac{1}{\mathcal{T}}\int_{0}^{\mathcal{T}}\Gamma^{2}(\tau)\,\mathrm{d}\tau\;.$
(12)

We will call $W$ the Wiener order parameter (WOP). One can easily see that
because the integral $\int_{0}^{\infty}\Gamma_{irr}^{2}(\tau)\mathrm{d}\tau$
converges, only the regular component contributes to $W$. If this regular
component vanishes, as it happens in an asynchronous regime, then the WOP
vanishes as well.

In the modern literature, expression (12) is often referred to as Wiener’s
lemma Kornfeld, Sinai, and Fomin (1982); Queffélec (2010). The following
simple calculation supports its validity. We are interested in the discrete
component of the power spectrum, which can be represented as a sum of delta-
functions

$S_{d}(\Omega)=\sum_{k}A_{k}\delta(\Omega-\Omega_{k}),\qquad\Gamma_{d}(\tau)=\sum_{k}A_{k}\cos(\Omega_{k}\tau)\;.$

Integrating the square of the ACF, we get the Wiener formula

$\lim_{\mathcal{T}\to\infty}\frac{1}{\mathcal{T}}\int_{0}^{\mathcal{T}}\Gamma_{d}^{2}(\tau)\,\mathrm{d}\tau=\frac{1}{2}\sum_{k}A^{2}_{k}\;.$

Some remarks about practical evaluation of the WOP are in order. First,
typically the ACF (10) is calculated from a long time series $0\leq t\leq T$
as an integral

$\Gamma[U_{M}](\tau)=\frac{1}{T-\tau}\int_{0}^{T-\tau}(U_{M}(t)-\bar{U}_{M})(U_{M}(t+\tau)-\bar{U}_{M})\,\mathrm{d}t\;.$
(13)

Thus, an additional parameter, the averaging time $T$ becomes relevant.
Second, formally Eq. (12) requires an integration over an infinite interval of
time lags, so that the contribution of the irregular component
$\Gamma_{irr}(\tau)$ vanishes due to the factor $1/\mathcal{T}$. Practically,
we calculate the WOP by starting integration not at zero, but at a time lag
$\mathcal{T}_{1}$:

$W(\Gamma)=\frac{1}{\mathcal{T}_{2}-\mathcal{T}_{1}}\int_{\mathcal{T}_{1}}^{\mathcal{T}_{2}}\Gamma^{2}(\tau)\,\mathrm{d}\tau\;.$
(14)

One should take $\mathcal{T}_{1}$ large enough so that in the interval
$\mathcal{T}_{1}<\tau<\mathcal{T}_{2}$, the irregular component $\Gamma_{irr}$
is very small. Additionally, one has to take the interval
$\mathcal{T}_{2}-\mathcal{T}_{1}$ much larger than the characteristic period
of the ACF (the corresponding error can be reduced by using an appropriate
windowing). We expect finite-observation errors for finite values of
$T,\mathcal{T}_{2}-\mathcal{T}_{1}$, hampering a perfect disentanglement of
regular and irregular components.

Summarizing, by virtue of the ACF calculations, we eliminate (not completely,
but significantly) the irregular component in the observed field $U_{M}(t)$.
Thus, we expect the Wiener order parameter $W$ value calculated according to
the suggested procedure to be small in the irregular state. Here, only the
finite-sample-induced irregular component is present in $U_{M}(t)$, which
leads to an ACF $\Gamma_{M}(\tau)$ that tends to zero at large time lags
$\tau$ (practically, one observes a decay to some level which depends on the
averaging time $T$ used for calculation of the ACF). Effectively, via the ACF
calculation, we perform an additional averaging of $U_{M}(t)$ in time, which
can be considered an equivalent of enlarging the sample size $M$.

In fact, dependence of the remnant fluctuations of the ACF on the averaging
time $T$ in (13) can be used for a better distinguishing of synchrony and
asynchrony (cf. Refs Ginzburg and Sompolinsky, 1994; Hansel and Sompolinsky,
1996; Golomb, Hansel, and Mato, 2001, which discuss in a similar spirit
exploration of the $M$-dependence). One can calculate the ACF for different
averaging times $T$. Because the remnant fluctuations are expected to be $\sim
T^{-1/2}$, the corresponding WOP value $W$ will scale as $W\sim T^{-1}$. If
one observes such a scaling, then “by extrapolation,” one can attribute $W=0$
and conclude that the observed regime is asynchronous. In contradistinction,
if the values of $W$ saturate starting from some $T$, this indicates a finite
value of $W$. One thus obtains a variant of a $0-1$ test for collective
synchrony.

We stress here that the case $M=1$ can also be included in the consideration.
Here, one does not have any sampling-induced irregularity, and the improvement
via the ACF analysis will be significant only if there is an internal
irregularity of the individual oscillators. We will show below in Sections
III,IV,V the examples of noisy and chaotic individual oscillators, where even
observations from one unit will provide a reasonable characterization of the
synchronization transition. However, observations of a few units deliver a
rather poor performance for regular oscillators in the deterministic Kuramoto
model (Section VI).

###  II.4 Finite-size effects

Above, we assumed that the mean fields are strictly regular, as in the
thermodynamic limit; see Eqs. (4-6). If, instead, a finite ensemble is
explored, then the finite-size fluctuations are inevitable Daido (1990);
Pikovsky and Ruffo (1999); Bertini, Giacomin, and Poquet (2014); Hong  _et
al._ (2015); Gottwald (2017); Peter and Pikovsky (2018); Yue and Gottwald
(2023), especially in the asynchronous state and close to the synchronization
threshold (a strongly synchronous regime can be perfectly regular even for
finite ensembles). If the fields $\vec{Y},\vec{Z}$ in Eqs. (1-3) weakly
fluctuate, the same can be expected for the complete observable $U(t)$. These
basic fluctuations will be complemented with finite-sample fluctuations for
all partial observables $U_{M}(t)$. For us most important is that the
regularity of $\vec{Y},\vec{Z},U$ is no more perfect: while these fields are
perfectly periodic in the thermodynamic limit, for finite $N$ one expects some
diffusion of the collective phase (diffusion constant $\sim N^{-1}$), which
means existence of a finite coherence time of the mean fields $t_{cor}\sim N$.
Thus, the peaks in the spectrum of the mean field fluctuations $U(t)$ are no
more delta-peaks but have finite width $\sim 1/t_{cor}$, and the ACF of $U(t)$
will decay at time $t_{cor}$.

Nevertheless, for large $N$, one still can separate the spectrum into a broad-
band component and narrow peaks; in terms of the ACF $\Gamma(\tau)$, one
should consider an interval of time lags $\mathcal{T}_{1},\mathcal{T}_{2}$
such that $\mathcal{T}_{1}$ is larger than the typical correlation time of the
broad-band component, and $\mathcal{T}_{2}$ is smaller than $t_{cor}$. Then,
we can still use the WOP (12), although the value of $W$ will depend on the
choice of $\mathcal{T}_{1},\mathcal{T}_{2}$. Suppose the main goal is to
explore the transition to synchrony. In that case, we suggest fixing
$\mathcal{T}_{1},\mathcal{T}_{2}$ and changing a parameter (e.g., coupling
strength) to see a difference in WOP between synchronous and asynchronous
regimes. All the arguments above assume large values of $N$; for small
population sizes $N$, the finite-size effects are rather prominent so that the
very notion of synchrony becomes fuzzy. We will illustrate this in Section
III.

###  II.5 Irregular mean fields

The described approach assumes that the mean fields in the thermodynamic limit
are regular (periodic or quasiperiodic) and relies on identifying the
corresponding regular component in the observation by virtue of the ACF
analysis. Here, we discuss the case of irregular mean fields that occur if a
solution of the mean-field equations (4-6) in the thermodynamic limit is
chaotic. Examples from the literature include globally coupled Stuart-Landau
oscillators Hakim and Rappel (1992); Nakagawa and Kuramoto (1994); Chabanol,
Hakim, and Rappel (1997); Ku, Girvan, and Ott (2015); Clusella and Politi
(2019); León and Pazó (2022) and coupled neurons Olmi, Politi, and Torcini
(2011); Pazó and Montbrió (2016); Ratas and Pyragas (2018).

We argue that the presented technique can still be used if the mean fields are
weakly irregular. The idea is basically the same as in the case of finite-size
fluctuations (Section II.4). Indeed, the ACF eventually decreases to zero for
an irregular process, so its power spectrum does not contain discrete
components, and the WOP (12) vanishes. However, suppose the process is close
to a regular one, with narrow peaks in the power spectrum. In that case, the
ACF decays rather slowly at large times (the characteristic time is the
diffusion time of the phase of the nearly-periodic oscillations). As a result,
the calculation of the WOP according to Eq. (14) provides a finite value. If
one chooses $\mathcal{T}_{1}$ in Eq. (14) to be larger than the time of an
initial drop of the ACF (due to a possible broadband component of the power
spectrum), and $\mathcal{T}_{2}$ in Eq. (14) is smaller than the
characteristic decay time of a long-living component of the ACF, then
expression (14) will still provide a reasonable estimate of the intensity of
the nearly regular component of the mean field variations.

If the irregularity of the mean fields is large, and there is no clear time
scale separation between different epochs of the ACF decay, the application of
the ACF analysis seemingly brings no advantage compared to the calculation of
the variance of the observed mean field.

###  II.6 Non-identical units

Above we concentrated on the case of identical globally coupled oscillators;
here we shortly argue, that the approach can also be applied to non-identical
units, provided they are irregular enough. A typical way to introduce
inhomogeneity in the population is to assume that the dynamics (1) depends
additionally on a parameter $b_{k}$, and this parameter is different for
different oscillators. In the thermodynamic limit, one describes the
inhomogeneity with a distribution $g(b)$. As a result, the probability density
contains $b$ as an additional variable, $w(\vec{x},t|b)$ and averaging in (6)
contains additional averaging over the distribution of $b$. With such a
modification, the basic conclusions of Sections II.1-II.3 are still valid.
Namely, at synchrony the mean fields oscillate regulary, while at asynchrony
they are constants. Each unit has a regular and irregular component, and the
former can be identified via the ACF analysis. The only modification is that
now sampled observations (Section II.2) are not sample-independent, because
the units are not identical. As a result, in Eq. (11), both components will
depend on the sample $\mathcal{S}_{M}$. Nevertheless, the separation of time
scales will be still valid, and the Wiener’s formula (12) will provide the
intensity of the regular component. This consideration shows that the WOP will
now depend on the sample $\mathcal{S}_{M}$. In particular, we expect to have
different values if just one unit is observed. Nevertheless, the contrast
between the asynchronous and synchronous states will be still present for any
sample $\mathcal{S}_{M}$.

Below, we apply the described approach to several examples. One of the main
questions is: how well can we characterize synchrony with an observation of
just one unit? We will start with strongly irregular units in Sections III,IV.
In Section V, we will explore identical noisy limit-cycle oscillators, also in
the case of a weak noise. In the last example (Section VI) we will apply our
approach to the Kuramoto model of non-identical phase oscillators.

##  III Noise-induced oscillations

A recent paper Klinshov  _et al._ (2021) studied the dynamics of globally
coupled active rotators. The equations for the units described by the phase-
like variables $\phi$ are formulated as

$\dot{\phi}_{k}=1-b\sin\phi_{k}+\frac{\varepsilon}{N}\sum_{j}\sin(\phi_{j}-\phi_{k})+\sqrt{2\sigma^{2}}\xi_{k}(t)\;.$
(15)

With the choice $b>1$, adopted in Ref. Klinshov  _et al._ , 2021, the
uncoupled noise-free systems have a stable steady state and do not oscillate.
In the presence of noise, the dynamics of uncoupled units are purely noise-
induced and thus highly irregular. The coupling term $\sim\varepsilon$ is of
Kuramoto type; it can be written in terms of mean fields

$C(t)=\frac{1}{N}\sum_{k}\cos\phi_{k}(t),\qquad
S(t)=\frac{1}{N}\sum_{k}\sin\phi_{k}(t)\;.$ (16)

In the thermodynamic limit, the dynamics of the probability density
$w(\phi,t)$ reduces to a nonlinear Fokker-Planck equation of type (4). The
regimes in this nonlinear PDE have been carefully studied in Ref. Klinshov
_et al._ , 2021, and a domain of parameters where the density $w(\phi,t)$
varies periodically in time, has been identified (together with other domains,
where the attractor in the nonlinear PDE is a steady state).

Below, we fix the parameters $\varepsilon=0.6$ and $b=1.025$. For these
values, the phase diagram elaborated in Ref. Klinshov  _et al._ , 2021
predicts a time-periodic (synchronous) solution for $\sigma^{2}=0.04$ and a
stationary (asynchronous) solution at $\sigma^{2}=0.06$.

###  III.1 Large number of active rotators

Figure 1: Synchronous and asynchronous states of noise-driven active rotators
(15). Presented are the “phase planes” of observables $C_{M}(t),S_{M}(t)$, see
Eq. (16), for different sample sizes $M$ (the total size of the ensemble is
$N=10^{4}$). Upper row (panels (a-e)): $\sigma^{2}=0.06$ (asynchrony), bottom
row (panels (f-k)): $\sigma^{2}=0.04$ (synchrony). The scales in all panels
are the same, $-1\leq C_{M},S_{M}\leq 1$.

Figure 1 illustrates how these solutions for an ensemble of $N=10^{4}$
oscillators are represented by finite samples. Here, we depict mean fields
$C_{M}$, $S_{M}$ calculated following Eqs. (16,7) over randomly chosen subsets
of $M$ units, starting with $M=1$ (just one active rotator is observed) up to
$M=N=10^{4}$, what corresponds to full mean fields (16). In the projections on
the partial mean field phase planes $(C_{M},S_{M})$ shown in Fig. 1, one can
see for $M\geq 100$ and $\sigma^{2}=0.04$ a “topological limit cycle”, where
the trajectory rotates with a void in the center. In contradistinction, for
$M=10$ and $\sigma^{2}=0.04$, and all $M$ and $\sigma^{2}=0.06$ one does not
see a topological limit cycle (such a cycle is, of course, present for $M=1$
because the observables defined in Eq. (16) fulfill $C_{1}^{2}+S_{1}^{2}=1$).
We note here that the existence of a void in a two-dimensional projection of
the mean fields has been suggested in Ref. Temirbayev  _et al._ , 2012 as a
practical criterion for synchrony in small populations (cf. a more detailed
analysis of this criterion for the Kuramoto model in Ref. Peter and Pikovsky,
2018).

Figure 2: Autocorrelation functions of different sampled observables
$C_{M}(t)$, calculated using different time intervals $T$. Left column (panels
(a,c,e)): small noise $\sigma^{2}=0.04$; right column (panels (b,d,f)): large
noise $\sigma^{2}=0.06$. Notice that the vertical scales are different in all
panels.

The autocorrelation functions of observables $C_{M}(t)$ calculated according
to Eq. (13) are shown in Fig. 2. Here we present only small time lags
$0\leq\tau<200$, to illustrate convergence with the averaging time $T$. One
can see that all the autocorrelation functions $\Gamma_{M}(t)$ averaged over
$T=10^{5}$ and $T=10^{6}$ practically coincide, while some deviations are seen
for the short averaging time $T=10^{4}$.

Figure 3: The ACFs for the synchronous case $\sigma^{2}=0.04$, $N=10^{4}$, for
different sample sizes $M$ at large time lags. The averaging time is
$T=10^{6}$. All the curves with $10\leq M\leq 10^{4}$ practically coincide.
Figure 4: The ACFs for the asynchronous case $\sigma^{2}=0.06$, $N=10^{4}$,
for different sample sizes $M$ at large time lags. The averaging time is
$T=10^{6}$. Panel (b) shows the ACFs for $M\geq 100$ in the appropriate scale.

Our primary interest is in the behavior of the ACFs at large time lags; we
show these data in Figs. 3,4 for $1500\leq\tau\leq 2000$. In the synchronous
case (Fig. 3), the ACFs for different sample sizes $M$ practically coincide if
the averaging time $T$ is large. This confirms the concept of section II that
the regularity of the mean fields can be extracted from samples with any $M$,
even with $M=1$. Another observation is that the amplitude of the ACF slowly
decreases. The reason for this are finite-size fluctuations for $N=10^{4}$, as
discussed in Section II.4.

The ACFs in the asynchronous case (Fig. 4) have relatively small values and
demonstrate strong dependence on the sample size $M$ (thus, we show two
panels, omitting cases $M=1,10$ in panel (b)). A remarkable observation is
that at a presented time interval of averaging $T=10^{6}$, the correlations of
the fields with $M=1000$ and $M=N=10^{4}$ are rather regular, although small
in amplitude. This issue needs further investigation.

Figure 5: The Wiener order parameter $W$ for different sample sizes $M$ and
ACF averaging times $T$ (other parameters: $\mathcal{T}_{1}=750$,
$\mathcal{T}_{2}=2000$). Panel (a): synchronous case $\sigma^{2}=0.04$,
$N=10^{4}$. Note that the vertical scale is linear. Panel (b): asynchronous
case $\sigma=0.06$, $N=10^{4}$; the verical scale is logarithmic. The black
dashed line shows the scaling $W\sim T^{-1}$.

Finally, we compare the WOPs $W$ for the two cases. In Fig. 5, we show $W$ for
different sample sizes $M$ and different averaging time intervals $T$. As is
expected from the discussion in Section II.3, dependence on $T$ is obviously
different for the two cases. In the synchronous case (panel (a)), $W$
saturates at some value for $T\geq 10^{5}$, so this value can be taken as a
correct estimate of $W$. Furthermore, this value is practically the same for
all $M$, as the coincidence of the ACFs in Fig. 3 suggests. On the contrary,
for the asynchronous case, the values of $W$ decay $\sim T^{-1}$ at large $T$.
While for purely fluctuating samples one would also expect $W\sim M^{-2}$;
such scaling is valid only for $M=1,10$; for $M\geq 100$, apparently, this
scaling is not applicable due to the relative regularity of fluctuations of
ACFs.

###  III.2 Finite size effects for small ensembles of active rotators

In Section III.1 above, we considered a rather large ensemble of irregular
oscillators with $N=10^{4}$. Here, we explore, for the same system, how the
discussed methods of synchrony characterization work for small population
sizes. In particular, we explore cases $N=10,\;100,\;1000$.

First, in Fig. 6 we present the phase portraits $(C_{N}(t),S_{N}(t))$ for the
small populations, for $\sigma^{2}=0.04$ and $\sigma^{2}=0.06$. One can see a
clear difference between two noise values for $N=1000$, while for $N=10,100$
for both noise levels, there is no void in the phase portraits around the
center of rotations. This allows for a preliminary conclusion that synchrony
detection with other methods will possibly work for $N=1000$, but might be
hard for $N=10,100$.

Figure 6: Phase portraits of the complete mean fields $C_{N}$ vs $S_{N}$ for
small system sizes $N=10,\;100,\;1000$ and two values of the noise level.
Figure 7: Characterization of the regimes depicted in Fig. 6 (c,f), $N=1000$,
with the ACFs. Both for small noise (synchrony in the thermodynamic limit),
$\sigma^{2}=0.04$, and for large noise (asynchrony in the thermodynamic
limit), $\sigma^{2}=0.06$, the ACFs from one element and from the complete
mean field coincide after a small initial interval of time lags. Figure 8: The
same as Fig. 5, but for the case $N=1000$. One can see that the WOP reasonably
contrasts synchronous and asynchronous states despite pronounced finite-size
irregularity. Parameters of the WOP calculation: $\mathcal{T}_{1}=200$,
$\mathcal{T}_{2}=600$.

In Fig. 7, we show different autocorrelation functions for $N=1000$. First, we
notice that the ACFs for $M=1$ and $M=1000$ practically coincide except for a
region of very small time lags (one cannot see the red and the green curves in
Fig. 7 because these curves are overlapped with the blue and the cyan ones).
Altogether, the different time behavior of the ACFs delivers a strong contrast
between synchrony and asynchrony for $N=1000$ (because the blue curve for
$\sigma^{2}=0.04$ demonstrates a much more ordered oscillatory tail compared
to the cyan curve for $\sigma^{2}=0.06$). On the other hand, since even in the
synchronous case, the ACF significantly decays after 10 periods, for the
calculation of the WOP $W$ we have to choose relatively small values
$\mathcal{T}_{1},\mathcal{T}_{2}$. With such a choice, this order parameter
allows for a reasonable synchrony characterization, as illustrated in Fig. 8.

Performing the same ACF analysis for small population sizes, we found that for
$N=100$ the difference in the correlation lengths for the two values on noise
strength is rather small, and for $N=10$ there is practically no difference.
These results confirm a general expectation in Section II that a transition to
synchrony for small population sizes becomes fuzzy.

##  IV Chaotic oscillators

Figure 9: Time series of the global mean field $X=N^{-1}\sum_{k}x_{k}(t)$ for
coupled chaotic Rössler oscillators (17) for different values of the coupling
strength $\varepsilon$. The ensemble size is $N=10^{4}$.

Here we describe synchronous states in coupled Rössler oscillators, first
reported in Ref. Pikovsky, Rosenblum, and Kurths, 1996:

$\displaystyle\dot{x}_{k}$
$\displaystyle=-y_{k}-z_{k}+\frac{\varepsilon}{N}\sum_{j}x_{j}\;,$ (17)
$\displaystyle\dot{y}_{k}$ $\displaystyle=x_{k}+ay_{k}\;,$
$\displaystyle\dot{z}_{k}$ $\displaystyle=b+(x_{k}-c)z_{k}\;.$

We use parameters $a=0.25$, $b=0.4$, $c=8.5$, at which each Rössler system has
a funnel-type attractor, i.e., the oscillators have an ill-defined phase (cf.
phase portraits with $M=1$ in Fig. 10 below). For this system we use the
observables

$X_{M}(t)=\frac{1}{M}\sum_{k\in\mathcal{S}_{M}}x_{k}(t)\;,\qquad
Y_{M}=\frac{1}{M}\sum_{k\in\mathcal{S}_{M}}y_{k}(t)\;.$

With the increase of the coupling strength $\varepsilon$, first periodic
oscillations of the mean field appear, then the amplitude of these
oscillations grows, and they become modulated. We illustrate this in Fig 9,
where we show the time series of the observable $X_{N}(t)$ for $N=10^{4}$ and
different values of $\varepsilon$. For $\varepsilon=0.01$, the mean field
fluctuates around a constant without any visible regularity; we attribute the
fluctuations to finite-size effects and denote this state as an asynchronous
one. In contrast, we observe regular macroscopic variations of the mean field
for $\varepsilon=0.04$ and $\varepsilon=0.2$ and correspondingly denote these
states as synchronous.

Figure 10: Phase portraits $(X_{M},Y_{M})$ for $M=1,10,100,10^{3},10^{4}$, for
synchronous and asynchronous regimes (from top to bottom: $\varepsilon=0.01$,
$\varepsilon=0.04$, $\varepsilon=0.2$).

In Fig. 10 we show trajectories of partial mean fields $(X_{M},Y_{M})$ for
$M=1,10,100,10^{3},10^{4}$, for synchronous and asynchronous regimes. One can
clearly see that these observables show very similar patterns for $M=1$ and
$M=10$ for all $\varepsilon$. At $M=100$ one can recognize a regularity in the
partial mean fields for $\varepsilon=0.2$, while a regularity for
$\varepsilon=0.04$ is only visible for $M\geq 10^{3}$.

Figure 11: Autocorrelation functions at small time lags (panels (a-c)) and at
large time lags (panels (d-f)) for different $M$ and $\varepsilon$ (averaging
time $T=10^{6}$). Panels (a,d): asynchronous state $\varepsilon=0.01$; here,
all the autocorrelation functions at large time lags fluctuate around zero.
Panels (b,e): $\varepsilon=0.04$; panels (c,f): $\varepsilon=0.2$. Here, the
ACFs for observables with different $M$ nearly coincide at large time lags,
where they demonstrate periodic behavior.

Finally, we apply the ACF analysis to reveal synchrony from observations with
small sample sizes $M$. In Fig. 11 we show the ACFs $\Gamma[X_{M}](t)$ for the
observables $X_{M}$, for all the three selected values of the coupling
strength $\varepsilon$. For $\varepsilon=0.01$ (panels (a,d)), the ACFs for
all $M$ fluctuate at a small level for large time lags, confirming the absence
of synchrony. For the synchronous cases $\varepsilon=0.04$ (panels (b,e)) and
$\varepsilon=0.2$ (panels (c,f)), at large time lags the autocorrelation
functions for all $1\leq M\leq N$ nearly coincide. This supports the
conclusion that observation of just one unit in a large population allows for
revealing synchrony by virtue of our method.

##  V Noisy limit-cycle oscillators

In this section, we present an example where we vary the irregularity level at
each unit to see the impact of this level on the performance of the method. We
take the system of $N$ Stuart-Landau oscillators with global nonlinear
coupling and add independent Gaussian white noise terms $\sim\sigma$ to each
unit:

$\dot{a}_{k}=(1+5\mathrm{i})a_{k}-|a_{k}|^{2}a_{k}+\varepsilon
A-\varepsilon_{nl}|A|^{2}A+\sigma\xi_{k}(t)\;.$ (18)

Here $a_{k}$ are complex variables, $k=1,2,\ldots,N$, and
$A=N^{-1}\sum_{k=1}^{N}a_{k}$ is the complex mean field. The terms
$\xi_{k}(t)$ are Gaussian white noises with zero mean and unit intensity
333Notice that the normalization of the noise strength differs from that used
in Eq. (15). Analysis of this system performed in Ref. Rosenblum and Pikovsky,
2015 without noisy perturbations shows that for specific parameter values the
ensemble exhibits partial synchronization, such that the mean field is
periodic, but its frequency differs from that of the oscillators and the
latter are thus quasiperiodic.

###  V.1 Linear coupling

Here, we fix $\varepsilon_{nl}=0$, $\sigma=0.2$, and vary real parameter
$\varepsilon$ to trace the synchronization transition. For this purpose, we
exploit two traditional measures (the mean-field variance and the Kuramoto
order parameter), as well as the WOP. We simulate the ensemble of $N=10^{4}$
units and take $U_{M}=M^{-1}\text{Re}\big{(}\sum_{k=1}^{M}a_{k}\big{)}$,
$M\leq N$, as our observable. (We use $10^{6}$ time points sampled with the
step $0.05$ to compute the autocorrelation function and fix the values
$\mathcal{T}_{1}=1.25\cdot 10^{4}$ and $\mathcal{T}_{2}=2\cdot 10^{4}$ for the
computation of $W$.) The results are illustrated in Fig. 12. Since $W$ is
proportional to the squared variance of the process, we plot
$[\text{Var}(U_{M})]^{2}$ in the panel (a) and $W$ in the panel (b). (We do
not show the KOP since its variation is, up to a nearly constant factor, very
similar to that of the variance measure.)

The main outcome is that the Wiener order parameter successfully traces the
synchronization transition from an observation of a single oscillator, while
the variance measure certainly fails here. Furthermore, for partial
observations with $M>1$, when both measures work, the WOP is more efficient:
the contrast between the asynchronous and synchronous states is much stronger.

Figure 12: The variance measure of synchrony (a) and the Wiener order
parameter (b) of the ensemble of noisy Stuart-Landau oscillators (18) with
global linear coupling as functions of the coupling strength $\varepsilon$.
Different curves in each panel correspond to observation of $M=1$,
$M=10,\ldots$, $M=N=10^{4}$ oscillators. Notice that $W$ successfully reveals
the synchronization transition even from an observation of a single
oscillator. Next, for the $M>1$ cases, the quantification with the WOP
provides an about two orders of magnitude stronger contrast between the states
before and after transition.

###  V.2 Nonlinear coupling

Here, we fix $\varepsilon=3+0.3\mathrm{i}$ and $\varepsilon_{nl}=8$. With this
nonlinear coupling, the noise-free ensemble in the thermodynamic limit
$N\to\infty$ exhibits the harmonic mean-field solution, while individual units
remain not locked to the field and thus dynamics of each unit is
quasiperiodic. This quasiperiodicity will be naturally seen in the one-unit
sampling $M=1$, but will be smeared for larger $M$.

Figure 13: The upper envelopes of the correlation functions for the
nonlinearly coupled ensemble of Stuart-Landau oscillators (18), for different
values of the noise intensity: $\sigma=0.02$ (a), $\sigma=0.1$ (b), and
$\sigma=0.5$ (c) and different values of $M$.

Figure 13 demonstrates the simulation results for $N=10^{4}$ and different
values of the noise intensity. Here, for a better visibility, we plot the
upper envelopes of the autocorrelation functions; thus, a nearly constant
value of $\bar{\Gamma}$ at large $\tau$ indicates for a purely periodic
correlation function with the corresponding amplitude of oscillations.

The main conclusion is that regularity of the oscillators is an obstacle for
observations with $M=1$: for a weak noise, $\sigma=0.02$ (panel (a)), the ACF
envelope for $M=1$ strongly oscillates, indicating modulation of the ACF, as
expected for a signal that is close to a quasiperiodic one. Correspondingly,
computation of $W$ from one oscillator will provide the sum of intensities of
periodic components on the oscillator’s frequency and that of the mean field
(and their harmonics). For this level of noise, oscillations of the ACF
envelope can be also seen for $M=10$, but they are much smaller. So, due to
relative regularity, the correlation decay is very slow; hence, inference of
synchrony from observation of one or several units requires an extremely long
observation.

The situation changes for an intermediate noise, $\sigma=0.1$ (panel (b)).
Now, we see that the oscillation of the ACF envelope for $M=10$ decays quickly
so that the partial mean field $U_{10}(t)$ suffices for the computation of
$W$. For stronger noise, $\sigma=0.5$ (panel (c)), even observation of one
oscillator (after a short decay time) yields the same ACF as the complete
observation of the mean field. This confirms the overall estimation that the
proposed method should work better for irregular oscillators and might not be
optimal when local units are highly regular.

Two remarks are in order. First, we plotted only the envelopes of the ACFs,
but we checked that when they coincided, the ACFs coincided as well. Second,
in panel (c), we see that the envelope is not constant but slightly decays so
that, strictly speaking, the ACFs are not stationary for the presented range
of $\tau$. This decay is because the mean field is not truly periodic due to
the large though finite ensemble size, as discussed in Section II.4.

###  V.3 Nonidentical noisy oscillators

To illustrate ideas of Section II.6 about the approach’s applicability to non-
identical units, we take a large ensemble of $N=10^{4}$ slightly inhomogeneous
noisy oscillators. The model now reads:

$\dot{a}_{k}=(1+\mathrm{i}\omega_{k})a_{k}-|a_{k}|^{2}a_{k}+\varepsilon
A+\sigma\xi_{k}\;,$ (19)

where frequencies $\omega_{k}$ are taken from a Gaussian distribution with the
mean value $\bar{\omega}=5$ and standard deviation $0.02$; the noise strength
is $\sigma=0.2$, and the coupling strength is varied. As before, we use
$U_{M}=M^{-1}\text{Re}\big{(}\sum_{k=1}^{M}a_{k}\big{)}$ to compute the WOP
$W$ as a function of $\varepsilon$. Figure 14 demonstrates the synchronization
transition revealed from the global mean field, $M=N=10^{4}$, and observations
of a small subpopulation, $M=20$. Since the units are non-identical, the
results definitely depend on the chosen subset of units. To check that this
dependence is weak, we performed calculations with $50$ random samples of
$M=20$ units and presented the results for all these samples as red dots. We
see that all the samples give very similar results. Thus, an observation of
only 2% of the oscillators reliably reveals the transition to synchrony.

Figure 14: The Wiener order parameter $W$ as a function of the coupling
parameter $\varepsilon$ for the system of $N=10^{4}$ inhomegeneous noisy
oscillators. Blue line shows the values obtained from the mean field of the
ensemble. Red dots present the results of partial observations, i.e., the
partial mean field is computed from 20 randomly chosen units. For each value
of $\varepsilon$, we perform 50 trials which provide generally different
results. Since the units are non-identical, the order parameter computed from
a partial observation depends on the chosen subset of units, what explain the
vertical scattering of the points. Notice that the scattering decreases with
$\varepsilon$ because of the tendency to synchrony.

##  VI Kuramoto model

Here we apply our approach to the standard Kuramoto model of globally coupled
phase oscillators

$\dot{\varphi}_{k}=\omega_{k}+\frac{\varepsilon}{N}\sum_{j=1}^{N}\sin(\varphi_{j}-\varphi_{k}),\quad
k=1,2,\ldots,N\;.$ (20)

The main differences to the cases above are: (i) the specific units, if
uncoupled, are not irregular (chaotic or noisy), but regular periodic
oscillators; (ii) the units are not identical but differ in natural
frequencies $\omega_{k}$. The latter property ensures the existence of an
asynchrony state for small coupling strengths $\varepsilon$.

In simulations below we take a Gaussian distribution of frequencies
$N(\bar{\omega},1)$ with unit variance centered at frequency
$\bar{\omega}=2\pi$. The theory Peter and Pikovsky (2018) predicts in the
thermodynamic limit a transition to periodic global oscillations at
$\varepsilon_{cr}=2\sqrt{2/\pi}\approx 1.596$; the dependence of the Kuramoto
order parameter (KOP) $R=|\langle e^{\mathrm{i}\varphi}\rangle|$ on
$\varepsilon$ can be represented in the parametric form with real-valued
parameter $p>0$ as

$R=e^{-p}\sqrt{\pi p/2}[I_{0}(p)+I_{1}(p)],\quad\varepsilon=2\sqrt{p}/R\;,$
(21)

where $I_{0},I_{1}$ are modified Bessel function of the first kind.

Below, we consider finite ensembles with moderate numbers of units $N$. We
will use regular sampling of natural frequencies $\omega_{k}$, like in Refs.
Carlu, Ginelli, and Politi, 2018; Peter and Pikovsky, 2018, calculating them
from the corresponding quantiles of the normal distribution. Because the
analytic solution (21) in the thermodynamic limit is available, we will
compare the numerical results with the theory. For this, we need to relate the
KOP to the WOP, see Eqs. (12,14). We will use the local observable
$u(\varphi)=\cos\varphi$ for the latter. Thus, the fully observed (i.e.,
calculated for the full ensemble size $N$) mean field $U(t)=\text{Re}(Z(t))$
is the real part of the complex Kuramoto order parameter
$Z(t)=N^{-1}\sum_{k}e^{\mathrm{i}\varphi_{k}(t)}$. In the thermodynamic limit,
the field $Z$ is periodic $Z=R\exp[\mathrm{i}\bar{\omega}t+\theta_{0}]$, where
$R$ is obtained according to (21). Thus,
$U(t)=R\cos(\bar{\omega}t+\theta_{0})$. The ACF function of this periodic
process is, according to Eq. (13), $\Gamma(t)=(R^{2}/2)\cos\bar{\omega}t$.
Thus, the WOP calculated according to Eq. (12), is $W=R^{4}/8$. This value
should be compared with the variance $\text{Var}(U)=\Gamma(0)=R^{2}/2$ of the
process $U(t)$, which is used in the traditional characterization of the
synchronization transition. Therefore, in the plots below, we use a function
of the variance $V=[\text{Var}(U)]^{2}/2$ because, in the thermodynamic limit,
it coincides with WOP $W$.

Figure 15: The variance-based order parameter $V$ (solid curves) and the WOP
$W$ (short-dashed curves with the corresponding colors), both for the complete
observable $U(t)$, for the Kuramoto model (20) with $N=20,50,100$. The long-
dashed black line shows the theoretical prediction in the thermodynamic limit
(21). One can see that the contrast between the synchronous and asynchronous
states increases, if WOP is used, by several orders for $N=100$, and by more
than one order for $N=20$. Parameters of the WOP calculation: $T=10^{5}$,
$\mathcal{T}_{1}=250$, $\mathcal{T}_{2}=500$.

We start with full global observables $U(t)$ of a population of Kuramoto
oscillators and report the order parameters close to the synchronization
transition in Fig. 15, for system sizes $N=20,\,50,\,100$. One can see that
the “sharpness” of the transition to synchrony is significantly improved if
the WOP is used instead of the variance-based characterization. Remarkably,
the WOP decreases prior to the transition threshold for a small number of
units $N=20$ (this effect is also present for $N=50$ but not so pronounced).
We attribute this to the chaoticity of the Kuramoto model, studied in Refs.
Popovych, Maistrenko, and Tass, 2005; Maistrenko, Popovych, and Tass, 2005;
Carlu, Ginelli, and Politi, 2018. In these papers, it has been shown that at
finite coupling strengths prior to the synchronization transition, the largest
Lyapunov exponent in the Kuramoto model is positive $\sim N^{-1}$ and reaches
a maximal value close to the transition. Because the chaoticity of the
population makes the observed mean field $U(t)$ irregular, its ACF decreases
faster for stronger chaos, and the WOP is smaller. This explains why, for
small $N$, the WOP has a minimum prior to the transition to synchrony.

Figure 16: Analysis of the synchronization transition for $N=50$ Kuramoto
oscillators via observations of $M=20$ (a), $M=10$ (b), $M=5$ (c), and $M=2$
(d) units. For each case, together with the variance-based parameter $V$
(based on the full observation; solid red lines) and the WOP $W$ (dashed blue
lines) calculated from the full mean field, we show WOPs for 50 random samples
of $M$ selected for observation units with grey lines. The black long-dashed
line is the theoretical prediction for the thermodynamic limit. Parameters of
the WOP calculation: $T=10^{5}$, $\mathcal{T}_{=}50$, $\mathcal{T}_{2}=100$.

Next, we explore in Fig. 16 how the synchronization transition can be
characterized with finite samples from the ensemble of $N=50$ units. Because
the oscillators are not identical (they have different natural frequencies),
their behavior at the transition is different. The oscillators with the
natural frequencies close to the middle one become locked by the mean field,
while those with frequencies at the tails of the distribution remain nearly
quasiperiodic (although they possess a periodic component from the mean field
driving). Thus, we considered not just one sample of $M$ units to calculate
the observed mean field $U_{M}(t)$ but performed calculations with a set of 50
randomly chosen samples. The results for all these 50 samples are shown in
Fig. 16 with grey lines to illustrate the diversity of the transition
characterizations due to the units’ non-identity. One can see that some, but
not all, cases with $M=2$ (panel (d)) show a reasonable contrast between
asynchrony and synchrony. On the contrary, already an observation of $M=5$
units (panel (c)) provides a reliable description of the transition. The
quality grows with $M$, and at $M=10$ (panel (b)), the quality of WOP is
almost the same as the quality of the variance-based characterization where
all units are used. For $M=20$ (panel (a)), the results from partial samplings
are close to those with the full observability $M=50$.

One remark is in order. Because the Kuramoto model (20) is invariant with
respect to time-dependent phase shifts $\varphi\to\varphi+\nu t$, one often
performs such a transformation with $\nu=\bar{\omega}$. As a result, in the
transformed phases, the mean fields are not periodic but time-independent. In
terms of the original oscillators, Eqs. (1-3), this corresponds to mean fields
in definitions of which (see Eq. (3)) an explicit periodic time dependence is
present. Thus, the example of the Kuramoto system does not contradict the
general considerations of Section II, where observables that are not
explicitly time-dependent are assumed.

The conclusion from this example is that although the individual oscillators
are regular, the mean fields and the fields from the samples with sufficiently
large size $M$ fluctuate in the asynchronous state and have a strong periodic
component in the synchronous regime. This property allows for a nice
separation of synchrony and asynchrony by virtue of the WOP. However, in
contradistinction to the case of local irregular units, the efficacy of the
WOP approach drops at small sample sizes $M$.

##  VII Conclusion

To summarize, we suggested an approach of synchrony quantification based on
the characterizing regularity of the mean fields in a population of globally
coupled units. The main difference to previous works is that we do not simply
identify the existence of macroscopic mean field oscillations by following the
variance but select the regular part of these oscillations by virtue of
Wiener’s formula applied to the auto-correlation function. We coin the term
“Wiener order parameter” for the corresponding quantity. As we show with
different examples, this approach strongly increases the contrast between the
synchronous and asynchronous states. Most successful is the application of
this approach to ensembles of highly irregular units, like excitable systems
exhibiting noise-induced oscillation (Section III) and chaotic oscillators
(Section IV). Furthermore, because in our method, we effectively separate the
regularity and irregularity of the observations by means of the ACF
calculations, there is, in fact, no need for ensemble averaging. This allows
for working with small samples, and in most cases, even observing one unit out
of the large ensemble allows for a reliable quantification of the
synchronization transition.

In our approach, we operate with a generic local observable that represents
the system’s state. We stress that this observable does not need to be the
phase of the oscillations, and there is no need to extract the phase from this
observable. Of course, in the cases where the phase-dynamics description is
adequate, knowing the phases would lead to an enhanced characterization of the
synchronization transition by including properties of frequency entrainment,
etc. The basic requirement for the observable is that it is suitable for
calculations of the ACF. This means that sparse in time observables, like
spikes resulting in a point process, require a special consideration to be
reported elsewhere.

While the description of synchrony/asynchrony contrast as that of
regularity/irregularity is valid in many circumstances, there are situations
where synchronous states show a large degree of irregularity (even in the
thermodynamic limit). Irregularity of the mean fields can also result from a
common noisy forcing. In such situations, the method should be significantly
modified; this is a work in progress.

###### Acknowledgements.

We thank O. Omelchenko and K. Krischer for fruitful discussions.

## Data availability

All numerical experiments are described in the paper and can be reproduced
without additional information.

## References

  * Winfree (1967) A. T. Winfree, “Biological rhythms and the behavior of populations of coupled oscillators,” J. Theor. Biol. 16, 15–42 (1967).
  * Kuramoto (1975) 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.
  * Strogatz and Stewart (1993) S. H. Strogatz and I. Stewart, “Coupled oscillators and biological synchronization,” Scientific American , 68–75 (1993).
  * Acebrón  _et al._ (1998) J. A. Acebrón, L. L. Bonilla, S. D. Leo, and R. Spigler, “Breaking the symmetry in bimodal frequency distributions of globally coupled oscillators,” Phys. Rev. E 57, 5287–5290 (1998).
  * Pikovsky, Rosenblum, and Kurths (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization. A Universal Concept in Nonlinear Sciences._ (Cambridge University Press, Cambridge, 2001).
  * 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).
  * Glass and Mackey (1988) L. Glass and M. C. Mackey,  _From Clocks to Chaos: The Rhythms of Life._ (Princeton Univ. Press, Princeton, NJ, 1988).
  * Buzsáki (2006) G. Buzsáki,  _Rhythms of the brain_ (Oxford UP, Oxford, 2006).
  * Lehnertz (2008) K. Lehnertz, “Epilepsy and nonlinear dynamics,” Journal of Biological Physics 34, 253–266 (2008).
  * Daffertshofer and Pietras (2020) A. Daffertshofer and B. Pietras, “Phase synchronization in neural systems,” Synergetics , 221–233 (2020).
  * Ginzburg and Sompolinsky (1994) I. Ginzburg and H. Sompolinsky, “Theory of correlations in stochastic neural networks,” Physical Review E 50, 3171 (1994).
  * Hansel and Sompolinsky (1996) D. Hansel and H. Sompolinsky, “Chaos and synchrony in a model of a hypercolumn in visual cortex,” Journal of Computational Neuroscience 3, 7–34 (1996).
  * Pikovsky, Rosenblum, and Kurths (1996) A. Pikovsky, M. Rosenblum, and J. Kurths, “Synchronization in a population of globally coupled chaotic oscillators,” Europhys. Lett. 34, 165–170 (1996).
  * Golomb, Hansel, and Mato (2001) D. Golomb, D. Hansel, and G. Mato, “Mechanisms of synchrony of neural activity in large networks,” in  _Handbook of Biological Physics_ , Vol. 4 (Elsevier, 2001) pp. 887–968.
  * Kapitaniak  _et al._ (2012) M. Kapitaniak, K. Czolczyński, P. Perlikowski, A. Stefański, and T. Kapitaniak, “Synchronization of clocks,” Physics Reports 517, 1–69 (2012).
  * Czolczyński  _et al._ (2013) K. Czolczyński, P. Perlikowski, A. Stefański, and T. Kapitaniak, “Synchronization of the self-excited pendula suspended on the vertically displacing beam,” Communications in Nonlinear Science and Numerical Simulation 18, 386 – 400 (2013).
  * Wiesenfeld, Colet, and Strogatz (1998) K. Wiesenfeld, P. Colet, and S. Strogatz, “Frequency locking in Josephson arrays: Connection with the Kuramoto model,” Physical Review E 57, 1563–1569 (1998).
  * Temirbayev  _et al._ (2012) A. A. Temirbayev, Z. Z. Zhanabaev, S. B. Tarasov, V. I. Ponomarenko, and M. Rosenblum, “Experiments on oscillator ensembles with global nonlinear coupling,” Phys. Rev. E 85, 015204 (2012).
  * Wenzel and Hamm (2021) M. Wenzel and J. P. Hamm, “Identification and quantification of neuronal ensembles in optical imaging experiments,” Journal of Neuroscience Methods 351, 109046 (2021).
  * Note (1) In this paper we assume ergodicity and stationarity, so practically the averaging is performed over time.
  * Note (2) In the calculation of the ACF it is usual to assume that the regular component possesses a random phase shift.
  * Wiener (1930) N. Wiener, “Generalized harmonic analysis,” Acta mathematica 55, 117–258 (1930).
  * Kornfeld, Sinai, and Fomin (1982) I. Kornfeld, Y. G. Sinai, and S. Fomin,  _Ergodic theory, Nauka, Moscow, 1980; English transl., Springer, New York_ (1982).
  * Queffélec (2010) M. Queffélec,  _Substitution dynamical systems-spectral analysis_ , Lecture Notes in Mathematics, Vol. 1294 (Springer, 2010).
  * Daido (1990) H. Daido, “Intrinsic fluctuations and a phase transition in a class of large population of interacting oscillators,” J. Stat. Phys. 60, 753–800 (1990).
  * Pikovsky and Ruffo (1999) A. Pikovsky and S. Ruffo, “Finite-size effects in a population of interacting oscillators,” Phys. Rev. E 59, 1633–1636 (1999).
  * Bertini, Giacomin, and Poquet (2014) L. Bertini, G. Giacomin, and C. Poquet, “Synchronization and random long time dynamics for mean-field plane rotators,” Probability Theory and Related Fields 160, 593–653 (2014).
  * Hong  _et al._ (2015) H. Hong, H. Chaté, L.-H. Tang, and H. Park, “Finite-size scaling, dynamic fluctuations, and hyperscaling relation in the Kuramoto model,” Phys. Rev. E 92, 022122 (2015).
  * Gottwald (2017) G. A. Gottwald, “Finite-size effects in a stochastic Kuramoto model,” Chaos: An Interdisciplinary Journal of Nonlinear Science 27 (2017).
  * Peter and Pikovsky (2018) F. Peter and A. Pikovsky, “Transition to collective oscillations in finite Kuramoto ensembles,” Phys. Rev. E 97, 032310 (2018).
  * Yue and Gottwald (2023) W. Yue and G. A. Gottwald, “A stochastic approximation for the finite-size Kuramoto-Sakaguchi model,” (2023), arXiv:2310.20048 [nlin.AO] .
  * Hakim and Rappel (1992) V. Hakim and W.-J. Rappel, “Dynamics of the globally coupled complex Ginzburg-Landau equation,” Physical Review A 46, R7347 (1992).
  * Nakagawa and Kuramoto (1994) N. Nakagawa and Y. Kuramoto, “From collective oscillations to collective chaos in a globally coupled oscillator system,” Physica D: Nonlinear Phenomena 75, 74–80 (1994).
  * Chabanol, Hakim, and Rappel (1997) M.-L. Chabanol, V. Hakim, and W.-J. Rappel, “Collective chaos and noise in the globally coupled complex Ginzburg-Landau equation,” Physica D: Nonlinear Phenomena 103, 273–293 (1997).
  * Ku, Girvan, and Ott (2015) W. L. Ku, M. Girvan, and E. Ott, “Dynamical transitions in large systems of mean field-coupled Landau-Stuart oscillators: Extensive chaos and cluster states,” Chaos: An Interdisciplinary Journal of Nonlinear Science 25 (2015).
  * Clusella and Politi (2019) P. Clusella and A. Politi, “Between phase and amplitude oscillators,” Physical Review E 99, 062201 (2019).
  * León and Pazó (2022) I. León and D. Pazó, “Enlarged Kuramoto model: Secondary instability and transition to collective chaos,” Physical Review E 105, L042201 (2022).
  * Olmi, Politi, and Torcini (2011) S. Olmi, A. Politi, and A. Torcini, “Collective chaos in pulse-coupled neural networks,” Europhysics Letters 92, 60007 (2011).
  * Pazó and Montbrió (2016) D. Pazó and E. Montbrió, “From quasiperiodic partial synchronization to collective chaos in populations of inhibitory neurons with delay,” Physical Review Letters 116, 238101 (2016).
  * Ratas and Pyragas (2018) I. Ratas and K. Pyragas, “Macroscopic oscillations of a quadratic integrate-and-fire neuron network with global distributed-delay coupling,” Physical Review E 98, 052224 (2018).
  * Klinshov  _et al._ (2021) V. V. Klinshov, S. Y. Kirillov, V. I. Nekorkin, and M. Wolfrum, “Noise-induced dynamical regimes in a system of globally coupled excitable units,” Chaos: An Interdisciplinary Journal of Nonlinear Science 31 (2021).
  * Note (3) Notice that the normalization of the noise strength differs from that used in Eq. (15).
  * Rosenblum and Pikovsky (2015) M. Rosenblum and A. Pikovsky, “Two types of quasiperiodic partial synchrony in oscillator ensembles,” Phys. Rev. E 92, 012919 (2015).
  * Carlu, Ginelli, and Politi (2018) M. Carlu, F. Ginelli, and A. Politi, “Origin and scaling of chaos in weakly coupled phase oscillators,” Physical Review E 97, 012203 (2018).
  * Popovych, Maistrenko, and Tass (2005) O. V. Popovych, Y. L. Maistrenko, and P. A. Tass, “Phase chaos in coupled oscillators,” Physical Review E 71, 065201 (2005).
  * Maistrenko, Popovych, and Tass (2005) Y. L. Maistrenko, O. V. Popovych, and P. A. Tass, “Chaotic attractor in the Kuramoto model,” International Journal of Bifurcation and Chaos 15, 3457–3466 (2005).