Jonas Ranft
email: jonas.ranft@ens.psl.eu

Benjamin Lindner

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

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 [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] , neural networks [START_REF] Ermentrout | Neural networks as spatio-temporal pattern-forming systems[END_REF][START_REF] Grytskyy | A unified view on weakly correlated recurrent networks[END_REF][START_REF] Di Volo | Transition from Asynchronous to Oscillatory Dynamics in Balanced Spiking Networks with Instantaneous Synapses[END_REF] , or power grids [START_REF] Filatrella | Analysis of a power grid using a Kuramoto-like model[END_REF] . 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 [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Pikovsky | Synchronization: A Universal Concept in Nonlinear Sciences[END_REF][START_REF] Acebrón | The Kuramoto model: A simple paradigm for synchronization phenomena[END_REF] . 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 [START_REF] Sompolinsky | Chaos in random neural networks[END_REF][START_REF] Brunel | Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons[END_REF][START_REF] Lerchner | Response variability in balanced cortical networks[END_REF][START_REF] Kriener | Correlations and population dynamics in cortical networks[END_REF][START_REF] Dummer | Self-consistent determination of the spike-train power spectrum in a neural network with sparse connectivity[END_REF][START_REF] Ostojic | Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons[END_REF][START_REF] Van Meegen | Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State[END_REF][START_REF] Pena | Self-Consistent Scheme for Spike-Train Power Spectra in Heterogeneous Sparse Networks[END_REF][START_REF] Vellmer | Theory of spike-train power spectra for multidimensional integrate-and-fire neurons[END_REF] . While such noisy dynamics can easily be observed e.g. in spiking network models where it manifests itself in the "asynchronous state" [START_REF] Renart | The asynchronous state in cortical circuits[END_REF] , 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 [START_REF] Van Meegen | Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State[END_REF] 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 internallygenerated 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 θ m , each driven by an intrinsic oscillation frequency ω m and receiving input from the other units in the network,

θm = ω m + ∑ n K mn f (θ 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 (•) that depends only on the 'input phases' θ n . To be specific, we will assume K mn to be a Gaussian random matrix with K mn = 0 and K mn K m n = K 2 N δ mm δ nn 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 θ m with an internally generated network noise ξ m ,

ξ m = ∑ n K mn f (θ 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 ξ (τ) ≡ ξ m (t)ξ m (t + τ) . To this end, following van Meegen & Lindner, we introduce the auxiliary function

Λ(τ) = τ 0 dt(τ -t)C ξ (t), (3)
which implies

C ξ (τ) = Λ(τ). (4)
In the case of the autonomous network described by Eq. ( 1), the function Λ(τ) is a solution of the ordinary differential equation

Λ(τ) = K 2 ∑ l |A l | 2 Φ(lτ) exp[-l 2 Λ(τ)]. (5)
Here, A l are the coefficients of a Fourier expansion of the coupling function, f (θ ) = ∑ l A l exp[ilθ ], and Φ(x) = exp[iω n x] ω n is the characteristic function of the intrinsic frequencies. The above equation for Λ follows from a selfconsistency condition for the noise autocorrelation, where the average is taken over realizations of the connectivity matrix K mn , intrinsic frequencies ω n , initial conditions, and the effective noise ξ m in a self-consistent manner.

We would like to emphasize that, while the relations of Λ(τ) in Eqs. ( 3) and ( 4) hold true throughout this manuscript, the dynamical equation that one obtains for Λ 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 Λ 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

C ξ (τ) = ∑ n,n K mn K mn f (θ n (t)) f (θ n (t + τ)) K,ω,θ 0 (6a) = ∑ n,n K mn K mn K f (θ n (t)) f (θ n (t + τ)) K,ω,θ 0 (6b) = K 2 N ∑ n f (θ n (t)) f (θ n (t + τ)) K,ω,θ 0 . ( 6c
)
Note that the disorder of the connectivity matrix in principle still contributes to the fluctuating dynamics of the θ n (t).

The self-consistency ansatz for statistics of the internallygenerated network fluctuations now consists in considering the effective stochastic dynamics θm = ω m + ξ m for the rotators in the network and to average over realizations of the network noise ξ m (t) in the following. Using the Fourier expansion of f (•) and θ m (t) = θ m (t 0 ) + t t 0 dt θm (t ), one thus obtains

C ξ (τ) K 2 /N = ∑ n ∑ l,l A l A l e ilθ n (t)+il θ n (t+τ) ξ,ω,θ 0 (7a) = ∑ n ∑ l,l A l A l e i(l+l )θ n (t 0 ) θ 0 × e il t t 0 dt θn (t )+il t+τ t 0 dt θ n (t ) ξ,ω (7b) = ∑ n ∑ l |A l | 2 e ilω n τ ω e il t+τ t dt ξ n (t ) ξ (7c) = ∑ n ∑ l |A l | 2 Φ(lτ) e ily n (τ;t) (7d) = ∑ n ∑ l |A l | 2 Φ(lτ)e -l 2 2 y n (τ;t) 2 (7e) = N ∑ l |A l | 2 Φ(lτ)e -l 2 τ 0 dt (τ-t )C ξ (t ) . (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 θ 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 ξ n . For uniformly distributed phases θ n (t 0 ), one immediately obtains e i(l+l )θ n (t 0 )

θ n (t 0 ) = δ l,l . We furthermore introduced the new variable y n (τ;t) = t+τ t dt ξ n (t ). In the limit of a large network (N 1), one can assume that the ξ n are Gaussian distributed, and thus that y n (τ;t) obeys the relation e ily n = e -l 2 y 2 n /2 . The second moment y 2 n can be found according to

y n (τ;t) 2 = t+τ t dt 1 t+τ t dt 2 ξ n (t 1 )ξ n (t 2 ) = 2 τ 0 dt(τ -t)C ξ (t) (8)
after a change of variables and using the symmetry of C ξ (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 Λ (Eq. ( 5)) is eventually obtained with the substitution Eq. (3), and the noise autocorrelation follows as

C ξ (τ) = K 2 ∑ l |A l | 2 Φ(lτ) exp[-l 2 Λ(τ)]. (9)
The autocorrelation function C x m (t) of the pointer x m (t) = e iθ m (t) for a given rotator m is then rather straightforwardly obtained as

x m (t) * x m (t + τ) = e i(θ m (t+τ)-θ m (t)) = e iω m τ e i t+τ t dt ξ m (t ) = e iω m τ-Λ(τ) .
Of note, C x m (t) explicitly depends on the particular intrinsic frequency ω m of the rotator m. Averaged over all rotators, the autocorrelation becomes C x (t) = Φ(τ)e -Λ(τ) .

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 [START_REF] Tuckwell | Stochastic Processes in the Neurosciences[END_REF][START_REF] Koch | Biophysics of Computation: Information Processing in Single Neurons[END_REF][START_REF] Rusakov | Noisy Synaptic Conductance: Bug or a Feature?[END_REF] .

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. [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF][START_REF] Pikovsky | Dynamics of globally coupled oscillators: Progress and perspectives[END_REF][START_REF] Schäfer | Non-Gaussian power grid frequency fluctuations characterized by Lévy-stable laws and superstatistics[END_REF] , but how noise affects the rotator dynamics considered here is not known. In this section, we will consider the case of intrinsic noise sources η m that are uncorrelated between different units (correlated inputs will be addressed in section IV). The model then reads

θm = ω m + ∑ n K mn f (θ n ) + η m , (10)
the recurrent input again being specified by the connectivity matrix K mn and coupling function f (•). We consider the intrinsic fluctuations to be independent Gaussian noise sources with η m = 0, η m (t)η n (t + τ) = δ mn C η (τ).

A. 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 (τ) of rotator m, the effective stochastic dynamics of which is given by

θm = ω m + ξ m + η m , (11)
where the statistics of the network noise ξ 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:

ξ m (t)η m (t ) = ∑ n K mn K f (θ n (t))η m (t ) = 0 (12)
since K mn 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 ξ (t) as follows. Along the lines of the derivation in the case without noise, Eq. (7c) now becomes

C ξ (τ) K 2 /N = ∑ n ∑ l |A l | 2 e ilω n τ ω e il t+τ t dt [ξ n (t )+η n (t )] ξ,η = ∑ n ∑ l |A l | 2 Φ(lτ) e ily n (τ;t) , (13)
where we redefined the variable

y n (τ;t) = t+τ t dt [ξ n (t ) + η n (t )]. (14)
Assuming that y n (τ;t) is Gaussian, we now find

e ily n (τ) = e -l 2 2 y n (τ) 2 (15a) = e -l 2 τ 0 dt(τ-t) [ξ n (t )+η n (t )][ξ n (t +t)+η n (t +t)] (15b) = e -l 2 τ 0 dt(τ-t)[C ξ (t)+C η (t)] . (15c)
Putting everything together, one obtains

C ξ (τ) = K 2 ∑ l |A l | 2 Φ(lτ) exp -l 2 [Λ(τ) + Λ η (τ)] . ( 16
)
Here Λ η (τ) = τ 0 dt(τt)C η (t) is given by the autocorrelation function C η (t) of the intrinsic noise sources [START_REF]Note that the integral can be either analytically calculated for simple processes or determined by numerically integrating the ordinary differential equation for Λ η (τ) with Λη (t) = C η (t)[END_REF] and the function Λ(τ) = τ 0 dt(τt)C ξ (t), defined as before, now obeys the ordinary differential equation

Λ(τ) = K 2 ∑ l |A l | 2 Φ(lτ) exp -l 2 [Λ(τ) + Λ η (τ)] . ( 17
)
Once we have calculated Λ(t), we know the autocorrelation function of the network noise C ξ (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θ m ):

C x m (τ) = exp[iω m τ -Λ(τ) -Λ η (τ)] . (18)
Similarly to the case without intrinsic noise, this can be derived using θ m (t) = θ m (t 0 ) + t t 0 dt θm (t ) 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 η m (t). For our numerical examples in the following, we consider for simplicity Gaussian white noise, for which

C η (τ) = 2D η δ (τ), implying Λ η (τ) = D η τ.

B. Effects of intrinsic white noise on the rotator autocorrelation

In order to test our theory, we simulated the rotator network with additional intrinsic white noise for various combinations of noise intensity D η 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 (τ) = exp[iω m τ -Λ(τ) -D η τ] , (19)
from which we can expect a damping of the autocorrelation function with increasing noise strength D η . 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 η n , n = m, on the autocorrelation function of rotator m enters the expression for C x m (t) via the modified network-noise autocorrelation function C ξ (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 η 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 η t , but where the network-noise autocorrelation C ξ (t) is the "intrinsic-noise free" version C ξ (t) = Λ(t) with Λ obeying Eq. ( 5).

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. 14 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 [START_REF] Tuckwell | Stochastic Processes in the Neurosciences[END_REF] , for which the assumption of independence seems clearly satisfied. Another possible interpretation of the noise term η 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. 14 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 η c is the additional input common to all units:

θm = ω m + ∑ n K mn f (θ n ) + η m + η c . ( 20
)
We assume that η c is a Gaussian noise which is uncorrelated with the intrinsic noise η m of each unit m, η c (t)η m (t + τ) ≡ 0, and that its temporal statistics are described by the autocorrelation function η c (t)η c (t + τ) = C c (τ). Furthermore, along the same lines as sketched for the intrinsic noise η m , Eq. ( 12), we note that η c is uncorrelated with the network noise ξ m , η c (t)ξ m (t + τ) = 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 ξ m (t)ξ m (t + τ) in terms of the connectivity matrix K mn , coupling function f (•), and rotator phases θ m . The additional average over independent realizations of the common noise η 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 ξ (τ) = K 2 N ∑ N ∑ l |A l | 2 e ilω n τ ω × e il t+τ t dt [ξ n (t )+η n (t )+η c (t )] ξ,η,η c . ( 21
)
By redefining the variable y n (τ;t) according to

y n (τ;t) = t+τ t dt [ξ n (t ) + η n (t ) + η c (t )] (22)
and noting that e ilω n τ ω = Φ(lτ), expression (21) can again be written more conveniently as

C ξ (τ) = K 2 N ∑ n ∑ l |A l | 2 Φ(lτ) e ily n (τ;t) ξ,η,η 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 (τ;t) were still Gaussian-distributed, it would follow that the respective expressions hold with the autocorrelation function C η (τ) of the intrinsic noise simply being replaced by the sum of the noise autocorrelation functions C η (τ) +C c (τ). 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 η (τ) = (1c)C tot (τ) and C c (τ) = cC tot (τ). 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 ξ (τ), we thus need to determine the only remaining unknown term e ily n (τ;t) ξ,η,η c in Eq. (23). Let us consider the general case of a stochastic variable y. By definition of the cumulant-generating function, we have

e ily y = exp ∞ ∑ k=1 (il) k k! κ k . (24)
Here, the κ k are the cumulants of y which we can calculate as a function of its moments y j y , j ≤ k. If y is zero-centered and Gaussian-distributed, one immediately finds e ily y = e -l 2 κ 2 /2 = e -l 2 y 2 y /2 , which led to the results obtained in the previous sections. We will see that, while y n (τ;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 e ily y by taking only the first few cumulants according to Eq. ( 24) into account, i.e.,

e ily y ≈ exp -l 2 κ 2 2 -il 3 κ 3 6 + l 4 κ 4 24 . (25)
Here, the higher-order cumulants κ 3 and κ 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 +η c and y n +y m . This strategy allows us to obtain a self-consistent equation for C ξ (τ) in the presence of common input that reproduces measured autocorrelation functions as well as higher-order cumulants of y n in stochastic network simulations.

A. Calculation of the cumulants of y n (τ;t)

The first cumulant is identical to the first moment and vanishes, (26) as ξ n = 0 because of the average over the K nn and η n = η c = 0 by definition. If the first moment vanishes, the second cumulant κ 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

y n (τ;t) = t+τ t dt [ ξ n (t ) + η n (t ) + η c (t ) ] = 0,
y n (τ) 2 = τ 0 dt(τ -t)[C ξ (t) +C η (t) +C 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

κ 3 = y 3 n , κ 4 = y 4 n -3 y 2 n 2 ,
to approximate the characteristic function in Eq. ( 23).

Third cumulant

Inserting the definition of y n , we obtain for the third cumulant

κ 3 (τ) = y 3 n (τ;t) = 3 ∏ j=1 t+τ t dt j [ξ n (t j ) + η n (t j ) + η c (t j )] . (28)
Note that κ 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 ξ n (t 1 )η n (t 2 )η c (t 3 ) or ξ n (t 1 )ξ n (t 2 )η c (t 3 ) . Because η n and η 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 ξ 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:

3 ∏ j=1 [ξ n (t j ) + η n (t j ) + η c (t j )] /3 = ξ n (t 1 )ξ n (t 2 )η n (t 3 ) + ξ n (t 1 )ξ n (t 2 )η c (t 3 ) , ( 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 ξ , C η , and C c , where the latter two are known and C ξ remains to be determined self-consistently. To be specific and in line with our treatment of η m in the previous section, we will consider both η m and η c to be Gaussian white noises with respective intensities D η and D c for the remainder of this section. The correlation functions are then given by C η (τ) = 2D η δ (τ), C c (τ) = 2D c δ (τ), which decidedly simplifies calculations. Using the definition of ξ m , we can expand the first contribution of Eq. ( 29) as

ξ m (t 1 )ξ m (t 2 )η m (t 3 ) = ∑ n,n K mn K mn f (θ n (t 1 )) f (θ n (t 2 ))η m (t 3 ) = ∑ n,n K mn K mn K f (θ n (t 1 )) f (θ n (t 2 ))η m (t 3 ) = K 2 N ∑ n f (θ n (t 1 )) f (θ n (t 2 ))η m (t 3 ) . ( 30
)
Again along the lines of Eqs. ( 6) and ( 7) but with θm = ω m + ξ m + η m + η c , we now obtain

ξ m (t 1 )ξ m (t 2 )η m (t 3 ) = K 2 N ∑ n ∑ l |A l | 2 Φ(l(t 2 -t 1 ))× e ily n (t 2 -t 1 ;t 1 ) η m (t 3 ) ξ,η,η c . ( 31
)
Instead of the characteristic function e ily n , we now need to find an expression for the average e ily n η m . We can distinguish the cases n = m and n = m. When n = m, y n (t 2t 1 )η m (t 3 ) = 0, since ξ n , η n , and η c are all uncorrelated with η m ; therefore, e ily n (t 2 -t 1 ;t 1 ) η m (t 3 ) = e ily n (t 2 -t 1 ;t 1 ) η m (t 3 ) = 0.

When n = m, the correlation e ily n (t 2 -t 1 ;t 1 ) η n (t 3 ) does not vanish but should remain bounded for large N. We thus expect that, according to Eq. ( 31), ξ n (t 1 )ξ n (t 2 )η n (t 3 ) scales like 1/N and vanishes for N → ∞. 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 ξ n (t 1 )ξ n (t 2 )η c (t 3 ) of Eq. ( 29), we need to find the average e ily n η c in analogy to Eq. ( 31) for η c . We exploit the following relation for stochastic variables a and b:

e ia b = -i d dr e i(a+rb) r=0 . (32)
In our case, a = ly n (t 2 -t 1 ), b = η c (t 3 ). In Eq. ( 32), the expectation e ia b is expressed as the derivative of a characteristic function e iuz z evaluated at u = 1, where z = a + rb. This characteristic function can again be expressed in terms of the cumulants of z,

e iuz z = exp ∞ ∑ k=1 (iu) k k! κk , (33)
where κk are the cumulants of z. The first cumulant κ1 of z vanishes. The second cumulant can be calculated from the variances and covariances of y n and η 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 higherorder cumulants of z, as these appear as the corrections to the corrections in Eq. (25). We then find

e iuz z ≈ exp - u 2 2 a 2 + 2r ab + r 2 b 2 , ( 34
)
and consequently for the required expression

e ily n η c ≈ il y n η c e -l 2 y 2 n /2 . ( 35
)
This term can now be inserted in the analogous version of Eq. ( 31) for η c ,

ξ m (t 1 )ξ m (t 2 )η c (t 3 ) = K 2 N ∑ n ∑ l |A l | 2 Φ(l(t 2 -t 1 ))× e ily n (t 2 -t 1 ;t 1 ) η c (t 3 ) ξ,η,η c ≈ K 2 ∑ l |A l | 2 ilΦ(l(t 2 -t 1 )) y n (t 2 -t 1 )η c (t 3 ) × e -l 2 y n (t 2 -t 1 ) 2 /2 . ( 36
)
Here, we furthermore used that e ily n η c does not depend on neuron index n and therefore all n terms of the sum contribute equally. The covariance of y n and η c can be calculated directly using the definition of y n in Eq. ( 22) and is given by

y n (t 2 -t 1 )η c (t 3 ) = 2D c Θ(t 3 -t 1 )Θ(t 2 -t 3 ). ( 37
)
We can now evaluate the triple integral of Eq. ( 28) to calculate the third cumulant κ 3 (τ), which based on the above can be expressed as

κ 3 (τ) ≈ 12D c K 2 ∑ l |A l | 2 il t+τ t dt 1 t+τ-t 1 0 dτ τ × Φ(lτ ) exp(-l 2 y n (τ ) 2 /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 ≥ t 1 , which gives an additional factor 2, and introduced the variable τ = t 2t 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 Λ(τ) = τ 0 dt(τt)C ξ (t), we eventually obtain

κ 3 (τ) ≈ 12D c K 2 ∑ l |A l | 2 il τ 0 dτ (τ -τ )τ × Φ(lτ ) exp(-l 2 [Λ(τ ) + (D η + D c )τ ]). (39)
In order to facilitate the numerical calculations and to avoid evaluating the integral over τ , we follow a similar strategy as in the cases without common noise, and compute the second time derivative of κ 3 (τ),

κ3 (τ) ≈ 12D c K 2 ∑ l |A l | 2 ilτΦ(lτ)e -l 2 [Λ(τ)+(D η +D c )τ] , (40)
which can be used to calculate κ 3 (τ) along with Λ(τ) by integrating numerically the coupled ordinary differential equations for both quantities. To be precise, if we stopped our approximation for e ily n at the third cumulant, Eq. ( 23) would become

Λ(τ) ≈ K 2 ∑ l |A l | 2 Φ(lτ)e -l 2 [Λ(τ)+(D η +D c )τ]-il 3 6 κ 3 (τ) . (41)
Eqs. ( 40) and (41) together with the initial conditions κ 3 (0) = κ3 (0) = Λ(0) = Λ(0) = 0 (found from Eqs. (39) and ( 3)) would then constitute the self-consistent theory for the autocorrelation function of the network noise C ξ (τ) = Λ(τ) up to third order.

Fourth cumulant

In order to improve our description of the non-Gaussian statistics of y n (τ), we also consider the fourth cumulant in the approximation of the average e ily n , see Eq. (25). The fourth cumulant of y n (τ) is given by

κ 4 (τ) = y 4 n (τ;t) -3 y 2 n (τ;t) 2 . ( 42
)
We first focus on the fourth moment

y 4 n (τ;t) = 4 ∏ j=1 t+τ t dt j [ξ n (t j ) + η n (t j ) + η c (t j )] . ( 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 ξ n have to vanish because of the averaging properties of the connectivity matrix K mn . In addition, all terms η p n η q c , p + q = 4, either vanish for odd powers p and q or are exactly compensated by the corresponding terms in the expression 3 y 2 n (τ) 2 . (Note that η n and η c are Gaussian variables for which the fourth cumulant vanishes.) The remaining terms that potentially contribute to the fourth cumulant therefore are

ξ n (t 1 ) • • • ξ n (t 4 ) , ξ n (t 1 )ξ n (t 2 )η n (t 3 )η n (t 4 ) , ξ n (t 1 )ξ n (t 2 )η c (t 3 )η c (t 4 ) , ξ n (t 1 )ξ n (t 2 )η n (t 3 )η c (t 4 ) . ( 44
)
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 ξ n can be expressed in terms of the rotator phases,

ξ n (t) = ∑ p,k K np |A k | 2 exp(ikθ p (t)), (45)
and that θ p (t) = θ p (t 0 ) + t t 0 dt (ω p + ξ p + η p + η c ). For all correlators containing the product ξ n (t i )ξ n (t j ), one eventually obtains expressions that involve the integrated input to the rotators y p (t jt i ), giving rise to averages of the type e iky p (t 2 -t 1 ) e ily q (t 4 -t 3 ) , e iky p (t 2 -t 1 ) η n (t 3 )η n (t 4 ) , e iky p (t 2 -t 1 ) η c (t 3 )η c (t 4 ) , and e iky p (t 2 -t 1 ) η n (t 3 )η c (t 4 ) , respectively. For the last three contributions, we next use a similar identity as used above for the third cumulant (cf. Eq. ( 32)):

e ia bc = - d dr d ds e i(a+rb+sc) r=0,s=0 . (46)
In all four cases of Eq. ( 44), we therefore have to evaluate averages e iz , where for the first case z = ky p (t 2 -t 1 ) + ly q (t 4t 3 ) for the first and z = a + rb + sc with a = ky p (t 2t 1 ) and b = η n/c (t 3 ), c = η n/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, e iz ≈ e -z 2 /2 , where only variances and covariances of y p , y q , η n , and η c appear in z 2 , depending on the contribution we consider. For the last three terms, the Gaussian approximation directly leads to e ia bc ≈ ( bcab ac ) e -a 2 /2 (47)

following relation (46), with appropriately chosen a, b, and c. Let us consider the first contribution, ξ n (t 1 ) • • • ξ n (t 4 ) , for which we have

z 2 = k 2 y 2 p (t 2 -t 1 ) + 2kl y p (t 2 -t 1 )y q (t 4 -t 3 ) + l 2 y 2 q (t 4 -t 3 ) . ( 48
)
For p = q, the covariance y p (t 2t 1 )y q (t 4t 3 ) scales with D c and vanishes in the absence of common input:

y p (t 2 -t 1 )y q (t 4 -t 3 ) = 2D c t 2 t 1 dt t 4 t 3 dt δ (t -t ) (49)
For p = q, one can show that the covariances y p (t 2t 1 )y q (t 4t 3 ) do not contribute to the final expression for ξ n (t 1 ) • • • ξ n (t 4 ) as the average e i(y p +y q ) 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

ξ n (t 1 ) • • • ξ n (t 4 ) ≈ 3K 4 ∑ k,l g k (t 2 -t 1 )g l (t 4 -t 3 ) × exp(-2klD c t 2 t 1 dt t 4 t 3 dt δ (t -t )), (50)
where we used the definition

g k (τ) = |A k | 2 Φ(kτ)e -k 2 [Λ(τ)+(D η +D c )τ] (51)
to simplify the expression. For the fourth cumulant, we are eventually interested in the (integrated) difference

κ 4,I = t+τ t dt 1 • • • t+τ t dt 4 ξ n (t 1 ) • • • ξ n (t 4 ) - 3 ξ n (t 1 )ξ n (t 2 ) ξ n (t 3 )ξ n (t 4 ) . ( 52
)
After a somewhat lengthy calculation (see appendix A), we obtain an expression for the second derivative of κ 4,I that can be used for the computation of the fourth cumulant,

d 2 dτ 2 κ 4,I = 24K 4 ∑ k,l τ 0 dt(τ -t)g k (τ)g l (t) e -2klD c t -1 + τ 0 dt a τ τ-t a dt b g k (t a )g l (t b ) e -2klD c (t a +t b -τ) -1 . ( 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 ξ n (t 1 )ξ n (t 2 )η n (t 3 )η n (t 4 ) , since y p η n = 0 for p = n (the covariances ab and ac 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µ [START_REF] Ermentrout | Neural networks as spatio-temporal pattern-forming systems[END_REF] 2 , see appendix A. The average ξ n (t 1 )ξ n (t 2 )η n (t 3 )η c (t 4 ) does not have a corresponding term in 3µ [START_REF] Ermentrout | Neural networks as spatio-temporal pattern-forming systems[END_REF] 2 , but vanishes nevertheless, as y p η n = 0 for p = n (the covariance ab of Eq. ( 47)) and in addition η n η c = 0 (the covariance bc of Eq. ( 47)). Thus, the only remaining contribution to the fourth cumulant comes from ξ n (t 1 )ξ n (t 2 )η c (t 3 )η c (t 4 ) , for which we have to sum over terms of the type

η c (t 3 )η c (t 4 ) -k 2 y p (t 2 -t 1 )η c (t 3 ) y p (t 2 -t 1 )η c (t 4 ) × e - k 2 y 2 p (t 2 -t 1 ) 2
, (54) see Eq. ( 47). While the first term with the covariance η c (t 3 )η c (t 4 ) is exactly cancelled by the corresponding product in 3µ 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 κ 4 is eventually found to be

d 2 dτ 2 κ 4,II = -48D 2 c K 2 ∑ k k 2 τ 2 g k (τ). ( 55
)
As expected, this contribution vanishes in the absence of common noise (D 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, κ 4,I (0) = 0 and κ 4,II (0) = 0 can be deduced from the defining quadruple integral; the first time derivatives of κ 4,I and κ 4,II vanish as well for τ = 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:

Λ(τ) = K 2 ∑ l g l (τ) exp - il 3 6 κ 3 (τ) + l 4 24 κ 4 (τ) (56) κ3 (τ) = 12D c K 2 ∑ l ilτg l (τ) (57) κ4 (τ) = 24K 4 ∑ k,l τ 0 dt(τ -t)g k (τ)g l (t) e -2klD c t -1 + τ 0 dt a τ τ-t a dt b g k (t a )g l (t b ) e -2klD c (t a +t b -τ) -1 -48D 2 c K 2 ∑ k k 2 τ 2 g k (τ), (58)
with the short-hand notation and initial conditions, respectively

g l (τ) = |A l | 2 Φ(lτ) exp -l 2 [Λ(τ) + (D η + D c )τ] , (59)
Λ(0) = Λ(0) = κ 3 (0) = κ3 (0) = κ 4 (0) = κ4 (0) = 0. ( 60
)
The network noise autocorrelation is then given as before by C ξ (τ) = Λ(τ), see Eq. ( 3).

B. 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 η + D c = D = const, we thus compare the cases of purely intrinsic noise (D η = D, D c = 0) and purely common noise (D c = D, D η = 0). We will restrict ourselves to a pure sine coupling f (•) = sin(•) and vary the coupling strength K and total noise intensity D.

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 τ, and are shown in Fig. 3 in a rescaled version that capture the deviation from Gaussianity relative to the (timedependent) standard deviation σ (τ) = κ 1/2 2 (τ) of the distribution of the y n . This corresponds to rewriting the expansion Eq. ( 24) as

e ily n (τ) = exp ∞ ∑ k=1 (ilσ (τ)) k s k , where s k = κ k /(κ k/2
2 k!) are the rescaled cumulants shown in Fig. 3. Because y n = 0 and by the definition of the variance, we have s 1 ≡ 0 and s 2 ≡ 1/2. For k > 2, the s k can be thought of as generalizations of the skewness κ 3 /κ 3/2 2 and excess kurtosis κ 4 /κ 2 2 to higher orders, weighted by 1/k! with which the kth cumulant enters the series Eq. ( 24). For common noise, s 3 (τ), s 4 (τ) and s 5 (τ) (solid lines) are different from zero, illustrating that the deviation from Gaussianity is significant. However, they remain always small compared to s 2 ≡ 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.

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 (τ)| and |s 4 (τ)| as functions of K and D, in Figs. 4 and5, 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 max ≈ 0.6 and D max ≈ 0.1 within our resolution for both |s 3 (τ)| and |s 4 (τ)| 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 → 0, y n becomes Gaussian because here the system approaches the previously studied purely autonomous network, in which fluctuations are Gaussian [START_REF] Van Meegen | Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State[END_REF] . For K → 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 → ∞ 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 → ∞ 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.

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 ξ 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 ξ 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 [START_REF] Doiron | Oscillatory activity in electrosensory neurons increases with the spatial correlation of the stochastic input stimulus[END_REF][START_REF] Lindner | Theory of oscillatory firing induced by spatially correlated noise and delayed inhibitory feedback[END_REF] .

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 ω 0 to a some- 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 num = 125. Other parameters as in Fig. 3. what 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 ξ (τ) of the network fluctuations and C x (τ) 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 κ 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. 14 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 [START_REF] Kaneko | Globally coupled chaos violates the law of large numbers but not the central-limit theorem[END_REF][START_REF] Kaneko | Mean field fluctuation of a network of chaotic elements: Remaining fluctuation and correlation in the large size limit[END_REF] ; an effect which is reproduced when independent maps are subject to a common noise [START_REF] Griniasty | Correlations and dynamics in ensembles of maps: Simple models[END_REF] . For the rotator network, the non-Gaussian statistics demand a serious revision of the theory developed by van Meegen & Lindner [START_REF] Van Meegen | Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State[END_REF] . 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 autocorrela- tion 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 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 c . Specifically, there are nonvanishing values K * and D * 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 sat-isfactory 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. 29 and 30, 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 [START_REF] Van Meegen | Self-Consistent Correlations of Randomly Coupled Rotators in the Asynchronous State[END_REF] . 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 [START_REF] Ostojic | Two types of asynchronous activity in networks of excitatory and inhibitory spiking neurons[END_REF][START_REF] Wieland | Slow fluctuations in recurrent networks of spiking neurons[END_REF] and one may speculate how common input may further increase low frequency power in those networks. On a different note, the increase in lowfrequency 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 [START_REF] Doiron | Oscillatory activity in electrosensory neurons increases with the spatial correlation of the stochastic input stimulus[END_REF][START_REF] Lindner | Theory of oscillatory firing induced by spatially correlated noise and delayed inhibitory feedback[END_REF] .

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 c (τ). 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 i j > 0 for all i) or inhibitory (K i j < 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.

ξ n (t 1 ) • • • ξ n (t 4 ) = ∑ m,o,p.q K nm K no K np K nq ∑ k,l,r,s A k A l A r A s e ik θ m (t 0 )+ t 1 t 0 dt θm (t ) • • • e is θ q (t 0 )+ t 4 t 0 dt θq (t ) (A7) = ∑ m,o,p,q ∑ k,l,r,s A k A l A r A s K nm K no K np K nq K e ik θ m (t 0 )+ t 1 t 0 dt θm (t ) • • • e is θ q (t 0 )+ t 4 t 0 dt θq (t )
.

(A8)

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

K nm K no K np K nq K = K 4 N 2 (δ mo δ pq + δ mp δ oq + δ mq δ op ) . (A9)
After eliminating the δ •• δ •• by summing over two of the four corresponding indices m, o, p and q and appropriate renaming of the remaining indices, one thus obtains .

ξ n (t 1 ) • • • ξ n (t 4 ) = K 4 N 2 ∑ p,q ∑ k,l,
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 We next average over random initial phases θ p (t 0 ), θ p (t 0 ) at the reference time t 0 , with e i(k+l)θ p (t 0 ) θ p (t 0 ) = δ k,-l and the analogous relation for θ q (t 0 ), leading to where care has to be taken in identifying the correct integration limits t low (τ,t a ,t b ,t c ) and t upp (τ,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 at b and 0, we eventually find t upp (τ,t a ,t b ,t c ) t low (τ,t a ,t b ,t c )

ξ n (t 1 ) • • • ξ n (t 4 ) = 3 K 4 N 2 ∑ p,...,s A k A l A r A s e ik θ p (t 0 )+ t 1 t 0 dt θp (t ) e il θ p (t 0 )+ t 2 t 0 dt θp (t ) • • • e is θ q (t 0 )+ t 4 t 0 dt θq (t ) = 3 K 4 N 2 ∑ p,...,
ξ m (t 1 ) • • • ξ m (t 4 ) = 3 K 4 N 2 ∑ p,q ∑ k,l |A k | 2 |A l | 2 e ik t
dt 1 =     
τt a + t c t aτ < t c < min(t at b , 0) τmax(t a ,t b ) min(t at b , 0) < t c < max(t at b , 0) τt bt c max(t at b , 0) < t c < τt b .

(A22)

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 τ. 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. (A22). The remaining time derivative of t upp t low dt 1 is surprisingly simple and equal to 1, as can be seen as well from Eq. (A22). Thus, The second time derivative is then given by as above at the following expression:

ξ n (t 1 )ξ n (t 2 )η n (t 3 )η c (t 4 ) = K 2 N ∑ p ∑ k |A k | 2 Φ(k(t 2 -t 1 )
)× e iky p (t 2 -t 1 ;t 1 ) η n (t 3 )η c (t 4 ) . (A31)

To evaluate the remaining average e iky p (t 2 -t 1 ;t 1 ) η n (t 3 )η c (t 4 ) , we first use relation (46) to express it as a derivative w.r.t. r and s of an average of the type e iz , where z = y p (t 2t 1 ;t 1 ) + rη n (t 3 ) + sη c (t 4 ). Here, we make again our second approximation and neglect higher-order cumulants of z when evaluating e iz , effectively treating z as Gaussian for the calculation of the higher-order cumulants of y n (τ). We can than express e iky p (t 2 -t 1 ;t 1 ) η n (t 3 )η c (t 4 ) in terms of the variances and covariances of y p (t 2t 1 ;t 1 ), η n (t 3 ), and η c (t 4 ), see Eq. (47), and obtain e iky p (t 2 -t 1 ;t 1 ) η n (t 3 )η c (t 4 ) = η n (t 3 )η c (t 4 )k 2 y p (t 2t 1 )η n (t 3 ) y p (t 2t 1 )η c (t 4 ) × e -k 2 y 2 p (t 2 -t 1 )

2

. (A32) However, upon closer inspection one finds that none of these terms contributes in the large-N limit, as η n (t 3 )η c (t 4 ) = 0 by definition and y p (t 2t 1 )η n (t 3 ) = δ pn 2D η Θ(t 3t 1 )Θ(t 2t 3 ). (A33)

The sum over p in Eq. (A31) thus reduces to a single term that scales with 1/N, This term also vanishes in the limit of N → ∞ and does not contribute to the fourth cumulant.

ξ

Figure 1 .

 1 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 (θ ) = sin(2θ ) + cos(3θ ) and parameters ω 0 = 1, σ ω = 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.

Figure 2 .

 2 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 η = 0.2; (b) strong noise level with K = 2, D η = 0.5. All curves are obtained from the theory, with dt = 0.001.

Figure 3 .

 3 Figure3. Rescaled cumulants s 3 (blue), s 4 (orange) and s 5 as function of time τ 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 ω 0 = 1 and σ ω = 0, and a simple sine coupling f (θ ) = sin θ was used. The size of the simulated networks was N = 200; further numerical parameters were as in Fig.1except that 30 networks were simulated for a duration of 10T 0 each.

Figure 4 .

 4 Figure 4. The maximum rescaled skewness max τ |s 3 (τ)| 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 .

 5 Figure 5. The maximum rescaled excess kurtosis max τ |s 4 (τ)| 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 6 .

 6 Figure 6. Power spectra S ξ and S x of the network noise ξ m (a,c) and the rotators x m (b,d), respectively, for purely common (c = 1, i.e. D c = D, D η = 0) and purely intrinsic (c = 0 i.e. D c = 0, D η = 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 num = 125. Other parameters as in Fig.3.

Figure 7 .

 7 Figure 7. Correlation function C ξ (τ) of the network noise ξ 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.

  k (t a )g l (t b ) e -2klD c G(t a ,t b ,t c ) -1 , (A23)where we furthermore introduced the short-hand notationG(t a ,t b ,t c ) = t a 0 dt t c +t b t c dt δ (tt ).(A24)

  s A k A l A r A s e i(k+l)θ p (t 0 ) e

	ik	t 1 t 0	dt θp (t ) e	il	t 2 t 0	dt θp (t ) e i(r+s)θ q (t 0 ) e ir	t 3 t 0	dt θq (t ) e is	t 4 t 0	dt θq (t ) .	(A11)

  Using θm = ω m + ξ m + η m + η c and averaging separately over the intrinsic frequencies, e ixω i ω = Φ(x) being the character-κ 4,I (τ) = dt 4 g k (t 2t 1 )g l (t 4t 3 )With the variable transformations t a = t 2t 1 , t b = t 4t 3 , and t c = t 3t 1 we can rewrite this expression asκ 4,I (τ) = 12K 4 ∑ upp (τ,t a ,t b ,t c )t low (τ,t a ,t b ,t c )

	t+τ		t+τ			4			2			4
	t	dt 1 • • •	t	dt 4		∏ j=1	ξ n (t j ) -3	∏ j=1	ξ n (t j )	j=3 ∏	ξ n (t j )
	= 12K 4 ∑ k,l	t t+τ dt 1	t 1 t+τ dt 2	t t+τ dt 3	t 3 t+τ					  e	-2klD c	t 2 t 1	dt	t 4 t 3	dt δ (t -t )	 -1   .	(A20)
				τ	dt a	τ	dt b	τ-t b dt c g k (ta)g l (t b )	  e -2klD c	ta 0	dt	tc+t b tc	dt δ (t -t )	-1	 	dt 1 ,	(A21)
		k,l	0	0		t a -τ
															2 t 1	dt θp (t ) e	il	t 4 t 3	dt θq (t ) .	(A12)

t

  n (t 1 )ξ n (t 2 )η n (t 3 )η c (t 4 ) = 4D η D c K 2 (t 2t 1 )k 2 Θ(t 3t 1 )Θ(t 2t 3 )Θ(t 4t 1 )Θ(t 2t 4 ),(A34)and vanishes as N → ∞.3. ξ n (t 1 )ξ n (t 2 )η n (t 3 )η n (t 4 )ξ n (t 1 )ξ n (t 2 ) η n (t 3 )η n (t 4 )We can directly apply the intermediate results from the previous contribution, in particular relation (47) also discussed in the main text, and writeξ n (t 1 )ξ n (t 2 )η n (t 3 )η n (t 4 ) = K 2 | 2 Φ(k(t 2t 1)) e iky p (t 2 -t 1 ;t 1 ) η n (t 3 )η n (t 4 ) (t 2t 1 ) η n (t 3 )η c (t 4 )k 2 y p (t 2t 1 ;t 1 )η n (t 3 ) y p (t 2t 1 ;t 1 )η n (t 4 ) . (t 1 )ξ n (t 2 ) η n (t 3 )η n (t 4 ) which we need to subtract to compute the cumulant. The averages lry p (t 2t 1 ;t 1 )η n (t 3 ) contribute only for p = n as discussed above, and we are left with ξ n (t 1 )ξ n (t 2 )η n (t 3 )η n (t 4 )ξ n (t 1 )ξ n (t 2 ) η n (t 3 )η n (t 4 ) k 2 g k (t 2t 1 )k 2 Θ(t 3t 1 )Θ(t 2t 3 )Θ(t 4t 1 )Θ(t 2t 4 ). (A37)

							(A35)
				=	K 2		(A36)
	The	first	term	is	exactly	canceled	by
	η ξ n = -4D 2	K 2

N ∑ k g k N ∑ p ∑ k |A k N ∑ p ∑ k g k N ∑ k

Here we aim to provide some details of the calculation of the fourth cumulant κ 4 (τ) of the stochastic variable y n (τ). Since y n (τ) = 0, the fourth cumulant is given by κ 4 = y [START_REF] Di Volo | Transition from Asynchronous to Oscillatory Dynamics in Balanced Spiking Networks with Instantaneous Synapses[END_REF] n -3 y

where we already used that odd powers of ξ n vanish when averaging over the Gaussian, zero-centered connectivity matrix K mn . (We also used that η 3 n η c = η 3 n η c = 0, η n η 3 c = η n η 3 c = 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 y 2 n 2 :

For a process y n (τ) with perfectly Gaussian statistics, this is exactly the expression the fourth moment would boil down to, and µ 4 -3µ 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 η n and η c are purely Gaussian by definition, and the differences of the corresponding terms vanish. We are left with four potentially finite contributions,

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 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.

By expressing the network noise again in terms of the couplings K mn and phases θ n based on the Fourier decomposition of the coupling function,

we obtain the following expression for the first term: istic function of the ω i , one eventually obtains

)) e i[ky p (t 2 -t 1 ;t 1 )+ly q (t 4 -t 3 ;t 3 )] , (A14)

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 e i[ky p (t 2 -t 1 ;t 1 )+ly q (t 4 -t 3 ;t 3 )] can again be expressed in terms of the cumulants of z = ky p (t 2t 1 ;t 1 ) + ly q (t 4t 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 e iz , 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 e iy n , i.e. the third and fourth cu-mulant of y n . This approximation allows one to express e i[ky p (t 2 -t 1 ;t 1 )+ly q (t 4 -t 3 ;t 3 )] only in terms of the variances and covariances of y p and y q , e iz = e -l 2 2 y 2 p (t 2 -t 1 ;t 1 ) e -k 2 2 y 2 q (t 4 -t 3 ;t 3 ) × e -kl y p (t 2 -t 1 ;t 1 )y q (t 4 -t 3 ;t 3 ) , (A15) which we can calculate within the self-consistent theory for the correlation function of the network noise. In particular, the covariance is given by y p (t 2t 1 ;t 1 )y q (t 4t 3 ;t 3 ) =

(Note that the crosscorrelation ξ p (t )ξ q (t ) vanishes for p = q, as can be seen from

The double sum over p and q of Eq. (A14) contains N terms with p = q and N(N -1) terms with p = 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. (A17). With C c (tt ) = 2D c δ (tt ) and y 2 {p,q} (τ,t) = 2Λ(τ) + 2(D η + D c )τ (see Eq. ( 27)), we thus obtain

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 y 2 n (τ) 2 is given by

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

Let us consider the four contributions separately.

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)

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 . . . τ and a double integral,

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

We can express the ξ n (t j ) in terms of the coupling matrix K n p and the phases θ p , and arrive along the very same lines

This last term follows again straightforwardly from the previous contribution, where we have simply to replace η n by η c . However, the averages y p (t 2t 1 ;t 1 )η n (t j ) contribute for each p in the sum which is the equivalent of Eq. (A36), and we eventually obtain

The contribution to the fourth cumulant is accordingly given by κ To further simplify the calculation of κ 4,IV (τ), we again take the second time derivative w.r.t. τ. After the variable substitution t a = t 2t 1 and changing the integration to t+τ t dt 1 t+τ-t 1 0 dt a , it is straightforward to obtain

where we furthermore used the substitution t = t +τ -t 1 . The second time derivative then follows immediately as

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