# A self-consistent analytical theory for rotator networks under stochastic
forcing: effects of intrinsic noise and common input

Jonas Ranft<EMAIL_ADDRESS>Institut de Biologie de l’ENS, Ecole
Normale Supérieure, CNRS, Inserm, Université PSL, 46 rue d’Ulm, 75005 Paris,
France  Benjamin Lindner  Bernstein Center for Computational Neuroscience
Berlin, Philippstraße 13, Haus 2, 10115 Berlin, Germany and Physics Department
of Humboldt University Berlin, Newtonstraß e 15, 12489 Berlin, Germany

(August 27, 2024)

###### Abstract

Despite the incredible complexity of our brains’ neural networks, theoretical
descriptions of neural dynamics have led to profound insights into possible
network states and dynamics. It remains challenging to develop theories that
apply to spiking networks and thus allow one to characterize the dynamic
properties of biologically more realistic networks. Here, we build on recent
work by van Meegen & Lindner who have shown that “rotator networks,” while
considerably simpler than real spiking networks and therefore more amenable to
mathematical analysis, still allow to capture dynamical properties of networks
of spiking neurons. This framework can be easily extended to the case where
individual units receive uncorrelated stochastic input which can be
interpreted as intrinsic noise. However, the assumptions of the theory do not
apply anymore when the input received by the single rotators is strongly
correlated among units. As we show, in this case the network fluctuations
become significantly non-Gaussian, which calls for a reworking of the theory.
Using a cumulant expansion, we develop a self-consistent analytical theory
that accounts for the observed non-Gaussian statistics. Our theory provides a
starting point for further studies of more general network setups and
information transmission properties of these networks.

> Not surprisingly, networks of recurrently connected neurons , such as e.g.
> found in the human cortex, are capable to support a very rich repertoire of
> dynamical states. A remarkable feature of such networks is that they can
> exhibit asynchronous, irregular activity even in the case of deterministic
> dynamics, where the disorder due to the random connectivity between neurons
> gives rise to temporal fluctuations of their activity. Understanding the
> statistics and dynamical properties of network states by mathematical
> analysis of models of recurrent networks remains a challenging problem
> especially for networks of spiking neurons, for which analytical expressions
> e.g. for the network noise autocorrelation function are hard to obtain. For
> these reasons, we investigate here a much simpler model class of recurrently
> connected oscillators that nevertheless has been shown to give a reasonable
> account of the statistics of spiking networks, and that allows one to obtain
> in comparison much simpler self-consistent equations for the network
> statistics. Using this model, we study the effects of intrinsic noise
> uncorrelated among oscillators and shared common noise. Interestingly, we
> find that in the presence of common noise the network fluctuations become
> non-Gaussian even if both the self-generated network noise in the absence of
> noise and the additional common noise are Gaussian. To account for the non-
> Gaussian statistics, we develop a self-consistent theory that includes
> higher cumulants of the network fluctuations. Our theory provides a
> satisfactory description of the observed power spectra and autocorrelation
> functions of the network fluctuations and rotator dynamics, and can serve as
> a starting point for further investigations of networks in which non-
> Gaussian statistics due to input correlations might play a role.

##  I Introduction

Networks of recurrently coupled oscillators have been studied in a wide
variety of contexts, representing spatially inhomogeneous chemical reactions
Kuramoto (1984), neural networks Ermentrout (1998); Grytskyy  _et al._ (2013);
di Volo and Torcini (2018), or power grids Filatrella  _et al._ (2008). A
classic paradigm has been and continues to be the Kuramoto model, which
allowed many analytical insights in e.g. the synchronization dynamics
depending on the disorder in the system Kuramoto (1984); Pikovsky  _et al._
(2001); Acebrón  _et al._ (2005). More generally, networks of recurrently
connected units can be in a regime where they exhibit apparently noisy
dynamics despite being entirely deterministic, a feature which has received a
lot of attention in the context of neural networks Sompolinsky  _et al._
(1988); Brunel (2000); Lerchner  _et al._ (2006); Kriener  _et al._ (2008);
Dummer  _et al._ (2014); Ostojic (2014); van Meegen and Lindner (2018); Pena
_et al._ (2018); Vellmer and Lindner (2019). While such noisy dynamics can
easily be observed e.g. in spiking network models where it manifests itself in
the “asynchronous state” Renart  _et al._ (2010), it is considerably harder to
analytically determine the effective noise characteristics for such models. In
order to thoroughly characterize the internally-generated noise as a function
of network parameters and understand its impact e.g. on the statistics of
individual units in the network, it seems important to have models at one’s
disposal that capture the phenomenon while still allowing one to gain
analytical insights.

Recently, van Meegen & Lindner van Meegen and Lindner (2018) introduced the
so-called “rotator network” of recurrently coupled rotators, in which the
coupling between the oscillating units is an arbitrary function of the
“presynaptic” oscillator phases. Despite being a strongly simplified
description of real spiking networks, the authors showed that the model still
allows to capture dynamical properties of model networks of spiking neurons.
Importantly, for the rotator network one can self-consistently calculate the
statistics of the internally-generated fluctuations (network noise). Here, we
will build upon this work to investigate the influence of additional intrinsic
noise sources at the level of the individual units as well as shared common
input. The additional input leads to modifications of the self-consistent
equation for the autocorrelation function of the network noise and impacts the
individual units’ statistics. When the external inputs are uncorrelated across
the units, the required modifications of the theory are straightforward.
However, if there are substantial correlations in the input across units,
basic assumptions of the existing theory break down and require a more
substantial revision with regard to the non-Gaussian statistics of the network
noise. Extending the theory accordingly, we aim to go a step further towards
increased biological realism in the context of neural networks while
maintaining the analytical tractability of the model.

After briefly introducing the basic model and reviewing previous results in
section II, we extend the rotator network model to the case where individual
units are subject to additional intrinsic noise in section III. Intrinsic
noise can e.g. be caused by stochastic ion channel openings in the case of
neural networks, or by fluctuating inputs or loads in power grids. In section
IV, we discuss the effect of shared external input on all units. Whereas the
internally-generated network noise remains Gaussian in the presence of
intrinsic noise only, it becomes non-Gaussian in the presence of correlated
inputs. This is captured by a modified self-consistent equation that takes
into account higher cumulants of the internally-generated noise. We conclude
in section V with a summary of our results and a brief discussion of their
implications.

##  II Previous results: network dynamics without external inputs

In a previous study, van Meegen & Lindner introduced a simple model of a
recurrent rotator network, where individual units are described as phase
variables $\theta_{m}$, each driven by an intrinsic oscillation frequency
$\omega_{m}$ and receiving input from the other units in the network,

$\dot{\theta}_{m}=\omega_{m}+\sum_{n}K_{mn}f(\theta_{n}),$ (1)

where the recurrent input is specified by the matrix of connection weights
$K_{mn}$ between units $n$ and $m$, and the coupling function $f(\cdot)$ that
depends only on the ‘input phases’ $\theta_{n}$. To be specific, we will
assume $K_{mn}$ to be a Gaussian random matrix with $\langle K_{mn}\rangle=0$
and $\langle
K_{mn}K_{m^{\prime}n^{\prime}}\rangle=\frac{K^{2}}{N}\delta_{mm^{\prime}}\delta_{nn^{\prime}}$
unless specified otherwise, but see section V for a discussion of other
choices for the connectivity. Based on dynamical mean-field theory, Eq. (1)
can be interpreted as a stochastic dynamics for the $\theta_{m}$ with an
internally generated network noise $\xi_{m}$,

$\xi_{m}=\sum_{n}K_{mn}f(\theta_{n}),$ (2)

the statistics of which has to be determined self-consistently. The relative
simplicity of the model allows one to derive an analytical solution for the
noise autocorrelation function
$C_{\xi}(\tau)\equiv\langle\xi_{m}(t)\xi_{m}(t+\tau)\rangle$. To this end,
following van Meegen & Lindner, we introduce the auxiliary function

$\Lambda(\tau)=\int_{0}^{\tau}dt(\tau-t)C_{\xi}(t),$ (3)

which implies

$C_{\xi}(\tau)=\ddot{\Lambda}(\tau).$ (4)

In the case of the autonomous network described by Eq. (1), the function
$\Lambda(\tau)$ is a solution of the ordinary differential equation

$\ddot{\Lambda}(\tau)=K^{2}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\exp[-l^{2}\Lambda(\tau)].$
(5)

Here, $A_{l}$ are the coefficients of a Fourier expansion of the coupling
function, $f(\theta)=\sum_{l}A_{l}\exp[il\theta]$, and
$\Phi(x)=\langle\exp[i\omega_{n}x]\rangle_{\omega_{n}}$ is the characteristic
function of the intrinsic frequencies. The above equation for $\Lambda$
follows from a self-consistency condition for the noise autocorrelation, where
the average is taken over realizations of the connectivity matrix $K_{mn}$,
intrinsic frequencies $\omega_{n}$, initial conditions, and the effective
noise $\xi_{m}$ in a self-consistent manner.

We would like to emphasize that, while the relations of $\Lambda(\tau)$ in
Eqs. (3) and (4) hold true throughout this manuscript, the dynamical equation
that one obtains for $\Lambda$ will depend on the situation considered
(absence or presence of intrinsic and/or common noise) and, for the most
involved case of common noise, on the level of approximation used. It would be
possible to add a label to the function $\Lambda$ that would distinguish these
different cases, but we abstain from doing so for the ease of notation.

We quickly recapitulate the derivation of van Meegen & Lindner for reference,
also because our new results will be obtained along similar lines. Using the
Gaussian statistics when averaging over the connectivity matrix, one first
obtains

$\displaystyle C_{\xi}(\tau)$ $\displaystyle=\sum_{n,n^{\prime}}\langle
K_{mn}K_{mn^{\prime}}f(\theta_{n}(t))f(\theta_{n^{\prime}}(t+\tau))\rangle_{K,\bm{\omega},\bm{\theta_{0}}}$
(6a) $\displaystyle=\sum_{n,n^{\prime}}\langle
K_{mn}K_{mn^{\prime}}\rangle_{K}\langle
f(\theta_{n}(t))f(\theta_{n^{\prime}}(t+\tau))\rangle_{K,\bm{\omega},\bm{\theta_{0}}}$
(6b) $\displaystyle=\frac{K^{2}}{N}\sum_{n}\langle
f(\theta_{n}(t))f(\theta_{n}(t+\tau))\rangle_{K,\bm{\omega},\bm{\theta_{0}}}.$
(6c)

Note that the disorder of the connectivity matrix in principle still
contributes to the fluctuating dynamics of the $\theta_{n}(t)$. The self-
consistency ansatz for statistics of the internally-generated network
fluctuations now consists in considering the effective stochastic dynamics
$\dot{\theta}_{m}=\omega_{m}+\xi_{m}$ for the rotators in the network and to
average over realizations of the network noise $\xi_{m}(t)$ in the following.
Using the Fourier expansion of $f(\cdot)$ and
$\theta_{m}(t)=\theta_{m}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}\dot{\theta}_{m}(t^{\prime})$,
one thus obtains

$\displaystyle\frac{C_{\xi}(\tau)}{K^{2}/N}$
$\displaystyle=\sum_{n}\sum_{l,l^{\prime}}A_{l}A_{l^{\prime}}\langle
e^{il\theta_{n}(t)+il^{\prime}\theta_{n}(t+\tau)}\rangle_{\bm{\xi},\bm{\omega},\bm{\theta_{0}}}$
(7a) $\displaystyle=\sum_{n}\sum_{l,l^{\prime}}A_{l}A_{l^{\prime}}\langle
e^{i(l+l^{\prime})\theta_{n}(t_{0})}\rangle_{\bm{\theta_{0}}}\times$
$\displaystyle\qquad\qquad\langle
e^{il\int_{t_{0}}^{t}dt^{\prime}\dot{\theta}_{n}(t^{\prime})+il^{\prime}\int_{t_{0}}^{t+\tau}dt^{\prime}\theta_{n}(t^{\prime})}\rangle_{\bm{\xi},\bm{\omega}}$
(7b) $\displaystyle=\sum_{n}\sum_{l}|A_{l}|^{2}\langle
e^{il\omega_{n}\tau}\rangle_{\bm{\omega}}\langle
e^{il\int_{t}^{t+\tau}dt^{\prime}\xi_{n}(t^{\prime})}\rangle_{\bm{\xi}}$ (7c)
$\displaystyle=\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\langle
e^{ily_{n}(\tau;t)}\rangle$ (7d)
$\displaystyle=\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l\tau)e^{-\frac{l^{2}}{2}\langle
y_{n}(\tau;t)^{2}\rangle}$ (7e)
$\displaystyle=N\sum_{l}|A_{l}|^{2}\Phi(l\tau)e^{-l^{2}\int_{0}^{\tau}dt^{\prime}(\tau-t^{\prime})C_{\xi}(t^{\prime})}.$
(7f)

A few remarks are in order: To ensure that the result does not depend on a
particular initialization, one can average over the absolute phases
$\theta_{n}(t_{0})$ at a given (arbitrary) reference time $t_{0}$ and require
that this average can be taken independently of the average over the
$\xi_{n}$. For uniformly distributed phases $\theta_{n}(t_{0})$, one
immediately obtains $\langle
e^{i(l+l^{\prime})\theta_{n}(t_{0})}\rangle_{\theta_{n}(t_{0})}=\delta_{l,l^{\prime}}$.
We furthermore introduced the new variable
$y_{n}(\tau;t)=\int_{t}^{t+\tau}dt^{\prime}\xi_{n}(t^{\prime})$. In the limit
of a large network ($N\gg 1$), one can assume that the $\xi_{n}$ are Gaussian
distributed, and thus that $y_{n}(\tau;t)$ obeys the relation $\langle
e^{ily_{n}}\rangle=e^{-l^{2}\langle y_{n}^{2}\rangle/2}$. The second moment
$\langle y_{n}^{2}\rangle$ can be found according to

$\displaystyle\langle y_{n}(\tau;t)^{2}\rangle$
$\displaystyle=\int_{t}^{t+\tau}\\!\\!\\!dt_{1}\int_{t}^{t+\tau}\\!\\!\\!dt_{2}\left\langle\xi_{n}(t_{1})\xi_{n}(t_{2})\right\rangle$
$\displaystyle=2\int_{0}^{\tau}\\!\\!\\!dt(\tau-t)C_{\xi}(t)$ (8)

after a change of variables and using the symmetry of $C_{\xi}(t)$ around the
origin. Since all units are statistically equivalent and their correlation
functions the same, the sum over $n$ then becomes trivial.

The differential equation for $\Lambda$ (Eq. (5)) is eventually obtained with
the substitution Eq. (3), and the noise autocorrelation follows as

$C_{\xi}(\tau)=K^{2}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\exp[-l^{2}\Lambda(\tau)].$
(9)

The autocorrelation function $C_{x_{m}}(t)$ of the pointer
$x_{m}(t)=e^{i\theta_{m}(t)}$ for a given rotator $m$ is then rather
straightforwardly obtained as

$\displaystyle\langle x_{m}(t)^{*}x_{m}(t+\tau)\rangle$ $\displaystyle=\langle
e^{i(\theta_{m}(t+\tau)-\theta_{m}(t))}\rangle$
$\displaystyle=e^{i\omega_{m}\tau}\langle
e^{i\int_{t}^{t+\tau}dt^{\prime}\xi_{m}(t^{\prime})}\rangle$
$\displaystyle=e^{i\omega_{m}\tau-\Lambda(\tau)}.$

Of note, $C_{x_{m}}(t)$ explicitly depends on the particular intrinsic
frequency $\omega_{m}$ of the rotator $m$. Averaged over all rotators, the
autocorrelation becomes $C_{x}(t)=\Phi(\tau)e^{-\Lambda(\tau)}$.

##  III Rotator network with individual noise sources

Focussing on recurrent input, van Meegen & Lindner did not address the effect
of external inputs or other noise sources acting on individual units.
Typically, however, in many applications, e.g. in biological neural networks,
the single units are subject to intrinsic noise, e.g. channel noise or
unreliable synaptic transmission in the case of neural networks Tuckwell
(1989); Koch (2004); Rusakov  _et al._ (2020).

It seems important to understand the degree to which these additional inputs
may shape the dynamics of the network and individual units. There is a rich
literature on the Kuramoto model driven by individual and common noise
beginning with Kuramoto himself, see e.g. Kuramoto (1984); Pikovsky and
Rosenblum (2015); Schäfer  _et al._ (2018), but how noise affects the rotator
dynamics considered here is not known. In this section, we will consider the
case of intrinsic noise sources $\eta_{m}$ that are uncorrelated between
different units (correlated inputs will be addressed in section IV). The model
then reads

$\dot{\theta}_{m}=\omega_{m}+\sum_{n}K_{mn}f(\theta_{n})+\eta_{m},$ (10)

the recurrent input again being specified by the connectivity matrix $K_{mn}$
and coupling function $f(\cdot)$. We consider the intrinsic fluctuations to be
independent Gaussian noise sources with $\langle\eta_{m}\rangle=0$,
$\langle\eta_{m}(t)\eta_{n}(t+\tau)\rangle=\delta_{mn}C_{\eta}(\tau)$.

###  III.1 Derivation of the self-consistent correlation functions in the
presence of intrinsic noise

We aim to determine how the presence of additional noise may shape the rotator
autocorrelation function $C_{x_{m}}(\tau)$ of rotator $m$, the effective
stochastic dynamics of which is given by

$\dot{\theta}_{m}=\omega_{m}+\xi_{m}+\eta_{m},$ (11)

where the statistics of the network noise $\xi_{m}$ is still to be determined.
First of all, it is not difficult to show that the network noise is not
correlated with the external noise, as seen as follows:

$\displaystyle\langle\xi_{m}(t)\eta_{m}(t^{\prime})\rangle$
$\displaystyle=\sum_{n}\langle K_{mn}\rangle_{K}\langle
f(\theta_{n}(t))\eta_{m}(t^{\prime})\rangle=0$ (12)

since $\langle K_{mn}\rangle_{K}=0$.

We now turn to the problem of calculating the autocorrelation function of the
network noise. When all rotators are subject to intrinsic noise, we find a
modified self-consistency condition for $C_{\xi}(t)$ as follows. Along the
lines of the derivation in the case without noise, Eq. (7c) now becomes

$\displaystyle\frac{C_{\xi}(\tau)}{K^{2}/N}$
$\displaystyle=\sum_{n}\sum_{l}|A_{l}|^{2}\langle
e^{il\omega_{n}\tau}\rangle_{\bm{\omega}}\langle
e^{il\int_{t}^{t+\tau}dt^{\prime}[\xi_{n}(t^{\prime})+\eta_{n}(t^{\prime})]}\rangle_{\bm{\xi},\bm{\eta}}$
$\displaystyle=\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\langle
e^{ily_{n}(\tau;t)}\rangle,$ (13)

where we redefined the variable

$y_{n}(\tau;t)=\int_{t}^{t+\tau}dt^{\prime}[\xi_{n}(t^{\prime})+\eta_{n}(t^{\prime})].$
(14)

Assuming that $y_{n}(\tau;t)$ is Gaussian, we now find

$\displaystyle\langle e^{ily_{n}(\tau)}\rangle$
$\displaystyle=e^{-\frac{l^{2}}{2}\langle y_{n}(\tau)^{2}\rangle}$ (15a)
$\displaystyle=e^{-l^{2}\int_{0}^{\tau}dt(\tau-t)\langle[\xi_{n}(t^{\prime})+\eta_{n}(t^{\prime})][\xi_{n}(t^{\prime}+t)+\eta_{n}(t^{\prime}+t)]\rangle}$
(15b)
$\displaystyle=e^{-l^{2}\int_{0}^{\tau}dt(\tau-t)[C_{\xi}(t)+C_{\eta}(t)]}.$
(15c)

Putting everything together, one obtains

$\displaystyle C_{\xi}(\tau)$
$\displaystyle=K^{2}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\exp\\!{\left(-l^{2}[\Lambda(\tau)+\Lambda_{\eta}(\tau)]\right)}.$
(16)

Here $\Lambda_{\eta}(\tau)=\int_{0}^{\tau}dt(\tau-t)C_{\eta}(t)$ is given by
the autocorrelation function $C_{\eta}(t)$ of the intrinsic noise
sources111Note that the integral can be either analytically calculated for
simple processes or determined by numerically integrating the ordinary
differential equation for $\Lambda_{\eta}(\tau)$ with
$\ddot{\Lambda}_{\eta}(t)=C_{\eta}(t)$, $\Lambda_{\eta}(0)=0$,
$\dot{\Lambda}_{\eta}(0)=0$. and the function
$\Lambda(\tau)=\int_{0}^{\tau}dt(\tau-t)C_{\xi}(t)$, defined as before, now
obeys the ordinary differential equation

$\ddot{\Lambda}(\tau)=K^{2}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\exp\\!{\left(-l^{2}[\Lambda(\tau)+\Lambda_{\eta}(\tau)]\right)}.$
(17)

Once we have calculated $\Lambda(t)$, we know the autocorrelation function of
the network noise $C_{\xi}(t)$ via Eq. (16), but can now also determine the
autocorrelation function $C_{x_{m}}(t)$ of the $m$-th oscillator (the complex
pointer $x_{m}=e^{i\theta_{m}}$):

$C_{x_{m}}(\tau)=\exp\\!\left[i\omega_{m}\tau-\Lambda(\tau)-\Lambda_{\eta}(\tau)\right].$
(18)

Similarly to the case without intrinsic noise, this can be derived using
$\theta_{m}(t)=\theta_{m}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}\dot{\theta}_{m}(t^{\prime})$
with Eq. (11) and using Eq. (15) for $l=1$. Eqs. (16-18) constitute the theory
for the self-consistent correlation functions for a rotator network with
individual intrinsic Gaussian noise.

So far we have not made any assumptions about the temporal correlations of the
intrinsic noise, i.e., our theory holds for a general colored noise
$\eta_{m}(t)$. For our numerical examples in the following, we consider for
simplicity Gaussian white noise, for which
$C_{\eta}(\tau)=2D_{\eta}\delta(\tau)$, implying
$\Lambda_{\eta}(\tau)=D_{\eta}\tau$.

###  III.2 Effects of intrinsic white noise on the rotator autocorrelation

Figure 1: Rotator autocorrelation function $C_{x_{m}}$ in the presence of
intrinsic noise, for weak (a) and strong (b) recurrent connection strengths.
Theoretical results (solid lines) and stochastic network simulations (crosses)
are shown for coupling function $f(\theta)=\sin(2\theta)+\cos(3\theta)$ and
parameters $\omega_{0}=1$, $\sigma_{\omega}=0$, and (a) $K=0.5$, respectively
(b) $K=2$. The size of the simulated networks was $N=100$. Simulations were
carried out with a time step $dt=0.01$, stationary initial conditions (i.e.
uniformly distributed phases), for a total duration of $T=25T_{0}$ with
$T_{0}=2500$, and $10$ different realizations of $K_{mn}$. Statistics were
computed for bouts of length $T_{0}$ of the phase trajectories, and averages
subsequently taken over bouts, network realizations, and rotators. The theory
was evaluated with a time step $dt=0.001$.

In order to test our theory, we simulated the rotator network with additional
intrinsic white noise for various combinations of noise intensity $D_{\eta}$
and coupling strength $K$. The corresponding theoretical predictions for
$C_{x_{m}}(t)$ obtained from Eqs. (17) and (18) are nicely confirmed for all
combinations used, see Fig. 1, although the number of units in the network
($N=100$) is not excessively large. For a white noise, Eq. (18) reads

$C_{x_{m}}(\tau)=\exp\\!\left[i\omega_{m}\tau-\Lambda(\tau)-D_{\eta}\tau\right],$
(19)

from which we can expect a damping of the autocorrelation function with
increasing noise strength $D_{\eta}$. This is indeed observed in the case of
both weak (Fig. 1a) and strong (Fig. 1b) network noise.

The network-mediated influence of the intrinsic noises $\eta_{n}$, $n\neq m$,
on the autocorrelation function of rotator $m$ enters the expression for
$C_{x_{m}}(t)$ via the modified network-noise autocorrelation function
$C_{\xi}(t)$. In order to assess this indirect effect of the intrinsic noise,
we compare our result to the case where only a specific, single unit $m$
receives additional intrinsic noise $\eta_{m}$ according to Eq. (10), while
all other units of the network evolve according to Eq. (1). For the particular
rotator $m$ that is subject to intrinsic noise, $C_{x_{m}}(t)$ is then still
given by Eq. (18) with the exponential damping term $e^{-D_{\eta}t}$, but
where the network-noise autocorrelation $C_{\xi}(t)$ is the “intrinsic-noise
free” version $C_{\xi}(t)=\ddot{\Lambda}(t)$ with $\Lambda$ obeying Eq. (5).

Figure 2: The network-mediated effect of noise in “background” units on
rotator autocorrelation function $C_{x_{m}}$ of noisy unit $m$ (solid lines)
can be assessed by comparison with the autocorrelation function $C_{x_{m}}$ of
unit $m$ in the case where only unit $m$ is subject to noise (dashed lines).
The autocorrelation in the absence of any noise is also shown (thin lines).
(a) Intermediate noise level with $K=0.5$, $D_{\eta}=0.2$; (b) strong noise
level with $K=2$, $D_{\eta}=0.5$. All curves are obtained from the theory,
with $dt=0.001$.

Stochastic network simulations and theoretical predictions for $C_{x}(t)$
again nicely agree (not shown) but look very similar to the case where all
units are subject to intrinsic noise shown in Fig. 1. We nevertheless identify
two distinct effects due to the network-mediated influence of the other
intrinsic noise sources, as shown in Fig. 2: For weak network noise (Fig. 2a),
we find that the intrinsic noise on the other units in the network further
decorrelates oscillator $m$ relative to the case of a single noisy unit with
“intrisic-noise free” recurrent input. Interestingly, in the case of strong
network noise (Fig. 2b), the opposite effect can be observed, as the
autocorrelation of the specific rotator $m$ decays faster in the case where it
is the only one subject to intrinsic noise while the rest of the network does
evolve without additional intrinsic noise.

##  IV Rotator network with common external input

So far we extended the rotator network model of Ref. van Meegen and Lindner,
2018 to the case where individual units are subject to intrinsic noise, and
assumed independence of the noise sources for different units. In the context
of neural networks, possible origins of such a noise term could be the
neurons’ stochastic ion channel kinetics and spontaneous release of
neurotransmitters Tuckwell (1989), for which the assumption of independence
seems clearly satisfied. Another possible interpretation of the noise term
$\eta_{m}$ is that it results from external (as opposed to recurrent) synaptic
inputs that are uncorrelated across units.

More generally however, neurons may also receive _correlated_ external inputs,
which could stem e.g. from other brain areas that project broadly to a local
recurrently connected network. Here, we will thus relax the assumption that
any additional inputs to the individual units are independent and study the
potential effects on the network dynamics. Regarding the theory, one might be
tempted to think that it can be extended to this new case without major
changes and that its main assumptions, e.g. the Gaussian distribution of the
network noise, are still valid. We will show however that the network noise
ceases to be Gaussian-distributed if the external noise is significantly
correlated across units. Consequently, the self-consistent theory of Ref. van
Meegen and Lindner, 2018 cannot easily be extended to this case. Instead, we
have to more substantially rework the theory. By taking higher-order
correlation functions of the network noise into account, we can arrive at an
approximate theory that captures the resulting non-Gaussian network statistics
surprisingly well.

To be specific, we consider the following model, where $\eta_{\mathrm{c}}$ is
the additional input common to all units:

$\dot{\theta}_{m}=\omega_{m}+\sum_{n}K_{mn}f(\theta_{n})+\eta_{m}+\eta_{\mathrm{c}}.$
(20)

We assume that $\eta_{\mathrm{c}}$ is a Gaussian noise which is uncorrelated
with the intrinsic noise $\eta_{m}$ of each unit $m$,
$\langle\eta_{\mathrm{c}}(t)\eta_{m}(t+\tau)\rangle\equiv 0$, and that its
temporal statistics are described by the autocorrelation function
$\langle\eta_{\mathrm{c}}(t)\eta_{\mathrm{c}}(t+\tau)\rangle=C_{\mathrm{c}}(\tau)$.
Furthermore, along the same lines as sketched for the intrinsic noise
$\eta_{m}$, Eq. (12), we note that $\eta_{\mathrm{c}}$ is uncorrelated with
the network noise $\xi_{m}$,
$\langle\eta_{\mathrm{c}}(t)\xi_{m}(t+\tau)\rangle=0$, by virtue of the
averaging properties of $K_{mn}$.

To develop the theory for the case with common input, we start again by
expressing the autocorrelation $\langle\xi_{m}(t)\xi_{m}(t+\tau)\rangle$ in
terms of the connectivity matrix $K_{mn}$, coupling function $f(\cdot)$, and
rotator phases $\theta_{m}$. The additional average over independent
realizations of the common noise $\eta_{\mathrm{c}}$ is the only departure
from the case considered in the previous section and does not interfere with
the initial steps of Eqs. (6) and (7), which now lead to

$C_{\xi}(\tau)=\frac{K^{2}}{N}\sum_{N}\sum_{l}|A_{l}|^{2}\langle
e^{il\omega_{n}\tau}\rangle_{\bm{\omega}}\times\\\ \langle
e^{il\int_{t}^{t+\tau}dt^{\prime}[\xi_{n}(t^{\prime})+\eta_{n}(t^{\prime})+\eta_{\mathrm{c}}(t^{\prime})]}\rangle_{\bm{\xi},\bm{\eta},\eta_{\mathrm{c}}}.$
(21)

By redefining the variable $y_{n}(\tau;t)$ according to

$y_{n}(\tau;t)=\int_{t}^{t+\tau}dt^{\prime}[\xi_{n}(t^{\prime})+\eta_{n}(t^{\prime})+\eta_{\mathrm{c}}(t^{\prime})]$
(22)

and noting that $\langle
e^{il\omega_{n}\tau}\rangle_{\bm{\omega}}=\Phi(l\tau)$, expression (21) can
again be written more conveniently as

$C_{\xi}(\tau)=\frac{K^{2}}{N}\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l\tau)\langle
e^{ily_{n}(\tau;t)}\rangle_{\bm{\xi},\bm{\eta},\eta_{\mathrm{c}}}.$ (23)

A naïve generalization of our previous approach would yield a simple
modification of Eqs. (16-18). Under the assumption that the $y_{n}(\tau;t)$
were still Gaussian-distributed, it would follow that the respective
expressions hold with the autocorrelation function $C_{\eta}(\tau)$ of the
intrinsic noise simply being replaced by the sum of the noise autocorrelation
functions $C_{\eta}(\tau)+C_{\mathrm{c}}(\tau)$. Consider the special
situation in which intrinsic noise and common noise share the same
autocorrelation function and we can describe by a parameter $c$ how much of
the input noise is common, i.e. $C_{\eta}(\tau)=(1-c)C_{\rm tot}(\tau)$ and
$C_{\mathrm{c}}(\tau)=cC_{\rm tot}(\tau)$. In this case, the output statistics
would not depend on $c$, implying that the network statistics is unaffected by
the presence or absence of correlations among the additional stochastic
forcings of the rotators. However, this lack of dependence is _not_ observed
in stochastic network simulations of Eq. (20), demonstrating that the naïve
generalization sketched above indeed fails to account for the effective
stochastic dynamics of the system.

In order to derive a correct self-consistent expression for the network-noise
autocorrelation function $C_{\xi}(\tau)$, we thus need to determine the only
remaining unknown term $\langle
e^{ily_{n}(\tau;t)}\rangle_{\bm{\xi},\bm{\eta},\eta_{\mathrm{c}}}$ in Eq.
(23). Let us consider the general case of a stochastic variable $y$. By
definition of the cumulant-generating function, we have

$\langle
e^{ily}\rangle_{y}=\exp\left(\sum_{k=1}^{\infty}\frac{(il)^{k}}{k!}\kappa_{k}\right).$
(24)

Here, the $\kappa_{k}$ are the cumulants of $y$ which we can calculate as a
function of its moments $\langle y^{j}\rangle_{y}$, $j\leq k$. If $y$ is zero-
centered and Gaussian-distributed, one immediately finds $\langle
e^{ily}\rangle_{y}=e^{-l^{2}\kappa_{2}/2}=e^{-l^{2}\langle
y^{2}\rangle_{y}/2}$, which led to the results obtained in the previous
sections. We will see that, while $y_{n}(\tau;t)$ as defined in Eq. (22) still
has a vanishing mean value in the presence of common input, it is no longer
Gaussian-distributed and higher cumulants need to be taken into account.

Assuming that the deviations from Gaussian statistics are small, our strategy
is the following: First, we approximate the average
$\langle{e^{ily}}\rangle_{y}$ by taking only the first few cumulants according
to Eq. (24) into account, i.e.,

$\langle
e^{ily}\rangle_{y}\approx\exp\\!\left(-l^{2}\frac{\kappa_{2}}{2}-il^{3}\frac{\kappa_{3}}{6}+l^{4}\frac{\kappa_{4}}{24}\right).$
(25)

Here, the higher-order cumulants $\kappa_{3}$ and $\kappa_{4}$ are finite only
in the presence of correlations between the inputs of different units, and
should vanish with vanishing correlations. As we will show, these higher-order
cumulants have in turn themselves Gaussian and non-Gaussian contributions in a
sense that will be made explicit below. In line with our assumption of a small
departure from Gaussian statistics, we thus, in a second approximation,
calculate each of these cumulants to zeroth order, that is, assuming purely
Gaussian statistics of composite variables of the type
$y_{n}+\eta_{\mathrm{c}}$ and $y_{n}+y_{m}$. This strategy allows us to obtain
a self-consistent equation for $C_{\xi}(\tau)$ in the presence of common input
that reproduces measured autocorrelation functions as well as higher-order
cumulants of $y_{n}$ in stochastic network simulations.

###  IV.1 Calculation of the cumulants of $y_{n}(\tau;t)$

The first cumulant is identical to the first moment and vanishes,

$\langle
y_{n}(\tau;t)\rangle=\int\limits_{t}^{t+\tau}dt^{\prime}[\langle{\xi_{n}(t^{\prime})}\rangle+\langle{\eta_{n}(t^{\prime})}\rangle+\langle{\eta_{\mathrm{c}}(t^{\prime})}\rangle]=0,$
(26)

as $\langle{\xi_{n}}\rangle=0$ because of the average over the
$K_{nn^{\prime}}$ and
$\langle{\eta_{n}}\rangle=\langle{\eta_{\mathrm{c}}}\rangle=0$ by definition.
If the first moment vanishes, the second cumulant $\kappa_{2}$, i.e. the
variance, is identical to the second moment. Its calculation is analogous to
the case of purely intrinsic noise discussed above, and given by

$\langle
y_{n}(\tau)^{2}\rangle=\int\limits_{0}^{\tau}dt(\tau-t)[C_{\xi}(t)+C_{\eta}(t)+C_{\mathrm{c}}(t)].$
(27)

The new non-trivial contributions in this theory will be the third and fourth
cumulants which can also be expressed in terms of the moments; in the
following we will calculate and use

$\displaystyle\kappa_{3}$ $\displaystyle=\langle{y_{n}^{3}}\rangle,$
$\displaystyle\kappa_{4}$
$\displaystyle=\langle{y_{n}^{4}}\rangle-3\langle{y_{n}^{2}}\rangle^{2},$

to approximate the characteristic function in Eq. (23).

####  IV.1.1 Third cumulant

Inserting the definition of $y_{n}$, we obtain for the third cumulant

$\displaystyle\kappa_{3}(\tau)$ $\displaystyle=\langle
y_{n}^{3}(\tau;t)\rangle$
$\displaystyle=\left\langle\prod_{j=1}^{3}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{j}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\right\rangle.$
(28)

Note that $\kappa_{3}$ does neither depend on $t$ nor on any neuron index $n$,
as these dependences disappear after averaging over uniform initial conditions
for the phases and the disorder of the connectivity, respectively.

Before evaluating the triple integral, we consider separately the different
contributions that arise from the correlator of the different noise sources,
such as
$\langle{\xi_{n}(t_{1})\eta_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})}\rangle$ or
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})}\rangle$.
Because $\eta_{n}$ and $\eta_{\mathrm{c}}$ are Gaussian and uncorrelated by
construction, all terms in which only these two appear vanish. Furthermore,
averaging over the zero-centered, Gaussian $K_{mn}$ implies that all terms
with odd powers of $\xi_{n}$ must also vanish, see e.g. Eq. (12). As all three
integrals in Eq. (28) cover the same domain, we can arbitrarily permute the
variables $t_{j}$, and retain the following remaining contributions:

$\big{\langle}\prod_{j=1}^{3}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\big{\rangle}/3=\\\
\langle\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\rangle+\langle\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\rangle,$
(29)

where the equality holds under the triple integral of Eq. (28).

We next aim to express both contributions in terms of the correlation
functions $C_{\xi}$, $C_{\eta}$, and $C_{\mathrm{c}}$, where the latter two
are known and $C_{\xi}$ remains to be determined self-consistently. To be
specific and in line with our treatment of $\eta_{m}$ in the previous section,
we will consider both $\eta_{m}$ and $\eta_{\mathrm{c}}$ to be Gaussian white
noises with respective intensities $D_{\eta}$ and $D_{\mathrm{c}}$ for the
remainder of this section. The correlation functions are then given by
$C_{\eta}(\tau)=2D_{\eta}\delta(\tau)$,
$C_{\mathrm{c}}(\tau)=2D_{\mathrm{c}}\delta(\tau)$, which decidedly simplifies
calculations. Using the definition of $\xi_{m}$, we can expand the first
contribution of Eq. (29) as

$\displaystyle\mkern-36.0mu\langle\xi_{m}(t_{1})\xi_{m}(t_{2})\eta_{m}(t_{3})\rangle$
$\displaystyle=\sum_{n,n^{\prime}}\langle
K_{mn}K_{mn^{\prime}}f(\theta_{n}(t_{1}))f(\theta_{n^{\prime}}(t_{2}))\eta_{m}(t_{3})\rangle$
$\displaystyle=\sum_{n,n^{\prime}}\langle
K_{mn}K_{mn^{\prime}}\rangle_{K}\langle
f(\theta_{n}(t_{1}))f(\theta_{n^{\prime}}(t_{2}))\eta_{m}(t_{3})\rangle$
$\displaystyle=\frac{K^{2}}{N}\sum_{n}\langle
f(\theta_{n}(t_{1}))f(\theta_{n}(t_{2}))\eta_{m}(t_{3})\rangle.$ (30)

Again along the lines of Eqs. (6) and (7) but with
$\dot{\theta}_{m}=\omega_{m}+\xi_{m}+\eta_{m}+\eta_{\mathrm{c}}$, we now
obtain

$\langle\xi_{m}(t_{1})\xi_{m}(t_{2})\eta_{m}(t_{3})\rangle=\frac{K^{2}}{N}\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l(t_{2}-t_{1}))\times\\\
\langle
e^{ily_{n}(t_{2}-t_{1};t_{1})}\eta_{m}(t_{3})\rangle_{\bm{\xi},\bm{\eta},\eta_{\mathrm{c}}}.$
(31)

Instead of the characteristic function $\langle e^{ily_{n}}\rangle$, we now
need to find an expression for the average $\langle
e^{ily_{n}}\eta_{m}\rangle$.

We can distinguish the cases $n\neq m$ and $n=m$. When $n\neq m$, $\langle
y_{n}(t_{2}-t_{1})\eta_{m}(t_{3})\rangle=0$, since $\xi_{n}$, $\eta_{n}$, and
$\eta_{\mathrm{c}}$ are all uncorrelated with $\eta_{m}$; therefore, $\langle
e^{ily_{n}(t_{2}-t_{1};t_{1})}\eta_{m}(t_{3})\rangle=\langle{e^{ily_{n}(t_{2}-t_{1};t_{1})}}\rangle\langle{\eta_{m}(t_{3})}\rangle=0$.
When $n=m$, the correlation $\langle
e^{ily_{n}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\rangle$ does not vanish but
should remain bounded for large $N$. We thus expect that, according to Eq.
(31), $\langle\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\rangle$ scales like
$1/N$ and vanishes for $N\to\infty$. Hence we are allowed to neglect the first
term on the r.h.s. of Eq. (29).

In order to calculate the last remaining term
$\langle\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\rangle$ of Eq.
(29), we need to find the average $\langle
e^{ily_{n}}\eta_{\mathrm{c}}\rangle$ in analogy to Eq. (31) for
$\eta_{\mathrm{c}}$. We exploit the following relation for stochastic
variables $a$ and $b$:

$\displaystyle\langle e^{ia}b\rangle$
$\displaystyle=-i\left(\frac{d}{dr}\left\langle
e^{i(a+rb)}\right\rangle\right)_{r=0}.$ (32)

In our case, $a=ly_{n}(t_{2}-t_{1})$, $b=\eta_{\mathrm{c}}(t_{3})$. In Eq.
(32), the expectation $\langle e^{ia}b\rangle$ is expressed as the derivative
of a characteristic function $\langle e^{iuz}\rangle_{z}$ evaluated at $u=1$,
where $z=a+rb$. This characteristic function can again be expressed in terms
of the cumulants of $z$,

$\langle
e^{iuz}\rangle_{z}=\exp\left(\sum_{k=1}^{\infty}\frac{(iu)^{k}}{k!}\hat{\kappa}_{k}\right),$
(33)

where $\hat{\kappa}_{k}$ are the cumulants of $z$. The first cumulant
$\hat{\kappa}_{1}$ of $z$ vanishes. The second cumulant can be calculated from
the variances and covariances of $y_{n}$ and $\eta_{\mathrm{c}}$. Above, we
assumed that $y_{n}$ is (weakly) non-Gaussian in the presence of common input,
and thus in turn $z$ cannot be Gaussian either. However, we make here a
Gaussian approximation and neglect higher-order cumulants of $z$, as these
appear as the corrections to the corrections in Eq. (25). We then find

$\langle
e^{iuz}\rangle_{z}\approx\exp\\!\left(-\frac{u^{2}}{2}\left(\langle{a^{2}}\rangle+2r\langle{ab}\rangle+r^{2}\langle{b^{2}}\rangle\right)\right),$
(34)

and consequently for the required expression

$\displaystyle\langle e^{ily_{n}}\eta_{\mathrm{c}}\rangle$
$\displaystyle\approx il\langle{y_{n}\eta_{\mathrm{c}}}\rangle
e^{-l^{2}\langle{y_{n}^{2}}\rangle/2}.$ (35)

This term can now be inserted in the analogous version of Eq. (31) for
$\eta_{\mathrm{c}}$,

$\displaystyle\mkern-24.0mu\langle\xi_{m}(t_{1})\xi_{m}(t_{2})\eta_{\mathrm{c}}(t_{3})\rangle$
$\displaystyle=\frac{K^{2}}{N}\sum_{n}\sum_{l}|A_{l}|^{2}\Phi(l(t_{2}-t_{1}))\times$
$\displaystyle\quad\quad\langle
e^{ily_{n}(t_{2}-t_{1};t_{1})}\eta_{\mathrm{c}}(t_{3})\rangle_{\bm{\xi},\bm{\eta},\eta_{\mathrm{c}}}$
$\displaystyle\approx
K^{2}\sum_{l}|A_{l}|^{2}il\Phi(l(t_{2}-t_{1}))\langle{y_{n}(t_{2}-t_{1})\eta_{\mathrm{c}}(t_{3})}\rangle\times$
$\displaystyle\qquad\qquad
e^{-l^{2}\langle{y_{n}(t_{2}-t_{1})^{2}}\rangle/2}.$ (36)

Here, we furthermore used that $\langle e^{ily_{n}}\eta_{\mathrm{c}}\rangle$
does not depend on neuron index $n$ and therefore all $n$ terms of the sum
contribute equally. The covariance of $y_{n}$ and $\eta_{\mathrm{c}}$ can be
calculated directly using the definition of $y_{n}$ in Eq. (22) and is given
by

$\displaystyle\langle{y_{n}(t_{2}-t_{1})\eta_{\mathrm{c}}(t_{3})}\rangle=2D_{\mathrm{c}}\Theta(t_{3}-t_{1})\Theta(t_{2}-t_{3}).$
(37)

We can now evaluate the triple integral of Eq. (28) to calculate the third
cumulant $\kappa_{3}(\tau)$, which based on the above can be expressed as

$\displaystyle\kappa_{3}(\tau)$ $\displaystyle\approx
12D_{\mathrm{c}}K^{2}\sum_{l}|A_{l}|^{2}il\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\\!\\!\int\limits_{0}^{t+\tau-
t_{1}}\\!\\!\\!d\tau^{\prime}\tau^{\prime}\times$
$\displaystyle\qquad\Phi(l\tau^{\prime})\exp(-l^{2}\langle{y_{n}(\tau^{\prime})^{2}}\rangle/2).$
(38)

Here, we furthermore used that $t_{1}$ and $t_{2}$ can be permuted to
constrain the integral over $t_{2}$ to $t_{2}\geq t_{1}$, which gives an
additional factor 2, and introduced the variable $\tau^{\prime}=t_{2}-t_{1}$.
By changing the order of integration, we can eliminate the integral over
$t_{1}$, and using also the result for the second moment of Eq. (27) with the
by now familiar definition
$\Lambda(\tau)=\int_{0}^{\tau}dt(\tau-t)C_{\xi}(t)$, we eventually obtain

$\displaystyle\kappa_{3}(\tau)$ $\displaystyle\approx
12D_{\mathrm{c}}K^{2}\sum_{l}|A_{l}|^{2}il\int\limits_{0}^{\tau}\\!\\!d\tau^{\prime}(\tau-\tau^{\prime})\tau^{\prime}\times$
$\displaystyle\qquad\Phi(l\tau^{\prime})\exp(-l^{2}[\Lambda(\tau^{\prime})+(D_{\eta}+D_{\mathrm{c}})\tau^{\prime}]).$
(39)

In order to facilitate the numerical calculations and to avoid evaluating the
integral over $\tau^{\prime}$, we follow a similar strategy as in the cases
without common noise, and compute the second time derivative of
$\kappa_{3}(\tau)$,

$\displaystyle\ddot{\kappa}_{3}(\tau)$ $\displaystyle\approx
12D_{\mathrm{c}}K^{2}\sum_{l}|A_{l}|^{2}il\tau\Phi(l\tau)e^{-l^{2}[\Lambda(\tau)+(D_{\eta}+D_{\mathrm{c}})\tau]},$
(40)

which can be used to calculate $\kappa_{3}(\tau)$ along with $\Lambda(\tau)$
by integrating numerically the coupled ordinary differential equations for
both quantities. To be precise, if we stopped our approximation for
$\langle{e^{ily_{n}}}\rangle$ at the third cumulant, Eq. (23) would become

$\displaystyle\ddot{\Lambda}(\tau)$ $\displaystyle\approx
K^{2}\sum_{l}|A_{l}|^{2}\Phi(l\tau)e^{-l^{2}[\Lambda(\tau)+(D_{\eta}+D_{\mathrm{c}})\tau]-\frac{il^{3}}{6}\kappa_{3}(\tau)}.$
(41)

Eqs. (40) and (41) together with the initial conditions
$\kappa_{3}(0)=\dot{\kappa}_{3}(0)=\Lambda(0)=\dot{\Lambda}(0)=0$ (found from
Eqs. (IV.1.1) and (3)) would then constitute the self-consistent theory for
the autocorrelation function of the network noise
$C_{\xi}(\tau)=\ddot{\Lambda}(\tau)$ up to third order.

Figure 3: Rescaled cumulants $s_{3}$ (blue), $s_{4}$ (orange) and $s_{5}$ as
function of time $\tau$ for the case of common noise (solid lines) and purely
individual noise (dotted lines) obtained from network simulations, shown
together with the theoretical prediction for $s_{3}$ and $s_{4}$ the case of
common noise (dashed black lines). Parameter values for $K$ and $D$ are
indicated in the respective panels; other parameter values are $\omega_{0}=1$
and $\sigma_{\omega}=0$, and a simple sine coupling $f(\theta)=\sin\theta$ was
used. The size of the simulated networks was $N=200$; further numerical
parameters were as in Fig. 1 except that $30$ networks were simulated for a
duration of $10T_{0}$ each.

####  IV.1.2 Fourth cumulant

In order to improve our description of the non-Gaussian statistics of
$y_{n}(\tau)$, we also consider the fourth cumulant in the approximation of
the average $\langle{e^{ily_{n}}}\rangle$, see Eq. (25). The fourth cumulant
of $y_{n}(\tau)$ is given by

$\displaystyle\kappa_{4}(\tau)$ $\displaystyle=\langle
y_{n}^{4}(\tau;t)\rangle-3\langle{y_{n}^{2}(\tau;t)}\rangle^{2}.$ (42)

We first focus on the fourth moment

$\displaystyle\langle y_{n}^{4}(\tau;t)\rangle$
$\displaystyle=\left\langle\prod_{j=1}^{4}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{j}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\right\rangle.$
(43)

Following the same reasoning as above for our calculation of the third
cumulant, several terms below the quadruple integral can be discarded. First,
all terms with odd powers of $\xi_{n}$ have to vanish because of the averaging
properties of the connectivity matrix $K_{mn}$. In addition, all terms
$\langle{\eta_{n}^{p}\eta_{\mathrm{c}}^{q}}\rangle$, $p+q=4$, either vanish
for odd powers $p$ and $q$ or are exactly compensated by the corresponding
terms in the expression $3\langle{y_{n}^{2}(\tau)}\rangle^{2}$. (Note that
$\eta_{n}$ and $\eta_{\mathrm{c}}$ are Gaussian variables for which the fourth
cumulant vanishes.) The remaining terms that potentially contribute to the
fourth cumulant therefore are

$\displaystyle\langle{\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})}\rangle,$ (44)
$\displaystyle\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle,$
$\displaystyle\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle,$
$\displaystyle\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle.$

We will briefly sketch the calculation for these terms but spare the reader
the rather laborious calculations, and refer to appendix A for the details.

First, we use again that the $\xi_{n}$ can be expressed in terms of the
rotator phases,

$\xi_{n}(t)=\sum_{p,k}K_{np}|A_{k}|^{2}\exp(ik\theta_{p}(t)),$ (45)

and that
$\theta_{p}(t)=\theta_{p}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}(\omega_{p}+\xi_{p}+\eta_{p}+\eta_{\mathrm{c}})$.
For all correlators containing the product $\xi_{n}(t_{i})\xi_{n}(t_{j})$, one
eventually obtains expressions that involve the integrated input to the
rotators $y_{p}(t_{j}-t_{i})$, giving rise to averages of the type
$\langle{e^{iky_{p}(t_{2}-t_{1})}e^{ily_{q}(t_{4}-t_{3})}}\rangle$,
$\langle{e^{iky_{p}(t_{2}-t_{1})}\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$,
$\langle{e^{iky_{p}(t_{2}-t_{1})}\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$,
and
$\langle{e^{iky_{p}(t_{2}-t_{1})}\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$,
respectively. For the last three contributions, we next use a similar identity
as used above for the third cumulant (cf. Eq. (32)):

$\displaystyle\langle e^{ia}bc\rangle$
$\displaystyle=-\left(\frac{d}{dr}\frac{d}{ds}\left\langle
e^{i(a+rb+sc)}\right\rangle\right)_{r=0,s=0}.$ (46)

In all four cases of Eq. (44), we therefore have to evaluate averages
$\langle{e^{iz}}\rangle$, where for the first case
$z=ky_{p}(t_{2}-t_{1})+ly_{q}(t_{4}-t_{3})$ for the first and $z=a+rb+sc$ with
$a=ky_{p}(t_{2}-t_{1})$ and $b=\eta_{n/{\rm c}}(t_{3})$, $c=\eta_{n/{\rm
c}}(t_{4})$ for the last three contributions. To compute these averages, we
then make the second simplifying approximation discussed at the beginning of
this section, effectively assuming Gaussian statistics for the stochastic
variable $z$. In this approximation, $\langle{e^{iz}}\rangle\approx
e^{-\langle{z^{2}}\rangle/2}$, where only variances and covariances of
$y_{p}$, $y_{q}$, $\eta_{n}$, and $\eta_{\mathrm{c}}$ appear in
$\langle{z^{2}}\rangle$, depending on the contribution we consider. For the
last three terms, the Gaussian approximation directly leads to

$\displaystyle\langle e^{ia}bc\rangle$
$\displaystyle\approx\left(\langle{bc}\rangle-\langle{ab}\rangle\langle{ac}\rangle\right)e^{-\langle{a^{2}}\rangle/2}$
(47)

following relation (46), with appropriately chosen $a$, $b$, and $c$.

Let us consider the first contribution,
$\langle{\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})}\rangle$, for which we have

$\langle{z^{2}}\rangle=k^{2}\langle{y_{p}^{2}(t_{2}-t_{1})}\rangle+2kl\langle{y_{p}(t_{2}-t_{1})y_{q}(t_{4}-t_{3})}\rangle\\\
+l^{2}\langle{y_{q}^{2}(t_{4}-t_{3})}\rangle.$ (48)

For $p\neq q$, the covariance
$\langle{y_{p}(t_{2}-t_{1})y_{q}(t_{4}-t_{3})}\rangle$ scales with
$D_{\mathrm{c}}$ and vanishes in the absence of common input:

$\langle{y_{p}(t_{2}-t_{1})y_{q}(t_{4}-t_{3})}\rangle=2D_{\mathrm{c}}\int\limits_{t_{1}}^{t_{2}}\\!\\!dt^{\prime}\\!\\!\\!\int\limits_{t_{3}}^{t_{4}}\\!\\!dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime})$
(49)

For $p=q$, one can show that the covariances
$\langle{y_{p}(t_{2}-t_{1})y_{q}(t_{4}-t_{3})}\rangle$ do not contribute to
the final expression for $\langle{\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})}\rangle$
as the average $\langle{e^{i(y_{p}+y_{q})}}\rangle$ appears in a double sum
over $p$ and $q$ and eventually scales with $1/N$, see appendix A.

Eventually, the first contribution can be written as

$\langle{\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})}\rangle\approx
3K^{4}\sum_{k,l}g_{k}(t_{2}-t_{1})g_{l}(t_{4}-t_{3})\\\
\times\exp(-2klD_{\mathrm{c}}\int\limits_{t_{1}}^{t_{2}}\\!\\!dt^{\prime}\\!\\!\\!\int\limits_{t_{3}}^{t_{4}}\\!\\!dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime})),$
(50)

where we used the definition

$g_{k}(\tau)=|A_{k}|^{2}\Phi(k\tau)e^{-k^{2}[\Lambda(\tau)+(D_{\eta}+D_{\mathrm{c}})\tau]}$
(51)

to simplify the expression. For the fourth cumulant, we are eventually
interested in the (integrated) difference

$\kappa_{4,{\rm
I}}=\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\big{(}\langle{\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})}\rangle-\\\
3\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\xi_{n}(t_{3})\xi_{n}(t_{4})}\rangle\big{)}.$
(52)

After a somewhat lengthy calculation (see appendix A), we obtain an expression
for the second derivative of $\kappa_{4,{\rm I}}$ that can be used for the
computation of the fourth cumulant,

$\frac{d^{2}}{d\tau^{2}}\kappa_{4,{\rm I}}=24K^{4}\sum_{k,l}\Bigg{\\{}\\\
\int\limits_{0}^{\tau}dt(\tau-t)g_{k}(\tau)g_{l}(t)\left[e^{-2klD_{\mathrm{c}}t}-1\right]\
+\\\ \int\limits_{0}^{\tau}dt_{a}\int\limits_{\tau-
t_{a}}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left[e^{-2klD_{\mathrm{c}}(t_{a}+t_{b}-\tau)}-1\right]\Bigg{\\}}.$
(53)

Although not as conveniently integrated as the expression for the third
cumulant, it still allows numerical solution.

Based on relation (47), we quickly discuss the three remaining contributions
in Eq. (44). For
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$,
since $\langle{y_{p}\eta_{n}}\rangle=0$ for $p\neq n$ (the covariances
$\langle{ab}\rangle$ and $\langle{ac}\rangle$ of Eq. (47)), one can show that
this term actually does not contribute to the fourth cumulant as it is exactly
canceled by its ‘counterpart’ in $3\mu_{2}^{2}$, see appendix A. The average
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$
does not have a corresponding term in $3\mu_{2}^{2}$, but vanishes
nevertheless, as $\langle{y_{p}\eta_{n}}\rangle=0$ for $p\neq n$ (the
covariance $\langle{ab}\rangle$ of Eq. (47)) and in addition
$\langle{\eta_{n}\eta_{\mathrm{c}}}\rangle=0$ (the covariance
$\langle{bc}\rangle$ of Eq. (47)). Thus, the only remaining contribution to
the fourth cumulant comes from
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$,
for which we have to sum over terms of the type

$\big{(}\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle-k^{2}\langle{y_{p}(t_{2}-t_{1})\eta_{\mathrm{c}}(t_{3})}\rangle\langle{y_{p}(t_{2}-t_{1})\eta_{\mathrm{c}}(t_{4})}\rangle\big{)}\\\
\times e^{-\frac{k^{2}\langle{y_{p}^{2}(t_{2}-t_{1})}\rangle}{2}},$ (54)

see Eq. (47). While the first term with the covariance
$\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$ is exactly
cancelled by the corresponding product in $3\mu_{2}^{2}$, the second term does
contribute to the fourth cumulant. After another somewhat lengthy calculation
(see appendix A), the second derivative of this contribution to $\kappa_{4}$
is eventually found to be

$\displaystyle\frac{d^{2}}{d\tau^{2}}\kappa_{4,{\rm II}}$
$\displaystyle=-48D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}\tau^{2}g_{k}(\tau).$
(55)

As expected, this contribution vanishes in the absence of common noise
($D_{\mathrm{c}}=0$), consistent with deviations from Gaussian statistics
being induced solely by input common to all units. The initial conditions for
each of the two contributions to the fourth cumulant, $\kappa_{4,{\rm
I}}(0)=0$ and $\kappa_{4,{\rm II}}(0)=0$ can be deduced from the defining
quadruple integral; the first time derivatives of $\kappa_{4,{\rm I}}$ and
$\kappa_{4,{\rm II}}$ vanish as well for $\tau=0$, see appendix A.

For the convenience of the reader, we provide below the final set of equations
that define the self-consistent theory for the network noise autocorrelation
function, including approximate corrections for non-Gaussian statistics up to
the fourth cumulant:

$\displaystyle\ddot{\Lambda}(\tau)$
$\displaystyle=K^{2}\sum_{l}g_{l}(\tau)\exp\left(-\frac{il^{3}}{6}\kappa_{3}(\tau)+\frac{l^{4}}{24}\kappa_{4}(\tau)\right)$
(56) $\displaystyle\ddot{\kappa}_{3}(\tau)$
$\displaystyle=12D_{\mathrm{c}}K^{2}\sum_{l}il\tau g_{l}(\tau)$ (57)
$\displaystyle\ddot{\kappa}_{4}(\tau)$
$\displaystyle=24K^{4}\sum_{k,l}\Bigg{\\{}\int\limits_{0}^{\tau}\\!\\!dt(\tau-t)g_{k}(\tau)g_{l}(t)\left[e^{-2klD_{\mathrm{c}}t}-1\right]\
+\int\limits_{0}^{\tau}\\!\\!dt_{a}\\!\\!\int\limits_{\tau-
t_{a}}^{\tau}\\!\\!dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left[e^{-2klD_{\mathrm{c}}(t_{a}+t_{b}-\tau)}-1\right]\Bigg{\\}}$
$\displaystyle\qquad-48D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}\tau^{2}g_{k}(\tau),$
(58)

with the short-hand notation and initial conditions, respectively

$\displaystyle g_{l}(\tau)$
$\displaystyle=|A_{l}|^{2}\Phi(l\tau)\exp\left(-l^{2}[\Lambda(\tau)+(D_{\eta}+D_{\mathrm{c}})\tau]\right),$
(59) $\displaystyle\Lambda(0)$
$\displaystyle=\dot{\Lambda}(0)=\kappa_{3}(0)=\dot{\kappa}_{3}(0)=\kappa_{4}(0)=\dot{\kappa}_{4}(0)=0.$
(60)

The network noise autocorrelation is then given as before by
$C_{\xi}(\tau)=\ddot{\Lambda}(\tau)$, see Eq. (3).

###  IV.2 Effects of common Gaussian white noise

We now turn to the effects that are specific to the common noise and are
captured with our modified theory. For a fixed total noise intensity
$D_{\eta}+D_{\mathrm{c}}=D=\text{const}$, we thus compare the cases of purely
intrinsic noise ($D_{\eta}=D$, $D_{\mathrm{c}}=0$) and purely common noise
($D_{\mathrm{c}}=D$, $D_{\eta}=0$). We will restrict ourselves to a pure sine
coupling $f(\cdot)=\sin(\cdot)$ and vary the coupling strength $K$ and total
noise intensity $D$.

####  IV.2.1 Common Gaussian white noise induces non-Gaussian network
fluctuations

We have already seen in our calculation that the network fluctuations cannot
be any longer assumed to be Gaussian in the presence of input shared among the
different units. In Fig. 3 we confirm this theoretical finding by numerical
simulations, and also illustrate that our approximations for the third and
fourth cumulants work reasonably well for the considered parameters. The
cumulants are functions of the time argument $\tau$, and are shown in Fig. 3
in a rescaled version that capture the deviation from Gaussianity relative to
the (time-dependent) standard deviation $\sigma(\tau)=\kappa_{2}^{1/2}(\tau)$
of the distribution of the $y_{n}$. This corresponds to rewriting the
expansion Eq. (24) as

$\langle
e^{ily_{n}(\tau)}\rangle=\exp\left(\sum_{k=1}^{\infty}(il\sigma(\tau))^{k}s_{k}\right),$

where $s_{k}=\kappa_{k}/(\kappa_{2}^{k/2}k!)$ are the rescaled cumulants shown
in Fig. 3. Because $\langle{y_{n}}\rangle=0$ and by the definition of the
variance, we have $s_{1}\equiv 0$ and $s_{2}\equiv 1/2$. For $k>2$, the
$s_{k}$ can be thought of as generalizations of the skewness
$\kappa_{3}/\kappa_{2}^{3/2}$ and excess kurtosis $\kappa_{4}/\kappa_{2}^{2}$
to higher orders, weighted by $1/k!$ with which the $k$th cumulant enters the
series Eq. (24). For common noise, $s_{3}(\tau)$, $s_{4}(\tau)$ and
$s_{5}(\tau)$ (solid lines) are different from zero, illustrating that the
deviation from Gaussianity is significant. However, they remain always small
compared to $s_{2}\equiv 1/2$, and especially $s_{5}$ is small. This lends
further support to the basic idea of our approximation, i.e., the expansion
around the Gaussian case and the neglect of even higher cumulants.

For purely intrinsic noise, the magnitudes of the rescaled higher-order
cumulants (dotted lines in Fig. 3) are extremely small – expected fluctuations
due to finite-size effects cannot even be resolved on the scale of this figure
and are zero for all practical purposes. This is a striking confirmation of
our assumption of Gaussian network noise statistics for the case of purely
individual noise.

Figure 4: The maximum rescaled skewness $\max{}_{\tau}|s_{3}(\tau)|$ obtained
from network simulations (a) and from our theory (b) as a function of coupling
strength $K$ and common noise strength $D$. Other parameters as in Fig. 3.
Crosses indicate the parameter values used in the different panels of Fig. 3.
Figure 5: The maximum rescaled excess kurtosis $\max{}_{\tau}|s_{4}(\tau)|$
obtained from network simulations (a) and from our theory (b) as a function of
coupling strength $K$ and common noise strength $D$. Other parameters as in
Fig. 3. Crosses indicate the parameter values used in the different panels of
Fig. 3.

Our theory that involves a number of approximations provides a quantitatively
reasonable description of the numerically measured cumulants for a variety of
parameters, irrespective of the intensity of the common noise $D$ and the
strength of the coupling $K$, Fig. 3. In particular, the theory reproduces the
pronounced oscillations at low total noise (e.g. $K=0.5$ and $D=0.05$ of Fig.
3e) and the respective maxima of the time-dependent functions at all
parameters considered. Small deviations are observed at intermediate times.

To provide a more comprehensive overview of the dependence on the parameters,
we plot the maxima of $|s_{3}(\tau)|$ and $|s_{4}(\tau)|$ as functions of $K$
and $D$, in Figs. 4 and 5, respectively. The quality of our approximation
becomes apparent from the comparison of the simulation results (panels a) and
the theoretical predictions (panels b). Remarkably, these functions display
global maxima at non-vanishing but finite values of $K$ and $D$, indicating a
maximal deviation from Gaussianity at this parameter combination ($K_{\rm
max}\approx 0.6$ and $D_{\rm max}\approx 0.1$ within our resolution for both
$|s_{3}(\tau)|$ and $|s_{4}(\tau)|$ and both simulations and theory).
Specifically, $s_{3}$ and $s_{4}$ seem to vanish for all limiting cases of
vanishing or infinite $K$ or $D$, which is consistent with the the integrated
input $y_{n}$ becoming effectively Gaussian. In the following, we argue why
this is so.

For $D\to 0$, $y_{n}$ becomes Gaussian because here the system approaches the
previously studied purely autonomous network, in which fluctuations are
Gaussian van Meegen and Lindner (2018). For $K\to 0$, the integrated inputs
$y_{n}$ is completely dominated by the common noise that is Gaussian by
assumption. Thus, it is also plausible that the cumulants vanish in this
limit. For $D\to\infty$ at fixed $K$, a similar argument holds true: The
common noise will dominate the integrated input, in which the network noise is
limited in amplitude, hence also in this limit the fluctuations will be
Gaussian. Finally, for $K\to\infty$ at fixed $D$, the dynamics of the network
approaches that of a purely autonomous network without external input, simply
because the network fluctuations are orders of magnitude stronger than the
common noise.

In summary, significantly non-vanishing values of $s_{3}$ and $s_{4}$ appear
as predicted by the theory when the network is subject to a common noise,
whereas a purely individual noise cannot evoke any non-Gaussian features of
the network fluctuations. Our theory captures the non-Gaussian statistics
reasonably well.

Figure 6: Power spectra $S_{\xi}$ and $S_{x}$ of the network noise $\xi_{m}$
(a,c) and the rotators $x_{m}$ (b,d), respectively, for purely common ($c=1$,
i.e. $D_{\mathrm{c}}=D,D_{\eta}=0$) and purely intrinsic ($c=0$ i.e.
$D_{\mathrm{c}}=0,D_{\eta}=D$) noise. (a,b) $K=0.8$, $D=0.2$, (c,d) $K=0.5$,
$D=0.1$. Network simulations (solid lines) are compared with the theoretical
prediction (dashed and dotted lines). The power spectra for the theory were
computed from the numerical solution until $t_{\rm num}=125$. Other parameters
as in Fig. 3.  Figure 7:  Correlation function $C_{\xi}(\tau)$ of the network
noise $\xi_{m}$ for strong (a) and weak (b) noise. Line style, color code, and
parameters as in Fig. 6. The solution of the simpler system obtained when
including only the third cumulant, Eqs. (40) and (41) is additionally shown
(thin grey line) for comparison.

####  IV.2.2 Common Gaussian white noise increases spectral low-frequency
power and shifts peaks to higher frequencies

We can use the computed cumulants to obtain the correlation functions and, by
an additional Fourier transform, the power spectra of the network noise
$\xi_{i}$ and the individual pointers $x_{i}$. We look at those statistics for
two selected points in the parameter plane ($K$,$D$), respectively (0.8, 0.2)
and (0.5,0.1), see Fig. 6. Two effects of the common noise can be observed if
we compare to the case of purely individual noise of the same strength.
Firstly, there is an increase of power around zero frequency, both for the
network fluctuations $\xi_{m}(t)$ (left) and the individual rotators $x_{m}$
(right). This increase is roughly a factor of two for $K=0.5$, $D=0.1$, but
smaller for the other set. Because the external noise has power for all
frequencies, this increase in low frequency power is a non-trivial prediction
of our theory that is confirmed in the simulations. Power at low frequencies
corresponds to slow fluctuations, which thus seem to become more important in
the presence of common noise. The described effect stands in marked contrast
to that observed in populations of spiking neurons with global inhibition, in
which common noise can induce a spectral peak at non-vanishing frequencies and
a reduction of power around zero frequency Doiron  _et al._ (2004); Lindner
_et al._ (2005).

There is a second effect of the common noise on the spectra. For $K=0.8$,
$D=0.2$, we observe that the main spectral peak is shifted from the intrinsic
rotator frequency $\omega_{0}$ to a somewhat higher frequency; the same effect
is in principle also present at the other parameter set but much weaker. The
increase in frequency corresponds to a slightly faster oscillation of the
units. It is also predicted by our theory and confirmed by the simulations.

For completeness, in Fig. 7 we show for the two parameter sets the self-
consistent autocorrelation functions $C_{\xi}(\tau)$ of the network
fluctuations and $C_{x}(\tau)$ of the rotators, which allow for another
comparison between theory and simulations. Indeed, for both parameter sets,
deviations of the simulations from the theory appear to be minor, and taking
into account only the first four cumulants yields a good agreement. We also
take the opportunity to show how the theory would work if we would only
include the third cumulant. Deviations from the simulations are stronger when
only $\kappa_{3}$ is used, and there is a clear benefit to include the fourth
cumulant.

##  V Summary & outlook

In this paper, we have studied a network of randomly coupled phase oscillators
which are driven by intrinsic and external (common) noise. We have generalized
the theory of self-consistent correlation functions developed in Ref. van
Meegen and Lindner, 2018 for the case of the autonomous network (i.e. in the
absence of stochastic forcing). It turned out that this generalization is
rather straightforward in the case of purely intrinsic (private) noise. Here,
the network fluctuations are still to a very good approximation Gaussian – in
fact, they approach Gaussian statistics in the strict thermodynamic limit of
an infinite network. Our theory works very well in describing the
autocorrelation functions of the network noise and of the individual rotators.
Furthermore, we also compared the situations in which all rotators receive a
private noise of a given intensity and in which only the observed rotator is
subject to noise of that intensity. The network-mediated effect of the private
noise on all rotators may have nontrivial consequences on the correlation
statistics of the observed rotator. The coherence of the oscillation may
decrease or increase by stochastic forcing of the rest of the network,
depending on the chosen parameters and the coupling function.

When the rotators are subject to common noise, the picture becomes more
complicated. We have shown that, somewhat paradoxically, the driving with a
common Gaussian noise turns a Gaussian network noise into a non-Gaussian one.
This is reminiscent of earlier observations made in the case of coupled
chaotic maps, where the mean-field ceases to be Gaussian in the thermodynamic
limit Kaneko (1990, 1992); an effect which is reproduced when independent maps
are subject to a common noise Griniasty and Hakim (1994). For the rotator
network, the non-Gaussian statistics demand a serious revision of the theory
developed by van Meegen & Lindner van Meegen and Lindner (2018). Here we put
forward a perturbative approach based on a cumulant expansion around the
Gaussian case that leads to a system of few equations for the network noise
autocorrelation and the third and fourth (time-dependent) cumulants of the
integrated network noise. The latter measure the deviations from Gaussianity
and vanish in the absence of common noise. We have evaluated these equations
for large ranges of coupling strength $K$ and noise intensity $D_{\mathrm{c}}$
for a simple sine coupling. We have found that our theory describes the
network statistics in most cases surprisingly well. Notably, looking at the
temporal maximum of the absolute values of third and fourth cumulants, we
found non-monotonic dependences of these statistics on $K$ and
$D_{\mathrm{c}}$. Specifically, there are non-vanishing values $K^{*}$ and
$D_{\mathrm{c}}^{*}$ for which the rescaled skewness $s_{3}$ and excess
kurtosis $s_{4}$ are maximal. Nonetheless, even in this case, the non-Gaussian
deviations are moderate and still captured by our theory. We recall that the
latter is based on two important assumptions: i) We consider higher-order
cumulants only up to the fourth cumulant, and ii) assume Gaussian statistics
in the calculations of the (non-Gaussian) third and fourth cumulants of the
integrated network noise. In principle, it is possible to relax these
assumptions and to calculate the fifth cumulant and/or consider non-Gaussian
contributions in the calculation of the cumulants. However, the calculations
would become even more involved and, at least for our system, the agreement
between predicted and measured correlation functions (the actual statistics of
interest) is already satisfactory with the approximations made. We note in
passing that the expansion in cumulants we propose here is somehow remindful
of the approach presented in Refs. Tyulkina  _et al._ , 2018; Goldobin  _et
al._ , 2018, where circular cumulants were used to improve the prediction of
the order parameter for Kuramoto oscillators with intrinsic noise and Lorentz-
distributed frequencies. These works were not concerned with the statistics of
the fluctuating “recurrent” input however, and the first Kuramoto order
parameter being finite, only the second cumulant was considered as a
correction in applications.

What are the effects of the non-Gaussian network fluctuations induced by
common noise at the level of the correlation functions and power spectra? We
found two effects on the spectra. For once, the common noise may shift the
main peak to higher frequencies. Secondly, the common Gaussian white, i.e.
temporally uncorrelated, noise gives rise to increased low-frequency power in
both network noise and individual rotator spectra. Both of these effects are
not easily explained, but are well captured by our self-consistent theory of
non-Gaussian fluctuations. It is still an open problem to extract these
features from our set of equations by analytical techniques, e.g. by an
appropriately simplified Fourier transformation of the equations in simple
limit cases as was done in van Meegen and Lindner (2018). We note that slow
fluctuations have also been observed to arise spontaneously in recurrent
networks of LIF neurons with sufficiently strong (current-based) synaptic
coupling in the absence of common noise Ostojic (2014); Wieland  _et al._
(2015) and one may speculate how common input may further increase low
frequency power in those networks. On a different note, the increase in low-
frequency power in our system stands in marked contrast to the reduction of
low-frequency power due to common noise in populations of integrate-and-fire
neurons with global inhibitory feedback Doiron  _et al._ (2004); Lindner  _et
al._ (2005).

Our theory may serve as a template for calculations of the self-consistent
correlation statistics for other systems of interest. A straightforward
generalization of our results is possible, for instance if the common noise is
temporally correlated, as characterized by a correlation function
$C_{\mathrm{c}}(\tau)$. More complicated will be the case of a coupling matrix
that respects Dale’s law, i.e. that a given unit $j$ is either excitatory
($K_{ij}>0$ for all $i$) or inhibitory ($K_{ij}<0$), possibly combined with a
given connection probability $p$ between units. Preliminary simulation results
show that in this case, the statistics of network fluctuations cease to be
Gaussian. Given the importance of this constraint in the neural context, a
generalization of our theory to this case is a particularly promising topic of
future research.

##  Appendix A Detailed calculation of the fourth cumulant

Here we aim to provide some details of the calculation of the fourth cumulant
$\kappa_{4}(\tau)$ of the stochastic variable $y_{n}(\tau)$. Since
$\langle{y_{n}(\tau)}\rangle=0$, the fourth cumulant is given by
$\kappa_{4}=\langle{y_{n}^{4}}\rangle-3\langle{y_{n}^{2}}\rangle^{2}$. For the
fourth moment, we have

$\displaystyle\langle y_{n}^{4}(\tau)\rangle$
$\displaystyle=\left\langle\prod_{j=1}^{4}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{j}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\right\rangle$
(61)
$\displaystyle=\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\big{(}\langle{\prod_{j=1}^{4}\xi_{n}(t_{j})}\rangle+\langle{\prod_{j=1}^{4}\eta_{n}(t_{j})}\rangle+\langle{\prod_{j=1}^{4}\eta_{\mathrm{c}}(t_{j})}\rangle+12\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$
$\displaystyle\qquad+6\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle+6\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle+6\langle{\eta_{n}(t_{1})\eta_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle\big{)},$
(62)

where we already used that odd powers of $\xi_{n}$ vanish when averaging over
the Gaussian, zero-centered connectivity matrix $K_{mn}$. (We also used that
$\langle{\eta_{n}^{3}\eta_{\mathrm{c}}}\rangle=\langle{\eta_{n}^{3}}\rangle\langle{\eta_{\mathrm{c}}}\rangle=0$,
$\langle{\eta_{n}\eta_{\mathrm{c}}^{3}}\rangle=\langle{\eta_{n}}\rangle\langle{\eta_{\mathrm{c}}^{3}}\rangle=0$
as both noise sources are uncorrelated by definition.) Of note, we can easily
match all terms but one to a corresponding term in
$3\langle{y_{n}^{2}}\rangle^{2}$:

$\displaystyle 3\langle y_{n}^{2}(\tau)\rangle^{2}$
$\displaystyle=3\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\left\langle\prod_{j=1}^{2}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\right\rangle\left\langle\prod_{j=3}^{4}[\xi_{n}(t_{j})+\eta_{n}(t_{j})+\eta_{\mathrm{c}}(t_{j})]\right\rangle$
(63)
$\displaystyle=\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\big{(}3\langle{\prod_{j=1}^{2}\xi_{n}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\xi_{n}(t_{j})}\rangle+3\langle{\prod_{j=1}^{2}\eta_{n}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\eta_{n}(t_{j})}\rangle+3\langle{\prod_{j=1}^{2}\eta_{\mathrm{c}}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\eta_{\mathrm{c}}(t_{j})}\rangle$
$\displaystyle\qquad+6\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle+6\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle+6\langle{\eta_{n}(t_{1})\eta_{n}(t_{2})}\rangle\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle\big{)}.$
(64)

For a process $y_{n}(\tau)$ with perfectly Gaussian statistics, this is
exactly the expression the fourth moment would boil down to, and
$\mu_{4}-3\mu_{2}^{2}=0$. Here, we are looking for a difference arising in the
fourth moment; this is slightly different from the third cumulant which—for a
centered process—is directly given by the third moment. All terms that only
involve powers of $\eta_{n}$ and $\eta_{\mathrm{c}}$ are purely Gaussian by
definition, and the differences of the corresponding terms vanish. We are left
with four potentially finite contributions,

$\langle y_{n}^{4}(\tau)\rangle-3\langle
y_{n}^{2}(\tau)\rangle^{2}=\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\Big{[}\underbrace{\langle{\prod_{j=1}^{4}\xi_{n}(t_{j})}\rangle-3\langle{\prod_{j=1}^{2}\xi_{n}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\xi_{n}(t_{j})}\rangle}_{\rm
I}+12\underbrace{\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle}_{\rm
II}\\\
+6\big{(}\underbrace{\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle}_{\rm
III}\big{)}+6\big{(}\underbrace{\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle}_{\rm
IV}\big{)}\Big{]},$ (65)

which we discuss separately in the following. We just mention that we can in
the following assume $t_{2}>t_{1}$ for all terms without loss of generality,
and consider the fourfold integral over the domain
$\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\int\limits_{t_{1}}^{t+\tau}\\!\\!\\!dt_{2}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{3}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}$
with the integrands being multiplied by two. For term I, we can and will
additionally assume $t_{4}>t_{3}$ without loss of generality for the same
symmetry reasons, and accordingly the multiply the integrand by four.

###  A.1
$\langle{\prod_{j=1}^{4}\xi_{n}(t_{j})}\rangle-3\langle{\prod_{j=1}^{2}\xi_{n}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\xi_{n}(t_{j})}\rangle$

By expressing the network noise again in terms of the couplings $K_{mn}$ and
phases $\theta_{n}$ based on the Fourier decomposition of the coupling
function,

$\xi_{n}(t)=\sum_{m}K_{nm}\sum_{l}A_{l}e^{ik\left(\theta_{p}(t_{0})+\int_{t_{0}}^{t}dt^{\prime}\dot{\theta}_{m}(t^{\prime})\right)},$
(66)

we obtain the following expression for the first term:

$\displaystyle\langle\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})\rangle\ $
$\displaystyle=\left\langle\sum_{m,o,p.q}K_{nm}K_{no}K_{np}K_{nq}\sum_{k,l,r,s}A_{k}A_{l}A_{r}A_{s}e^{ik\left[\theta_{m}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{m}(t^{\prime})\right]}\cdots
e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\right\rangle$
(67)
$\displaystyle=\sum_{m,o,p,q}\sum_{k,l,r,s}A_{k}A_{l}A_{r}A_{s}\left\langle
K_{nm}K_{no}K_{np}K_{nq}\right\rangle_{K}\left\langle
e^{ik\left[\theta_{m}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{m}(t^{\prime})\right]}\cdots
e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\right\rangle.$
(68)

We separately average over the disorder, which in the case of Gaussian
couplings gives

$\displaystyle\left\langle K_{nm}K_{no}K_{np}K_{nq}\right\rangle_{K}$
$\displaystyle=\frac{K^{4}}{N^{2}}\left(\delta_{mo}\delta_{pq}+\delta_{mp}\delta_{oq}+\delta_{mq}\delta_{op}\right).$
(69)

After eliminating the $\delta_{\cdot\cdot}\delta_{\cdot\cdot}$ by summing over
two of the four corresponding indices $m$, $o$, $p$ and $q$ and appropriate
renaming of the remaining indices, one thus obtains

$\displaystyle\langle\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})\rangle=\frac{K^{4}}{N^{2}}\sum_{p,q}\sum_{k,l,r,s}A_{k}A_{l}A_{r}A_{s}\Bigg{[}$
$\displaystyle\Bigg{\langle}e^{ik\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{il\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{2}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{ir\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{3}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\Bigg{\rangle}+$
$\displaystyle\Bigg{\langle}e^{ik\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{il\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{3}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{ir\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{2}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\Bigg{\rangle}+$
$\displaystyle\Bigg{\langle}e^{ik\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{il\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{ir\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{2}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{3}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\Bigg{\rangle}\Bigg{]}.$
(70)

The three terms differ in their combinations of the time arguments $t_{1}$,
$t_{2}$, $t_{3}$ and $t_{4}$, and in principle cannot be subsumed into a
single term. As they appear below a quadruple integral over all four time
arguments, we can for the ease of notation and without loss of generality
permute time indices separately for the three terms however, and consider

$\displaystyle\langle\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})\rangle$
$\displaystyle=3\frac{K^{4}}{N^{2}}\sum_{p,\ldots,s}A_{k}A_{l}A_{r}A_{s}\Bigg{\langle}e^{ik\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}e^{il\left[\theta_{p}(t_{0})+\int_{t_{0}}^{t_{2}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})\right]}\cdots
e^{is\left[\theta_{q}(t_{0})+\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})\right]}\Bigg{\rangle}$
$\displaystyle=3\frac{K^{4}}{N^{2}}\sum_{p,\ldots,s}A_{k}A_{l}A_{r}A_{s}\Bigg{\langle}e^{i(k+l)\theta_{p}(t_{0})}e^{ik\int_{t_{0}}^{t_{1}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})}e^{il\int_{t_{0}}^{t_{2}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})}e^{i(r+s)\theta_{q}(t_{0})}e^{ir\int_{t_{0}}^{t_{3}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})}e^{is\int_{t_{0}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})}\Bigg{\rangle}.$
(71)

We next average over random initial phases $\theta_{p}(t_{0})$,
$\theta_{p}(t_{0})$ at the reference time $t_{0}$, with $\left\langle
e^{i(k+l)\theta_{p}(t_{0})}\right\rangle_{\theta_{p}(t_{0})}=\delta_{k,-l}$
and the analogous relation for $\theta_{q}(t_{0})$, leading to

$\displaystyle\langle\xi_{m}(t_{1})\cdots\xi_{m}(t_{4})\rangle$
$\displaystyle=3\frac{K^{4}}{N^{2}}\sum_{p,q}\sum_{k,l}|A_{k}|^{2}|A_{l}|^{2}\left\langle
e^{ik\int_{t_{1}}^{t_{2}}dt^{\prime}\dot{\theta}_{p}(t^{\prime})}e^{il\int_{t_{3}}^{t_{4}}dt^{\prime}\dot{\theta}_{q}(t^{\prime})}\right\rangle.$
(72)

Using $\dot{\theta}_{m}=\omega_{m}+\xi_{m}+\eta_{m}+\eta_{\mathrm{c}}$ and
averaging separately over the intrinsic frequencies,
$\langle{e^{ix\omega_{i}}}\rangle_{\bm{\omega}}=\Phi(x)$ being the
characteristic function of the $\omega_{i}$, one eventually obtains

$\displaystyle\langle\xi_{m}(t_{1})\cdots\xi_{m}(t_{4})\rangle$
$\displaystyle=3\frac{K^{4}}{N^{2}}\sum_{p,q}\sum_{k,l}|A_{k}|^{2}|A_{l}|^{2}\left\langle
e^{ik\omega_{p}(t_{2}-t_{1})}\right\rangle\left\langle
e^{il\omega_{q}(t_{4}-t_{3})}\right\rangle\times$
$\displaystyle\qquad\qquad\qquad\qquad\left\langle
e^{ik\int_{t_{1}}^{t_{2}}dt^{\prime}[\xi_{p}(t^{\prime})+\eta_{p}(t^{\prime})+\eta_{\mathrm{c}}(t^{\prime})]}e^{il\int_{t_{3}}^{t_{4}}dt^{\prime}[\xi_{q}(t^{\prime})+\eta_{q}(t^{\prime})+\eta_{\mathrm{c}}(t^{\prime})]}\right\rangle$
(73)
$\displaystyle=3\frac{K^{4}}{N^{2}}\sum_{p,q}\sum_{k,l}|A_{k}|^{2}|A_{l}|^{2}\Phi(k(t_{2}-t_{1}))\Phi(l(t_{4}-t_{3}))\langle{e^{i[ky_{p}(t_{2}-t_{1};t_{1})+ly_{q}(t_{4}-t_{3};t_{3})]}}\rangle,$
(74)

where we furthermore used the definition of Eq. (22) for the $y_{p}$, $y_{q}$
that appear in the average in the last equation.

In principle, the average
$\langle{e^{i[ky_{p}(t_{2}-t_{1};t_{1})+ly_{q}(t_{4}-t_{3};t_{3})]}}\rangle$
can again be expressed in terms of the cumulants of
$z=ky_{p}(t_{2}-t_{1};t_{1})+ly_{q}(t_{4}-t_{3};t_{3})$, in perfect analogy to
Eq. (24). However, we will make here again the second approximation of our
theory and consider only the second cumulant of $z$ when evaluating the
average $\langle{e^{iz}}\rangle$, effectively assuming that the deviations
from Gaussianity of the $y_{p}$ and $y_{q}$ are small and can be neglected in
the calculation of the higher-order corrections to the average
$\langle{e^{iy_{n}}}\rangle$, i.e. the third and fourth cumulant of $y_{n}$.
This approximation allows one to express
$\langle{e^{i[ky_{p}(t_{2}-t_{1};t_{1})+ly_{q}(t_{4}-t_{3};t_{3})]}}\rangle$
only in terms of the variances and covariances of $y_{p}$ and $y_{q}$,

$\langle{e^{iz}}\rangle=e^{-\frac{l^{2}}{2}\langle{y_{p}^{2}(t_{2}-t_{1};t_{1})}\rangle}e^{-\frac{k^{2}}{2}\langle{y_{q}^{2}(t_{4}-t_{3};t_{3})}\rangle}\times\\\
e^{-kl\langle{y_{p}(t_{2}-t_{1};t_{1})y_{q}(t_{4}-t_{3};t_{3})}\rangle},$ (75)

which we can calculate within the self-consistent theory for the correlation
function of the network noise. In particular, the covariance is given by

$\displaystyle\langle{y_{p}(t_{2}-t_{1};t_{1})y_{q}(t_{4}-t_{3};t_{3})}\rangle$
$\displaystyle=\int\limits_{t_{1}}^{t_{2}}\\!\\!dt^{\prime}\int\limits_{t_{3}}^{t_{4}}\\!\\!dt^{\prime\prime}\langle{[\xi_{p}(t^{\prime})+\eta_{p}(t^{\prime})+\eta_{\mathrm{c}}(t^{\prime})][\xi_{q}(t^{\prime\prime})+\eta_{q}(t^{\prime\prime})+\eta_{\mathrm{c}}(t^{\prime\prime})]}\rangle$
(76)
$\displaystyle=\int\limits_{t_{1}}^{t_{2}}\\!\\!dt^{\prime}\int\limits_{t_{3}}^{t_{4}}\\!\\!dt^{\prime\prime}\left[\delta_{pq}\left(C_{\xi}(t^{\prime\prime}-t^{\prime})+C_{\eta}(t^{\prime\prime}-t^{\prime})\right)+C_{\mathrm{c}}(t^{\prime\prime}-t^{\prime})\right].$
(77)

(Note that the crosscorrelation
$\langle{\xi_{p}(t^{\prime})\xi_{q}(t^{\prime\prime})}\rangle$ vanishes for
$p\neq q$, as can be seen from
$\sum_{m,n}\langle{K_{pm}K_{qn}f(\theta_{m}(t^{\prime}))f(\theta_{n}(t^{\prime\prime}))}\rangle=\sum_{m,n}\langle{K_{pm}K_{qn}}\rangle\langle{f(\theta_{m}(t^{\prime}))f(\theta_{n}(t^{\prime\prime}))}\rangle$
and $\langle{K_{pm}K_{qn}}\rangle=0$ when $p\neq q$.)

The double sum over $p$ and $q$ of Eq. (74) contains $N$ terms with $p=q$ and
$N(N-1)$ terms with $p\neq q$, while the sum has a prefactor $1/N^{2}$. In the
large-$N$ limit, we neglect all contributions that scale with $1/N$ and keep
thus only the last term in Eq. (77). With
$C_{\mathrm{c}}(t^{\prime\prime}-t^{\prime})=2D_{\mathrm{c}}\delta(t^{\prime\prime}-t^{\prime})$
and
$\langle{y_{\\{p,q\\}}^{2}(\tau,t)}\rangle=2\Lambda(\tau)+2(D_{\eta}+D_{\mathrm{c}})\tau$
(see Eq. (27)), we thus obtain

$\langle\xi_{n}(t_{1})\cdots\xi_{n}(t_{4})\rangle=3K^{4}\sum_{k,l}g_{k}(t_{2}-t_{1})g_{l}(t_{4}-t_{3})\times\\\
e^{-2klD_{\mathrm{c}}\int\limits_{t_{1}}^{t_{2}}dt^{\prime}\int\limits_{t_{3}}^{t_{4}}dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime})},$
(78)

where we used the definition for the $g_{i}(t)$ given by Eq. (51).

Along the same lines, it is straightforward to show that the corresponding
term from $\langle{y_{n}^{2}(\tau)}\rangle^{2}$ is given by

$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\xi_{n}(t_{3})\xi_{n}(t_{4})}\rangle=K^{4}\sum_{k,l}g_{k}(t_{2}-t_{1})g_{l}(t_{4}-t_{3}),$
(79)

so that we eventually obtain the following expression for the first
contribution to the fourth cumulant:

$\displaystyle\kappa_{4,{\rm I}}(\tau)$
$\displaystyle=\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\left(\langle{\prod_{j=1}^{4}\xi_{n}(t_{j})}\rangle-3\langle{\prod_{j=1}^{2}\xi_{n}(t_{j})}\rangle\langle{\prod_{j=3}^{4}\xi_{n}(t_{j})}\rangle\right)$
$\displaystyle=12K^{4}\sum_{k,l}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\int\limits_{t_{1}}^{t+\tau}\\!\\!\\!dt_{2}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{3}\int\limits_{t_{3}}^{t+\tau}\\!\\!\\!dt_{4}g_{k}(t_{2}-t_{1})g_{l}(t_{4}-t_{3})\left(e^{-2klD_{\mathrm{c}}\int\limits_{t_{1}}^{t_{2}}dt^{\prime}\int\limits_{t_{3}}^{t_{4}}dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime})}-1\right).$
(80)

With the variable transformations $t_{a}=t_{2}-t_{1}$, $t_{b}=t_{4}-t_{3}$,
and $t_{c}=t_{3}-t_{1}$ we can rewrite this expression as

$\displaystyle\kappa_{4,{\rm I}}(\tau)$
$\displaystyle=12K^{4}\sum_{k,l}\int\limits_{0}^{\tau}\\!\\!dt_{a}\int\limits_{0}^{\tau}\\!\\!dt_{b}\\!\int\limits_{t_{a}-\tau}^{\tau-
t_{b}}\\!\\!\\!dt_{c}\,g_{k}(ta)g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}\int\limits_{0}^{t_{a}}dt^{\prime}\int\limits_{t_{c}}^{t_{c}+t_{b}}dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime})}-1\right)\int\limits_{t_{\rm
low}(\tau,t_{a},t_{b},t_{c})}^{t_{\rm
upp}(\tau,t_{a},t_{b},t_{c})}\\!\\!\\!dt_{1},$ (81)

where care has to be taken in identifying the correct integration limits
$t_{\rm low}(\tau,t_{a},t_{b},t_{c})$ and $t_{\rm
upp}(\tau,t_{a},t_{b},t_{c})$ for $t_{1}$ when changing the order of
integration. Making first a case distinction for $t_{a}>t_{b}$ and
$t_{b}>t_{a}$, respectively, and then considering $t_{c}$ relative to
$t_{a}-t_{b}$ and 0, we eventually find

$\displaystyle\int\limits_{t_{\rm low}(\tau,t_{a},t_{b},t_{c})}^{t_{\rm
upp}(\tau,t_{a},t_{b},t_{c})}dt_{1}$ $\displaystyle=\begin{cases}\tau-
t_{a}+t_{c}&t_{a}-\tau<t_{c}<\min(t_{a}-t_{b},0)\\\
\tau-\max(t_{a},t_{b})&\min(t_{a}-t_{b},0)<t_{c}<\max(t_{a}-t_{b},0)\\\ \tau-
t_{b}-t_{c}&\max(t_{a}-t_{b},0)<t_{c}<\tau-t_{b}\end{cases}.$ (82)

While the numerical evaluation of the resulting triple integral is possible _a
priori_ , we try and simplify the expression further by taking the second time
derivative with respect to $\tau$. For the first time derivatives w.r.t. the
integration limits of the triple integral over $t_{a}$, $t_{b}$, and $t_{c}$,
one can see that they vanish when we use the respective values for $t_{a}$,
$t_{b}$, and $t_{c}$ with Eq. (82). The remaining time derivative of
$\int_{t_{\rm low}}^{t_{\rm upp}}dt_{1}$ is surprisingly simple and equal to
1, as can be seen as well from Eq. (82). Thus,

$\displaystyle\frac{d}{d\tau}\kappa_{4,{\rm I}}(\tau)$
$\displaystyle=12K^{4}\sum_{k,l}\int\limits_{0}^{\tau}\\!\\!dt_{a}\int\limits_{0}^{\tau}\\!\\!dt_{b}\\!\int\limits_{t_{a}-\tau}^{\tau-
t_{b}}\\!\\!\\!dt_{c}\,g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(t_{a},t_{b},t_{c})}-1\right),$
(83)

where we furthermore introduced the short-hand notation

$G(t_{a},t_{b},t_{c})=\int_{0}^{t_{a}}dt^{\prime}\int_{t_{c}}^{t_{c}+t_{b}}dt^{\prime\prime}\delta(t^{\prime\prime}-t^{\prime}).$
(84)

The second time derivative is then given by

$\displaystyle\frac{d^{2}}{d\tau^{2}}\kappa_{4,{\rm I}}(\tau)$
$\displaystyle=12K^{4}\sum_{k,l}\Big{[}\int_{0}^{\tau}dt_{b}\int_{0}^{\tau-
t_{b}}dt_{c}g_{k}(\tau)g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(\tau,t_{b},t_{c})}-1\right)+$
$\displaystyle\qquad\qquad\int_{0}^{\tau}dt_{a}\int_{t_{a}-\tau}^{0}dt_{c}g_{k}(t_{a})g_{l}(\tau)\left(e^{-2klD_{\mathrm{c}}G(t_{a},\tau,t_{c})}-1\right)+$
$\displaystyle\qquad\qquad\int_{0}^{\tau}dt_{a}\int_{0}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(t_{a},t_{b},\tau-
t_{b})}-1\right)+$
$\displaystyle\qquad\qquad\int_{0}^{\tau}dt_{a}\int_{0}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(t_{a},t_{b},t_{a}-\tau)}-1\right)\Big{]}.$
(85)

Let us consider the four contributions separately.

1) $t_{a}=\tau:\quad t_{c}>0,\ t_{c}+t_{b}<\tau\ \Rightarrow\
G(\tau,t_{b},t_{c})=t_{b}$

$\int_{0}^{\tau}dt_{b}\int_{0}^{\tau-
t_{b}}dt_{c}g_{k}(\tau)g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(\tau,t_{b},t_{c})}-1\right)=\int_{0}^{\tau}dt_{b}(\tau-
t_{b})g_{k}(\tau)g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}t_{b}}-1\right)$ (86)

2) $t_{b}=\tau:\quad t_{c}<0,\ t_{c}+\tau>t_{a}\ \Rightarrow\
G(t_{a},\tau,t_{c})=t_{a}$

$\int_{0}^{\tau}dt_{a}\int_{t_{a}-\tau}^{0}dt_{c}g_{k}(t_{a})g_{l}(\tau)\left(e^{-2klD_{\mathrm{c}}G(t_{a},\tau,t_{c})}-1\right)=\int_{0}^{\tau}dt_{a}(\tau-
t_{a})g_{k}(t_{a})g_{l}(\tau)\left(e^{-2klD_{\mathrm{c}}t_{a}}-1\right)$ (87)

Note that this term is completely analogous to the previous one and because of
the $k,l$ symmetry both can be combined by simply computing one of them and
taking it twice.

3) $t_{c}=\tau-t_{b}:\quad t_{c}>0,\ t_{c}+t_{b}>t_{a}\ \Rightarrow\
G(t_{a},t_{b},\tau-t_{b})=\begin{cases}0&t_{a}+t_{b}<\tau\ (t_{c}>t_{a})\\\
t_{a}+t_{b}-\tau&t_{a}+t_{b}>\tau\end{cases}$

$\displaystyle\int_{0}^{\tau}dt_{a}\int_{0}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(t_{a},t_{b},\tau-
t_{b})}-1\right)=\int_{0}^{\tau}dt_{a}\int_{\tau-
t_{a}}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}(t_{a}+t_{b}-\tau)}-1\right)$
(88)

4) $t_{c}=t_{a}-\tau:\quad t_{c}<0,\ t_{c}+t_{b}<t_{a}\ \Rightarrow\
G(t_{a},t_{b},t_{a}-\tau)=\begin{cases}0&t_{a}+t_{b}<\tau\ (t_{c}+t_{b}<0)\\\
t_{a}+t_{b}-\tau&t_{a}+t_{b}>\tau\end{cases}$

$\displaystyle\int_{0}^{\tau}dt_{a}\int_{0}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}G(t_{a},t_{b},t_{a}-\tau)}-1\right)=\int_{0}^{\tau}dt_{a}\int_{\tau-
t_{a}}^{\tau}dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}(t_{a}+t_{b}-\tau)}-1\right)$
(89)

This term is obviously identical to the previous one.

Based on the above, the four terms can be regrouped into a simple integral
from $0\ldots\tau$ and a double integral,

$\displaystyle\frac{d^{2}}{d\tau^{2}}\kappa_{4,{\rm I}}(\tau)$
$\displaystyle=24K^{4}\sum_{k,l}\left[\int\limits_{0}^{\tau}dt_{b}(\tau-
t_{b})g_{k}(\tau)g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}t_{b}}-1\right)+\int\limits_{0}^{\tau}dt_{a}\int\limits_{\tau-
t_{a}}^{\tau}\\!dt_{b}g_{k}(t_{a})g_{l}(t_{b})\left(e^{-2klD_{\mathrm{c}}(t_{a}+t_{b}-\tau)}-1\right)\right],$
(90)

which is the first contribution to the expression given in Eq. (IV.1.2).

###  A.2
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$

We can express the $\xi_{n}(t_{j})$ in terms of the coupling matrix $K_{n}p$
and the phases $\theta_{p}$, and arrive along the very same lines as above at
the following expression:

$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle=\\\
\frac{K^{2}}{N}\sum_{p}\sum_{k}|A_{k}|^{2}\Phi(k(t_{2}-t_{1}))\times\\\
\langle{e^{iky_{p}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle.$
(91)

To evaluate the remaining average
$\langle{e^{iky_{p}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$,
we first use relation (46) to express it as a derivative w.r.t. $r$ and $s$ of
an average of the type $\langle{e^{iz}}\rangle$, where
$z=y_{p}(t_{2}-t_{1};t_{1})+r\eta_{n}(t_{3})+s\eta_{\mathrm{c}}(t_{4})$. Here,
we make again our second approximation and neglect higher-order cumulants of
$z$ when evaluating $\langle{e^{iz}}\rangle$, effectively treating $z$ as
Gaussian for the calculation of the higher-order cumulants of $y_{n}(\tau)$.
We can than express
$\langle{e^{iky_{p}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle$
in terms of the variances and covariances of $y_{p}(t_{2}-t_{1};t_{1})$,
$\eta_{n}(t_{3})$, and $\eta_{\mathrm{c}}(t_{4})$, see Eq. (47), and obtain

$\langle{e^{iky_{p}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle=\big{(}\langle{\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle-\\\
k^{2}\langle{y_{p}(t_{2}-t_{1})\eta_{n}(t_{3})}\rangle\langle{y_{p}(t_{2}-t_{1})\eta_{\mathrm{c}}(t_{4})}\rangle\big{)}\\\
\times e^{-\frac{k^{2}\langle{y_{p}^{2}(t_{2}-t_{1})}\rangle}{2}}.$ (92)

However, upon closer inspection one finds that none of these terms contributes
in the large-$N$ limit, as
$\langle{\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle=0$ by definition and

$\langle{y_{p}(t_{2}-t_{1})\eta_{n}(t_{3})}\rangle=\delta_{pn}2D_{\eta}\Theta(t_{3}-t_{1})\Theta(t_{2}-t_{3}).$
(93)

The sum over $p$ in Eq. (91) thus reduces to a single term that scales with
$1/N$,

$\displaystyle\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle=4D_{\eta}D_{\mathrm{c}}\frac{K^{2}}{N}\sum_{k}g_{k}(t_{2}-t_{1})k^{2}\Theta(t_{3}-t_{1})\Theta(t_{2}-t_{3})\Theta(t_{4}-t_{1})\Theta(t_{2}-t_{4}),$
(94)

and vanishes as $N\to\infty$.

###  A.3
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$

We can directly apply the intermediate results from the previous contribution,
in particular relation (47) also discussed in the main text, and write

$\displaystyle\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$
$\displaystyle=\frac{K^{2}}{N}\sum_{p}\sum_{k}|A_{k}|^{2}\Phi(k(t_{2}-t_{1}))\langle{e^{iky_{p}(t_{2}-t_{1};t_{1})}\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$
(95)
$\displaystyle=\frac{K^{2}}{N}\sum_{p}\sum_{k}g_{k}(t_{2}-t_{1})\big{(}\langle{\eta_{n}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle-k^{2}\langle{y_{p}(t_{2}-t_{1};t_{1})\eta_{n}(t_{3})}\rangle\langle{y_{p}(t_{2}-t_{1};t_{1})\eta_{n}(t_{4})}\rangle\big{)}.$
(96)

The first term is exactly canceled by
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle$
which we need to subtract to compute the cumulant. The averages
$lr{y_{p}(t_{2}-t_{1};t_{1})\eta_{n}(t_{3})}$ contribute only for $p=n$ as
discussed above, and we are left with

$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{n}(t_{3})\eta_{n}(t_{4})}\rangle\\\
=-4D_{\eta}^{2}\frac{K^{2}}{N}\sum_{k}k^{2}g_{k}(t_{2}-t_{1})k^{2}\Theta(t_{3}-t_{1})\Theta(t_{2}-t_{3})\Theta(t_{4}-t_{1})\Theta(t_{2}-t_{4}).$
(97)

This term also vanishes in the limit of $N\to\infty$ and does not contribute
to the fourth cumulant.

###  A.4
$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{c}(t_{3})\eta_{c}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{c}(t_{3})\eta_{c}(t_{4})}\rangle$

This last term follows again straightforwardly from the previous contribution,
where we have simply to replace $\eta_{n}$ by $\eta_{\mathrm{c}}$. However,
the averages $\langle{y_{p}(t_{2}-t_{1};t_{1})\eta_{n}(t_{j})}\rangle$
contribute for each $p$ in the sum which is the equivalent of Eq. (96), and we
eventually obtain

$\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{\mathrm{c}}(t_{3})\eta_{\mathrm{c}}(t_{4})}\rangle\\\
=-4D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}g_{k}(t_{2}-t_{1})\Theta(t_{3}-t_{1})\Theta(t_{2}-t_{3})\Theta(t_{4}-t_{1})\Theta(t_{2}-t_{4}).$
(98)

The contribution to the fourth cumulant is accordingly given by

$\displaystyle\kappa_{4,{\rm IV}}(\tau)$
$\displaystyle=6\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\cdots\\!\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{4}\left(\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})\eta_{c}(t_{3})\eta_{c}(t_{4})}\rangle-\langle{\xi_{n}(t_{1})\xi_{n}(t_{2})}\rangle\langle{\eta_{c}(t_{3})\eta_{c}(t_{4})}\rangle\right)$
$\displaystyle=-48D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}\int\limits_{t}^{t+\tau}\\!\\!\\!dt_{1}\int\limits_{t_{1}}^{t+\tau}\\!\\!\\!dt_{2}g_{k}(t_{2}-t_{1})(t_{2}-t_{1})^{2}.$
(99)

To further simplify the calculation of $\kappa_{4,{\rm IV}}(\tau)$, we again
take the second time derivative w.r.t. $\tau$. After the variable substitution
$t_{a}=t_{2}-t_{1}$ and changing the integration to
$\int_{t}^{t+\tau}dt_{1}\int_{0}^{t+\tau-t_{1}}dt_{a}$, it is straightforward
to obtain

$\displaystyle\frac{d}{d\tau}\kappa_{4,{\rm IV}}(\tau)$
$\displaystyle=-48D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}\int\limits_{0}^{\tau}\\!dt^{\prime}g_{k}(t^{\prime})t^{\prime
2},$ (100)

where we furthermore used the substitution $t^{\prime}=t+\tau-t_{1}$. The
second time derivative then follows immediately as

$\displaystyle\frac{d^{2}}{d\tau^{2}}\kappa_{4,{\rm IV}}(\tau)$
$\displaystyle=-48D_{\mathrm{c}}^{2}K^{2}\sum_{k}k^{2}g_{k}(\tau)\tau^{2},$
(101)

which is the second contribution to the expression given in Eq. (IV.1.2).

## References

  * Kuramoto (1984) Y. Kuramoto,  _Chemical Oscillations, Waves, and Turbulence_ (Springer, Berlin, 1984).
  * Ermentrout (1998) B. Ermentrout, “Neural networks as spatio-temporal pattern-forming systems,” Reports on Progress in Physics, 61, 353–430 (1998).
  * Grytskyy  _et al._ (2013) D. Grytskyy, T. Tetzlaff, M. Diesmann, and M. Helias, “A unified view on weakly correlated recurrent networks.” Front Comput Neurosci, 7, 131 (2013).
  * di Volo and Torcini (2018) M. di Volo and A. Torcini, “Transition from Asynchronous to Oscillatory Dynamics in Balanced Spiking Networks with Instantaneous Synapses,” PRL, 121, 128301 (2018).
  * Filatrella  _et al._ (2008) G. Filatrella, A. H. Nielsen, and N. F. Pedersen, “Analysis of a power grid using a Kuramoto-like model,” Eur Phys J B, 61, 485–491 (2008).
  * Pikovsky  _et al._ (2001) A. Pikovsky, M. Rosenblum, and J. Kurths,  _Synchronization: A Universal Concept in Nonlinear Sciences_ (Cambridge University Press, 2001).
  * 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–185 (2005).
  * Sompolinsky  _et al._ (1988) H. Sompolinsky, A. Crisanti, and H. J. Sommers, “Chaos in random neural networks,” Phys Rev Lett, 61, 259–262 (1988).
  * Brunel (2000) N. Brunel, “Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons.” J Comput Neurosci, 8, 183–208 (2000).
  * Lerchner  _et al._ (2006) A. Lerchner, U. Cristina, J. Hertz, M. Ahmadi, P. Ruffiot, and S. Enemark, “Response variability in balanced cortical networks,” Neural Comput, 18, 634–659 (2006).
  * Kriener  _et al._ (2008) B. Kriener, T. Tetzlaff, A. Aertsen, M. Diesmann, and S. Rotter, “Correlations and population dynamics in cortical networks.” Neural Comput, 20, 2185–2226 (2008).
  * Dummer  _et al._ (2014) B. Dummer, S. Wieland, and B. Lindner, “Self-consistent determination of the spike-train power spectrum in a neural network with sparse connectivity.” Front Comput Neurosci, 8, 104 (2014).
  * Ostojic (2014) S. Ostojic, “Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons,” Nat Neurosci, 17, 594–600 (2014).
  * van Meegen and Lindner (2018) A. van Meegen and B. Lindner, “Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State,” Phys Rev Lett, 121, 258302 (2018).
  * Pena  _et al._ (2018) R. F. O. Pena, S. Vellmer, D. Bernardi, A. C. Roque, and B. Lindner, “Self-Consistent Scheme for Spike-Train Power Spectra in Heterogeneous Sparse Networks.” Front Comput Neurosci, 12, 9 (2018).
  * Vellmer and Lindner (2019) S. Vellmer and B. Lindner, “Theory of spike-train power spectra for multidimensional integrate-and-fire neurons,” Physical Review Research, 1, 023024 (2019).
  * Renart  _et al._ (2010) A. Renart, J. de la Rocha, P. BarthO, L. Hollender, N. Parga, A. Reyes, and K. D. Harris, “The asynchronous state in cortical circuits.” Science, 327, 587–590 (2010).
  * Tuckwell (1989) H. C. Tuckwell,  _Stochastic Processes in the Neurosciences_ , CBMS-NSF Regional conference Series in Applied Mathematics (Society for Industrial and Applied Mathematics, 1989).
  * Koch (2004) C. Koch,  _Biophysics of Computation: Information Processing in Single Neurons_ (Oxford University Press, 2004).
  * Rusakov  _et al._ (2020) D. A. Rusakov, L. P. Savtchenko, and P. E. Latham, “Noisy Synaptic Conductance: Bug or a Feature?” Trends in Neurosciences, 43, 363–372 (2020).
  * Pikovsky and Rosenblum (2015) A. Pikovsky and M. Rosenblum, “Dynamics of globally coupled oscillators: Progress and perspectives,” Chaos, 25, 097616 (2015).
  * Schäfer  _et al._ (2018) B. Schäfer, C. Beck, K. Aihara, D. Witthaut, and M. Timme, “Non-Gaussian power grid frequency fluctuations characterized by Lévy-stable laws and superstatistics,” Nature Energy, 1–8 (2018).
  * Note (1) Note that the integral can be either analytically calculated for simple processes or determined by numerically integrating the ordinary differential equation for $\Lambda_{\eta}(\tau)$ with $\mathaccentV{ddot}07F\Lambda_{\eta}(t)=C_{\eta}(t)$, $\Lambda_{\eta}(0)=0$, $\mathaccentV{dot}05F\Lambda_{\eta}(0)=0$.
  * Doiron  _et al._ (2004) B. Doiron, B. Lindner, A. Longtin, L. Maler, and J. Bastian, “Oscillatory activity in electrosensory neurons increases with the spatial correlation of the stochastic input stimulus.” PRL, 93, 048101 (2004).
  * Lindner  _et al._ (2005) B. Lindner, B. Doiron, and A. Longtin, “Theory of oscillatory firing induced by spatially correlated noise and delayed inhibitory feedback,” Phys Rev E, 72, 061919 (2005).
  * Kaneko (1990) K. Kaneko, “Globally coupled chaos violates the law of large numbers but not the central-limit theorem.” PRL, 65, 1391–1394 (1990).
  * Kaneko (1992) K. Kaneko, “Mean field fluctuation of a network of chaotic elements: Remaining fluctuation and correlation in the large size limit,” Physica D, 55, 368–384 (1992).
  * Griniasty and Hakim (1994) M. Griniasty and V. Hakim, “Correlations and dynamics in ensembles of maps: Simple models.” Phys Rev E, 49, 2661–2667 (1994).
  * Tyulkina  _et al._ (2018) I. V. Tyulkina, D. S. Goldobin, L. S. Klimenko, and A. Pikovsky, “Dynamics of Noisy Oscillator Populations beyond the Ott-Antonsen Ansatz,” PRL, 120, 264101 (2018).
  * Goldobin  _et al._ (2018) D. S. Goldobin, I. V. Tyulkina, L. S. Klimenko, and A. Pikovsky, “Collective mode reductions for populations of coupled noisy oscillators,” Chaos, 28, 101101 (2018).
  * Wieland  _et al._ (2015) S. Wieland, D. Bernardi, T. Schwalger, and B. Lindner, “Slow fluctuations in recurrent networks of spiking neurons,” Phys Rev E Stat Nonlin Soft Matter Phys, 92, 2870–5 (2015).