Network desynchronization by non-Gaussian fluctuations

TL;DR

This paper develops a general theory for how non-Gaussian noise causes desynchronization in oscillator networks, revealing how noise statistics influence the rate and providing a simplified modeling approach.

Contribution

It introduces a novel theoretical framework for analyzing network desynchronization under non-Gaussian fluctuations, including a reduction technique for complex networks.

Findings

Higher moments of noise affect desynchronization rates exponentially.

The theory predicts when noise accelerates or decelerates desynchronization.

A reduction method simplifies network analysis to a single oscillator in some cases.

Abstract

Many networks must maintain synchrony despite the fact that they operate in noisy environments. Important examples are stochastic inertial oscillators, which are known to exhibit fluctuations with broad tails in many applications, including electric power networks with renewable energy sources. Such non-Gaussian fluctuations can result in rare network desynchronization. Here we build a general theory for inertial oscillator network desynchronization by non-Gaussian noise. We compute the rate of desynchronization and show that higher-moments of noise enter at specific powers of coupling: either speeding up or slowing down the rate exponentially depending on how noise statistics match the statistics of a network's slowest mode. Finally, we use our theory to introduce a technique that drastically reduces the effective description of network desynchronization. Most interestingly, when instability is associated with a single edge, the reduction is to one stochastic oscillator.

Equations

$$M{\dot{v}}_i+\gamma v_i=P_i+\sum _j{K}_{ij}\sin ({\varphi }_j-{\varphi }_i),$$

$$P_i(t)={\bar{P}}_i+p_i(t),$$

$$\dot{p_i}=-\alpha p_i+{\xi }_i(t),$$

$${\xi }_i(t)=\sum _{bn}{g}_{ib}\delta (t-t_{ib}[n]).$$

$$\displaystyle\frac{\partial\rho}{\partial t}=\sum_{i}\!\Bigg{[}\!-\frac{\partial}{\partial{\phi_{i}}}\!\Big{[}v_{i}\rho\Big{]}+\frac{\partial}{\partial{p_{i}}}\!\Big{[}\alpha p_{i}\rho\Big{]}$$

$$\displaystyle-\frac{\partial}{\partial{v_{i}}}\!\Big{[}\Big{(}\!\!-\frac{\gamma v_{i}}{M}+\frac{1}{M}\big{(}P_{i}+\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\big{)}\!\Big{)}\rho\Big{]}$$

$$\displaystyle+\sum_{b}\nu_{ib}\Big{[}\rho(\bm{\phi},\mathbb{v},\mathbb{p}-g_{ib}\mathbb{1}_{i})-\rho(\bm{\phi},\mathbb{v},\mathbb{p})\Big{]}\!\Bigg{]},$$

$$\displaystyle H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})=\sum_{i}\!\!\Bigg{[}\;\lambda^{\phi}_{i}v_{i}-\alpha p_{i}\lambda^{p}_{i}$$

$$\displaystyle+\;\sum_{b}\nu_{ib}\Big{(}\!\exp\{g_{ib}\lambda_{i}^{p}\}-1\Big{)}\;\;+$$

$$\displaystyle\frac{\lambda_{i}^{v}}{M}\Big{(}\!\!-\gamma v_{i}+\bar{P}_{i}+p_{i}+\!\sum_{j}K_{ij}\sin(\phi_{j}-\phi_{i})\!\Big{)}\Bigg{]}.$$

$$\displaystyle S(\bm{\phi},\mathbb{v},\mathbb{p})=\sum_{i}\!\Bigg{[}\!\int\!\lambda_{i}^{\phi}d\phi_{i}\;+\!\int\!\lambda_{i}^{v}dv_{i}\;+\!\int\!\lambda_{i}^{p}dp_{i}\Bigg{]}.$$

$$\ln ⟨T⟩\approx S({\mathbf{\varphi }}^s,0,0)+\text{constant}.$$

$$\displaystyle C=\sqrt{\frac{2\big{|}\sum_{ij}A_{ij}\sin(\phi_{j}^{SN}\!\!-\!\phi_{i}^{SN})r_{i}\big{|}}{\big{|}\sum_{ij}A_{ij}\sin(\phi_{j}^{SN}\!\!-\!\phi_{i}^{SN}){[r_{j}\!-r_{i}]}^{2}r_{i}\big{|}}}.$$

$$\frac{v_i}{r_i}$$

$$\frac{{\lambda }_i^{\varphi }}{r_i}$$

$$\frac{{\lambda }_i^p}{r_i}$$

$$S({\mathbf{\varphi }}^s,0,0)\approx \frac{8\sqrt{2}{K}_{\text{SN}}{\kappa }^{3/2}\gamma {\alpha }^2{R}_0^{3/2}}{3R_2^{1/2}{\mu }_2}.$$

$${\Delta }^{(n)}S=$$

$$\frac{2^{n+\frac{1}{2}}\gamma {\alpha }^n{R}_0^{n-\frac{1}{2}}{\int }_{-1}^1(1-x^2{)}^{n-1}dx}{R_2^{1/2}n!}.$$

$$\dot{V}$$

$$\dot{P}$$

$${\dot{\Lambda }}_{\Psi }$$

$${\dot{\Lambda }}_V$$

$$S(\Psi ,V,P)=\int {\Lambda }_{\Psi }d\Psi +\int {\Lambda }_{V}dV+\int {\Lambda }_{P}dP.$$

$${\dot{\varphi }}_i=$$

$$M{\dot{v}}_i=$$

$${\dot{p}}_i=$$

$${\dot{\lambda }}_i^{\varphi }=$$

$${\dot{\lambda }}_i^v=$$

$${\dot{\lambda }}_i^p=$$

$$0=$$

$$0_i=\sum _j{A}_{ij}\cos ({\varphi }_j^{SN}-{\varphi }_i^{SN})[q_j-q_i],$$

$$0_i=$$

$$+\sum _j{A}_{ij}\sin ({\varphi }_j^{SN}-{\varphi }_i^{SN})$$

$$-\frac{C^2}{2}\sum _j{A}_{ij}\sin ({\varphi }_j^{SN}-{\varphi }_i^{SN})[r_j-r_i{]}^2.$$

$$z=-\sqrt{2R_0{R}_2}{\kappa }^{1/2},$$

$$0=\sum _{i\in {\mathcal{S}}_1}{P}_i+K_{\text{SN}}\sum _{i\in {\mathcal{S}}_1,j}\sin ({\varphi }_j^{SN}-{\varphi }_i^{SN}).$$

$$K_{\text{SN}}\stackrel{\text{SCSN}}{=}\sum _{i\in {\mathcal{S}}_1}{\bar{P}}_i,$$

Full text

I INTRODUCTION

Networks of coupled oscillators form the basis for complex physical, biological and technological systemsRN909 , such as smart gridsRN914 , Josephson junction arraysISI:000258975100054 ; ISI:A1982NL11700063 , optical networksRN915 , biological networksRN910 ; RN919 ; RN926 ; RN925 ; RN920 , and coupled mechanical devicesRN916 ; RN918 . From a deterministic viewpoint, most networks operate in stable attractor regimes, such as synchronized oscillations. However, in reality uncertainties and noise produce fluctuations from an attractor, which over long time scales may bring a network into dynamically unstable states and result in large deviations. Consequently, much recent attention has been given to the effects of noise on networked oscillators, e.g., stochastic escapeTimme2 ; ISI:000428996900014 ; Tyloo ; HindesKM ; DeVille , noise cancellationRonellenfitsch , noise propagationISI:000351565600008 ; Haehne ; Zhang , and synchronizationISI:000383052300003 ; ISI:000417337800032 .

A motivating application for stochastic networked oscillators is that of fluctuating power grids driven by renewable energy sources, such as wind and solar. Such fluctuations are highly non-GaussianMilan ; Anvari ; Schafer , and may significantly impact power-grid stabilityKamps ; RN935 ; Kurths2 ; Nesti . Non-Gaussian noise is understood to exponentially alter the rates for large, rare fluctuations in simple oscillator systems, including Josephson junctions and micro-mechanical oscillatorsRN937 ; RN939 ; RN940 ; RN936 ; RN938 ; RN941 . Yet, predicting escape from synchrony in complex oscillator networks, subjected to general noise patterns, remains an outstanding problem.

This work takes a first step in this direction by analyzing desynchronization in networks of inertial oscillators driven by broadly distributed Poisson noise. Our approach is the first to connect desynchronization events to general noise statistics and network modes and subgraphs. We explicitly show how higher fluctuation moments control desynchronization rates near bifurcation points, both speeding up and slowing down rates depending on whether noise and network statistics are aligned. Our analytical methods allow us to predict desynchronization rates from power-fluctuation data, and explain how general noise tends to effectively desynchronize only certain network subgraphs.

Consider a model for $N$ coupled phase oscillators with inertia. We assume that the acceleration of the $i$th oscillator’s phase, $\phi_{i}$, is determined by the velocity $v_{i}\equiv\dot{\phi_{i}}$, natural frequencies (input-power) $P_{i}$, and coupling between oscillators $\sum_{j}\!K_{ij}\sin(\phi_{j}-\phi_{i})$. The oscillator dynamics satisfy the second-order differential equation,

[TABLE]

where $M$ and $\gamma$ are inertial and damping constantsFN5 . Eq. (1) gives an approximation to the swing equations describing the transient dynamics of high-voltage electric power grids Nishikawa ; Machowski , which we discuss as an illustration. In what follows we take $K_{ij}=KA_{ij}$, where $A_{ij}$ is a symmetric adjacency matrixTimme2 .

In order to study the effects of non-Gaussian fluctuations on Eq.(1), we take the input power to be

[TABLE]

where the fluctuation $p_{i}(t)$ has over-damped dynamics

[TABLE]

with a damping rate $\alpha$ and stochastic drive $\xi_{i}(t)$Ronellenfitsch ; Timme2 ; Jacquod1 ; HindesKM ; Tyloo . Without loss of generality, we take the average natural frequency to vanish, $\sum_{i}\!\bar{P}_{i}\!=\!0$. For electric power grids this reflects a balance between production and consumption.

In modeling the noisy drive $\xi_{i}(t)$ we are interested in cases where $p_{i}(t)$ exhibits large intermittent fluctuations, as seen in the output from wind and solar sourcesKurths2 ; Milan ; Anvari ; Schafer . The non-Gaussian features of such sources can be captured by tracking the change in $p_{i}(t)$ (or power increments) over fixed time-intervals, $\tau$ Anvari . In practice, $\tau$ may be set by the time resolution of data. For example, in Fig.1(a) we show a histogram of power increments, $g(\tau)\!\equiv\!p(t+\tau)-p(t)$, from 12 wind turbines in northern Germany as inHaehne . The probability distribution, $\text{pr}(g)$, is highly non-Gaussian, e.g., the fourth moment is much larger than for a Gaussian with the same variance.

In order to build a flexible noise model that approximates the measured increment distribution, we assign an independent Poisson pulse for every bin, $b$, in the histogram. Each Poisson pulse has an amplitude, $g_{b}$, which is equal to the bin average, and an occurrence rate $\nu_{b}\!=\!\text{pr}(g_{b})/\tau$. Given this choice, a pulse occurs on average every $\tau$ units of time. More generally, let there be $\mathcal{M}$ arbitrary power increments such that the amplitude for the $b$th increment on the $i$th oscillator is $g_{ib}$, where $b\!\in\!\{1,2,...,\mathcal{M}\}$. By denoting the time at which the $n$th such increment occurs as $t_{ib}[n]$, the stochastic drive $\xi_{i}(t)$ can be represented by a sum of Dirac delta functionsFeynman

[TABLE]

Because the noise $\xi_{i}(t)$ is built from Poisson pulses, each $t_{ib}[n+1]-t_{ib}[n]$ is a stochastic variable with an exponential distribution whose rate is $\nu_{ib}$. For simplicity, we take the time average of the power fluctuations for each oscillator to be zero, i.e., $\sum_{b}g_{ib}\nu_{ib}\!=\!0$. An example time series given our model is shown in Fig.1(b).

II LARGE FLUCTUATION PICTURE OF DESYNCHRONIZATION

When the coupling constant $K$ is sufficiently large a synchronized state is a stable fixed point of Eq.(1). This stable phase-locked state (PLS) depends on network topology and the distribution of $\bar{P}_{i}$ FN5 . In general, the PLS emerges through a saddle-node bifurcation as $K$ is increasedMirollo ; Manik , implying the existence of unstable, saddle phase-locked states. Noise on Eq.(1) can cause networks to fluctuate to these saddles HindesKM ; Timme2 ; Tyloo . Once a saddle is reached the network can desynchronize by either undergoing a large phase slip upon returning to the PLS modulo $2\pi$, or exiting its basin of attraction altogether. An example is shown in Fig.1(b)(lower). Such noise-induced desynchronization whereby fluctuations drive oscillator networks to saddle points, are examples of the general phenomenon of basin-escapeKramers ; Mark1 .

Given the non-Gaussian noise discussed in Sec.I, our strategy is to construct the most-likely (optimal) path of noise and network dynamics that maximizes the probability of reaching a saddle. When desynchronization is rare, the optimal path is describable using analytical mechanics toolsKramers ; Mark1 ; WoillezPRL . Our approach is a generalization of Kramer’s theory for escape, and is valid as long as typical fluctuations are small compared to the distances to saddlesKramers .

We begin our analysis with the network probability distribution, $\rho(\bm{\phi},\mathbb{v},\mathbb{p},t)$. It’s dynamics satisfy a generalized master equation

[TABLE]

where the vector $\mathbb{1}_{i}\!=\!\left<0\;{}_{1},0\;{}_{2},..,1\;{}_{i},..,0\;{}_{N}\right>$ LoraPRL2010 . Note that Eq.(5) is similar to a Fokker-Planck equation, except that, instead of the typical diffusive term, there is a sum over discrete increments to $\mathbb{p}$, given in the last line of Eq.(5), as in a master equation for Poisson processes.

To analyze rare events encoded in the exponential tail of $\rho(\bm{\phi},\mathbb{v},\mathbb{p},t)$, we substitute a WKB ansatz, $\rho(\bm{\phi},\mathbb{v},\mathbb{p},t)\!\cong\!B\exp\{-S(\bm{\phi},\mathbb{v},\mathbb{p},t)\}$, into Eq.(5), assuming $S(\bm{\phi},\mathbb{v},\mathbb{p},t)\gg 1$, and keep the leading order terms in $\partial_{\bm{\phi}}S$, $\partial_{\bm{v}}S$, and $\partial_{\bm{p}}S$ LoraPRL2010 ; Mark2 . This approximation converts Eq.(5) to a Hamilton-Jacobi equation (HJE) for the probability exponent, $S(\bm{\phi},\mathbb{v},\mathbb{p},t)$, called the action, in terms of $\bm{\phi}$, $\mathbb{v}$, $\mathbb{p}$, and their conjugate momenta: $\bm{\lambda}^{\phi}\equiv\partial_{\bm{\phi}}S$, $\bm{\lambda}^{v}\equiv\partial_{\bm{v}}S$, and $\bm{\lambda}^{p}\equiv\partial_{\bm{p}}S$. The network Hamiltonian is

[TABLE]

For the given function $H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})$, the optimal-path dynamics satisfy Hamilton’s equationsFN2 . Once Hamilton’s equations are solved, the action can be calculated as

[TABLE]

In general, solutions are computable numerically subject to boundary conditionsLindleyIAMM . Once the action is known, so is the expected waiting time (or inverse rate) for desynchronizationKramers ; Mark1 ; Mark ; HindesKM

[TABLE]

II.1 Optimal paths near bifurcation

Using our desynchronization mechanics, let us first consider the optimal path (OP) from PLS, $\bm{\phi}^{*}$, to saddles, $\bm{\phi}^{s}$, near the saddle-node bifurcation (SN). We denote the critical coupling $K_{\text{SN}}$, where $K\!=\!K_{\text{SN}}[1+\kappa]$. When $\kappa\!\ll\!1$, the dynamics slows onto a one-dimensional manifold with universal properties. As we will show, statistical moments of the noise first contribute to the action at specific powers of $\kappa$, from which we can calculate their effects on desynchronization rates. In order to simplify the analysis, we assume that the noise for all nodes is independent and identically distributed. Hence, we drop the subscript $i$ in $\nu$, $g$ and $\mu$ from now on.

First, we construct the lowest-order solution and thereby demonstrate our natural expansion in $\kappa^{1/2}$. Higher order terms and further calculation details are given in Apps.C-D. From the fixed-point boundary conditions we expand $\bm{\phi}^{*}$ and $\bm{\phi}^{s}$ around the SN value, $\bm{\phi}^{SN}$, in powers of $\kappa$. At the saddle-node bifurcation, the Fiedler modeFiedler of the network Laplacian, $L_{ij}(\bm{\phi}^{*})\!=A_{ij}\cos(\phi_{j}^{*}-\phi_{i}^{*})-\delta_{ij}\sum_{k}A_{ik}\cos(\phi_{k}^{*}-\phi_{i}^{*})$, has zero eigenvalue. In general, the Fiedler mode is the slowest mode of $L_{ij}(\bm{\phi}^{*})$, and we denote its components $r_{i}$. The normalized Fiedler mode is defined up to an overall sign, and we choose the convention $r_{i}\!>\!0$ if $\phi_{i}^{s}\!-\!\phi_{i}^{*}\!>\!0$, so that $\mathbb{r}$ points from $\bm{\phi}^{*}$ to $\bm{\phi}^{s}$. Quite remarkably, the fixed points close to bifurcation can be expressed in terms of $r_{i}$, $\phi_{i}^{*}\!=\!\phi_{i}^{SN}\!\!-\!C\kappa^{1/2}r_{i}$ and $\phi_{i}^{s}\!=\!\phi_{i}^{SN}\!\!+\!C\kappa^{1/2}r_{i}$, where:

[TABLE]

Because the sums in Eq.(9) appear when considering the noise moments, we write $C\!\equiv\!\sqrt{2R_{0}/R_{2}}$. Next, we write the phases emerging from the SN in terms of a coordinate $x(t)$, $\phi_{i}(t)\!=\!\phi_{i}^{SN}\!+\!C\kappa^{1/2}r_{i}x(t)$, where $x\!\in\![-1,1]$, which is valid near bifurcationKuznetsov1 . Substituting this form into Hamilton’s equations, and collecting terms at order $\kappa$, we find:

[TABLE]

where $\mu_{2}\!\equiv\!\sum_{b}\nu_{b}g_{b}^{2}$ is the noise variance. From Eqs.(7) and (II.1), we get

[TABLE]

The structure of Eq.(11) is interesting. The action at lowest order in $\kappa$ is proportional to the damping rate for power-fluctuations squared, implying that doubling the rate, increases the time-scale for desynchronization as $\left<T\right>\!\rightarrow\!\left<T\right>^{4}$. What’s more, each node’s contribution is proportional to $\mu_{2}^{\;-1}$, and hence noise distributions with the same variance produce the same rate of desynchronization at lowest order. Therefore, the effects of higher fluctuation moments must appear at (and be suppressed by) higher powers of $\kappa$.

Continuing the OP construction at higher-powers in $\kappa$, we notice that the lowest-order contribution to the action from the $n$th moment of the noise, $\mu_{n}\!\equiv\!\sum_{b}\nu_{b}g_{b}^{n}$, is determined by energy-conservation for large fluctuations, $H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})\!=\!0$, at $\mathcal{O}(\kappa^{n})$ FN2 . Interestingly, such contributions, $\Delta^{(n)}S$, only depend on Eq.(II.1) and can be calculated (see App.D for derivation):

[TABLE]

Equation (12) is very useful for comparing non-Gaussian and Gaussian white noise (GWN) effects. In the latter, the stochastic drive in Eq.(3) is replaced by a Gaussian process with time-correlation $\left<\xi_{i}(t)\xi_{j}(t^{\prime})\right>\!=\!\mu_{2}\delta_{ij}\delta(t-t^{\prime})$ $\forall\{i,j\}$, and Eq.(5) becomes a Fokker-Planck equation. Simple GWN is considered in most works on stochastic oscillatorsTimme2 ; Tyloo ; ISI:000383052300003 . Because the action for GWN is equivalent to keeping only the variance in the Hamiltonian, Eq.(12) represents the correction to the GWN action from fluctuations with a first non-zero moment $\mu_{n}$ ($n\!>\!2$).

Interestingly, we see that desynchronization rates exhibit a critical behavior characterized by a spectrum of exponents $n-\frac{1}{2}$ depending on the $n$th moment of the noise distribution. Most importantly, the sign of the first non-Gaussian correction in Eq.(12) is given by the product of the $n$th noise moment times the $n$th moment of the Fiedler components. In the particular case when $n\!=\!3$, the sign of $\Delta S$ is determined by the product of the skewness of the Fiedler components with that of the noise. When the two skewnesses are aligned, $\Delta S^{(3)}\!<\!0$, and desynchronization occurs at an exponentially faster rate. When they have opposite signs, $\Delta S^{(3)}\!>\!0$, and desynchronization occurs at an exponentially slower rate. Generally speaking, the Fiedler modes in high voltage power grids mostly reside on peripheral nodes - but not all peripheral nodesPhilippePlosOne . Hence, such networks are highly skewed, by our definition, and will show exponential sensitivity in desynchronization rates to noise skewness.

Two examples are shown in Fig.2, where we compare desynchronization times for GWN (red), positively skewed $\mu_{3}\!>\!0$ (green), and wind turbine (blue) fluctuations. All power fluctuations have the same variance. The skewed noise distribution is a simple two-pulse model. Solid lines denote Eq.(7) computations, while points indicate Monte-Carlo simulations. The dashed lines in green and blue denote the action for GWN plus Eq.(12), which is in good agreement (in most cases hard to distinguish). The networks are drawn in each subpanel. Nodes are blue if $r_{i}\!>\!0$ and red if $r_{i}\!<\!0$; sizes are proportional to $|r_{i}|$. In Fig.2(a), the network has a negatively skewed Fiedler mode, $\sum_{i}\!{r_{i}}^{3}\!<\!0$, and hence desynchronization occurs at an exponentially slower rate, just as predicted. In contrast in Fig.2(b) the network has a positively skewed Fiedler mode, and therefore desynchronization occurs at an exponentially faster rate. Note that for both networks, blue lines are always under the red, meaning symmetric noise ($\mu_{3}\!\approx\!0$) produces an increase in desynchronization rates.

III SYNCHRONIZED SUBGRAPH APPROXIMATION

Since we have shown that rare desynchronization occurs along the Fiedler mode near bifurcation, we expect that nodes that are topologically nearby in the network and have similar Fiedler-mode values, do not tend to desynchronize during a large fluctuation – as a first approximation for general coupling strengths. Partitioning the network according to the Fiedler modeFiedler at bifurcation is a useful approximation for finding much lower-dimensional desynchronization pathways, which we call the synchronized subgraph approximation (SSA). An algorithm for constructing Hamilton’s equations for a SSA is given in the App.E.

First, the SSA can be made exact for networks where the SN corresponds to a single (overloaded) edge with a phase-difference $\pi/2$. This condition is satisfied for tree topologies and frequently satisfied for sparse networksManik . For such single-cut saddle-nodes (SCSN), we can always construct a SSA where the network splits into exactly two subgraphs at bifurcation and the nodes within each remain synchronized on average. If we denote the two subgraphs divided by the overloaded edge $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$, then we can reduce Hamilton’s equations to a single noisy oscillator system in relative phase-space coordinates:

[TABLE]

In Eqs.(III) $\Psi$ is the phase difference between the two subgraphs, $V$ is the difference in their velocities, $P$ is the difference in the power fluctuations, and the $\Lambda$’s are conjugate momenta. The system is closed since $K_{\text{SN}}\!=\!|\sum_{i\in\mathcal{S}_{1}}\bar{P}_{i}|$ for SCSN Dekker . Remarkably, we can see that the dynamics is parameterized only by $\kappa$ and the sizes of the two subgraphs, $|\mathcal{S}_{1}|$ and $|\mathcal{S}_{2}|$, but is otherwise independent of topology. The action corresponds to a single stochastic oscillator:

[TABLE]

Comparisons between desynchronization times and SSA predictions are shown in Fig.3(a) for three networks with a SCSN: red, blue, and green. In each case $N\!\sim\!\mathcal{O}(100)$. In particular, the UK power grid is shown in green, and its two subgraphs are drawn in the top of Fig.3(b); the red and blue examples in Fig.3(a) correspond to two randomly generated 100-node trees with power-law degree distributions, and an exponent equal to three. Theory computations are shown with solid lines (Eq.14) and Monte-Carlo simulations with points. Examples show quantitative agreement despite two-orders of magnitude reduction in the effective dynamics. Moreover, in App.F we show that Eqs. (12) and (14) are accurate even in networks with significant variation in $M$, $\gamma$, and $\alpha$.

Similarly, the SSA can be applied to networks without a SCSN, but where there is an exact symmetry in the Fiedler mode at bifurcation. An example is shown in Fig.3(a) in cyan for a block network. In this case the network splits into exactly four subgraphs at bifurcation. The Fiedler-mode has only four unique component values– one for each subgraph. The subgraphs are drawn in the lower panel of Fig.3(b). As a consequence, the SSA consists of four effectively coupled oscillators.

Lastly, the SSA can be applied to networks whose Fiedler modes do not have exact symmetry at bifurcation, by partitioning the network into subgraphs with approximately uniform Fiedler-mode values. An example is shown in Fig.3(a) (magenta), again for the UK power grid, but with a different distribution of $\bar{P}_{i}$. In this case, there is no symmetry at bifurcation and each node has a unique Fiedler component value. The solid line represents a SSA assuming 20 subgraphs.

IV CONCLUSION

In this work, we analyzed desynchronization in complex oscillator networks by non-Gaussian noise. It is often thought that broad-tailed fluctuations perturb dynamical systems more strongly. However, we found that such noise does not always increase the rate of network desynchronization; the latter depends on whether higher-moments of fluctuations and a network’s slowest mode have the same sign, and is therefore a topological as well as a noise effect. Our approach was based on arbitrarily distributed Poisson fluctuations, which we fit to power-increment data from renewable energy sources. In addition, we developed a reduction technique for predicting desynchronization, based on the Fiedler mode at bifurcation, and the observation that noise tends to effectively desynchronize only certain network subgraphs. Such a reduction should be valuable for studying rare processes in high-dimensional oscillator networks more broadly, where predicting rare events is both analytically and computationally difficult.

Lastly, our approach revealed a spectrum of scaling exponents that determine at what powers in the coupling the $n$th moment of noise contributes to desynchronization rates. Our results are general for escape through a saddle. However, our methods can be further generalized to rare events induced by non-Gaussian noise in other dynamical processes in networks including: extinctionHindesPRL2016 ; HindesPRL2019 , switchingHavlin2017 , and more general oscillator transitionsLoos ; HindesChaos .

JH was supported through the U.S Naval Research Laboratory Karle Fellowship. PJ was supported by the Swiss National Science Foundation grants (200020_182050) and (PYAPP2_154275). IBS was supported by the U.S. Naval Research Laboratory funding (N0001419WX00055) and the Office of Naval Research (N0001419WX01166) and (N0001419WX01322).

V. APPENDIX

Appendix A Computing optimal paths

Optimal paths for rare desynchronization satisfy Hamilton’s equations:

[TABLE]

Numerical solutions of Eqs.(15-20) were found using the Iterative-Action-Minimization-Method (B. S. Lindley and I. B. Schwartz, Physica D 255, 22 (2013)) with fixed-point boundary conditions. Example Matlab code can be found in the supplementary material of (J. Hindes and I. B. Schwartz, EPL 120, 56004 (2017)), for instance, and is available upon request. The method requires a trial solution. For small $\kappa$ we used Eqs.(II.1) from the main text, and then bootstrapped to other regions of parameter space.

Appendix B Fixed points

The coupling at which the saddle-node occurs, $K_{\text{SN}}$, can be computed numerically by solving the following $N+1$ equations with a quasi-Newton method:

[TABLE]

where $z$ is the second smallest (in magnitude) eigenvalue of the Laplacian matrix $L_{ij}(\bm{\phi}^{SN})\!=A_{ij}\cos(\phi_{j}^{SN}-\phi_{i}^{SN})-\delta_{ij}\sum_{k}A_{ik}\cos(\phi_{j}^{SN}-\phi_{i}^{SN})$.

In order to calculate optimal paths as a function of the distance to bifurcation $\kappa$, where $K\!=\!K_{\text{SN}}[1+\kappa]$, we first calculate the fixed points of Eq.(1) in powers of $\kappa$. Let us substitute $\phi_{i}^{*}\!=\!\phi_{i}^{SN}\!+\!\kappa^{1/2}q_{i}\!+\!\kappa w_{i}\!+\!...$ into Eq.(1), given $\dot{v}_{i}\!=\!0$ and $v_{i}\!=\!0$ $\forall i$. Note, the sub(super)-script SN implies evaluation at the saddle-node. Our goal is to find $q_{i}$. At $\mathcal{O}(\kappa^{1/2})$ we find the equation

[TABLE]

which expresses the saddle-node condition that the network Laplacian has a Fiedler mode with eigenvalue zero. Hence, we may write $q_{i}\!=\!-Cr_{i}$, where $r_{i}$ is the Fiedler mode at bifurcation and $C$ is a constant that we wish to determine. Continuing to $\mathcal{O}(\kappa)$, gives

[TABLE]

If we take the product of Eq.(B) with $r_{i}$, sum over $i$, and solve for $C$, we get Eq.(9) from the main text. This is easy to see by expanding $w_{i}$ in the eigenmodes of the Laplacian. In Eq.(B), the component parallel to $r_{i}$ vanishes from the saddle-node condition, while all other components vanish due to orthonormality. Note: there are two possible solutions for $C$: a stable phase-locked state (positive), and a saddle (negative). Given the definitions in the main text, the Fiedler value can be shown to be

[TABLE]

by an analogous expansion.

B.1 Single-cut saddle node

The procedure outlined above can be carried out indefinitely. Here, we restrict ourselves to a particular class of saddle-node bifurcations, which we call single-cut saddle-nodes (SCSN). In this special case, a single edge becomes overloaded at bifurcation (with a $\pi/2$ phase-difference between the nodes), and is a cut edge of the network specified by the adjacency matrix, $A$. The edge partitions the network into exactly two subgraphs (denoted $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$), whose nodes only share one edge (the cut edge) in common. For example, tree networks always have SCSN.

Networks with SCSN have useful properties. For instance, $K_{\text{SN}}$ can be determined by summing Eq.(1) over all nodes in $\mathcal{S}_{1}$ (or $\mathcal{S}_{2}$), given $\dot{v}_{i}\!=\!0$ and $v_{i}\!=\!0$ $\forall i$:

[TABLE]

Only one term survives in the second sum in Eq.(26): corresponding to the cut edge connecting $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$– because sine is an odd function. Moreover the phase difference for the cut edge is $\phi_{j}^{SN}\!-\!\phi_{i}^{SN}\!=\!-\pi/2$. This latter property can be shown from the saddle-node condition and the fact that for SCSN $r_{i}\!=\!\sqrt{|\mathcal{S}_{2}|/[N|\mathcal{S}_{1}|]}$ if $i\in\mathcal{S}_{1}$, and $r_{i}\!=\!-\sqrt{|\mathcal{S}_{1}|/[N|\mathcal{S}_{2}|]}$ if $i\in\mathcal{S}_{2}$. Note: the specifications for $\mathcal{S}_{1}$ and $\mathcal{S}_{2}$ are given by the convention that the node in $\mathcal{S}_{1}$ connected to $\mathcal{S}_{2}$ along the cut edge has a larger phase than its counterpart in $\mathcal{S}_{2}$. Therefore,

[TABLE]

as first noted for trees in Dekker .

Using the stated SCSN properties, it is straightforward to show that

[TABLE]

by expanding the fixed-point conditions to higher orders in $\kappa$. In Fig.4(a) we plot the mean (blue) and standard deviation (red) of the error-vector, $\text{Error}_{i}\!=\![\phi_{i}^{s}-\phi_{i}^{*}-2Cr_{i}\kappa^{1/2}(1-\frac{5}{12}\kappa+\frac{43}{160}\kappa^{2})]/r_{i}$ for an example SCSN shown in Fig.4(b). We can see that the fixed-point expression Eq.(28) is accurate to $\mathcal{O}(\kappa^{7/2})$ and is parallel to $r_{i}$. The fixed points were computed numerically with the choice of zero-average-phase $\left<\phi_{i}^{s}\right>\!=\!\left<\phi_{i}^{*}\right>\!=\!0$.

Appendix C Near bifurcation paths

Using the fixed-point boundaries expressed in powers of $\kappa$, we can construct the optimal paths for sufficiently small $\kappa$. To make the calculation simpler, at this point we restrict ourselves to $iid$ power fluctuations, and drop the node-subscript $i$ in the noise moments, $\mu$. As mentioned in the main text, at lowest order in $\kappa$, we substitute the ansatz $\phi_{i}(t)\!=\!\phi_{i}^{SN}\!+\!C\kappa^{1/2}r_{i}x(t)+...$ into Eqs.(15-20). For the other phase-space variables, we have the general expansions:

[TABLE]

where $m\!=\!2,3,...$ When $m\!=\!2$, we assume that all variables are parallel to the Fiedler mode according to $v_{i,2}\!=\!a_{v}(x)r_{i}$, $p_{i,2}\!=\!a_{p}(x)r_{i}$, $\lambda_{i,2}^{\phi}\!=\!a_{\lambda}^{\phi}(x)r_{i}/\mu_{2}$ , $\lambda_{i,2}^{v}\!=\!a_{\lambda}^{v}(x)r_{i}/\mu_{2}$ , and $\lambda_{i,2}^{p}\!=\!a_{\lambda}^{p}(x)r_{i}/\mu_{2}$.

We note that since $\dot{x}\!=\!a_{v}(x)\kappa^{\frac{1}{2}}/C$, time-derivatives for phase-space variables other than $\phi$ are $\mathcal{O}(\kappa^{\frac{3}{2}})$, and so the LHS of Eqs.(16-20) can be ignored at $\mathcal{O}(\kappa)$. By substituting the expansion into Eqs.(16-17) and Eqs.(19-20), we find the following relations for $\mathcal{O}(\kappa)$:

[TABLE]

One more equation is needed to close the system. We use the zero-energy condition $H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})\!=\!0$, at $\mathcal{O}(\kappa^{2})$, which results in

[TABLE]

Solving for the non-zero solution gives Eqs.(II.1) from the main text. Note: we have assumed that the product, $\sum_{ij}A_{ij}\sin(\phi_{j}^{SN}\!\!-\!\phi_{i}^{SN})r_{i}$, is negative in writing Eqs.(II.1), given our sign convention. This property appears to be general for the SN. Hence, $R_{0}\!=-\!\sum_{ij}A_{ij}\sin(\phi_{j}^{SN}\!\!-\!\phi_{i}^{SN})r_{i}\!=\!|\sum_{ij}A_{ij}\sin(\phi_{j}^{SN}\!\!-\!\phi_{i}^{SN})r_{i}|$.

C.1 Single-cut saddle node

As with the fixed-point boundaries in Sec.B, we restrict ourselves to SCSN for higher-order OPs, since the results are comparatively simple. First, we start with the ansatz $\phi_{i}(x)\!=\!\phi_{i}^{SN}\!+Cx(t)r_{i}\kappa^{1/2}[1-\frac{5}{12}\kappa]$, from Eq.(28). The next order in $\kappa$ requires substituting this ansatz, as well as the general expansion Eqs.(C) into Eqs.(16-17) and Eqs.(19-20) and collecting terms of $\mathcal{O}(\kappa^{3/2})$. As noted in the previous section, the relevant time derivatives are at lowest order $\mathcal{O}(\kappa^{3/2})$:

[TABLE]

When $m\!=\!3$, again, all variables are parallel to the Fiedler mode according to $v_{i,3}\!=\!b_{v}(x)r_{i}$, $p_{i,3}\!=\!b_{p}(x)r_{i}$, $\lambda_{i,3}^{\phi}\!=\!b_{\lambda}^{\phi}(x)r_{i}/\mu_{2}$ , $\lambda_{i,3}^{v}\!=\!b_{\lambda}^{v}(x)r_{i}/\mu_{2}$ , and $\lambda_{i,3}^{p}\!=\!b_{\lambda}^{p}(x)r_{i}/\mu_{2}$. Substituting into Eqs.(C.1), Eqs.(16-17) and Eqs.(19-20), we find:

[TABLE]

Equation (36) requires elaboration. For SCSN, the coupling term is

[TABLE]

Note: all cosine terms vanish, since if $i$ and $j$ are in the same subgraph, then $r_{i}\!=\!r_{j}$, and if they are not, then $|\phi_{j}^{SN}-\phi_{i}^{SN}|=\frac{\pi}{2}$. This means, for instance, that $\bar{P}_{i}+K_{\text{SN}}[1+\kappa]\sum_{j}A_{ij}\sin(\phi_{j}-\phi_{i})$ only contains integer powers of $\kappa$.

To close the system, we again use the zero-energy condition at $\mathcal{O}(\kappa^{5/2})$. The resulting equation is:

[TABLE]

Substituting Eqs.(30-34) into Eqs.(36-39) and Eq.(C.1)

[TABLE]

In order to calculate the action at next order, $\mathcal{O}(\kappa^{5/2})$, we need to know $\lambda_{i,4}^{\phi}$. This can be achieved using Eqs.(42-46) and the zero-energy condition at $\mathcal{O}(\kappa^{3})$. If we assume $\lambda_{i,4}^{\phi}\sim r_{i}$, then the calculation only requires a fair bit of algebra. The result is:

[TABLE]

Note, that all constants in Eq.(C.1) can be given a closed-form expression for SCSN, and depend only on $|\mathcal{S}_{1}|$ and $|\mathcal{S}_{2}|$. In Fig.5 we show the accuracy of our expansion for the UK grid with a SCSN, Fig.3(b)(top). In Fig.5 panel (a), solutions of Eqs.(III) are shown in blue for $\kappa\!=\!0.001$, and the $\mathcal{O}(\kappa)$ solution, Eqs.(II.1), is shown in red. In panel (b), we plot the difference between the two curves in panel (a) in blue, and compare to the $\mathcal{O}(\kappa^{2})$ solution Eq.(C.1) shown in red. The agreement is excellent at both orders of $\kappa$.

Integrating our expressions for the phase-space variables to $\mathcal{O}(\kappa^{2})$ in $\lambda_{i}^{\phi}$ and $\mathcal{O}(\kappa^{\frac{3}{2}})$ for all others, according to Eq.(7), gives the action for SCSN at $\mathcal{O}(\kappa^{\frac{5}{2}})$.

Appendix D Non-Gaussian effects

We would like to know how the higher-moments of the noise distribution contribute to the action in general. As in the main text, we restrict ourselves to $iid$ power fluctuations. A clue comes from our solution for the OP at lowest order in $\kappa$. Recall from Sec.C, that $\lambda_{i,2}^{\phi}$ is determined by the zero-energy condition at $\mathcal{O}(\kappa^{2})$ – where the noise variance $\mu_{2}$ first appears in $H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})$. The general pattern is the following: the contribution proportional to $\mu_{n}$ in $\lambda_{i,2(n-1)}^{\phi}$ is determined by the zero-energy condition at $\mathcal{O}(\kappa^{n})$. First, we can check that this pattern holds for $\mu_{3}$, by looking at $H(\bm{\phi},\mathbb{v},\mathbb{p},\bm{\lambda}^{\phi},\bm{\lambda}^{v},\bm{\lambda}^{p})\!=\!0$ at $\mathcal{O}(\kappa^{3})$. Using our general expansion, we get an equation:

[TABLE]

Note, we have introduced the notation $T_{i}\!\equiv\!-\gamma v_{i}+p_{i}+\bar{P}_{i}+K\sum_{j}A_{ij}\sin(\phi_{j}-\phi_{i})\!$ and $T_{i}\!=\!\sum_{m}T_{i,m}\kappa^{\frac{m}{2}}$.

It is important to realize a few properties of the phase-space coordinates. First, only the conjugate momenta, (the $\lambda$’s), depend explicitly on the noise statistics. Second, since $\mathcal{O}(\kappa^{n})$ is the lowest order at which $\mu_{n}$ enters the Hamiltonian, $\lambda_{i,2(n-1)-j}^{\phi}$, $\lambda_{i,2(n-1)-j}^{v}$, and $\lambda_{i,2(n-1)-j}^{p}$ $\forall j\in\{1,2,...,2(n-2)\}$ depend on $\mu_{n-1}$, $\mu_{n-2}$,… $\mu_{2}$, but not on $\mu_{n}$. Therefore, if $\lambda_{i,4}^{v}$ and $\lambda_{i,4}^{p}$ vanish from Eq.(D), then only two terms depend on $\mu_{3}$, i.e., $\lambda_{i,4}^{\phi}v_{i,2}$ and $\frac{\mu_{3}}{3!}[\lambda_{i,2}^{p}]^{3}$, and we can solve explicitly for the contribution to $\lambda_{i,4}^{\phi}$ from $\mu_{3}$.

Luckily, the terms proportional to $\lambda_{i,4}^{v}$ and $\lambda_{i,4}^{p}$ in Eq.(D) depend on the lowest order solution, which we have explicit expressions for. In fact, $T_{i,2}\!=\!0$ and $-\alpha p_{i,2}+\mu_{2}\lambda_{i,2}^{p}\!=\!0$. Finally, we add an additional assumption that $\lambda_{i,4}^{\phi}\sim r_{i}$, as with the SCSN, Eq.(C.1). With this assumption, the lowest-order (in $\kappa$) contribution to $\lambda_{i}^{\phi}$ from $\mu_{3}$, denoted $\Delta^{(3)}\lambda_{i}^{\phi}$ is found from Eq.(D):

[TABLE]

Exactly the same argument can be used to calculate $\Delta^{(4)}\lambda_{i}^{\phi}$, etc. The general expression is

[TABLE]

Finally, the lowest-order contribution to the action from $\mu_{n}$ is found from Eq.(7) – namely, integrating $\Delta^{(n)}\lambda_{i}^{\phi}(x)$ over $x$. The result is Eq.(12) in the main text. Note: the integral contributions from $\lambda_{i}^{v}$ and $\lambda_{i}^{p}$ can be ignored, since $v_{i}$ and $p_{i}$ are at lowest order $\mathcal{O}(\kappa)$, while $\phi_{i}^{s}-\phi_{i}^{*}$ is $\mathcal{O}(\kappa^{1/2})$.

Furthermore, Eq. (12) can be used to compare the actions of two increment distributions that differ in the $n$th moment. For instance, the difference in the action between a symmetric increment distribution ($\mu_{3}\!=\!0$) and a non-symmetric distribution ($\mu_{3}\!\neq\!0$), each with the same variance $\mu_{2}$, is given by Eq.(12) with $n=3$. Moreover, as noted in the main text, if the power fluctuations are assumed to be Gaussian white noise with a variance $\mu_{2}$, then Eq. (12) gives the exponential correction to the rate of desynchronization for non-symmetric $(n=3)$, and symmetric $(n=4)$ increment distributions. Our expansion gives very accurate results for desynchronization rates near bifurcation, as demonstrated in Fig.2.

For completeness, we note that for Gaussian white noise, Hamilton’s equations are identical, except for Eq.(17):

[TABLE]

and the Hamiltonian Eq.(6):

[TABLE]

Appendix E Synchronized subgraph approximation

As mentioned in the main text, in many cases the SN bifurcation occurs with exact symmetry in the Fiedler mode; namely, the network is partitioned into $\mathcal{N}$ subgraphs, $\mathcal{S}_{1},\mathcal{S}_{2},...,\mathcal{S}_{\mathcal{N}}$, where two nodes $i$ and $j$ in the same subgraph $\mathcal{S}_{n}$ have $r_{i}\!=\!r_{j}=r^{(n)}$ at bifurcation. As already noted, SCSN cases have $\mathcal{N}=2$, including trees. Larger values of $\mathcal{N}$ occur for block-networks (or clique trees), such as Fig.3(b)–the lower panel.

Since, the Fiedler mode is the weakest stable mode of the network and the optimal desynchronization path is parallel to $r_{i}$ over several orders in $\kappa$, as demonstrated in previous sections, we simply assume (as an approximation) that all nodes within the same subgraph, e.g., $n$, are synchronized: $\phi_{i}-\phi_{i}^{SN}=\Phi_{n}$, $v_{i}=\mathcal{V}_{n}$, $p_{i}=\mathcal{P}_{n}$, $\lambda_{i}^{\phi}=l_{n}^{\phi}$, $\lambda_{i}^{v}=l_{n}^{v}$, and $\lambda_{i}^{p}=l_{n}^{p}$ $\forall i\in\mathcal{S}_{n}$. We can find an approximate set of Hamilton’s equations for such OPs by simply averaging over all nodes within a subgraph:

[TABLE]

Performing these averages over the original Eqs.(15-20) results in network coupling terms

[TABLE]

Let us specify an index function which maps the node number $i$ to its subgraph number $n$, i.e., $F(i)\!=\!n$. Using this notation, the coupling terms become:

[TABLE]

[TABLE]

Now, since all phase-space variables are assumed to be synchronized within subgraphs, we can define the following coupling matrices between subgraphs $n$ and $n^{\prime}$

[TABLE]

Therefore, the synchronized subgraph equations become

[TABLE]

In general, the matrices Eqs.(58-59) have to be computed numerically at the saddle-node bifurcation in order to solve Eqs.(60-65).

However, for SCSN the results simplify significantly. In this case $n\!=\!1,2$, and $C_{2,1}=C_{1,2}=0$ and $S_{2,1}=-S_{1,2}=1$. In addition, we notice that the effective two-oscillator system for SCSN can be reduced to one by introducing the relative coordinates: $\Psi\!=\!\Phi_{1}\!-\!\Phi_{2}$, $V\!=\!\dot{\Psi}$, $P\!=\!\mathcal{P}_{1}\!-\!\mathcal{P}_{2}$, $\Lambda_{\Psi}\!=\!|\mathcal{S}_{1}|l_{1}^{\phi}\!=\!-|\mathcal{S}_{2}|l_{2}^{\phi}$, $\Lambda_{V}\!=\!|\mathcal{S}_{1}|l_{1}^{v}\!=\!-|\mathcal{S}_{2}|l_{2}^{v}$, and $\Lambda_{P}\!=\!|\mathcal{S}_{1}|l_{1}^{p}\!=\!-|\mathcal{S}_{2}|l_{2}^{p}$. The result is Eqs.(III) in the main text.

In addition, Eq.(60-65) can be used as a coarse-grained approximation for networks without exact symmetry in the Fiedler mode at bifurcation. One simply groups together nodes in a subgraph if they have similar Fiedler-mode values, where similar means within some tolerance $\epsilon$, e.g., \big{|}\!\sum_{l\in\mathcal{S}_{n}}\!r_{l}/|\mathcal{S}_{n}|\!-\!r_{i}\big{|}\!<\!\epsilon. An example is shown in Fig.6, for the UK power-grid example from Fig.3 without symmetry at bifurcation (magenta). Plotted is the Fiedler mode at bifurcation versus the node number, $i$. Blue, magenta, green, and cyan nodes were placed in the same subgraphs, since node Fiedler-mode values were within $5\%$ of one another (for similarly colored nodes). Each of the red nodes was treated as a subgraph (of size one). The total number of subgraphs was twenty given this partition. The coupling matrices Eqs.(58-59) were computed at bifurcation, and the OPs were solved from Eqs.(60-65) with twenty subgraphs.

Appendix F Parameter heterogeneity

In this section we show that our expansion in $\kappa$ and the SSA are robust to parameter variation. Instead of being homogeneous, we let $M$, $\gamma$, and $\alpha$ be drawn from uniform distributions, independently for each node. Given a heterogeneity scale parameter $\eta$, we let $M_{i}\!=\!M_{0}[1.0+\eta(\mathcal{R}_{i,M}-0.5)]$, $\gamma_{i}\!=\!\gamma_{0}[1.0+\eta(\mathcal{R}_{i,\gamma}-0.5)]$, and $\alpha_{i}\!=\!\alpha_{0}[1.0+\eta(\mathcal{R}_{i,\alpha}-0.5)]$, where $\mathcal{R}_{i,M}$, $\mathcal{R}_{i,\gamma}$, and $\mathcal{R}_{i,\alpha}$ are independent, random numbers drawn uniformly over $[0,1]$. In this section , we take $M_{0}\!=\!0.02546s^{2}$, $\gamma_{0}\!=\!0.10053s$, and $\alpha_{0}\!=\!1s^{-1}$.

Figure 8 compares predictions of Eq.(12) with the computed difference in the actions between Gaussian and non-Gaussian noise. The computed actions were found by solving Eqs.(15-20). Note: in Eqs.(15-20) $M\!\rightarrow\!M_{i}$, $\alpha\!\rightarrow\!\alpha_{i}$ and $\gamma\!\rightarrow\!\gamma_{i}$. Results are shown for several levels of the heterogeneity parameter $\eta$ and are labeled by the average coefficient of variation, i.e, \text{cv}\!\equiv\!\big{[}\frac{\sigma_{M}}{\left<M\right>}+\frac{\sigma{\gamma}}{\left<\gamma\right>}+\frac{\sigma_{\alpha}}{\left<\alpha\right>}\big{]}/3, where $\sigma_{q}$ and $\left<q\right>$ denote the standard deviation and average for parameter $q$. The underlying network corresponds to Fig.2(b). We see good agreement between numerics and Eq.(12) despite significant heterogeneity, particularly for skewed noise distributions.

Similarly, we perform computations using the SSA for single-cut-saddle-nodes with homogeneous parameters, Eqs.(III-14), and compare with heterogeneous parameters. For example, Figure 8 shows such a compare for the network drawn in Fig.4(b). The computed actions are shown in blue and the SSA results are shown with a black-line for skewed noise (see Sec.G). Again, we see good agreement despite the heterogeneity.

Appendix G Simulations

Monte-Carlo simulations were performed using a combination of Gillespie’s algorithm for power fluctuations, $p_{i}$, and Euler’s method for integrating Eqs.(1-3) between reaction times. For illustration, let us assume that the noise is independent and identically distributed (iid) for each node. The stochastic rate for the next reaction is $\mathcal{R}\!=\!\sum_{i,m}\nu_{im}=N\sum_{m}\nu_{m}$. The next reaction time is stochastically selected $\Delta T\!=-\!\ln(r_{1})/\mathcal{R},$ where $r_{1}$ is a uniformly distributed random number over the unit interval. Since the noise is $iid$, the node which receives the increment is selected uniformly at random from the $N$ nodes. Another random number is generated, $\mathcal{P}\!=\!r_{2}\sum_{m}\nu_{m}$, where $r_{2}$ is a uniformly distributed random number over the unit interval. The $n$th increment is chosen if $\sum_{m=1}^{n-1}\nu_{m}<\mathcal{P}<\sum_{m=1}^{n}\nu_{m}$. Equations (1-3) are integrated with Euler’s method from $t$ to $t+\Delta T$ with time steps $dt\!=\!4*10^{-5}$, at which time the selected node, e.g,. $i$, has its power incremented: $p_{i}\!\rightarrow\!p_{i}\!+\!g_{n}$. On the other hand, for Gaussian noise the Euler-Milstein method was used with $dt\!=\!4*10^{-5}$.

The wind-turbine data was taken from Haehne . The data consists of power measurements at 1-second resolution for 12 turbines, each rated at 2MW. We averaged over the available data at each time step and histogrammed the power increments $p(t+1)\!-\!p(t)$ using $60$ uniformly spaced bins ($\mathcal{M}\!=\!60$). The result is shown Fig.1(a). In order to make desynchronization less rare, we artificially increased the occurrence rate for the Poisson pulses, such that $\sum_{m}\nu_{m}\!=\!10s^{-1}$ (not $1s^{-1}$). For reference, the noise variance for $iid$ wind-turbine noise given our model is $\mu_{2}\!=\!5.8317*10^{-4}s^{-1}$. On the other hand for skewed noise comparisons, we chose a simple two-pulse model: $\nu_{m}\!\in\!\{10,3.33333\}$ and $g_{m}\!\in\!\{-0.0038183,0.0114549\}$, which has the same variance as the wind-turbine model.

Lastly, each Monte-Carlo point in Figs.2-3 represent the log of the average of 200 slip times – defined as the time it takes to see a phase difference greater than $2\pi$ develop between any connected oscillators. Each of the slips were generated from the phase-locked state initial conditions and with different random number seeds. All simulations in Fig.3 were done with skewed noise except for the block network (cyan), for which we used the wind-turbine noise.