Basin sizes depend on stable eigenvalues in the Kuramoto model

Following Wiley2006 , we consider here a system of $N$ identical Kuramoto oscillators on a regular ring where each oscillator is coupled with equal strength to its $R$ nearest neighbors on either side. In the co-rotating reference frame the time evolution of this system is governed by the following set of ODEs

$$\dot{\theta_{j}}=\sum_{k=j-R}^{j+R}G\,\sin(\theta_{k}-\theta_{j})\,,\,\,j=1,2,...,N;$$

(1)

where the index $k$ is periodic mod $N$, the coupling constant $G$ is positive and can be completely removed from the equations by rescaling the time.

As pointed out in Wiley2006 , the set of equations (1) is a gradient system that can be recast as $\dot{\theta}=-\nabla V$ with e.g. $V(\theta_{1},\ldots,\theta_{N})=G\sum_{i,j}\cos(\theta_{i}-\theta_{j})$ so that $V$ is bounded both from below and from above: $-N^{2}G\leq V\leq N^{2}G$. Therefore $dV/dt=-(\nabla V)^{2}$ and all trajectories, except the points of equilibrium and their stable manifolds, are flowing “monotonically downhill” and asymptotically tend to those of the equilibria that correspond to the local minima of $V$. For this reason, “we need not concern ourselves with the possibility of more complicated long-term behavior, such as limit cycles, attracting tori, or strange attractors” for (1) Wiley2006 .

The system (1) has a family of equilibrium states which can be characterized by an integer $q$

$$\theta_{j}=\frac{2\pi q}{N}j+C\,,$$

(2)

the “winding number”, which measures the number of full twists in phase as one goes around the ring once and can assume the values $q=0,1,2,...,N-1$; and $C$ is a real constant. The state with $q=0$ corresponds to $\theta_{j}=C,\forall j$, i.e. with all oscillators synchronized in phase.

In the states with $q\neq 0$ all oscillators are synchronized in frequency but with a constant phase difference between two successive units: $\Delta_{j}\equiv\theta_{j}-\theta_{j-1}=2\pi q/N,\forall\,j$. Such states are also known as “$q$-twisted” states. Recalling that $\Delta_{j}$ is taken mod $2\pi$, we notice that for a state with $q=N-n,\,1\leq n\leq N/2$, the phase difference between two successive oscillators is

$$\Delta_{j}=\frac{2\pi}{N}(N-n)=2\pi-\frac{2\pi}{N}n=-\frac{2\pi}{N}n<0\,,$$

so that the phases of oscillators are distributed (in the order: $\theta_{1},\theta_{2},...$) clockwise around the circle, in such a way that the winding number of this state can also be represented by $q=-n$, denoting $n$ clockwise twists. So one can change the range of winding numbers to $q=-m,...,-1,0,1,...,m$ with $m=N/2$ (for odd $N$, $m=(N-1)/2$).

In turn, the eigenvalues of the Jacobian matrices near those equilibrium states (parametrized by $q$) are real and obey the expression  Mihara2019

	$\displaystyle\gamma_{\ell}(q)=-4\,G\,\,\sum_{k=1}^{R}\cos\left(k\frac{2\pi}{N}q\right)\sin^{2}\left(k\frac{\pi}{N}\ell\right)\,,$
	$\displaystyle\ell=1,2,...,N-1.$		(3)

Clearly, the linear stability of these equilibria is determined by the eigenvalues, and it was observed that the stability of each $q$-twisted state depends on the ratio $R/N$ ($\propto$ connectivity) : (i) for $R/N\gtrsim 0.34$, $q=0$ (all oscillators in phase) is the only stable equilibrium state; (ii) below this threshold, and as the ratio decreases, more twisted states (with $q\neq 0$) become stable Wiley2006 ; Maistrenko2012 ; Mihara2019 . As mentioned before, the authors of Ref. Wiley2006 concluded that the distribution of volumes of basins of attraction of the $q$-twisted states has Gaussian shape with standard deviation $\sigma\sim 0.2\sqrt{N/R}$.

The Koopman (or composition) operator approach can provide a global description of dynamical systems in terms of the time evolution of observables (functions) of the state space. In this approach a nonlinear dynamical system is represented in terms of an infinite–dimensional (but linear) operator acting on a Hilbert space of functions of the system states. The spectral decomposition of the Koopman operator provides complete description of the nonlinear system. Despite being an infinite–dimensional operator, there are several numerical methods (such as DMD and its variations) capable of obtaining finite–dimensional approximations for Koopman eigenvalues / modes and they have applications in various real-world problems such as fluids dynamics, power grids, epidemiology, climatology, etc. , see Ref.koopmanbook2020 and references therein. In turn the Perron-Frobenius (or transfer) operator evolves densities of trajectories in the state space. It is also linear and is dual to the Koopman operator. Both operators share the same spectral properties and can provide global descriptions of a dynamical system koopmanbook2020 .

Recently it was demonstrated that for a Morse-Smale gradient flow acting on a smooth, compact and oriented manifold with no boundary, the spectrum of the transfer operator is given by linear combinations of the Lyapunov exponents at the critical points of the Morse function (i.e., the eigenvalues of Jacobian at the fixed points) DangRiviere2019 and it holds globally on the manifold. This result agrees with the observation that for a $d$-dimensional autonomous system with a hyperbolic fixed point $x^{*}\in X$, where $X$ is a compact, connected and forward–invariant subset of $\mathds{R}^{d}$, the spectrum of the Koopman operator is given by the eigenvalues of the Jacobian matrix evaluated at $x^{*}$ mauroy2016 and also with a relevant property of this operator: if $\phi_{1}$ and $\phi_{2}$ are Koopman eigenfunctions associated with the eigenvalues $\mu_{1}$ and $\mu_{2}$, then $\Phi=\phi_{1}^{a}\phi_{2}^{b}$, with $a,b\in\mathds{R}$, is also a Koopman eigenfunction with eigenvalue $a\mu_{1}+b\mu_{2}$.

On the other hand one should notice that the set of Kuramoto equations (1) above is a Morse-Smale system with Morse function given by the “potential” $V(\theta_{1},...,\theta_{N})$, the critical points of $V$ are the $q$–twisted states and the Lyapunov exponents (at critical points) are the eigenvalues $\gamma_{\ell}(q)$. Then the mathematical results DangRiviere2019 ; mauroy2016 above ensure that the $\gamma_{\ell}(q)$, despite being obtained by local methods, somehow contain global information of system (1), and we shall use them to explore the basins of attraction of the $q$-twisted states. In order to represent the stability in all directions of the phase space, we consider the sum of the eigenvalues of a $q$-twisted state,

$$\hat{\gamma}_{q}\equiv\sum_{\ell=1}^{N-1}\,\gamma_{\ell}(q)$$

(4)

which resembles the entropy functional for Morse-Smale diffeomorphisms in the framework of (Gibbs) variational principle for dynamical systems takahashi84 .

Then we define $\Gamma_{q}$, the equilibrium stability measure, as the sum of the eigenvalues of a $q$-twisted state, normalized by the most negative $\hat{\gamma}_{q}$ which in this case is $\hat{\gamma}_{0}$:

$$\Gamma_{q}\equiv\frac{\hat{\gamma}_{q}}{\hat{\gamma}_{0}}\,.$$

(5)

We present in Fig. 1(a) the plot of $\Gamma_{q}$ with respect to $q$ for different network sizes with local coupling ($R=1$). (hyperbolic) stable states, i.e. those with $\gamma_{\ell}<0,\forall\,\ell$, were taken into account. It is remarkable the similarity between the plot of $\Gamma_{q}$ and the plot of $\alpha_{\tau}(q)$, the typical linear size of the basin of attraction presented in Fig. 3 (inset) of Ref. Delabays2017 . Therefore, if $\Gamma_{q}$ has a behavior similar to the linear size of a basin of attraction, it is reasonable to expect that for a network with $N$ oscillators the volume of the basin of attraction of a given $q$-state will behave as $\sim\Gamma_{q}^{N}$. In Fig. 1(b) we plot $\Gamma_{q}^{N}$, for $N=60$ and $R=2$; that can be well approximated by a Gaussian curve. To establish this result, in the next subsection we show explicitly for low values of $R/N$ that $\Gamma_{q}^{N}$ can be approximated by a Gaussian function with respect to the winding number $q$, .

Our goal here is to derive an explicit expression for the standard deviation $\sigma$ in terms of the number of nodes $N$ and the connection $R$ of the network. Considering $G=1/N$, Eq. (3) returns

$$\gamma_{\ell}(q)=-4\,\sum_{k=1}^{R}\,\frac{1}{N}\,\cos(k\,q\,\delta)\,\sin^{2}(k\,\ell\,\delta/2)\,,$$

(6)

where $\ell=1,2,\dots,N-1$ and $\delta=2\pi/N$.

By using trigonometric identities, Eq. (6) can be rewritten as

	$\displaystyle\gamma_{\ell}(q)=$	$\displaystyle-\frac{1}{N}\,\sum_{k=1}^{R}\,\{2\cos(k\,q\,\delta)\,-\cos(k\,\delta\,(q+\ell))$
		$\displaystyle-\cos(k\,\delta\,(q-\ell))\}\,,$		(7)

and on performing the summations one obtains

	$\displaystyle\gamma_{\ell}(q)=$	$\displaystyle-\frac{1}{N}\,\left\{\frac{\sin\left(M\,q\,\pi/N\right)}{\sin\left(q\,\pi/N\right)}-\frac{1}{2}\left[\frac{\sin\left(M\,(q+\ell)\,\pi/N\right)}{\sin\left((q+\ell)\,\pi/N\right)}\right.\right.$
		$\displaystyle\left.\left.+\frac{\sin\left(M\,(q-\ell)\,\pi/N\right)}{\sin\left((q-\ell)\,\pi/N\right)}\right]\right\}\,,$		(8)

where $M\equiv 2R+1$.

Now let us evaluate $\hat{\gamma}_{q}=\sum_{\ell}\gamma_{\ell}(q)$. The summation of the first term in the RHS of Eq. (8) is trivial, but for the second and the third terms, the sums are approximated by integrals since we are regarding $N\gg R,\,q$ :

	$\displaystyle B_{\pm}=$	$\displaystyle\sum_{\ell=1}^{N-1}\,\frac{1}{N}\,\frac{\sin\left(M\,(q\pm\ell)\,\pi/N\right)}{\sin\left((q\pm\ell)\,\pi/N\right)}$
	$\displaystyle\approx$	$\displaystyle\int_{0}^{1-\varepsilon}\frac{\sin\left(M\,(y\pm x)\,\pi\right)}{\sin\left((y\pm x)\,\pi\right)}dx\,,$		(9)

where $\varepsilon=1/N$, $x=\ell/N$, $y=q/N$ and $M=2R+1$. Then one arrives at

$$\hat{\gamma}_{q}=1-\frac{\sin\left(M\,y\,\pi\right)}{\sin\left(y\,\pi\right)}+\mathcal{O}(\varepsilon^{2})\,.$$

(10)

In turn, the normalized sum of eigenvalues $\Gamma_{q}$, can be expanded around $y\sim 0$

	$\displaystyle\Gamma_{q}$	$\displaystyle=$	$\displaystyle 1-\frac{\pi^{2}}{3}(2R+1)(R+1)\,y^{2}+\dots\,$		(11)
		$\displaystyle\approx$	$\displaystyle 1-\beta\frac{q^{2}}{N}\,,$

where

$$\beta=\frac{\pi^{2}}{3}\frac{(2R+1)(R+1)}{N}\,.$$

(12)

For $N\gg R,q$, the quantity $\Gamma_{q}^{N}$ can be approximated by a Gaussian function

$$\Gamma_{q}^{N}\approx\left(1-\beta\frac{q^{2}}{N}\right)^{N}\approx e^{\displaystyle-\beta q^{2}}\,,$$

(13)

with standard deviation equal to $1/{\sqrt{2\beta}}$. As a result, the standard deviation can be written as

$$\sigma_{T}=\sqrt{N}\,\mathcal{F}(R),\quad\mathcal{F}(R)=\left[\frac{3}{2\pi^{2}(2R+1)(R+1)}\right]^{1/2}.$$

(14)