The Size of the Sync Basin Revisited

TL;DR

This paper introduces an efficient numerical method to estimate the size of basins of attraction in high-dimensional dynamical systems, and applies it to the Kuramoto model, revealing new scaling laws for basin volumes.

Contribution

The authors develop a novel algorithm for accurately estimating basin volumes and demonstrate its effectiveness on complex oscillator models, revisiting previous assumptions about basin size distributions.

Findings

Basin volumes scale as (1-4q/n)^n, differing from previous Gaussian assumptions.

The method successfully applies to models with diverse natural frequencies and network structures.

Revisits and refines understanding of basin size distributions in the Kuramoto model.

Abstract

In dynamical systems, the full stability of fixed point solutions is determined by their basin of attraction. Characterizing the structure of these basins is, in general, a complicated task, especially in high dimensionality. Recent works have advocated to quantify the non-linear stability of fixed points of dynamical systems through the relative volumes of the associated basins of attraction [D. A. Wiley {\it et al.} Chaos {\bf 16}, 015103 (2006), P. J. Menck {\it et al.} Nat. Phys. {\bf 9}, 89 (2013)]. Here we revisit this issue and propose an efficient numerical method to estimate these volumes. The algorithm first identifies stable fixed points. Second, a set of initial conditions is considered that are randomly distributed at the surface of hypercubes centered on each fixed point. These initial conditions are dynamically evolved. The linear size of each basin of attraction is finally determined by the proportion of initial conditions which converge back to the fixed point. Armed with this algorithm, we revisit the problem considered by Wiley et al. in a seminal paper [D. A. Wiley {\it et al.} Chaos {\bf 16}, 015103 (2006)] that inspired the title of the present manuscript, and consider the equal-frequency Kuramoto model on a cycle. Fixed points of this model are characterized by an integer winding number $q$ and the number $n$ of oscillators. We find that the basin volumes scale as $(1-4q/n)^n$, contrasting with the Gaussian behavior postulated in Wiley et al.'s paper. Finally, we show the applicability of our method to complex models of coupled oscillators with different natural frequencies and on meshed networks.

Tables

Table 1. Table 1: Maximal values of q 𝑞 q with respect to β 𝛽 \beta , estimated by Eq. ( B.6 ) and obtained with 5 random frequency distributions following the simulation process of Sec.  II.4 .

$\beta$	$q_{\max }\left(\beta \right)$: Eq. (B.6)	$q_{\max }\left(\beta \right)$: Sec. II.4
0	20.75	$\{20,20,20,20,20\}$
0.01	17.14	$\{17,17,18,18,18\}$
0.02	15.62	$\{16,16,17,17,17\}$
0.05	12.56	$\{13,14,14,15,15\}$
0.1	8.97	$\{10,11,12,12,13\}$

Equations

$${\dot{\theta }}_i$$

$$p(q)=(\sqrt{2\pi }\sigma {)}^{-1}\exp \left(-q^2/2{\sigma }^2\right).$$

$${\dot{\theta }}_i$$

$${\Delta }_{i,i+1}$$

$$q(\vec{\theta })$$

$$\Delta =2\pi q/n$$

$${\theta }_i^{(q)}$$

$$(n-1)\cdot {\Delta }^{\prime }+\pi -{\Delta }^{\prime }=2\pi q^{\prime }$$

$${\phi }_i^{(q^{\prime })}$$

$$T_i^{(k)}$$

$${\theta }_i^{(q)}-{\phi }_i^{(q)}$$

$$\Vert {\vec{\theta }}^{(q)}-{\vec{\phi }}^{(q)}{\Vert }_{\infty }$$

$$\underset{n\to \infty }{\lim }\Vert {\vec{\theta }}^{(q)}-{\vec{\phi }}^{(q)}{\Vert }_{\infty }$$

$$V_q$$

$${\Delta }^{\prime }=-\frac{\pi }{n-2}$$

$${\vec{\eta }}_{q,j,\alpha }:={\vec{\theta }}^{(q)}+\pi \alpha \text{leavevmode}{\overrightarrow{ϵ}}_j,$$

$$p_q(\alpha )$$

$${\alpha }_{\tau }(q)$$

$${\dot{\theta }}_i$$

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

$$-\frac{\partial \mathcal{V}}{\partial {\theta }_i}$$

$${\vec{q}}_G(\vec{\theta })=(q_1,q_2,...,q_m),$$

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

$${\theta }_i=q\arctan \left(\frac{y_i-y_0}{x_i-x_0}\right),$$

$${\vec{q}}_{\mathrm{sh}}=(q_1,q_2,q_3,q_4,q_5).$$

$$\sin (\Delta )=\sin ({\Delta }^{\prime })+\sin ({\Delta }^{\prime \prime })\text{leavevmode},$$

$$(m-\ell )\cdot {\Delta }^{\prime \prime }-\ell \cdot {\Delta }^{\prime }=2\pi q^{(m)}=0.$$

$${\phi }_i^{(q^{\prime })}$$

$$S_k:=n^{-1}\left[\sum _{j=0}^{k-1}j{\Delta }^{\prime }+\sum _{j=k}^{n-1}\left((j-2){\Delta }^{\prime }+\pi \right)\right].$$

$${\phi }_i^{(q^{\prime })}$$

$$T_i^{(k)}$$

Full text

Many natural systems of coupled elements, such as fireflies, pacemaker cells or electrical grids, exhibit synchronization phenomena. When a system synchronizes, each of its components behaves coherently with respect to the others, due to the coupling between them. A central issue in various fields of science and engineering is to understand how robust is this synchronized state against external perturbations or imperfections in the system. Mathematically, this problem is usually hard. Here we propose a tractable numerical approach which first identifies the synchronization states and second evaluates the magnitude of the largest perturbation such that the system converges back to its initial synchronized state. Our method allows to tackle large complex systems in a reasonable computation time with good resolution, which was not the case for numerical methods proposed so far.

I Introduction

Models of coupled dynamical systems are widely used to investigate collective behaviors in complex systems. One particularly puzzling phenomena is that of synchrony, where different individual dynamical systems start to behave coherently when sufficiently strongly coupled. Ott (2002); Strogatz (2004, 2000); Acebrón et al. (2005); Dörfler and Bullo (2014) The Kuramoto model was introduced Kuramoto (1975, 1984a) to describe such synchronizing behaviors. The model considers a set of $n$ coupled harmonic oscillators, with angle coordinate $\theta_{i}$ and natural frequency $P_{i}$,

[TABLE]

where $K_{ij}\in\mathbb{R}$ is the coupling constant between oscillators $i$ and $j$. In its original formulation, the Kuramoto model considers identical all-to-all coupling, $K_{ij}\equiv K/n$. Kuramoto (1975, 1984a) It was found that for a coupling constant $K$ exceeding a critical value $K_{c}$, a finite non-empty set $F\subset\{1,...,n\}$ of oscillators synchronizes, i.e. $\dot{\theta}_{i}-\dot{\theta}_{j}=0$, for $i,j\in F$. This type of synchrony, where oscillators rotate at the same frequency, $\dot{\theta}_{i}=\dot{\theta}_{j}$, for $i,j\in F$, but not necessarily the same phase, is called frequency synchronization. Phase synchronization, where additionally $\theta_{i}=\theta_{j}$, $\forall i,j$, is in general not achievable for heterogeneous natural frequencies. In this manuscript, “synchronization” refers to “frequency synchronization”.

The value of $K_{c}$ can be computed by solving an implicit equation. Dörfler and Bullo (2014); Aeyels and Rogge (2004); Mirollo and Strogatz (2005); Verwoerd and Mason (2008) The critical coupling $K_{c}$ depends on the distribution of the natural frequencies $g(P)$. In particular, if the support of the distribution $g(P)$ is compact, full-synchrony, i.e. $\dot{\theta}_{i}-\dot{\theta}_{j}=0$ for all $i,j$, is reached for large enough $K$. Ermentrout (1985); van Hemmen and Wreskinski (1993)

The Kuramoto model in its various versions has evolved into a paradigm for investigating synchronizing behaviors. Its popularity stems from its simple formulation, which allows a tractable analytical treatment while still capturing the essence of the synchronizing behavior of many real systems, in fields as diverse as physics, Wiesenfeld et al. (1996) chemistry, Kuramoto (1984b) biology Sompolinsky et al. (1991) or electrical engineering. Dörfler et al. (2013) These two advantages explain the popularity of the Kuramoto model in science and engineering, which led to many generalizations reviewed for instance in Refs. Acebrón et al., 2005; Dörfler and Bullo, 2014.

Full frequency synchrony in the Kuramoto model is reached at a fixed point of Eq. (1). A natural question is then to assess the stability of these fixed points. Most works on the stability of fixed points of dynamical systems rely on the seminal work of Lyapunov Lyapunov (1892, 1992) at the end of the XIXth century. Lyapunov first investigated linear stability, by linearizing the dynamics around a given fixed point. For small deviations, the dynamics is determined by a stability matrix whose eigenvalues are called Lyapunov exponents. The fixed point is linearly stable if all its Lyapunov exponents are non-positive. This guarantees that small enough deviations go exponentially fast to zero. Linear stability is however a local concept and gives no information about the stability of the system against large perturbations.

Lyapunov went beyond linear stability with his second method, which assesses stability based on the existence of what is now called a Lyapunov function. Lyapunov (1892, 1992) The latter generalizes the concept of energy for the states of a dynamical system. The Lyapunov function of a system can be used to determine the basin of attraction of a given fixed point, Pai (1981); Chiang and Thorp (1989) which is the set of all initial conditions converging dynamically to this fixed point. A global measure of the stability of a fixed point is given by the volume of its basin of attraction - this has been called basin stability. Wiley et al. (2006); Menck et al. (2013); Schultz et al. (2017) Clearly, the larger the basin of attraction, the more likely it is to reach the corresponding fixed point dynamically. This gives a global measure of the stability of a fixed point.

In networks of all-to-all coupled oscillators, tight estimates of the volume of the basin of attraction of the sychronous state are known. Dörfler and Bullo (2011) Much less is known about the basins of attraction of cycle networks. As pointed out by Korsak Korsak (1972) already in 1972 in the context of electrical networks, the Kuramoto model on a cycle network admits several stable fixed points, characterized by their winding numbers (to be defined in Sec. II). In Ref. Wiley et al., 2006, Wiley et al. considered such network topology, with identical frequencies and investigated how the volume of the basin of attraction of a stable fixed point is related to its winding number. In particular, they are interested in the likelihood of the system to reach the phase synchronous state. This likelihood is directly related to the volume of the basin of attraction of the phase synchronous state, which they call “sync basin”. Starting from random initial conditions, they numerically evolved the system until it converged to a stable fixed point. The volume of the basin of attraction of every fixed point was then estimated by the proportion of initial conditions that converged to it. It was found that the volume of the basins of attraction follow a Gaussian distribution with respect to the winding numbers $q$, as shown in Fig. 1 (red dots).

One issue with that procedure is that the winding numbers of randomly chosen initial conditions also follow a Gaussian distribution,

[TABLE]

From the data shown in Fig. 1, we obtain a standard deviations of $\sigma=2.63$ for the distribution of winding numbers of the initial conditions and a narrower distribution with $\sigma=1.63$ for the converged fixed points (the latter value in agreement with Ref. Wiley et al., 2006). The Gaussian distribution for the initial conditions can easily be understood once one realizes that picking an initial condition is similar to a random walk. Pathria and Beale The node index along the cycle corresponds to a time step index and the angle on each node gives by how much and in what direction the random walk progresses. Large winding numbers correspond then to random walks with large excursions. This analogy explains the obtained Gaussian distribution for initial winding numbers.

We also observe that the winding number of the initial conditions and of the converged fixed points are significantly correlated, with a correlation coefficient of 0.47. Therefore, if one does not have enough resolution for the initial conditions, the distribution of the winding number of the final states may, at least partially, reflect the initial distribution of $q$ instead of the volume of the basins of attraction. Due to the high dimensionality of the state space ($n=83$ in Fig. 1 and $n=80$ in Ref. Wiley et al., 2006), simulations with random initial conditions would need an unfeasible number of runs to representatively cover the whole state space.

To the best of our knowledge, the only paper, beside Ref. Wiley et al., 2006, focusing specifically on the basins of attraction of the Kuramoto model on cycle network is Ref. Ha and Kang, 2012, which analytically obtains lower bounds on the volume of the basins of attraction for a cycle of Kuramoto oscillators with unidirectional coupling.

We therefore revisit this issue by constructing a new systematic numerical method. Our approach is first to find all the stable fixed points of Eq. (1) describing all the possible frequency-synchronous states, and second to perturb them in random directions with an increasing magnitude to assess the volume of their basins of attraction. The volume is estimated from the magnitude of the largest perturbation still converging to the initial fixed point. In the case of a single cycle with identical frequencies, one can analytically identify all stable fixed points and the problem is sufficiently tractable to obtain an analytical estimate of the volume of the basins of attraction, which we confirm numerically. We show that for $n\gg 1$ and $q$ not too small, the volume of the basins of attraction scales as $V_{q}\sim(1-4q/n)^{n}$ instead of the Gaussian law of Wiley et al. [Eq. (2)]. We then extend our perturbation procedure to cycle networks with non-identical frequencies and to meshed networks with identical frequencies. Our numerical method guarantees that we investigate every basin of attraction with a representative number of initial conditions. It is based on (i) a numerical procedure to systematically find stable fixed points of Eq. (1) on any meshed network and (ii) the perturbation procedure described above. We believe that our method could be applied to other fields of research such as planar spin glasses Mézard et al. (1987) or disordered Josephson junction arrays Fazio and van der Zant (2001) in condensed matter physics.

II Single cycle

We first revisit the model of Ref. Wiley et al., 2006 and consider the Kuramoto model on a cycle with $n$ nodes, equal frequencies $P_{i}\equiv P_{0}$, for all $i$, and identical coupling $K_{ij}\equiv K$, for all connected nodes $i$ and $j$. In a frame rotating with angular frequency $P_{0}$, after the change of variables $\theta_{i}\to\theta_{i}+P_{0}t$, Eq. (1) reduces to

[TABLE]

where indices are taken modulo $n$.

We define the angle difference $\Delta_{ij}\coloneqq\theta_{i}-\theta_{j}$ taken modulo $2\pi$ in the interval $(-\pi,\pi]$. Fixed points of Eq. (3) satisfy either

[TABLE]

where the sign in front of $\pi$ in the right-hand side is chosen to ensure that $\Delta_{i,i+1}\in(-\pi,\pi]$.

Given an angle vector $\vec{\theta}=(\theta_{1},...,\theta_{n})$ we define the integer winding number on the cycle

[TABLE]

When summing the angle differences around the cycle, we have to end at an integer multiple of $2\pi$ to guarantee single-valuedness of angles. Therefore, $q(\vec{\theta})\in\mathbb{Z}$.

Remark**.**

For any graph topology, we can always define a winding number on every cycle. A fixed point is then characterized by a winding vector, whose components are the winding numbers on each cycle.

According to Ref. Delabays et al., 2016, any stable fixed point of Eq. (3) must have all angle differences between neighboring oscillators in $[-\pi/2,\pi/2]$. Eq. (4) then implies that there is a unique stable fixed point, $\vec{\theta}^{(q)}$, with winding number $q$,

[TABLE]

for all $i=1,...,n$, where $\theta_{0}^{(q)}$ is an arbitrary uniform angle shift. This implies that for a fixed point to be stable, the winding number cannot be larger than $q_{\max}\coloneqq{\rm Int}(n/4)$, which would imply angle differences larger than $\pi/2$ otherwise. Winding numbers for stable fixed points then range from $-q_{\max}$ to $+q_{\max}$, which gives $2q_{\max}+1$ stable fixed points, Delabays et al. (2016) with the caveat that if $n$ is a multiple of $4$, the fixed point with winding number $q_{\max}$ has all Lyapunov exponents equal to zero, and its basin of attraction has measure zero. Manik et al. (2016)

A fixed point of Eq. (3) is unstable if one Delabays et al. (2016) or more Taylor (2012) angle differences are larger than $\pi/2$. In this case, a fixed point has $n-j$ angle differences $\Delta_{k+1,k}\equiv\Delta\in[-\pi/2,\pi/2]$ and $j$ angle differences $\Delta_{k+1,k}\equiv\pm\pi-\Delta$, with $j>0$, and $\sum_{k}\Delta_{k+1,k}=2\pi q$.

II.1 Identical frequencies: Analytical approach

For a cycle network of length $n$ with identical frequencies, we derive an analytical expression for the volume of the basins of attraction. Our approach is to approximate the basin of attraction of a given stable fixed point by the hypercube centered at the fixed point and whose radius is the distance to the closest unstable fixed point.

As Eq. (3) is invariant under a constant shift of all angles $\theta_{0}^{(q)}$, we will work in the hyperplane ${\cal H}_{n-1}$ orthogonal to the vector $(1,...,1)$. The angle vector $\vec{\theta}^{(q)}$ of Eq. (6) projected on ${\cal H}_{n-1}$ has components

[TABLE]

According to Ref. DeVille, 2012, a fixed point of Eq. (3) has a unique unstable direction in angle space if and only if it has a single angle difference between neighboring oscillators which is larger than $\pi/2$. Such an unstable fixed point is called a 1-saddle point. Consider then a 1-saddle point with winding number $q^{\prime}$, where the $k^{\text{th}}$ angle difference is larger than $\pi/2$. Combining Eqs. (4) and (5), its winding number $q^{\prime}$ is given by

[TABLE]

This allows to compute the components of the 1-saddle angle vector, $\vec{\varphi}^{(q^{\prime})}$, projected on ${\cal H}_{n-1}$ (see Appendix A),

[TABLE]

where

[TABLE]

We have found numerically that stable fixed points and 1-saddles are closest when they have the same winding number (see Appendix C). We thus investigate the case $q=q^{\prime}$. The difference between angles is easily obtained as

[TABLE]

which gives the distance between the stable fixed point and the unstable 1-saddle point

[TABLE]

In particular, the large $n$ limit is

[TABLE]

We remark that we lost the dependence on $k$, meaning that the stable fixed point $\vec{\theta}^{(q)}$ is equidistant to all 1-saddles with the same winding number. This indicates that the basins of attraction are isotropic in the directions of the 1-saddle points.

Since $\vec{\theta}^{(q)}$ is equidistant to all 1-saddles, one expect, for $n\gg 1$, that the volume $V_{q}$ of the basin of attraction is well approximated by an hypercube of side $\|\vec{\theta}^{(q)}-\vec{\varphi}^{(q)}\|_{\infty}$, up to a constant factor, i.e.

[TABLE]

As $q_{\max}={\rm Int}(n/4)$, $V_{q}\to e^{-4q}$ in the limit $n\to\infty$ for fixed $q$.

We further found numerically that 1-saddles are closer to the stable fixed points than $p$-saddles, $p>1$ (see Appendix C). This suggests that Eq. (18) underestimates the volume of basins of attraction.

Remark**.**

In our convention, angle differences are taken in the interval $(-\pi,\pi]$. Trying to construct a 1-saddle with $q^{\prime}=0$, one obtains

[TABLE]

The angle vector obtained in this way has winding number $q^{\prime}=-1$. This means that there is no 1-saddle with winding number zero. Eq. (15) then applies to $q>0$. It can be checked that there are no $p$-saddles with winding number $q=0$ for any $p\geq 1$.

II.2 Identical frequencies: Numerical approach

To validate the scaling of Eq. (18), we numerically estimate the volume of the basin of attraction of each stable fixed point of Eq. (3). For a cycle network of length $n$ with identical frequencies, all stable fixed points are known and given by Eq. (7). To estimate the volume of the basin of attraction of each $\vec{\theta}^{(q)}$, we randomly choose $d$ normalized perturbation vectors $\vec{\epsilon}_{j}\in{\cal H}_{n-1}\subset\mathbb{R}^{n}$, for $j=1,...,d$, $\|\vec{\epsilon}_{j}\|_{\infty}=1$. We then consider perturbed states

[TABLE]

as initial conditions for the dynamics of Eq. (3), with $\alpha\geq 0$. The parameter $\alpha$ is increased from zero to $\alpha_{q,j}$ which we define as the largest value such that $\vec{\eta}_{q,j,\alpha}$ converges back to $\vec{\theta}^{(q)}$ under the dynamics of Eq. (3). The distance between the stable fixed point $\vec{\theta}^{(q)}$ and the boundary of its basin of attraction in the direction $\vec{\epsilon}_{j}$ is given by $\pi\alpha_{q,j}$.

We performed $4^{\text{th}}$-order Runge-Kutta simulations of the dynamics of Eq. (3) for $n=23,43,83,163,323$. These values are chosen to maximize the volume of the basin of attraction for the largest winding number $q_{\max}={\rm Int}(n/4)$, which, as mentioned above, vanishes when $n$ is a multiple of $4$. Manik et al. (2016) We took $d=1000$ randomly chosen perturbation directions and increased $\alpha$ by steps of $0.01$. For each $q$, we can then estimate the proportion of the hypercube of side $\alpha$ centered at $\vec{\theta}^{(q)}$ which belongs to its basin of attraction as

[TABLE]

where ${\rm Card}$ stands for the cardinality of the ensemble. In Fig. 2 we see that this proportion stays close to 1 for small values of $\alpha$ and quickly drops to zero around some $q$-dependent value of $\alpha$.

Given a threshold $\tau\in[0,1]$ we can then define

[TABLE]

as a typical linear size of the basin of attraction. The abrupt drop of the curves in Fig. 2 implies that the precise value of $\tau$ is not too significant to understand the behavior of $\alpha_{q,j}$ with respect to $q$, provided that $\tau$ is neither too close to 1, nor to 0. We arbitrarily chose $\tau=0.7$, but checked that similar conclusions follow for $\tau=0.6$ and $0.8$.

In Fig. 3, we plot $\alpha_{\tau}(q)$ for various system sizes.

Except a saturation for small $q$’s, we observe a linear behavior of $\alpha_{\tau}$ with respect to $q$. Furthermore, curves for different values of $n$ varying by more than one order of magnitude are rescaled almost on top of one another when plotting them against $q/q_{\max}$. Both findings corroborate Eq. (16). Fig. 4, shows for each $q$, the quartiles of the values of $\alpha_{q,j}$ [defined by $\alpha_{\tau}(q)$ for $\tau=0.25,0.5,0.75$] and the extreme values $\min_{j}\alpha_{q,j}$ and $\max_{j}\alpha_{q,j}$, for $n=323$, as well as the distance between the stable fixed point $\vec{\theta}^{(q)}$ and the 1-saddle $\vec{\varphi}^{(q)}$ given by Eq. (16) (dashed line).

All curves have linear behavior, except for small $q$. The discrepancy between numerics and Eq. (16) comes from the fact that the random perturbations are not aligned with the direction of shortest distance to a 1-saddle. In other words, Eq. (16) is a lower bound on the distance between the stable fixed point $\vec{\theta}^{(q)}$ and the boundary of its basin of attraction. As mentioned above, 1-saddles are closer to the stable fixed points than other saddle points, which are also on the boundary of the basins of attraction. Our method of using the 1-saddles to evaluate the volume of the basins of attraction underestimates it, but clearly gives its right parametric dependence in $n$ and $q$.

II.3 The size of the sync basin revisited

The scaling obtained in Eq. (18), for large values of $n$, is different from the Gaussian scaling postulated in Ref. Wiley et al., 2006. The numerical method used there took initial conditions at random in the angle space $(-\pi,\pi]^{n}$, which would need a huge number of runs to reach a resolution allowing a fair estimate of the volumes of the basins of attraction. Even for a moderate resolution of $0.5$ in each angle direction, one would need approximately $(2\pi/0.5)^{83}\approx 10^{91}$ different initial condition, which is obviously unfeasible numerically. Hence, estimates based on brute-force numerical methods cannot catch the scaling behavior in dynamical systems with large dimensionality, especially for large winding numbers, which have very small basins of attractions.

Our approach overcomes this difficulty. Taking advantage of our knowledge of the stable fixed points, we are able to restrict the exploration of the basins of attraction to the neighborhood of the stable fixed points. We avoid scanning the whole angle space which significantly reduces the computation time and increase the accuracy of the method.

II.4 Non-identical frequencies

We introduced our method in the simplest case of a cycle network with identical frequencies. To generalize our understanding of the problem, we now add non-identical frequencies to the same cycle network. Even if we cannot obtain the stable fixed points analytically, we can find them numerically and then apply the same numerical procedure as in the identical frequency case. Instead of Eq. (3), our single-cycle model is now defined by

[TABLE]

with $P_{i}$ randomly and homogeneously taken in $[-\beta,\beta]$, satisfying $\sum_{i}P_{i}=0$. For small values of $\beta$, the non-identical frequencies almost always lead to small variations of the fixed points, Gilmore (1981) and thus the volume of the basins of attraction should not change much. To find the stable fixed points of Eq. (23), we start with the fixed points for $\beta=0$ given in Eq. (7) and follow them with a $4^{\rm th}$-order Runge-Kutta implementation of Eq. (23), while gradually increasing $\beta$ to the desired value. This allows to identify and follow numerically the location of the stable fixed point $\vec{\theta}^{(q)}(\{P_{i}\})$, which is not anymore given by Eq. (7), but is still characterized by its winding number $q$. We then perturb this stable fixed point in $1000$ random directions with increasing magnitude as in Eq. (20) and apply the same procedure as in Sec. II.2 to evaluate the volume of the basins of attraction.

Results are shown in Fig. 5 for $n=83$ and $\beta=0,0.01,0.02,0.05,0.1$.

For values of $q$ which are neither too small nor too large, the linear behavior of $\alpha_{\tau}$ is preserved, especially for small $\beta$. As can be expected, Coletta et al. (2016) as soon as we add some finite natural frequencies, the fixed points with large $q$ lose stability. More surprising, at first glance, is the abrupt drop of $\alpha_{\tau}$ for large $q$ with little change at small $q$. We offer an explanation for this behavior.

The dynamics of Eq. (1) is given by the gradient of the Lyapunov function

[TABLE]

Increasing $\beta$ modifies $\mathcal{V}$ [the first term on the right-hand side of Eq. (24)] and makes the fixed points move in angle space. Eventually, a stable fixed point $\vec{\theta}^{(q)}$ will meet an unstable fixed point and then lose stability through a saddle-node bifurcation. In Fig. 6, we give a schematic illustration of $\mathcal{V}$ on a cycle, projected on an appropriate direction in angle space, such that stable and unstable fixed points are aligned in one angle dimension. As long as a fixed point remains stable, the volume of its basin of attraction does not change much [compare the green segments in Figs. 6(a) and 6(b)]. The fixed point then abruptly vanishes when $\beta$ becomes too large. Since $\alpha_{\tau}$ is an average over many randomly chosen directions, its value abruptly drops when the stable fixed point vanishes.

In Appendix B, we furthermore estimate the maximal winding number possible for a given width $\beta$ of frequencies distribution and find that it agrees qualitatively with the numerically observed maximal winding numbers.

III Meshed Networks

We finally extend the perturbation method described above to more complicated, meshed networks. It is a two-stage method where we first numerically identify fixed points of Eq. (1) on complex graphs and second perturb the obtained stable fixed points in the same way as in Sec. II.

III.1 Identifying stable fixed points

Stable fixed points are much harder to find on complex graphs. Except for the $\vec{\theta}^{(0)}=(0,...,0)$ fixed point for equal frequencies, they are usually impossible to find analytically. To tackle this problem, we construct a numerical algorithm similar to but different from the one proposed in Ref. Manik et al., 2016.

From Refs. Delabays et al., 2016; Dörfler et al., 2013, we know that two fixed points of Eq. (1) differ only by a collection of loop flows quantized by their winding numbers (similar to vortices in superconductors). We define a vector composed of the winding numbers on each cycle of a graph ${{G}}$ as

[TABLE]

where $m$ is the total number of cycles in ${{G}}$. Assuming $|\theta_{i}-\theta_{j}|<\pi/2$ for all connected nodes $i$, $j$, each stable fixed point $\vec{\theta}^{*}$ can be uniquely labelled by its winding vector $\vec{q}_{{{G}}}(\vec{\theta}^{*})$. Delabays et al. (2016); Dörfler et al. (2013); Manik et al. (2016) For $P_{i}\equiv 0$, the Lyapunov function, Eq. (24), reduces to

[TABLE]

which is the Hamiltonian of a XY model, describing the interaction of planar classical spins.Fazio and van der Zant (2001) Stable fixed points are the local minima of this energy function and it is known that they correspond to vortex-carrying states, i.e. states with non-zero winding vector $\vec{q}_{G}(\vec{\theta})$. We therefore search for stable fixed points via an iterative process starting from vortex-carrying initial states described by

[TABLE]

where $(x_{0},y_{0})$ are the coordinates of the center of the vortex with charge/vorticity $q\in\mathbb{Z}$, and $(x_{i},y_{i})$ the coordinates of the position of $\theta_{i}$. To find stable fixed points, the algorithm reads:

Define a two-dimensional embedding of the network. Use this to superimpose a regular lattice of coordinates on the network. This is shown in Fig. 7. 2. 2.

Set $(x_{0},y_{0})$ to a node of the regular lattice.

1. (a)

Using Eq. (28), define a new initial state. 2. (b)

Follow numerically Eq. (1) on the considered meshed graph until a stable fixed point is reached. 3. (c)

Each stable fixed point can be unambiguously identified by its winding vector. Use this to determine if the fixed point just found is a new one. If yes, store it. 3. 3.

Go back to step 2.

As complex meshed network, we consider the UK high voltage grid which is composed of 120 nodes and 165 lines. To illustrate the algorithm just described, Figs. 7 and 7 show an initial condition and the stable state toward which it dynamically converges respectively. Stable fixed points with many vortices are obtained by setting $|q|$ to large values, in our case $q\in\{-50,...,50\}$. The dynamics will then split this initial vortex into several vortices with smaller $q$’s, located on different cycles of the network. This method can be used on any network whether complex or regular. Time evolving Eq. (1) on the UK grid with this initial condition returns only stable fixed points. In this way we found more than 4000 different stable fixed points of Eq. (1) with $P_{i}=0$.

III.2 Estimating the volume of basins of attraction

Having identified stable fixed points $\vec{\theta}^{(\vec{q})}$ of Eq. (1) on the UK grid, we next follow the same procedure as in Sec. II.2, and measure the volume of their basins of attraction.

We focus on stable fixed points with non-zero winding number only on the five cycles in red in Fig. 8. We introduce a shorthand notation with the winding numbers of these cycles only

[TABLE]

Taking each cycle independently, the maximum winding numbers are ${q}_{1}^{\rm max}=4$, ${q}_{2}^{\rm max}=2$, ${q}_{3}^{\rm max}=1$, ${q}_{4}^{\rm max}=2$, ${q}_{5}^{\rm max}=2$.

In Fig. 9, we show $\alpha_{\tau}$ of Eq. (22) for $\tau=0.5$ and the interquartile values for various stable fixed points identified by their unique combination of winding numbers. For $P_{i}=0$, Eq. (1) is symmetric under $\theta_{i}\rightarrow-\theta_{i}$, $\forall i$. This implies that the volume of each basin of attraction is invariant under $\vec{q}_{G}\rightarrow-\vec{q}_{G}$. This symmetry can be seen on Fig. 9(a) which is symmetric under $q_{1}\rightarrow-q_{1}$, as there is only one cycle with a non-zero winding number. Alternatively, $\vec{q}_{G}\to-\vec{q}_{G}$ interchanges panels (b) and (c) with $q_{1}\to-q_{1}$. For fixed points with more than one cycle carrying a non-zero winding number, there is no symmetry for $q_{1}\rightarrow-q_{1}$ [see Figs. 9(b)–(f)] and the asymmetry is even more significant if the vortices are close to each other [see Figs. 9(b) and 9(d)], because there are few intermediate nodes which can screen the effect of one vortex on the other and so the cycles interact strongly.

In Fig. 9(d), we did not find any fixed point with winding vector $\vec{q}_{\rm sh}=(1,0,1,0,0)$ but we found one with $\vec{q}_{\rm sh}=(-1,0,1,0,0)$. To understand this, we consider the simplified situation of two connected cycles of size $m$ and $n$, sharing $\ell$ edges as depicted on Fig. 10. We consider the case when $q^{(n)}=1$ and $q^{(m)}=0$. At a fixed point, Eq. (1) implies

[TABLE]

and Eq. (5) gives

[TABLE]

From Ref. Delabays et al., 2017, if $\ell>1$, we have $|\Delta|$,$|\Delta^{\prime}|$,$|\Delta^{\prime\prime}|\leq\pi/2$. Eq. (30) and Eq. (31) imply that $\Delta>\Delta^{\prime}$. Therefore, to have $q^{(n)}=1$, we must have $\Delta>{2\pi}/{n}$. Thus, if we add edge-sharing cycles with zero winding number to a main cycle with a non-zero winding number, some of the angle differences must increase. When $\Delta$ is large, this can bring $\Delta>{\pi}/{2}$ where stability is lost.Delabays et al. (2017) Adding a cycle carrying a non-vanishing winding number makes the situation even more critical. If we isolate cycles 1 and 3, which correspond to Fig. 10 with $n=16$, $m=7$ and $\ell=2$, an easy calculation shows that there exist stable fixed points with $({q^{(16)}},q^{(7)})=(1,1)$ and $({q^{(16)}},q^{(7)})=(1,-1)$. However, when we consider the complete network, only the solution with $({q^{(16)}},q^{(7)})=(1,-1)$ remains stable. This comes from the fact that, when both winding numbers have the same signs, the angle differences on the shared edges benefit only to one of the cycles. The other cycle then has to make a winding number out of a reduced number of edges, implying larger angle differences. Finally, when we take the complete network, we put cycles next to the two initial ones and make the angle differences even larger, until stability is lost. When the winding numbers have opposite signs, both cycles benefit from the angle differences on the shared edges, which leads to smaller angle differences than in the previous case. This explains why the fixed point $\vec{q}_{\rm sh}=(-1,0,1,0,0)$ is stable while $\vec{q}_{\rm sh}=(1,0,1,0,0)$ is not.

To conclude this section, we note that a meshed network has an effect similar to the case considered in Sec. II.4 with $P_{i}\neq 0$ in that, compared to the single-cycle network, (i) there are fewer stable fixed points with large winding numbers and (ii) the volume of basins of attraction of fixed points with small winding numbers seem to be unaffected.

IV Conclusion

We have developed a numerical method to investigate the volume of basins of attraction of fixed points in dynamical systems. Our method first locates the stable fixed points of the dynamical system, using an algorithm based on the concept of loop flows. Delabays et al. (2016) Second, it pertubs them in random directions with increasing magnitude. The proportion of pertubed states that converge back to the initial fixed point allows to evaluate the radius of the basin of attraction and then its volume.

We then used our method to investigate the Kuramoto model on a cycle with identical frequencies. We obtained that the volume of the basin of attraction is proportional to $(1-4q/n)^{n}$, contrasting with the Gaussian distribution suggested in Ref. Wiley et al., 2006. We then extended the application of our method to the Kuramoto model on a cycle with non-identical frequencies and to the Kuramoto model on meshed networks. These two generalizations render the investigations of the basins of attraction much less tractable, which imposes to rely on numerics. We believe that our method significantly speeds up these investigations.

Compared to other existing methods to investigate basins of attraction, our method has three main advantages:

• •

It does not require a Lyapunov function of the dynamical system considered, Genesio et al. (1985); Najafi et al. (2016) which is complicated to find in general;

• •

It is not limited to quadratic or polynomial systems; Amato et al. (2007)

• •

The investigation is guided by our knowledge of the system and avoids to randomly pick initial conditions in the state space. Wiley et al. (2006); Menck et al. (2013)

These advantages come with the drawback that we limit our investigations to the volume of the basins of attraction and have no indications about their shape. In particular, our method is probably not adapted to the investigation of fractal basins of attraction. Grebogi et al. (1987) We think that our method may be useful in many other contexts including finding local energy minima in planar spin glasses Mézard et al. (1987) and disordered Josephson junction arrays Fazio and van der Zant (2001) among others.

Acknowledgment

This work has been supported by the Swiss National Science Foundation under an AP Energy Grant. We thank Jimmy Dubuisson for discussions at the early stage of this project.

Appendix A Angle vector of the 1-saddle points

We give here details of the computation of the components of the angle vector $\vec{\varphi}^{(q^{\prime})}\in{\cal H}_{n-1}$ defined in Sec. II.1. The vector $\vec{\varphi}^{(q^{\prime})}$ is the 1-saddles on a cycle of length $n$, with winding number $q^{\prime}$ and the only angle difference exceeding $\pi/2$ located on the edge between vertices $k-1$ and $k$. Its components are given by

[TABLE]

where $S_{k}$ is a constant angle shift guaranteeing that the sum of components is zero,

[TABLE]

Some algebra gives

[TABLE]

where

[TABLE]

Appendix B Maximal winding number on a cycle with random natural frequencies

According to Ref. Delabays et al., 2016, the angle difference on the edge between vertices $k$ and $k+1$ is given by

[TABLE]

where $P_{k,k+1}^{*}\coloneqq\sum_{j=1}^{k}P_{j}$ is a reference flow and $\varepsilon_{q}$ is the loop flow parameter determining the winding number of the fixed point. For the sake of simplicity we take $K=1$. The frequencies $P_{k}$ are taken randomly and homogeneously in the interval $[-\beta,\beta]$. Their expectation and variance are

[TABLE]

Expectation and variance for the $P_{k,k+1}^{*}$ are then

[TABLE]

For $n$ sufficiently large, we then expect typical excursions of magnitude $\sqrt{k}\beta/\sqrt{3}$ of $P_{k,k+1}^{*}$ away from its average $\mathbb{E}=0$. It is known Delabays et al. (2016) that on cycles with finite natural frequencies, stable fixed points may have one angle diffence slightly larger than $\pi/2$ before losing stability at $\pi/2+\delta$. As $\delta$ is always small, we will approximate the loss of stability to happen when $\Delta_{k,k+1}=\pi/2$, i.e. when the argument of the arcsine in Eq. (B.1) is equal to one. Finally, we approximate $\varepsilon_{q}$ by its value when $\beta=0$,

[TABLE]

Putting everything together and taking $k$ to be the average value of the indices, i.e. $k=n/2$, Eq. (B.1) gives

[TABLE]

which gives a maximal possible value of $q$ before losing stability, with respect to $\beta$,

[TABLE]

Simulated and estimated values of $q_{\max}$ are given in Table 1 for various values of $\beta$. Even if the simplifications assumed to obtain Eq. (B.6) underestimates $q_{\max}$, it is in fair agreement with numerical obtained values.

Appendix C Distance between stable fixed points and saddle points

We justify numerically the two statements of Sec. II.1, that the 1-saddles are the closest unstable fixed points to the stable fixed points, and that stable fixed points and 1-saddles are the closest if they have the same winding number. We write $\vec{\varphi}^{(q^{\prime},\ell)}_{i}$ for the unstable fixed point with winding number $q^{\prime}$ and $\ell$ angle differences larger than $\pi/2$, where the index $i$ labels the different fixed points with same $q^{\prime}$ and $\ell$.

Fig. 11 shows that the 1-saddles are the closest unstable fixed points to stable fixed points. Fig. 12 shows that stable fixed points and 1-saddles are the closest if they have the same winding number. As we remarked in Sec. II.1, the case of $q=0$ is special because there are no unstable fixed points $\vec{\varphi}^{(0,\ell\neq 0)}_{i}$, with winding number zero.