Finite-time Correlations Boost Large Voltage-Angle Fluctuations in Electric Power Grids

Decarbonization in the energy sector has been accompanied by an increased penetration of new renewable energy sources in electric power systems. Such sources differ from traditional productions in that, first, they induce larger, undispatchable fluctuations in power generation and second, they lack inertia. Therefore, substituting new renewables for traditional generation induces stronger and more frequent disturbances and modifies the way disturbances propagate across AC electric power grids. Recent measurements have indeed reported long, non-Gaussian tails in the distribution of local grid-frequency data. Large frequency deviations may induce grid instabilities, leading in worst-case scenarios to cascading failures and large-scale blackouts. In this manuscript, we investigate how correlated noise disturbances, characterized by the cumulants of their distribution, propagate through meshed, high-voltage power grids. We show that for a single source of fluctuations, non-Gaussianities in the form of finite skewness and positive kurtosis of the noise distribution propagate over the entire network when the noise correlation time is larger than the network's intrinsic time scales, but that they vanish over short distances if the noise fluctuates rapidly. We furthermore show that a Berry-Esseen theorem leads to the vanishing of non-Gaussianities as the number of uncorrelated noise sources increases. Our results show that the persistence of non-Gaussian fluctuations of feed-in power have a global impact on power-grid dynamics when they fluctuate over time scales larger than the intrinsic time scales of the system, which, we argue, is the relevant regime in real power grids. Our predictions are corroborated by numerical simulations on realistic models of power grids.


I. INTRODUCTION
The fight against climate change is one of the biggest challenges currently facing humankind [1]. Globally increasing atmospheric and oceanic temperatures have been directly related to the emission of greenhouse gases [2]. Therefore, key to mitigating climate changes is our ability to reduce emissions of such gases. Of particular interest is carbon dioxide, because of its large emission volumes and century-long lifetime in the atmosphere. Decarbonization, i.e. the reduction of carbon dioxide emissions from human activities, requires a fast, fundamental shift to renewable, low-carbon energy sources and a systematic electrification of the energy sector. This will affect the operation of AC electric power grids as both productions and consumptions will change [3]. In particular, higher penetration of renewable energy sources means fluctuating and uncertain power productions [4], as well as reduced electromechanical inertia [5,6], which impacts dynamic properties of power systems. It is expected -and in fact already observed -that future power grids will be subjected more often to stronger external perturbations to which they may respond more strongly [7,8].
AC power grids are technological entities that can be modeled as network-coupled dynamical systems. Since they operate according to market rules, it is often difficult to obtain true, reliable data on their operational state. Only recently has it been possible to get access to suffi-ciently large, statistically significant frequency datasets. Analyses of these datasets have emphasized the non-Gaussian nature of frequency fluctuations in AC power grids [9][10][11][12], with distributions exhibiting long tails and large increments. The source of these large deviations is often attributed to the presence of new renewable sources of energy [4,9,13]. Large frequency deviations are an important risk factor for the stability and hence for the safety of operation of present and future AC electric power systems. It is therefore of the utmost importance to understand how non-Gaussian disturbances propagate through electric networks. Many recent papers have investigated the propagation of disturbances originating from noisy power feed-in into power grids [9,[13][14][15][16][17][18][19]. However, most of them considered Gaussian noise, with two notable exceptions. First, Ref. [9] showed analytically that the variance of noise-induced frequency fluctuations decays exponentially away from its feed-in source, while the frequency kurtosis was numerically observed to exhibit a slower, possibly power-law decay. Second, Ref. [13] conjectured that the structure of power grids amplifies non-Gaussianities in power feed-in. Different noise probability distributions have been considered, yet the influence of fundamental noise characteristics such as correlation time, or the presence of multiple, independent noise sources has been neglected so far. Below we show that these characteristics are indeed key to understanding how non-Gaussian fluctuations of voltage angles propagate through high-voltage power grids.
The short-time dynamics of power grids is commonly modelled by the swing equations [5], which are nonlinear, damped wave equations on discrete networks. Source terms, representing fluctuating power feed-in, generate voltage angle and frequency waves that spread through the system. In this manuscript we investigate the propagation of such waves through realistic power grids and characterize the noisy source terms by the cumulants of their distributions and by their correlation time τ 0 . Given a single, or several sources of non-Gaussian noise, we calculate cumulants δθ p i c , p ≤ 4 of the voltage angle distributions at any node i on the power grid, over the distribution of the noise injected at one or several nodes. Non-Gaussianities are quantified by nonzero third and fourth cumulants -skewness and kurtosis. We are particularly interested in finding how far they propagate away from the noise sources, and what is their fate in the presence of multiple independent sources of noise. First, we find that non-Gaussianities in noise disturbances propagate differently, depending on the relation between τ 0 and the intrinsic network timescales. When τ 0 is the shortest time scale, non-Gaussianities disappear over short distances relative to Gaussian fluctuations, while when τ 0 is the longest time scale, they propagate through the whole system just like Gaussian fluctuations do, independently of the distribution of inertia. This is what happens when a single noise source is present. Second, from a Berry-Esseen theorem we show that, for identically but independently distributed sources of noise, non-Gaussianities disappear with the number of noise sources in both asymptotics of short and long τ 0 . Our analytical results are corroborated by numerical simulations on realistic power grids. Compared to earlier works on noise propagation in complex synchronous networks of oscillators and power grids [9,[13][14][15][16][17][19][20][21][22][23][24][25], our manuscript (i) goes beyond the white-noise limit, and includes in particular regimes of long noise correlation time that are particularly relevant for high-voltage power networks, (ii) is based on analytical calculations, and (iii) considers the case of multiple sources of power feed-in noise. Our approach relies on a single restrictive assumption, that the non-Gaussianities can be modelled by the first few cumulants of their distribution. While this excludes Lorentzian and power-law distributions with small expo-nents, it is not an important restriction, however, since the frequency fluctuations that have been reported in power systems so far exhibit close to exponential tails [9][10][11][12].
The manuscript is organized as follows. Following this introduction, we construct our dynamical model for the dynamics of high-voltage AC power grids in Section II. Analytical results are derived and presented in Section III. We confirm them numericallly in Section IV and discuss their importance and relevance in Section V.

A. The Swing Equations
The operational state of an AC power grid is determined by complex voltages V i = |V i | exp[iφ i ] at each of the i = 1, . . . N nodes of the grid. In normal operation, voltage amplitudes are fixed not far from their rated value, and voltage angles rotate close to the rated frequency, φ i (t) = Ω 0 t + θ i (t), with Ω 0 /2π = 50 or 60Hz. Over time intervals ranging roughly from seconds to several tens of seconds, the transient dynamics of highvoltage power grids is given by the swing equations [5]. They govern the time-evolution of voltage angles θ i (t) in a frame rotating at the rated frequency. In high-voltage power grids, a standard approximation is the lossless line approximation, which neglects Ohmic losses. The swing equations then read with the inertia (m i ) and damping (d i ) parameters. The active power P i is positive for generators and negative for loads, and B ij = b ij |V i ||V j | denotes the product of the voltage magnitudes at nodes i and j with the line susceptance. In the lossless line approximation, line conductances are neglected. At equilibrium, electric power grids are synchronized network systems [26]. They lie close to an operational synchronous state where voltage angles are solutions to a set of transcendental equations called the power flow equations. Under our assumptions, they read Their solution {θ ..N corresponds to the instantaneous, synchronous operational state of the power grid.

B. Wave Propagation
We want to investigate how a local perturbation about the solution to Eq. (2) propagates across the system and influences voltage angles far away from it. Such an occurence is illustrated for the PanTaGruEl [27] model of the synchronous grid of continental Europe in Fig. 1. Initially, the system is in a steady-state solution of Eq. (2). An abrupt power loss P i → 0, corresponding to the disconnection of a large power plant in Spain, brings the system out of equilibrium. Following that perturbation, a voltage-angle wave propagates across the grid, which is represented in five consecutive color-coded snapshots in Fig. 1.
In large power grids, even the loss of large power plants is a relatively weak perturbation in a mathematical sense. For instance the European Network of Transmission System Operators for Electricity (ENTSO-E) reference incident considers the simultaneous tripping of two of the largest power plants, connected to the same bus [28]. This corresponds to less than one percent of the total power injected in the synchronous grid of continental Europe. For the case plotted in Fig. 1 of a power loss of ∆P = 900MW, corresponding to a large power plant, frequency deviations never exceed 2π∆ω = 0.12 Hz, i.e. a fraction of a percent of the rated frequency [16]. Power feed-in noises being by nature smaller than the rated power on which they are superimposed, it is therefore legitimate to investigate them through the linearization of Eq. (1) about the operational synchronous state. With where we grouped the voltage angle deviations into a vector δθ, and introduced the diagonal inertia and damping matrices, M = diag(m i ) (with m i = 0 on load nodes), D = diag(d i ) as well as the weighted network Laplacian matrix L, The perturbation generating a wave of voltage angle and frequency disturbances is encoded in the source term vector δP, whose components are non-zero at nodes where the perturbation is active. Below we consider cases of (i) a single noisy perturbation and (ii) a collection of independent, geographically distributed noisy perturbations.
Power grids have two types of nodes, corresponding to power plants and loads. They have very different dynamical inertia and damping parameters. Most loads as well as inverter-connected, new renewable sources of energy have no inertia, m i = 0, while traditional power plants have an inertia roughly proportional to their rated power output [5]. Furthermore, loads have a damping parameter significantly smaller than generators [29]. While it is crucial to incorporate these dynamic inhomogeneities in any analysis of realistic power grids, they render analytical approaches intractable. Recent works, took a perturbation theory approach to incorporate small deviations about the homogeneous case [30,31], however most are based on homogeneity assumptions [32][33][34][35]. To justify it, one often invokes a Kron-reduction of the network [36] into an effective network with modified line susceptances connecting only inertiaful, generator nodes. This transformation is based on Schur's complement formula [37], and since the reduced load nodes have no inertia and a much smaller damping term, this reduction modifies the dynamics on the generators only marginally. Once the reduction is performed, one furthermore argues that considering uniform damping and inertia, d i = d, m i = m is justified, because, only large plants, all with with large rated power are connected to the high-voltage grids we are focusing on here. Additionally, machine measurements indicate that the ratio of damping over inertia does not vary by much from one machine to another [38]. Hence, in our analytical treatment, we consider noise propagation from Eq. (3) for a Kron-reduced network with homogeneous dynamic parameters, d i = d, m i = m. However, our numerical treatment is entirely free from this assumption and is based on a realistic, inhomogeneous power-grid model. Our numerical results validate our analytical results in the particularly relevant regime of long noise correlation time.

III. DISTURBANCE WAVE PROPAGATION : ANALYTICAL APPROACH
A. Linearized swing equations and modal decomposition Eq. (3) is a damped wave equation with a source term. It is defined on a discrete, meshed complex network encoded in the Laplacian matrix L, which accordingly replaces the Laplace operator of continuous wave equations. Given a source term, we compute the moments µ p ≡ δθ p i (t → ∞) , p ≤ 4 of the distribution of angle deviations at any node i on the network. Here, t → ∞ means that the observation takes place long after the onset of the noisy perturbation, to avoid transient behaviors. To do so, we use a modal expansion of Eq. (3) over the set of eigenmodes {u α } of the Laplacian matrix L. Writing δθ i (t) = α c α (t) u α,i , Eq. (3) becomes where Lu α = λ α u α , with λ α ≥ 0, α = 1, . . . N . Eq. (5) is the differential equation for a damped, driven harmonic oscillator. It is easily solved by means of Laplace transforms. The general solution reads [27] with Γ α = γ 2 − 4λ α /m and γ = d/m . Moments µ p of voltage-angle deviations are calculated as averages over the noise distribution. From Eq. (6), µ p contains an aver-age δP i1 (t 1 )δP i2 (t 2 ) . . . δP ip (t p ) , over the product of p sources of noise inside exponential integrals. One therefore needs to specify the moments of the noise distribution. We start from a geographically uncorrelated feed-in noise on nodes labeled i 0 , whose first two moments are given by to which we add non-Gaussianities in the skewness and kurtosis of the noise distribution as where . . . c explicitely refers to the cumulant. This in particular substracts all disconnected averages such . The parameters a 3,4 characterize non-Gaussianities in the noise distribution. They correspond to skewed distributions (a 3 = 0), with tails longer (a 4 > 0) or shorter (a 4 < 0) than the normal distribution. The moments µ p are given by exponential integrals and are straightforwardly calculated. However, their exact expressions are somewhat complicated. We give them for the variance only and, for the third and fourth cumulants, discuss limiting cases of long and short correlation time.  [17]. For the synchronous grid of continental Europe, a detailed analysis based on realistic line admittances and dynamic parameters in largescale power grids gave estimates for these time scales as γ −1 2.5s, T α < 1s and γT 2 α < 0.4s ∀α. Therefore the regime of long noise correlation time is already reached for τ 0 5 − 10s, while the short correlation time regime requires τ 0 1µs [16,17].
While circuit breakers and other switches may disconnect power lines and put power plant off-line in a fraction of a second, disconnection-reconnection sequences may occur at most two to three times consecutively by design of these switches. It is hard to think of a significant noise perturbation fluctuating persistently on a time scale shorter than a few seconds. Moreover, in the Sup-plemental Material we show several examples of feed-in power fluctuations from renewables with τ 0 2 − 5 mins. Hence, we conclude that the long correlation time regime is the relevant one for high-voltage transmission grids we are interested in.
We first consider the case of a single source of noise and discuss multiple noisy nodes in paragraph III F.

C. Voltage angle variance
In the limit of large observation time, the voltage angle variance is given in Eq. (S8) of the Supplemental Material. The two limiting cases of long and short noise correlation time τ 0 give key insights on noise propagation.
First, when τ 0 is the largest time scale, the voltage angle variance at node i reads The quantity squared inside the parenthesis is the Green's function for the linear operator L, from the noise source to the observation node i. For optical or electronic waves propagating through disordered mesoscopic systems, quantities similar to δθ 2 i in Eq. (9) decay as power laws with the distance between i 0 and i, when averaged over a relatively narrow but high-lying spectral interval [39]. Eq. (9) instead corresponds to a "zeroenergy" Green's function, indicating that fluctuations with long correlation times are transmitted by a few lowlying, long-wavelength eigenmodes of L, for which the perturbative approaches of Ref.
[39] cannot be directly applied.
Second, when τ 0 is the shortest time scale, one obtains (10) In the inertialess limit, m = 0, the variance is given by a two-particle Green's function.
D. Higher voltage angle cumulants : Long correlation time regime In the limit of long correlation time, it is straightforward to show that higher cumulants behave similarly to the variance, Eq. (9), namely with the parameters a p giving the deviation from Gaussianity in the feed-in fluctuations, Eqs. (8). The most remarkable thing is that, from Eqs. (9) and (11), standardized higher cumulants are given by δθ p i c / δθ 2 i p/2 = a p , regardless of the distance between the measured node and the noise source. This is one of the main results of this paper: long correlation times boost non-Gaussian fluctuations from a single noise source so much, that they propagate and persist over the whole network. This result is in particular independent of inertia, suggesting, as previously found in Ref. [40], that disturbances with long characteristic times are affected only marginally by inertia, even when the latter is inhomogeneous. Below, we confirm the validity of this conjecture with numerical simulations on realisticially inhomogeneous power grids corroborating our theoretical predictions.

E. Higher voltage angle cumulants : Short correlation time regime
The result for the third moment is given in Eq. (S9) in the Supplemental Material. In the limit of vanishing inertia, m = 0, it gives, which, together with Eq. (10), reflects the fact that, for Kuramoto, i.e. inertialess oscillators, non-Gaussianities in the p th cumulants propagate as a p−particle Green's function in the white-noise limit. As mentioned previously, Green's functions decay exponentially with distance in disordered mesoscopic systems, however it is not clear whether this behavior applies to the "zero-energy" case considered here, nor to p-particle Green's functions. To understand better the propagation of non-Gaussianities in the short correlation time regime, we therefore numerically evaluate the expressions in Eqs. (10) and (12). Fig. 2 shows the theoretically predicted standardized skewness and kurtosis (normalized by a 3 and a 4 respectively) according to Eqs. (10) 3/2 2 of the 3−particle Green's function and the 3/2−power of the twoparticle Green's function. Right panel: Normalized ratio G4/a4G 2 2 of the 4−particle Green's function and the square of the two-particle Green's function. In the inertialess case, these ratios give the standardized skewness and kurtosis of the voltage angle fluctuations in the case of short noise correlation time. Black symbols correspond to the full PanTaGruEl model of the synchronous grid of continental Europe [16,17], red ones to a connected subsection with N=1000 nodes of PanTaGruEl and green ones to the SciGRID model of the high voltage power grid of Germany [41]. PanTaGruEl data for G4 are missing because they require prohibitively large computation times.
and (12). We consider inertialess networks, and plot the results as a function of the geodesic distance to the source of noise, for various power-grid models. In contrast to the the long correlation time prediction, both skewness and kurtosis decay fast away from the noisy node (seemingly exponentially fast) before they saturate at a constant, small value. We conclude that short noise correlation times suppress the propagation of non-Gaussianities through meshed networks over a network-dependent distance.

F. Multiple sources of noise
We finally consider the case with M distinct, independently but identically distributed sources of power feed-in fluctuations. In that case, there are M contributions similar to that in Eq. (11) to the cumulant, but M !(M − p/2)!/(p/2)! pairings of the noise sources for the moment of even order p. These latter contributions are much more numerous and they result in a Gaussian p th moment -this is the standard mechanism behind central limit and Berry-Esseen theorems [42]. For instance for the p = 4 moment in the large correlation time limit, one obtains lim τ0→0 where the factor 3 in the second line accounts for all possible pairings between the product of four noise sources, δP i l , l = 1, . . . 4. The pairing mechanism giving the second term on the right-hand side of Eq. (13) leads to the convergence of the voltage angle distribution to a Gaussian distribution, with δθ 4 i / δθ 2 i 2 → 3 [to see this, sum over the noisy nodes i 0 = 1, . . . M in Eq. (11) and compare the result with (13)]. The convergence is the same as in the Berry-Esseen theorem [42]. When M is large, this second term dominates over the first one by a factor (M − 1)/2, so that the ratio of the fourth cumulant -a measure of non-Gaussianity -to the fourth moment becomes ∝ M −1 . With the standard definition [42], non-Gaussianities disappear at a rate ∝ M −1/2 .

IV. NUMERICAL RESULTS
Our two main theoretical predictions are that, (i) non-Gaussianities propagate over the entire network when they originate from a noisy source with long correlation time, while they disappear exponentially with the distance from the source for short noise correlation time, (ii) non-Gaussianities become smaller with the number M of uncorrelated sources of noise. We confirm these two predictions by numerical integration of Eq. (1) for various networks with single or multiple sources of noise with short and long correlation times. Fig. 3 first illustrates how voltage-angle fluctuations behave as more uncorrelated sources of noise are added. Blue and red crosses correspond to long and short noise correlation times respectively, and three situations of a single (left column), 40 (middle) and 120 (right) uncorrelated sources of noise are shown. When a single noise source is present, skewness and kurtosis of voltage-angles directly reflect their value for the noise source for long noise correlation time, while they are significantly reduced for short correlation time. In the latter case, a finite skewness persists over more than half of the network, and fluctuates about zero for the rest of the network nodes. As the number of noise sources increases, both skewness and kurtosis are suppressed as voltageangles become normally distributed following the action of the Berry-Esseen theorem. Numerical data still fluctuate due to the discreteness of time steps and the finiteness of the integration time. Error bars indicate the largest data variation upon doubling of the integration time. A remarkable feature is the sign change in the longcorrelation-time skewness in the middle-left panel. It is easily understood when re-expressing the single-particle Green's function in terms of graph theoretic indicators as where Ω i0,i is the resistance distance between node i 0 and i [43], C 1 (i) = (n −1 j Ω ij ) −1 is the resistance cen- trality of node i and Kf 1 = i<j Ω i,j is the so-called Kirchhoff index. The sign of all odd-p cumulants in the long correlation time limit, see Eq. (11), is given by the p th power of the Green's function. It is therefore determined by a trade-off between the centralities of the input and measured nodes on the one hand, and the resistance distance between them on the other hand. As but one consequence, the skewness changes sign as the measurement point i is taken further and further away from i 0 , when the resistance distance Ω i0i becomes larger than the sum of the inverse node centralities in Eq. (14). Theory-simulation agreement for the UK model with homogeneous damping and inertia parameters is excellent, and well within the error bars of finite-time integration.
We next turn our attention to a larger-scale, more realistic model of a high-voltage power grid and consider the PanTaGruEl model of the synchronous grid of continental Europe [16,17]. As discussed in Section III A, intrinsic time scales in such large-scale power grids are such that the short correlation time limit corresponds to τ 0 1µs while the long correlation time regime corresponds to τ 0 5 − 10s. Persistent sources of noise therefore correspond to the long correlation time regime, which we focus on. map, and voltage-angle fluctuations are measured at the five other colored nodes. One sees first, that all voltageangle distributions are the same, up to a sign inversion δθ i → −δθ i . This corroborates our prediction of Eq. (11), according to which all standardized cumulants are the same, up to possible sign changes in odd cumulants, in the case of noise with long correlation time. The observed sign change is consistent with Eq. (14), where the blue node has the same normalized voltage-angle distribution as the source, pink node, because it is close to it and the right-hand side in Eq. (14) is dominated by the sum of the inverse centralities, C −1 1 (i 0 ) + C −1 1 (i) > Ω i0,i . All other nodes are further away and correspond to a regime where the inequality is reversed, C −1 1 (i 0 )+C −1 1 (i) < Ω i0,i and odd cumulants undergo a sign change.
In the regime of long correlation time, we saw in Section III D that voltage-angle cumulants depend neither on inertia, nor on damping, and following Ref. [40] we conjectured that the prediction of Eq. (14) also applies to cases with inhomogeneous inertia and damping. We confirm this conjecture in Fig. 5, where the normalized voltage-angle distributions also keep their  Fig. 4, but for an inhomogeneous, realistic configuration of inertia and damping in the PanTaGruEl model of the synchronous grid of continental Europe [16,17].
non-Gaussianities all over the network, regardless of the distance between source and measurement nodes, in the PanTaGruEl model with realistically inhomogeneous dynamic parameters m i and d i .
Finally, we investigate the case when multiple uncorrelated sources of noise are present. Fig. 6 confirms that, for 381 sources of power feed-in fluctuations, non-Gaussianities disappear and voltage-angle deviations become Gaussian distributed. In summary, our numerical simulations with realistic power-grid models fully confirm the theoretical predictions presented in Section III.

V. CONCLUSION
The theory we just presented has uncovered two previously neglected, yet crucial characteristics that determine how voltage-angle disturbances propagate through power grids: (i) the correlation time τ 0 , i.e. the characteristic time over which sources fluctuate, and (ii) the number of sources of fluctuations. In the white-noise limit of short τ 0 , non-Gaussianities decay with the distance from their source and saturate at small values. These values are determined by the relevant multi-particle Green's function in the limit of small inertia. In the other limit of long correlation times, non-Gaussian noise propagates through the whole network, regardless of the distance to the source and independently of inertia, leading to voltage angle fluctuations with the same non-Gaussian distribution as the feed-in power noise. These non-Gaussianities are, however, averaged out in the presence of multiple, uncorrelated sources of noise.
Modern power grids are rather resilient and in particular able, in a normal operational mode to absorb not too strong voltage angle fluctuations. Future grids will be subjected to more disturbances, in particular from new renewable energy sources. A major planification and operational concern is that electro-mechanical inertia may be significantly reduced in the future, in particular at times of large renewable power production. Our results indicate that, from the point of view of disturbance propagation, one should not be concerned by this reduction of inertia.
We finally point out that out theory applies generically to diffusively coupled agents, well beyond the power-grid models and similar second-order Kuramoto models considered in this manuscript.

ACKNOWLEDGMENTS
PJ and MT were supported by the Swiss National Science Foundation under grant No. 200020 − 182050. JH was supported through the U.S Naval Research Laboratory Karles Fellowship. We thank Julian Fritzsch for discussions on the manuscript, Laurent Pagnier for Fig. 1 and general discussions, and Will Holmgren and Tucson Electric Power for data on photovoltaic and wind turbine productions.
[36] G. Kron, Tensor Analysis of Networks (Wiley, 1939 Cumulants are then straightforwardly calculated from the noise cumulants given in Eqs. (7) and (8) in the main text.

A. Angle variance
We first calculate δθ 2 i (t) . It is straightforward to obtain in terms of the components u α,i of the eigenvectors of L. Neglecting terms decaying exponentially with time, one gets The first term in this expression makes it clear that, at the level of the variance of the angle deviation, the perturbation propagates as the square of the Green's function of the Laplacian, in the limit of long noise correlation times, τ 0 → ∞. In the other limit of short noise correlation time, τ 0 → 0, the perturbation propagates as the two-particle Green's function

B. Angle skewness
A similar calculation gives the third cumulant as The two asymptotic limits of large and small correlation time are interesting. First, for τ 0 → ∞, one gets i.e. the Green's function G 1 to the third power, up to a prefactor. Taking this result together with that for the second moment, Eq. (S3), it therefore seems to be a general property that cumulants are dominantly given by the appropriate power of Green's functions in the limit of long noise correlation times. We will confirm this conjecture below by calculating the kurtosis.
In the other limit of short noise correlation times, τ 0 → 0, one obtains i.e. a three-particle Green's function.

C. Angle kurtosis
Finally the fourth cumulant of the angle deviations is given by The structure of the expression is very similar to Eq. (S4). In particular it is easily checked that, here again, the τ 0 → ∞ limit is dominated by the fourth power of the Green's function G 1 (i 0 , i), while the other limit τ 0 → 0 gives which up to a prefactor is a 4-particle Green's function. This confirms the above conjecture that p th cumulants of angle fluctuations are given by p−particle Green's functions when the noise correlation time is vanishingly small.

VII. PROPAGATION OF NON-GAUSSIANITIES WITH INERTIA
The same approach allows one to calculate cumulants of angle deviations in the presence of homogeneous inertia terms, with Eq. (6) in the main text instead of Eq. (S1). Instead of Eq. (S3), one obtains for the variance of angle deviations, The validity of Eq. (S8) is confirmed numerically in Fig. S1, for a general case where the noise correlation is neither long, nor short, λ 2 < τ −1 0 < λ N . In the long correlation time limit, it can furthermore be shown that the p th cumulant is given by the p th power of the Green's function. Higher moments can also be calculated in the finite and short correlation time limit, however the expressions become rather long and do not bring much direct information. We here only discuss the skewness in the short correlation time limit. It reads with We easily check that, in the m = 0 limit, we recover the result for Kuramoto oscillators, u α,i0 u β,i0 u γ,i0 u α,j u β,j u γ,j d 2 (λ α + λ β + λ γ ) . (S10)

VIII. NUMERICAL SIMULATIONS AND NOISE GENERATION
All numerical simulations are done by time-evolving Eq. (1) in the main text using standard fourth-order Runge-Kutta algorithms. More interesting is the generation of time-correlated non-Gaussian noise. To do that we follow the standard procedure to generate time-correlated noise as an Ornstein-Uhlenbeck process, δṖ (t) = −τ −1 0 δP (t) + δP 0 ξ(t) . (S11) However, instead of taking ξ as a normal random variable, we take it as a non-Gaussian process. In the numerical simulations presented in the main text, we take it as a partially de-symmetrized log-normal process. This simple trick allows to generate noise sequences with finite correlation time and non-Gaussian properties. Figure S2 shows a few example of noise sequences generated using Eq. (S11). The top panel has a longer correlation time than the bottom one.

IX. CORRELATION TIME IN NEW RENEWABLE POWER FEED-IN NOISE
In Fig. S3 we show power feed-in time series for photovoltaic and wind turbine productions. Fluctuations magnitude and correlation time depend on production type, geographical location and time, i.e. seasonal and daily weather fluctuations. Here, we are interested in getting qualitative information on typical correlation time. From the data of Fig. S3 we extract the normalized correlator C(t) = δP i (t 0 )δP i (t 0 + t) / δP 2 i (t 0 ) . For the photovoltaic data, the correlator is calculated over production periods of 6 hours, individually for each of the seven days shown in the left panel of Fig. S3, while for the wind turbines, the correlator is calculated for the full duration shown in Fig. S3.
The rightmost panel of Fig. S3 shows the correlators for the three most strongly fluctuating days of PV production and for the full wind turbine production sequences. It is seen that the correlator decays over a scale of 1-2 minutes or more. Given that intrinsic time scales of large-scale power grids lie in the subsecond range, we conclude that new renewable power productions are in the long correlation time regime. [S1] M. Tyloo, T. Coletta, and Ph. Jacquod, Phys. Rev. Lett. 120 084101 (2018).