Noise-Induced Desynchronization and Stochastic Escape from Equilibrium in Complex Networks

TL;DR

This paper investigates how noise can cause complex physical systems, modeled by Kuramoto-like equations, to escape from stable states, revealing that inertia can accelerate this process under certain noise conditions.

Contribution

It derives conditions under which noise induces stochastic escape in complex networks and uncovers the counterintuitive role of inertia in this phenomenon.

Findings

Inertia can lead to faster escape from stable states under noise.

Noise strength influences the escape dynamics significantly.

Counterintuitive effects of inertia on stochastic escape are identified.

Abstract

Complex physical systems are unavoidably subjected to external environments not accounted for in the set of differential equations that models them. The resulting perturbations are standardly represented by noise terms. We derive conditions under which such noise terms perturb the dynamics strongly enough that they lead to stochastic escape from the initial basin of attraction of an initial stable equilibrium state of the unperturbed system. Focusing on Kuramoto-like models we find in particular that, quite counterintuitively, systems with inertia leave their initial basin faster than or at the same time as systems without inertia, except for strong white-noise perturbations.

Tables

Table 1. Table 1: Final winding number q ( 1 ) superscript 𝑞 1 q^{(1)} and number of iterations before the escape for m = 0 𝑚 0 m=0 (simulations 1-3) and finite inertia (simulations 4-6). Each triplet is obtained by integrating Eq. (1) of the main text with the same noise sequence.

$\mathrm{Simulation}$		1			2			3			4			5			6
$\mathrm{Method}$	1	2	3	1	2	3	1	2	3	1	2	3	1	2	3	1	2	3
$q^{\left(1\right)}$	-1	-1	-1	-1	-1	-1	-1	-1	-1	1	1	-1	-1	-1	-1	1	1	-1
#iterations	400	400	400	685	685	685	558	558	550	1609	1609	950	1664	1664	1249	1887	1887	1151

Equations

$$m{\ddot{\theta }}_i+d{\dot{\theta }}_i$$

$$m\delta \ddot{\mathbf{\theta }}+d\delta \dot{\mathbf{\theta }}$$

$${\mathbb{L}}_{ij}$$

$$⟨\delta P_i(t)\delta P_j(t^{\prime })⟩={\delta }_{ij}\delta P_0^2\exp [-|t-t^{\prime }|/{\tau }_0].$$

$$\begin{array}{rl}\underset{t\to \infty }{\lim }⟨\delta {\mathbf{\theta }}^2(t)⟩& =\delta P_0^2\sum _{\alpha \ge 2}\frac{{\tau }_0+m/d}{{\lambda }_{\alpha }({\lambda }_{\alpha }{\tau }_0+d+m/{\tau }_0)}.\end{array}$$

$$\underset{t\to \infty }{\lim }⟨\delta {\mathbf{\theta }}^2(t)⟩$$

$$\mathcal{V}(\mathbf{\theta },t)$$

$$\delta P_0^2\sum _{\alpha \ge 2}\frac{{\tau }_0+m/d}{{\lambda }_{\alpha }({\lambda }_{\alpha }{\tau }_0+d+m/{\tau }_0)}$$

$$T_{\mathrm{esc}}\propto {\left[2{\int }_{\beta \Delta }^{\infty }P(\overline{\delta \theta })\mathrm{d}(\overline{\delta \theta })\right]}^{-1}$$

$$m{\ddot{c}}_{\alpha }(t)+d{\dot{c}}_{\alpha }(t)$$

$$c_{\alpha }(t)$$

$$⟨\delta {\mathbf{\theta }}^2(t)⟩$$

$${\lambda }_{\alpha }$$

$$\delta P_0^2\le \frac{{\pi }^{2}dn}{{\tau }_0(n-2{)}^2},$$

$$\delta P_0^2\le \frac{60{\pi }^{2}n}{(n-2{)}^2(n^2+11)},$$

$${\lambda }_{\alpha }$$

$$m_i{\ddot{\theta }}_i+d_i{\dot{\theta }}_i$$

$${\Delta }^2$$

$$\mathcal{V}(\mathbf{\theta })$$

$$P_i^{(0)}$$

Full text

Introduction. Complex physical systems are mathematically modelled as dynamical systems. Equilibrium states, if they exist, are determined and characterized by fixed points, limit cycles and tori, or even strange attractors of the corresponding differential equations Ott (2002). In principle the latter should be complemented by stochastic terms to account for unavoidable perturbations from unaccountable environmental degrees of freedom van Kampen (1976). A central question of broad interest is to determine the magnitude and statistical properties of the relevant stochastic terms that could lead to the loss of equilibrium or induce transitions between different equilibria. Some physically important situations where such stochastic escape phenomena may occur are electric power grids with high penetration of fluctuating renewable energy sources Machowski et al. (2008); Auer et al. (2017); Schäfer et al. (2017), superconducting rings Gough et al. (1987) and Josephson junction arrays Il’ichev and Omelyanchouk (2008) subjected to noisy magnetic fields, as well as neuronal systems subjected to synaptic, ion-channel, neurotransmitter or membrane potential noise Braun et al. (1994); Liu et al. (2018).

Despite decades of investigations, theoretical studies of problems related to stochastic escape are generally extensions of the pioneering work of Kramers Kramers (1940), which relates chemical reaction rates to action integrals between different potential minima. The problem is analytically tractable in low dimensions only, and several recent works considered noise-induced large fluctuations in the dynamical behavior of higher-dimensional network-coupled systems through the numerical determination of action minimizing paths DeVille (2012); Schäfer et al. (2017); Hindes and Schwartz (2016, 2018). A better analytical understanding of the interplay of noise characteristics with the network topology is clearly desirable. In this manuscript we propose a resolutely different approach to stochastic escape from stable equilibria in complex, network-coupled dynamical systems, incorporating noise characteristics as well as network dynamics and topology.

For sufficiently weak, bounded noise, fluctuations are small and there is no stochastic escape Lee et al. (2018). Noise makes the system fluctuate about its equilibrium, and typical deviation amplitudes can be evaluated from a linearized dynamics about the equilibrium Bamieh et al. (2012); Tyloo et al. (2018); Haehne et al. (2018). The situation becomes fundamentally different for stronger noise. This is illustrated in Fig. 1, which shows the time-evolution of the winding number $q$ labelling different equilibrium fixed points of Kuramoto-like models, Eq. (1), with additive Ornstein-Uhlenbeck noise. Changes in $q$ indicate that the system visits other basins of attraction, surrounding different equilibrium states. Depending on the oscillator inertia and the noise amplitude and correlation time, this happens more or less quickly and for longer or shorter periods of time. Due to the high dimensionality of the state space and the nonlinear coupling between oscillators, the exact shape and size of the basins are impossible to capture Wiley et al. (2006); Menck et al. (2013); Delabays et al. (2017a), consequently, the escape time from one basin is hard to predict. For the Kuramoto model with cyclic interactions, DeVille DeVille (2012) showed that the escape time scales as the exponential of the potential barrier height between the initial and final equilibrium states. In the spirit of Kramers Kramers (1940), Hindes and Schwartz Hindes and Schwartz (2016, 2018) further relate the escape time to the numerically computed action on the action-minimizing trajectory between the two equilibria. It is hard to see how these numerical approaches could give analytical estimates for stochastic escape in higher dimension.

In this manuscript, we follow an altogether different approach. We specify to synchronous fixed points of Kuramoto-like models, but stress that the approach is applicable to more general systems. We subject the initial, synchronous state to additive Ornstein-Uhlenbeck noise. Linearizing the dynamics about the synchronous state, we calculate the standard deviation of the noise-induced fluctuations about that state. The linearized dynamics is no longer accurate when the standard deviation exceeds some threshold distance $D_{c}$. Clearly, $D_{c}$ is bounded from above by the distance $\Delta$ between the stable synchronous state and the closest saddle point to the next basin of attraction. We postulate that $D_{c}$ is parametrically proportional to $\Delta$, so that the breakdown of linear response coincides with the occurence of stochastic escapes. This postulate allows us to derive a criterion for stochastic escape based on the distance $\Delta$ between the initial stable synchronous fixed point and the nearest saddle point and not as in Kramers’ and other approaches Kramers (1940); DeVille (2012); Schäfer et al. (2017); Hindes and Schwartz (2016, 2018) on their potential height difference. We validate numerically our postulate that $D_{c}\sim\Delta$ for four, very different networks and furthermore show it gives precise estimates for the first stochastic escape time.

The Model. We consider generic, Kuramoto-like models of nonlinearly coupled oscillators on complex graphs defined by the differential equations Kuramoto (1975)

[TABLE]

Oscillators with inertia $m$ and damping parameter $d$ are described by compact angle coordinates $\theta_{i}\in(-\pi,\pi]$ and natural frequencies $P_{i}\in\mathbb{R}$. They are located on nodes $i=1,...,n$ of a connected coupling network defined by the adjacency matrix, $b_{ij}\geq 0$. Without loss of generality, we consider $\sum_{i}P_{i}=0$, which is equivalent to considering the system in a rotating frame, because Eq. (1) is invariant under $\theta_{i}(t)\rightarrow\theta_{i}(t)+\Omega t$ $P_{i}\rightarrow P_{i}+d\,\Omega$. For bounded distributions of natural frequencies on small enough intervals, synchronous states exist with $\dot{\theta}_{i}\equiv 0$, $\forall i$.

We consider a stable synchronous state ${\bm{\theta}}^{(0)}=(\theta_{1}^{(0)},\ldots,\theta_{n}^{(0)})$ corresponding to natural frequencies $\bm{P}^{(0)}$. We subject this state to a time-dependent perturbation $\bm{P}(t)=\bm{P}^{(0)}+\delta\bm{P}(t)$. Linearizing the dynamics defined by Eq. (1) with $\bm{\theta}(t)=\bm{\theta}^{(0)}+\delta\bm{\theta}(t)$, one obtains

[TABLE]

with the weighted Laplacian ${\mathbb{L}}(\{\theta_{i}^{(0)}\})$ defined by

[TABLE]

This matrix is positive semidefinite, with a single eigenvalue $\lambda_{1}=0$ and associated eigenvector ${\bf u}_{1}=(1,1,1,...1)/\sqrt{n}$, while $\lambda_{\alpha}>0$, $\alpha=2,3,...n$.

The dynamics of Eq. (2) is characterized by different times scales. The first one characterizes the noisy perturbations. We consider spatially uncorrelated noise with vanishing average and Ornstein-Uhlenbeck correlator

[TABLE]

Thus, the perturbation is characterized by its variance, $\delta P_{0}^{2}$ and its correlation time, $\tau_{0}>0$. The second time scale is $m/d$. It gives the typical time over which local excitations are damped by $d$, neglecting the network dynamics. Finally, one has a set of time scales $d/\lambda_{\alpha}$, $\alpha=2,...n$, each of them defined by the ratio of the damping parameter and an eigenvalue of the Laplacian. For $m/d>d/4\lambda_{\alpha}$ these correspond to oscillation time scales of the Laplacian modes, while for $m/d<d/4\lambda_{\alpha}$ they give network-dynamical corrections to the damping time scale. We consider $\tau_{0}$ as a tunable parameter allowing us to explore different regimes depending on its relation with $m/d$ and $d/\lambda_{\alpha}$.

We measure the distance between the state of the system and the initial synchronous state as the square root of the variance $\langle\delta{\bm{\theta}}^{2}(t)\rangle=\sum_{i}\langle[\delta{\theta_{i}}(t)-\delta\overline{\theta}(t)]^{2}\rangle$ with $\delta\overline{\theta}(t)=n^{-1}\sum_{i}\delta\theta_{i}(t)$ and brackets indicating an average over different realizations of noise with the same first two moments. It appropriately gives the standard deviation of the angle deviations in the subspace orthogonal to ${\bf u}_{1}$, because displacements in that subspace do not change the state. To calculate $\langle\delta{\bm{\theta}}^{2}(t)\rangle$, we expand angle deviations over the eigenbasis of ${\mathbb{L}}$ and solve Eq. (2) for the coefficients of that expansion (See Supplemental Material SM ). We obtain the long-time limit

[TABLE]

In the two limits of long and short $\tau_{0}$, one has

[TABLE]

with $Kf_{p}=n\sum_{\alpha\geq 2}\lambda_{\alpha}^{-p}$ Klein and Randić (1993); Tyloo et al. (2018). Interestingly, none of these asymptotics depend on inertia.

Escape from the basin. The dynamics of Eq. (1) is described by a vector function $\bm{\theta}(t)$ following the gradient of the potential

[TABLE]

starting from $\bm{\theta}(t=0)=\bm{\theta}^{(0)}$. When the noisy perturbation tilts this potential strongly enough, $\bm{\theta}$ can escape the basin of attraction of $\bm{\theta}^{(0)}$. DeVille showed that, for not too large $\delta P_{0}$, the system almost surely escapes the basin in a neighborhood of a $1$-saddle DeVille (2012). Comparing the typical distance between $\bm{\theta}$ and $\bm{\theta}^{(0)}$ of Eq. (7) with the distance $\Delta$ between $\bm{\theta}^{(0)}$ and its closest $1$-saddle $\bm{\varphi}$ gives us a parametric condition for noise-induced stochastic escape

[TABLE]

Our task is therefore to identify the position of the $1$-saddles. This is in general no trivial task because the geometry of basins of attraction in such high-dimensional problems is impossible to fully capture. For single-cycle networks with identical frequencies, $1$-saddles can be identified analytically DeVille (2012); SM . For more general networks, we construct a numerical algorithm which locates $1$-saddles $\bm{\varphi}$ and constructs the distribution of their distance to $\bm{\theta}^{(0)}$ (See Supplemental Material SM ).

Numerical simulations. We first check Eq. (12) against numerical simulations of the Kuramoto model of Eq. (1) with $m=0$. We consider four different networks with constant couplings $b_{0}=1$ and identical frequencies, which are a single-cycle network with nearest-neighbor coupling, a single-cycle with nearest- and $3^{\rm rd}$-neighbor coupling, a model of the UK transmission network Delabays et al. (2017a) and a realization of a small-world network Watts and Strogatz (1998). Details about these networks are given in the Supplemental Material SM . At each node, natural frequencies are perturbed by spatially uncorrelated Gaussian noisy sequences $\delta P_{i}(t)$ satisfying Eq. (6). We integrate the dynamics of Eq. (1) during an observation time $T_{\rm obs}$ and check for a stochastic escape at every time step. Our method for detecting such occurences is based on Refs. Dörfler et al. (2013); Delabays et al. (2017b); Manik et al. (2017) which showed that on meshed networks, different fixed-point solutions of Eq. (1) correspond to a vector of winding numbers ${\bm{q}}$, each component corresponding to one of the cycles of the network. Refs. DeVille (2012); Hindes and Schwartz (2018) observed that transitions between different such equilibrium states occur by phase slips of few oscillators, and we show in the Supplemental Material SM that these slips can be detected by recording the time evolution of ${\bm{q}}$, as illustrated on Fig. 1. We therefore detect desynchronizing events through variations of winding numbers. Details of the method and comments on its accuracy are presented in the Supplemental Material SM . For each set of noise parameters $\delta P_{0}$ and $\tau_{0}$ we perform several calculations corresponding to different noise realizations.

Fig. 2 shows the fraction $\mathbb{P}$ of runs that remain in the initial basin for $t\leq T_{\rm obs}$. The parameter space is sharply divided into (a) the red region (denoted $U$ for ”unstable”) where all runs left the basin of attraction before $T_{\rm obs}$, (b) the blue region (denoted $S$ for ”stable”), where none of the runs left the initial basin of attraction and (c) a rather narrow intermediate region between $U$ and $S$ where some runs left and some runes stayed in the initial basin.

It is quite remarkable that the intermediate region (c) is qualitatively if not quantitatively identifed by Eq. (12) with a network-dependent $\Delta$. As discussed above, $\Delta$ is given by a typical distance between the initial stable fixed point $\bm{\theta}^{(0)}$ and the nearest saddle point $\bm{\varphi}$ roughly giving the smallest linear size of the basin of attraction. For the single-cycle network, all $1$-saddles are located at the same distance from $\bm{\theta}^{(0)}$, which can be obtained analytically DeVille (2012) (See Supplemental Material SM ). For the other three networks, many, though likely not all 1-saddles are identified numerically (See the Supplemental Material for details of the method SM ). For the single-cycle network with nearest- and $3^{\rm rd}$-neighbor coupling, all the $1$-saddles we find are located at the same distance $\Delta$ from $\bm{\theta}^{(0)}$. For the UK and small-world networks, on the other hand, we find a distribution of $\Delta\in[\Delta_{\rm min},\Delta_{\rm max}]$, which is likely due to the complexity of those meshed networks. The yellow dashed lines in Fig. 2 then indicate our theoretical prediction Eq. (12) for the obtained value $\Delta$ for the two single-cycle networks and for values of $\Delta$ corresponding to the 25th, the 50th and the 75th precentiles of the distribution of $\Delta$ for the UK and small-world networks. In all cases, the shape of the boundary is well predicted. For the more complex UK transmission network, Fig. 2(c), there is a horizontal shift between theory and numerics, presumably due to to stronger anisotropies of the basins of attraction in this more complex network, effectively requiring a larger $T_{\rm obs}$.

In the case of bounded noise, we expect an inertialess system to remain in its initial basin for weak enough noise Lee et al. (2018). However, the noise considered in our case is Gaussian and arbitrarily large excursion will occur if one waits long enough. As a matter of fact, we found that increasing $T_{\rm obs}$ shifts the boundary between stable and unstable regions to lower $\delta P_{0}$ (see Supplemental Material SM ). Fig. 3 further shows the stochastic escape time as a function of $\delta P_{0}$. A superexponential behavior is observed which can be understood as follows. The noise generates a distribution of angle deviations which we expect to be Gaussian with a variance given by Eq. (7). The escape time is then inversely proportional to the probability to have such a deviation exceeding $\Delta$, i.e.

[TABLE]

with a free parameter $\beta$ of order $1$. Fig. 3 validates this argument using a Gaussian distribution of single-angle deviation $P(\overline{\delta\theta})$ with variance $\langle\delta{\bm{\theta}}^{2}(t)\rangle/n$, see Eq. (7). We have found, but do not show, that $T_{\mathrm{esc}}$ diverges at a finite value of $\delta P_{0}$ for a box-distributed, bounded noise.

We finally consider Eq. (1) with nonzero inertia. We focus on the single-cycle network with nearest- and $3^{\rm rd}$-neighbor coupling, and tune the inertia parameter $m$ to explore different regimes defined by the different time scales of Eq. (1). Fig. 4 shows the difference in survival probabilities with and without inertia in the regimes (a) ${d}/{\lambda_{\alpha}}\gtrsim{m}/{d}$, (b) ${d}/{\lambda_{\alpha}}\lesssim{m}/{d}$ and (c) $d/\lambda_{\alpha}\ll m/d$. Deep in the stable (unstable) regions, both inertialess and inertiaful models have $\mathbb{P}=0$ ($\mathbb{P}=1$) and the difference $\mathbb{P}(m=0)-\mathbb{P}(m)=0$. Somehow counterintuitively, however, there is an intermediate region where the presence of inertia facilitates stochastic escape compared to the inertialess case, $\mathbb{P}(m=0)-\mathbb{P}(m)>0$. The boundary of that region are in excellent agreement with the prediction of Eq. (12), giving the two dashed yellow lines for $m=0$ and $m\neq 0$.

For large $\tau_{0}$, the faster escape of the system with finite inertia is easily understood. With long correlation time, the noise tends to push the system in the same direction for long sequences. This is sufficient to have the inertiaful system accumulate a significant kinetic energy. The system keeps then moving, even if, after some time, the noise starts pushing the other way and allows it to move above a saddle point with inertia, whereas the inertialess system is immediately stopped by noise reversal.

For smaller $\tau_{0}$, on the other hand, inertia resists short sequences of pushes in rapidly varying directions and accordingly, we found that inertia stabilizes the system in that case (See Supplemental Material SM ). This is not predicted by Eq. (12) and is probably due to contributions beyond our linear response theory, because discrepancies appear for values of $\delta P_{0}$ comparable to the coupling strength $b_{0}$. The influence of inertia on stochastic escapes is perhaps best illustrated in Fig. 1, where the presence of inertia stabilizes the system under short-correlated noise [panel (a)] but leads to more frequent stochastic escapes for long-correlated noise [panel (b)].

Conclusion. We have constructed a novel approach to stochastic escape, based on a spectral calculation of typical distances of stochastic excursions about equilibrium states and the evaluation of the distance between this equilibrium state and 1-saddles. The method provides analytical results with a single, model-dependent free parameter of order one [$\beta$ in Eq. (13)]. It gives remarkably accurate estimates for stochastic escape times, as is illustrated in Fig. 3. Interestingly, we found that the presence of inertia leads to faster, more frequent escapes for long noise coherence times, while the effect is reversed for short noise coherence times. This is illustrated in Fig. 1. Further studies should consider the effect of spatially correlated noise and non-Gaussian, long-tailed noise distributions Haehne et al. (2018).

This work has been supported by the Swiss National Science Foundation under grants 200020_182050 and PYAPP2_154275.

I Details of Calculations for the variance of the angle displacements

We give some details of the calculation that leads to Eq. (5) of the main text. Expanding the angle deviations over the eigenmodes of the Laplacian Eq. (3) of the main text, i.e., $\delta{\bm{\theta}}(t)=\sum_{\alpha}c_{\alpha}(t){\bf u}_{\alpha}$, Eq. (2) of the main text becomes,

[TABLE]

With the help of a Laplace transform, the solution of Eq. (S1) is given by

[TABLE]

with $\Gamma_{\alpha}=\sqrt{(d/m)^{2}-4\lambda_{\alpha}/m}$. Taking advantage of the orthogonality between eigenmodes of the Laplacian we have,

[TABLE]

with $\delta\overline{\theta}(t)=n^{-1}\sum_{i}\delta\theta_{i}(t)$. Inserting Eq. (S1) into Eq. (S3), using the time correlator of $\delta{\bm{P}}$ Eq. (4) of the main text, and finally taking the long time limit one obtains, after some algebra, Eq. (5) of the main text.

II Method to determine escape time

Various methods can be used to determine, at any iteration step of the simulation, if the system under consideration has escaped its initial basin of attraction. We compared three of them, which we detail here.

Method 1. As stated in the main text, stable equilibria of Eq. (S8) can be unambiguously distinguished by their winding vector $\bm{q}$. The method that we used for the numerical simulations in the main text proceeds as:

At each time step, compute $\bm{q}$; 2. 2.

If ${\bm{q}}\neq{\bm{q}^{(0)}}$ the winding vector of the initial basin of attraction, check if the system is still in the initial basin. To do so, simulate the dynamics without noise, taking the current state of the system as initial conditions. Once synchrony is reached, compute the winding vector $\bm{q}^{(1)}$; 3. 3.

If ${\bm{q}^{(1)}}\neq{\bm{q}^{(0)}}$, then the system was out of the initial basin. Otherwise, if ${\bm{q}^{(1)}}={\bm{q}^{(0)}}$, the system was still in the basin and thus the simulation can move to the next time step.

Method 2. This method is based on DeVille’s observation DeVille (2012) that escapes from basins of attraction occur on a short time interval and can be identified by a fast slip of a small group of angles. It proceeds as:

At each time step, check if some angles made a large excursion, i.e., $\|{\bm{\theta}}(t)-{\bm{\theta}^{(0)}}\|_{\infty}>2\pi$; 2. 2.

If so, then simulate the dynamics without noise, taking the current state of the system as initial conditions, until it synchronizes to the state ${\bm{\theta}}^{(1)}$; 3. 3.

If ${\bm{\theta}^{(1)}}\neq{\bm{\theta}^{(0)}}$, then the system was out of the initial basin. Otherwise, if ${\bm{\theta}^{(1)}}={\bm{\theta}^{(0)}}$, the system was still in the basin and thus the simulation can move to the next time step.

Method 3. Finally, we tested the method in which we check at every time step whether the system returns to the initial basin or not. This method guarantees to find the best estimate of the escape time, at least for the Kuramoto model ($m=0$), but is very time-consuming.

Table 1 compares escape times and final winding numbers for a single-cycle of $n=83$ nodes. For the Kuramoto model ($m=0$) the three methods give very similar results. For the case with inertia, the first two give larger escape times compared to the last method. We explain this as follows. When the noise is removed, the system may have accumulated some kinetic energy that will drive it out of the basin of attraction. And this can happen before the winding number changes or a large angle excursion occurs. Furthermore, if the perturbation was still active, it could have pushed the system back towards the stable fixed point before it leaves the basin of attraction, increasing the escape time.

III The four networks

We briefly describe the networks used for the numerical simulations of the main text.

III.1 Cycle with nearest neighbors coupling

We consider a cycle network of size $n$, with identical natural frequencies. The eigenvalues of its weighted Laplacian, Eq. (3) of the main text, can be obtained analytically,

[TABLE]

where $\delta$ is the angle difference between neighboring sites (which are identical at a stable equilibrium Delabays et al. (2016)) and $k_{\alpha}=2\pi(\alpha-1)n^{-1}$. For $n=83$ we have $\lambda_{\alpha}\in[0,4\cos(\delta)]$ and $\lambda_{2}=0.0057$.

Eq. (6) in the main paper can be explicitly calculated for cyclic networks as functions of the number of nodes $n$

[TABLE]

Fig. S1 shows the maximum values of $\delta P_{0}$ satisfying Eqs. (S5), (S6). One remarks that, while increasing the size of the cycle, the stable region gets smaller and even vanishes for $n\rightarrow\infty$ similarly to fluctuations that destroy long-range order in $1$ dimensional locally interacting quantum magnets Giamarchi (2004).

III.2 Cycle with nearest- and $3^{\mathrm{rd}}$-neighbors coupling

We consider a cycle network of size $n$, where each vertex is connected to its nearest- and $3^{\rm rd}$-neighbors [see Fig. S2(b)]. With identical natural frequencies, the eigenvalues of its weighted Laplacian, Eq. (3) of the main text, can be obtained analytically,

[TABLE]

where $\delta$ is the angle difference between neighboring sites (which are identical at a stable steady-state Delabays et al. (2016)) and $k_{\alpha}=2\pi(\alpha-1)n^{-1}$. For $n=83$ we have $\lambda_{\alpha}\in[0,8\cos(\delta)]$ and $\lambda_{2}=0.057$.

III.3 UK transmission grid

Model of the electrical transmission grid of UK depicted in Fig. S2(c). It is composed of $120$ nodes and $165$ edges making $44$ cycles. During the numerical simulations, to check whether the system has left the initial basin of attraction or not, we check the winding number on each cycle, i.e., the winding vector ${\bm{q}}=(q_{1},...,q_{44})$. The second eigenvalue of its Laplacian matrix is $\lambda_{2}\approx 0.013$.

III.4 Small world

A small world network is constructed from an initial network, where some edges are randomly rewired (see Watts and Strogatz (1998)). In our case, the initial network ${\cal G}_{0}$ is a cycle with $n=200$ vertices and where each vertex is connected to its first and second neighbors [see Fig. S2(a)]. Each edge $(i,j)$ is then replaced with probability $p=0.05$ by the edge $(i,k)$, where $k$ is chosen at random among the vertices not already connected to $i$. The network obtained ${\cal G}_{p}$ is illustrated in Fig. S2(d). It is a small world as it has a large relative clustering coefficient $C({\cal G}_{p})/C({\cal G}_{0})\approx 0.89$ and a small relative characteristic path length $L({\cal G}_{p})/L({\cal G}_{0})\approx 0.32$ (see Watts and Strogatz (1998) for more details). The second eigenvalue of its Laplacian matrix is $\lambda_{2}\approx 0.046$.

IV Finding 1-saddles

We detail our methods for finding 1-saddles (equilibria with a unique unstable direction) of the dynamical system

[TABLE]

for arbitrary coupling graph.

IV.1 Cycle Networks

For cycle networks with nearest neighbor coupling and identical natural frequencies, the distance between the stable equilibrium $\bm{\theta}^{(0)}=(0,...,0)$, and the 1-saddle $\bm{\varphi}$, can be computed analytically as Delabays et al. (2017a)

[TABLE]

IV.2 General Networks

For general networks, the anisotropy of the basins of attraction renders the 1-saddles complicated to identify analytically. We propose a numercial method to locate 1-saddles, which is based on two results of DeVille DeVille (2012):

• •

Escapes from basins of attraction almost always occur in a neighborhood of a 1-saddle of the potential

[TABLE]

• •

Transitions from a basin to another occur on a short time interval compared to the time the system remains in a basin of attraction.

We numercially integrate Eq. (S8), where $\delta P_{i}$ is a noise with small variance, and keep track of the angles in order to identify iterations where the system is close to a 1-saddle. As observed in DeVille (2012), when the system is driven (by the noise) to another basin of attraction, its trajectory goes close to a 1-saddle, and this can be seen in the time-evolution of the angles as a fast jump of a set of angles of amplitude $2\pi$ (see Fig. S3). The state $\bm{\varphi}^{(0)}$ of the system in the middle of this jump will be a candidate for a 1-saddle. This state is probably not exactly a 1-saddle, but according to DeVille (2012), it should be close to one. We then solve the steady-state equations

[TABLE]

using a Newton-Raphson method with initial conditions $\bm{\varphi}^{(0)}$. This gives an equilibrium $\bm{\varphi}^{*}$ of Eq. (S8), which we expect to be close to $\bm{\theta}^{(0)}$. Computing the eigenvalues of the Jacobian of Eq. (S8), the equilibrium $\bm{\varphi}^{*}$ is a $p$-saddle if and only if it has $p$ positive eigenvalues. Note that one eigenvalue is always zero due to invariance of Eqs. (S8) and (S10) under a constant shift of all angles.

Running this simulation for a long enough time, we identified:

• •

$284$ 1-saddles for the cycle with nearest- and $3^{\rm rd}$-neighbor. The distribution of their distance to the stable equilibrium $\bm{\theta}^{(0)}$ is given in Fig. S4(a). Looking more into details, we observe that each value in Fig. S4(a) corresponds to a unique 1-saddle, up to an index shift or the angles’ sign reversal. The 1-saddles with the two smallest norm, $\bm{\varphi}^{(1)}$ and $\bm{\varphi}^{(2)}$, are represented in Fig. S5. The first one [Fig. S5(a)] has the smallest 2-norm, but its configuration with $n-1$ equal angles and one angle $\pi$ apart from all others is, in our opinion, unlikely to occur. As we consider noisy perturbation at all nodes, a configuration with a single large angle excursion and no excursion for all other nodes seems less likely than a configuration where all angles are slightly displaced from their neighbors. In the main text, we performed our study using $\bm{\varphi}^{(2)}$ as 1-saddle for the cycle with nearest- and $3^{\rm rd}$-neighbor.

• •

$788$ 1-saddles for the UK network, whose distribution of the distances to the stable equilibrium is given in Fig. S4(b). Distances cover a large range of value, due to the anisotropy of the basin of attraction;

• •

$4956$ 1-saddles for the small-world network. The distribution of the distances to $\bm{\theta}^{(0)}$ is given in Fig. S4(c). Most of the 1-saddles are at similar distance.

V Superexponential Escape Time

To evaluate the influence of the observation time $T_{\rm obs}$ on Fig. (2) of the main text, we performed the simulation for the cycle, increasing the observation time. Fig. S6 shows the fraction on simulations that stay in the initial basin of attraction after an observation time satisfying $\lambda_{2}T_{\rm obs}/d=14.2$ [Fig. S6(a)], $142.4$ [Fig. S6(b)], $569$ [Fig. S6(c)], for a cycle network with $n=83$ nodes. As $T_{\rm obs}$ increases exponentially, we observe the boundary between region $U$ and $S$ drifting to the left due to the escape time that is superexponential as $\delta P_{0}$ decreases.

VI Linearization Break-Down

In the main text, we show that, according to our theory, inertia always destabilizes the system compared to the inertialess case. However, for the cycle network, we found that for small $\tau_{0}$ and large $\delta P_{0}$, inertia stabilizes the system, as illustrated on Fig. S7. The blue area where inertia stabilized the system is not predicted by our theory, Eq. (8) of the main text. This can be explained by the breakdown of the linear approximation. Indeed, the blue region on Fig. S7 starts for value of the order of the coupling $\delta P_{0}\cong b_{0}\equiv 1$.