Chaos in networks of coupled oscillators with multimodal natural frequency distributions

TL;DR

This paper investigates the emergence of chaos in networks of coupled oscillators modeled by the Kuramoto model with multimodal natural frequency distributions, identifying conditions under which chaos occurs or is prevented.

Contribution

It introduces a collective coordinate reduction approach to analyze chaos in multimodal Kuramoto models and clarifies the role of cluster synchronization and desynchronization.

Findings

Chaos requires at least four clusters for three degrees of freedom.

Intermittent desynchronization can induce chaos in trimodal distributions.

Chaos cannot occur in bimodal distributions, even with asymmetry.

Abstract

We explore chaos in the Kuramoto model with multimodal distributions of the natural frequencies of oscillators and provide a comprehensive description under what conditions chaos occurs. For a natural frequency distribution with $M$ peaks it is typical that there is a range of coupling strengths such that oscillators belonging to each peak form a synchronized cluster, but the clusters do not globally synchronize. We use collective coordinates to describe the inter- and intra-cluster dynamics, which reduces the Kuramoto model to $2M-1$ degrees of freedom. We show that under some assumptions, there is a time-scale splitting between the slow intracluster dynamics and fast intercluster dynamics, which reduces the collective coordinate model to an $M-1$ degree of freedom rescaled Kuramoto model. Therefore, four or more clusters are required to yield the three degrees of freedom necessary for chaos. However, the time-scale splitting breaks down if a cluster intermittently desynchronizes. We show that this intermittent desynchronization provides a mechanism for chaos for trimodal natural frequency distributions. In addition, we use collective coordinates to show analytically that chaos cannot occur for bimodal frequency distributions, even if they are asymmetric and if intermittent desynchronization occurs.

Equations

$$\dot{{\varphi }_i}={\omega }_i+\frac{K}{N}\sum _{j=1}^N{A}_{ij}\sin ({\varphi }_j-{\varphi }_i),$$

$$g(\omega )=\sum _{m=1}^M{\gamma }_m{g}_m(\omega ;{\Omega }_m,{\sigma }_m^2),$$

$$r(t)e^{i\psi (t)}=\frac{1}{N}\sum _{j=1}^N{e}^{i{\varphi }_j(t)},$$

$$\bar{r}=\frac{1}{T}{\int }_{t_0}^{t_0+T}r(t)dt,$$

$$r_m(t)e^{i{\psi }_m(t)}=\frac{1}{N_m}\sum _{j\in {\mathcal{C}}_m}{e}^{i{\varphi }_j(t)},$$

$$\mathbf{\varphi }(t)\approx \sum _{m=1}^M{\alpha }_m(t){\mathbf{\Phi }}^{(m)}+f_m(t)1_{N_m},$$

$$\dot{\mathbf{\varphi }}=\mathbf{\omega }-\frac{K}{N}L\mathbf{\varphi },$$

$$\hat{\mathbf{\varphi }}=\frac{N}{K}{L}^{+}\mathbf{\omega },$$

$$\mathbf{\Phi }=\alpha (t)\hat{\mathbf{\varphi }},$$

$$\hat{\mathbf{\varphi }}=\frac{1}{K}(\mathbf{\omega }-\Omega 1_N).$$

$${\mathcal{E}}_i=\dot{\alpha }{\hat{\varphi }}_i-{\omega }_i-\frac{K}{N}\sum _{j=1}^N{A}_{ij}\sin \left(\alpha \left({\hat{\varphi }}_j-{\hat{\varphi }}_i\right)\right).$$

$$\dot{\alpha }=\frac{{\hat{\mathbf{\varphi }}}^T\mathbf{\omega }}{{\hat{\mathbf{\varphi }}}^T\hat{\mathbf{\varphi }}}+\frac{1}{{\hat{\mathbf{\varphi }}}^T\hat{\mathbf{\varphi }}}\frac{K}{N}\sum _{i,j=1}^N\widehat{{\varphi }_i}{A}_{ij}\sin \left(\alpha \left({\hat{\varphi }}_j-{\hat{\varphi }}_i\right)\right).$$

$$\dot{\alpha }=K+\frac{K^2}{{\Sigma }^2{N}^2}\sum _{i,j=1}^N{\omega }_i\sin \left(\frac{\alpha }{K}\left({\omega }_j-{\omega }_i\right)\right),$$

$$\dot{\beta }=1+\frac{K}{{\Sigma }^2{N}^2}\sum _{i,j=1}^N{\omega }_i\sin \left(\beta \left({\omega }_j-{\omega }_i\right)\right).$$

$$\dot{\beta }=1+\frac{K}{{\Sigma }^2}\iint \omega \sin \left(\beta \left(\eta -\omega \right)\right)g(\omega )g(\eta )d\omega d\eta =\mathcal{I}(\beta ).$$

$$\mathcal{I}(\beta )=1-K\beta \exp \left(-{\sigma }^2{\beta }^2\right).$$

$$K\ge \sigma \sqrt{2e}.$$

$$r(t)=\exp \left(-\frac{{\sigma }^2{\beta }^2}{2}\right).$$

$$\mathcal{I}(\beta )=$$

$$\left[{\Omega }_j\sin \left(\beta \left({\Omega }_i-{\Omega }_j\right)\right)-\beta {\sigma }_j^2\cos \left(\beta \left({\Omega }_i-{\Omega }_j\right)\right)\right].$$

$$\underset{\beta }{\min }\mathcal{I}(\beta )<0.$$

$${\dot{\varphi }}_i^{(m)}={\omega }_i^{(m)}+\frac{K}{N}\sum _{j\in {\mathcal{C}}_m}{A}_{ij}\sin \left({\varphi }_j^{(m)}-{\varphi }_i^{(m)}\right),$$

$${\hat{\mathbf{\varphi }}}^{(m)}=\frac{N}{K}{L}_m^{+}{\mathbf{\omega }}^{(m)},$$

$${\hat{\mathbf{\varphi }}}^{(m)}=\frac{N}{K{N}_m}\left({\mathbf{\omega }}^{(m)}-{\Omega }_m{1}_{N_m}\right),$$

$${\mathbf{\Phi }}^{(m)}=\alpha {\pi }^{(m)}\hat{\mathbf{\varphi }}+{\alpha }_m{\hat{\mathbf{\varphi }}}^{(m)}+f_m{1}_{N_m},$$

$${\hat{\mathbf{\varphi }}}^T\hat{\mathbf{\varphi }}\dot{\alpha }+\sum _{m=1}^M({\pi }^{(m)}\hat{\mathbf{\varphi }}{)}^T\left({\hat{\mathbf{\varphi }}}^{(m)}{\dot{\alpha }}_m+1_{N_m}{\dot{f}}_m\right)={\hat{\mathbf{\varphi }}}^T\mathbf{G}(\mathbf{\Phi }),$$

$$\mathbf{G}(\mathbf{\Phi })=\mathbf{\omega }+{\left(\frac{K}{N}\sum _{j=1}^N\sin ({\Phi }_j-{\Phi }_i)\right)}_{i=1,\dots ,N}$$

$$({\hat{\mathbf{\varphi }}}^{(m)}{)}^T{\pi }^{(m)}\hat{\mathbf{\varphi }}\dot{\alpha }+({\hat{\mathbf{\varphi }}}^{(m)}{)}^T{\hat{\mathbf{\varphi }}}^{(m)}{\dot{\alpha }}_m=({\hat{\mathbf{\varphi }}}^{(m)}{)}^T{\pi }^{(m)}\mathbf{G}(\mathbf{\Phi }).$$

$$1_{N_m}^T{\pi }^{(m)}\hat{\mathbf{\varphi }}\dot{\alpha }+N_m{\dot{f}}_m=1_{N_m}^T{\pi }^{(m)}\mathbf{G}(\mathbf{\Phi }).$$

$$\mathcal{A}\dot{\mathbf{x}}=\mathbf{b}(\mathbf{x}),$$

$${\Phi }_i^{(m)}$$

Full text

Synchronization of coupled oscillators occurs in many natural processes and engineering applications. The dynamics of the globally synchronized state is regular and the phases typically rotate with a constant mean frequency. In the case of multimodal distributions of natural frequencies of the oscillators, one can observe more complex dynamics including chaos. Under which conditions the synchronized state may exhibit chaos has not been fully addressed. Distinct peaks in a multimodal natural frequency distribution correspond to synchronized clusters for a range of coupling strengths and network parameters. We study the intercluster and intracluster dynamics using a collective coordinate approach, which reduces the dimension of the full Kuramoto model to a small number of active degrees of freedom. We find necessary conditions for chaos to occur. In particular, at least four peaks in the natural frequency distribution are required to produce phase chaos, and chaos can also occur for three peaks via intermittent desynchronization of clusters.

I Introduction

Synchronization in networks of coupled oscillators occurs in many natural systems, including the activity of the brain Sheeba, Stefanovska, and McClintock (2008); Bhowmik and Shanahan (2012) and synchronous firefly flashing, Mirollo and Strogatz (1990) as well as many engineering applications, such as power grids, Filatrella, Nielsen, and Pedersen (2008) and Josephson junction arrays. Watanabe and Strogatz (1994); Wiesenfeld, Colet, and Strogatz (1998)

In typical models of synchronization, the dynamics is either incoherent, partially synchronized, or fully synchronized. In the case of a unimodal frequency distribution, the dynamics transitions upon increasing the coupling strength from the incoherent state at low coupling strength, to a partially synchronized state where a collection of oscillators synchronize (those with native frequency closest to the mean frequency), to the fully synchronized state at high coupling strengths. For multimodal frequency distributions, however, several synchronized clusters may emerge in the partially-synchronized regime. That is, there are clusters of oscillators that remain synchronized within themselves, but the oscillators do not form a single synchronized cluster. These clusters may have complex interactions, both inter-cluster and intra-cluster, producing complex dynamics, including chaos.

Chaos in coupled oscillator networks has been previously studied. For the Kuramoto model, Kuramoto (1984); Strogatz (2000); Pikovsky, Rosenblum, and Kurths (2001); Acebrón et al. (2005); Osipov, Kurths, and Zhou (2007); Arenas et al. (2008); Dörfler and Bullo (2014); Rodrigues et al. (2016) which is the model focused on here, chaos has been observed in the incoherent state, termed phase chaos, Popovych, Maistrenko, and Tass (2005); Maistrenko, Popovych, and Tass (2005); Miritello, Pluchino, and Rapisarda (2009) provided there are at least four oscillators. This type of phase chaos occurs at the microscopic level and is associated with the chaotic dynamics of individual phase oscillators. For such microscopic phase chaos, the Lyapunov exponent was found to scale inversely proportionally to the number of oscillators. Carlu, Ginelli, and Politi (2018) In particular, this implies that in the thermodynamic limit of infinitely many oscillators the Lyapunov exponent is zero, i.e., no chaos. Here we focus on collective chaotic behavior of synchronized subpopulations of phase oscillators. Such collective chaos has been studied for systems with symmetric bimodal natural frequency distributions which were subjected to a time-periodic coupling strength, So and Barreto (2011) or for different inter- and intra-cluster coupling strengths as well as a phase lag. Bick, Panaggio, and Martens (2018) However, for the classical Kuramoto model, it has been shown that in the thermodynamic limit with bimodal natural frequency distributions chaos is impossible. So and Barreto (2011) For trimodal frequency distributions, which yield three synchronized subpopulations, chaos has been observed for superposed Lorentzian natural frequency distributions, but only in the partially synchronized state, which involves microscopic chaos of incoherent oscillators. Cheng et al. (2017)

Here we present and analytically study generic situations of collective chaos in which the dynamics of synchronized subpopulations of coupled oscillators, termed clusters, can be chaotic. We distinguish between two types of chaotic dynamics, one akin to phase chaos and the other due to intermittent desynchronization. Here we refer to collective phase chaos when each of the synchronized clusters preserves their shape while the phases of the clusters show chaotic behavior. In this case, the possibility of chaos is determined by the number of synchronized clusters, which determines the number of active degrees of freedom. We shall see that to obtain phase chaos at least four synchronized clusters are necessary. This is analogous to needing at least four oscillators to generate microscopic phase chaos in the incoherent state of the Kuramoto model. Popovych, Maistrenko, and Tass (2005); Maistrenko, Popovych, and Tass (2005); Miritello, Pluchino, and Rapisarda (2009)

A different type of chaos is observed when clusters intermittently desynchronize through their mutual interactions. In this case, as we will show, chaos may occur even for trimodal natural frequency distributions.

The key underlying reason for both types of chaos is that chaos can only occur when there are at least three degrees of freedom. Each synchronized cluster can be characterized by a time-varying shape variable and a mean phase variable, which are the active degrees of freedom, and the interaction of these collective coordinates can lead to chaos. We reduce the full Kuramoto model to the evolution equations for these collective coordinates. Gottwald (2015, 2017); Hancock and Gottwald (2018) We demonstrate a time-scale splitting between the (slow) shape and the (fast) phase variables, that enables further reduction. Under this reduction, the full Kuramoto model with $M$ clusters reduces to a renormalized Kuramoto model with $M$ oscillators, which has $M-1$ degrees of freedom, implying that $M\geq 4$ is necessary for phase chaos to occur. However, when a cluster intermittently desynchronizes, the time-scale splitting is invalid, yielding additional active degrees of freedom, and the potential for chaos with three clusters.

The paper is organized as follows; in §II we describe the Kuramoto model. Then in §III we present the collective coordinate ansatz and derive the evolution equations for the collective coordinates. In §IV we show that phase chaos occurs for four clusters, and that there is quantitative agreement between the leading Lyapunov exponent for the full Kuramoto model and the collective coordinate reduction. In §V we show that chaos can occur for three clusters via intermittent desynchronization of a cluster, and provide a detailed description of this mechanism. Again, there is quantitative agreement between the leading Lyapunov exponent for the full Kuramoto model and the collective coordinate reduction. In §VI we show that chaos is not possible for two clusters in the thermodynamic limit of infinitely many oscillators. Lastly, in §VII we summarize our results and provide an outlook for future studies.

II The model

The Kuramoto model has been widely used to model networks of coupled oscillators, Kuramoto (1984); Strogatz (2000); Pikovsky, Rosenblum, and Kurths (2001); Acebrón et al. (2005); Osipov, Kurths, and Zhou (2007); Arenas et al. (2008); Dörfler and Bullo (2014); Rodrigues et al. (2016) in large part due to its analytical tractability. For a network of $N$ coupled oscillators, each with phase $\phi_{i}$, the dynamics is given by

[TABLE]

where the natural frequencies $\omega_{i}$ are drawn from a distribution $g(\omega)$, $A$ is the adjacency matrix of the network, i.e., $A_{ij}=1$ if nodes $i$ and $j$ are connected, otherwise $A_{ij}=0$, and $K$ is the coupling strength. We shall restrict our study of collective chaos to an all-to-all coupling topology with $A_{ij}=1-\delta_{ij}$. For the exposition of the model reduction technique presented in §III, however, we choose to present the Kuramoto model (1) with a general topology. It is widely known that if the coupling strength is sufficiently large, then the oscillators spontaneously synchronize, all oscillating at the same frequency, even though their natural frequencies are different. Furthermore, below the global synchronization threshold, synchronized clusters may emerge due to either clusters in the network topology, or distinct modes in the natural frequency distribution, or both.

We consider multimodal natural frequency distributions $g(\omega)$ of the form

[TABLE]

such that each $g_{m}$ is a normal distribution with mean $\Omega_{m}$ and variance $\sigma_{m}^{2}$, and the weights $0\leq\gamma_{m}\leq 1$ satisfy $\sum\gamma_{m}=1$. In particular, we primarily consider the case of well-separated peaks, such as the the example shown in Fig. 1. The distribution (2) has $M$ peaks, which typically correspond to $M$ clusters of synchronized oscillators for a range of coupling strengths. Note that the Kuramoto model is invariant under uniform phase shifts. Therefore, we may assume without loss of generality that the mean natural frequency is zero, i.e., $\sum_{m}\gamma_{m}\Omega_{m}=0$.

A characterization of the state of the system is the instantaneous order parameter $r(t)$ which is defined as

[TABLE]

and describes the mean position of all oscillators in the complex plane. The long term dynamics can be characterized by the time-averaged order parameter

[TABLE]

which is independent of $t_{0}$ and $T$ for sufficiently long transient times $t_{0}$ and averaging times $T$. If $\bar{r}$ is close to $1$, the oscillators are globally synchronized. If $\bar{r}\approx 1/\sqrt{N}$, the oscillators are in the incoherent state. In addition, for cases with multiple synchronized clusters, we can define analogous instantaneous and time-averaged order parameters for each cluster. For example, the instantaneous order parameter for the $m$-th cluster is

[TABLE]

where $\mathcal{C}_{m}$ is the set of oscillators in cluster $m$, and $N_{m}$ is the number of oscillators in $\mathcal{C}_{m}$.

Consider, for example, the frequency distribution shown in Fig. 1, with four peaks. The modes $\Omega_{m}$ are equally spaced between $-1$ and $1$, and the standard deviations are all the same with $\sigma_{m}=0.05$. For $N=100$ equiprobably111We find the values $\omega_{i}$, $i=1,\dots,N$, such that $G(\omega_{i})=(i+1)/(N+1)$, where $G(\omega)$ is the cumulative density function of the natural frequencies. By sampling equiprobably we eliminate the need to average over many randomly drawn realizations. drawn oscillators from this distribution, the time-averaged order parameter, $\bar{r}$, is shown for $0<K<2.5$ in Fig. 2. For $K<0.3$, the oscillators are incoherent, and $\bar{r}$ is of the order $1/\sqrt{N}$. For $K>1.6$, the oscillators globally synchronize, forming a single cluster, and $\bar{r}\approx 1$. For intermediate values, i.e., $0.3<K<1.6$, the oscillators corresponding to each peak in $g(\omega)$ synchronize to form a cluster, but they do not globally synchronize, resulting in $\bar{r}\approx 0.45$. In this study we are mostly interested in these intermediate values, where there can be complex interactions within and in between clusters. Note that $\bar{r}$ exhibits unusual non-monotonic behavior around $1<K<1.25$, which, as we shall see, is the region where chaotic dynamics occurs.

Synchronization of clusters is shown by the snapshots of oscillators in the complex plane in Fig. 3 for four different values of $K$. The oscillators of each color (corresponding to the same colored peak in Fig. 1) are synchronized, but there are clearly four distinct clusters. These clusters have both their own internal dynamics and interact with other clusters. For $K=0.9$, the dynamics is quasiperiodic, demonstrated by the trajectory of the complex order parameter $r(t)e^{i\psi(t)}$ shown as the blue curve inside the circle in Fig. 3(a) (Multimedia view). Increasing $K$, at a critical coupling strength $K_{c}$ the dynamics becomes chaotic. For example, with $K=1.2$, shown in Fig. 3(b) (Multimedia view), the dynamics is chaotic (which is confirmed by computing the leading Lyapunov exponent, $\lambda=6.18\times 10^{-2}$). The dynamics then becomes regular again, for example with $K=1.22$ and $K=1.3$ the trajectory of the complex order parameter is periodic (cf. Fig. 3(c) (Multimedia view) and Fig. 3(d) (Multimedia view), respectively). For $K=1.3$, the trajectory is confined to a straight line due to the existence of an attracting symmetric manifold. Four cluster cases such as these will be discussed in more detail in §IV.

III Model reduction via collective coordinates

Since we are primarily interested in the macroscopic inter- and intracluster dynamics, we use model reduction to reduce the high dimensional full Kuramoto model (1) to a small number of active degrees of freedom. One frequently used method is the Ott-Antonsen approach, Ott and Antonsen (2008) which assumes infinitely many oscillators. Recently, an alternative approach for model reduction has been proposed, termed collective coordinate reduction,Gottwald (2015, 2017); Hancock and Gottwald (2018) which can be readily applied to finite networks of coupled oscillators.

The idea of the collective coordinate reductionGottwald (2015, 2017); Hancock and Gottwald (2018) is to express the $N$-dimensional phase vector $\bm{\phi}$ as a linear combination of a small number of dynamically relevant modes. Intuitively, the reduction is motivated by the fact that synchronization is characterized by oscillators forming a collective entity which is described by its mean phase and its shape. The time-varying coefficients of the linear combination are coined collective coordinates, and encode the temporal evolution of the modes. Identification of the relevant modes is situation-dependent. In the case of a single synchronized cluster of oscillators, where the global phase is not relevant, a single mode $\bm{\Phi}$ describing the shape suffices, and we approximate $\bm{\phi}(t)\approx\alpha(t)\bm{\Phi}$. When multiple clusters interact, phase variables need to be accounted for. We will denote the shape modes by $\bm{\Phi}^{(m)}$ and the phase modes by $\bm{1}_{N_{m}}$ (the vector consisting of $N_{m}$ $1$’s, where $N_{m}$ is the size of the $m$-th cluster $\mathcal{C}_{m}$), with associated collective coordinates $\alpha_{m}$ and $f_{m}$, respectively, such that

[TABLE]

where typically $2M\ll N$.

The method of collective coordinates Gottwald (2015, 2017); Hancock and Gottwald (2018) is in effect a Galerkin approximation, where the residual error made by the ansatz (3) is minimized and the minimization leads to a system of evolution equations for the collective coordinates $\alpha_{m}(t)$ and $f_{m}(t)$.

The choice of basis functions is crucial. The shape mode $\bm{\Phi}^{(m)}$ can be chosen via linearization of the Kuramoto model (1), restricted to oscillators in $\mathcal{C}_{m}$. For sufficiently large coupling strengths $K$, $\bm{\Phi}^{(m)}$ will solve the Kuramoto model to good accuracy (ignoring the interactions with any oscillators outside of $\mathcal{C}_{m}$).

We follow the methods outlined previously Gottwald (2015, 2017); Hancock and Gottwald (2018) and derive a collective coordinate reduction for multimodal natural frequency distributions of the form (2). We first present the reduction for a single synchronized cluster of oscillators, and then present results for several interacting clusters.

III.1 Single cluster ansatz

Linearizing the full Kuramoto model (1) around $\phi_{i}-\phi_{j}=0$ for all $i,j$ results in

[TABLE]

where $L=D-A$ is the graph Laplacian, and $D$ is the diagonal degree matrix, i.e., $D_{ii}$ is the degree of node $i$. Note that $L$ has a nontrivial kernel with $L\bm{1}_{N}=0$, associated with the invariance to a global constant phase shift. Global synchronization corresponds to all oscillators rotating at the mean natural frequency $\Omega=(1/N)\sum_{i}\omega_{i}$. Substituting $\dot{\bm{\phi}}=\Omega\bm{1}_{N}$ into (4) we obtain the global synchronization mode

[TABLE]

where $L^{+}$ denotes the pseudoinverse of $L$, and we note that $L^{+}\bm{1}_{n}=\bm{0}$. Therefore, the single cluster ansatz function is

[TABLE]

with collective coordinate $\alpha(t)$. For all-to-all coupling, $L=NI_{N}-\bm{1}_{N}\bm{1}_{N}^{T}$. Therefore, $L^{+}=\frac{1}{N}(I_{N}-\frac{1}{N}\bm{1}_{N}\bm{1}_{N}^{T})$ 222$L$ has a zero eigenvalue with multiplicity one, and a repeated eigenvalue $\lambda=N$ with multiplicity $N-1$. Diagonalizing $L$, we have $L=VDV^{T}$ where $D=\text{diag}(N,\dots,N,0)$ and so $L^{+}=VD^{+}V^{T}=\frac{1}{N^{2}}L$. and

[TABLE]

Note that for a single synchronized cluster, as a result of the phase shift invariance, we may assume, without loss of generality, that $\Omega=0$. For multiple synchronized clusters, the different mean natural frequencies of each cluster must be accounted for, which we show in §III.2.

The evolution equation for the collective coordinate $\alpha(t)$ can be found as a Galerkin approximation using the same approach as in previous studies. Gottwald (2015, 2017); Hancock and Gottwald (2018) The ansatz (6) is substituted into the Kuramoto model (1), yielding a residual error

[TABLE]

This residual error, which is a two-dimensional manifold parametrized by $\alpha$ and $\dot{\alpha}$, is minimized when it is orthogonal to the one-dimensional line $\alpha\hat{\bm{\phi}}$ that we are restricting the solution to. Setting $\bm{\mathcal{E}}\cdot\hat{\bm{\phi}}=0$ we obtain an evolution equation for the collective coordinate $\alpha$

[TABLE]

For all-to-all coupling with mean frequency $\Omega=0$, this simplifies to

[TABLE]

where $\Sigma^{2}=(1/N)\sum_{i}\omega_{i}^{2}$ is the variance of the natural frequencies. Setting $\beta=\alpha/K$, so that $\bm{\phi}\approx\beta\bm{\omega}$, yields

[TABLE]

Stationary points of (7) correspond to synchronized states for the Kuramoto model.

In the thermodynamic limit, $N\to\infty$, (7) becomes

[TABLE]

For normally distributed natural frequencies, with mean zero and variance $\sigma^{2}$, we obtain

[TABLE]

Since $\mathcal{I}(0)=1$, it follows that $\beta$ has a stationary point if and only if $\mathcal{I}$ has a negative local minimum. Solving $\frac{d\mathcal{I}}{d\beta}=0$ and $\frac{d^{2}\mathcal{I}}{d\beta^{2}}>0$ yields $\beta=(\sigma\sqrt{2})^{-1}$. Therefore, $\beta$ has a stationary point if and only if $\mathcal{I}\left(\left(\sigma\sqrt{2}\right)^{-1}\right)\leq 0$, which is equivalent to

[TABLE]

If the condition (10) is satisfied, the oscillators synchronize and form a single cluster.

The instantaneous order parameter for the collective coordinates can be calculated as

[TABLE]

This relation shows that $\beta$ measures the spread of the oscillators. Large values of $\beta$ correspond to small $r$, meaning the oscillators are evenly distributed on the circle, whereas small values of $\beta$ for which $|\bm{\Phi}|\ll 1$ correspond to $r\approx 1$, corresponding to tightly clustered oscillators.

For the multimodal natural frequency distribution (2) with $M$ peaks we obtain

[TABLE]

As for a unimodal distribution, $\mathcal{I}(0)=1$, and a stable stationary solution of (8), corresponding to global synchronization of oscillators, exists if and only if the minimum of $\mathcal{I}(\beta)$ (obtained numerically) is negative. Therefore, the condition for global synchronization is

[TABLE]

III.2 Multiple cluster ansatz

For multimodal frequency distributions, there is generally a range of $K$ values which are sufficiently large that oscillators form synchronized clusters, $\mathcal{C}_{1},\dots,\mathcal{C}_{M}$, corresponding to each peak in the distribution, but which are not sufficiently large to allow for global synchronization. In such a case, we use a modified ansatz, which accounts for intracluster and intercluster dynamics. Note that while we are primarily concerned with clusters originating from a multimodal natural frequency distribution, the same analysis can be performed for topological clusters. Hancock and Gottwald (2018) For oscillators in cluster $\mathcal{C}_{m}$, the intracluster dynamics is given by the restricted Kuramoto model

[TABLE]

where for now we ignore the influence of oscillators belonging to different clusters $k\neq m$. Following the same linearization procedure as for the full Kuramoto model yields the intracluster mode

[TABLE]

where $L_{m}^{+}$ is the pseudo-inverse of the graph Laplacian of the subgraph obtained by restricting to nodes in cluster $\mathcal{C}_{m}$. In the case of all-to-all coupling we obtain

[TABLE]

where $N_{m}$ is the number of oscillators in cluster $m$, and $\Omega_{m}$ is the mean frequency of cluster $m$. Note that for well separated peaks in the frequency distribution, such as the example in Fig. 1, $N_{m}/N\approx\gamma_{m}$, where $\gamma_{m}$ is the weighting of peak $m$ in the natural frequency distribution (2).

The intracluster mode $\hat{\bm{\phi}}^{(m)}$ does not account for interactions with oscillators not belonging to cluster $m$. Therefore, $\hat{\bm{\phi}}^{(m)}$ does not capture the asymptotic dynamics of the system for large $K$, where the oscillators will globally synchronize and form a single cluster. For global synchronization, the single cluster ansatz $\hat{\bm{\phi}}$ in (5) is a more appropriate mode. We remark that one can perform a Galerkin approximation valid for all coupling strengths by considering a linear superposition of the single cluster mode (5), the superposition of all possible synchronized clusters (13), as well as all possible mergings of synchronized clusters, each of these equipped with their own collective coordinate. However, since the advantage of employing the collective coordinate reduction is simplicity, which allows us to study the dynamics of the $N$-dimensional Kuramoto model, we prefer to use Galerkin approximations tailored for a particular dynamical range, parametrized by the coupling strength $K$.

When studying intercluster dynamics between cluster modes (13), the Galerkin approximation needs to account for the mean phases of each cluster, denoted $f_{m}$. These phases vary in time due to interactions between clusters. Accounting for these phase interactions, and the possibility of all clusters merging into a single cluster, for oscillators in cluster $m$ we propose the ansatz

[TABLE]

where $\pi^{(m)}$ denotes projection onto the nodes in cluster $m$, i.e. $\pi^{(m)}(v_{i})=v_{i}$ if $i\in\mathcal{C}_{m}$ and $\pi^{(m)}(v_{i})=0$ if $i\notin\mathcal{C}_{m}$. Here $\alpha,\alpha_{m}$ and $f_{m}$, $m=1,\dots,M$, are the collective coordinates. As for the single cluster ansatz, the dynamics for the collective coordinates are obtained by substituting the ansatz (14) into the Kuramoto model (1) to determine the residual error. Then, to ensure errors are minimized, we require the error to be orthogonal to the restricted solution hyperplane, spanned by $\hat{\bm{\phi}}$, $\hat{\bm{\phi}}^{(m)}$ and $\bm{1}_{N_{m}}$. The condition that the residual error is orthogonal to $\hat{\bm{\phi}}$ is given by

[TABLE]

where

[TABLE]

is the right hand side of the Kuramoto model (1) in vector form. The condition that the residual error is orthogonal to $\hat{\bm{\phi}}^{(m)}$ is given by

[TABLE]

(We note that since $\hat{\bm{\phi}}^{(m)}$ is orthogonal to $\bm{1}_{N_{m}}$ there is no $\dot{f}_{m}$ term). Lastly, the condition that the residual error is orthogonal to $\bm{1}_{N_{m}}$ is given by

[TABLE]

Equations (15)–(17) form a system of linear equations

[TABLE]

where $\bm{x}=(\alpha,\alpha_{1},\dots,\alpha_{m},f_{1},\dots,f_{m})^{T}$ is the vector comprised of the collective coordinates. This linear system can be solved to find the evolution equations for each of the collective coordinates.

In the case of all-to-all coupling, the projection $\pi^{(m)}\hat{\bm{\phi}}$ (5), the cluster modes $\hat{\bm{\phi}}^{(m)}$ (13) and the constant vectors $\bm{1}_{N_{m}}$ are linearly dependent, and so the ansatz (14) simplifies to

[TABLE]

where $\beta_{m}=\frac{1}{K}\left(\frac{N}{N_{m}}\alpha_{m}+\alpha\right)$ and $\tilde{f}_{m}=f_{m}+\alpha\frac{\Omega_{m}}{K}$. This means that the global synchronization ansatz (5) can be fully described by the cluster modes (13) with suitable mean phases of each mode, and so the collective coordinate $\alpha$ associated with global synchronization can effectively be ignored. 333For Erdős–Renyi networks with large number of oscillators, we observe that the projection $\pi^{(m)}\hat{\bm{\phi}}$ (5), the cluster mode $\hat{\bm{\phi}}^{(m)}$ (13) and the constant vector $\bm{1}_{N_{m}}$ are almost linearly dependent, and so $\alpha$ can be ignored from the ansatz to leading order. In essence, $\beta_{m}$ measures the spread of the oscillators within cluster $m$, and $\tilde{f}_{m}$ determines the collective phase of the cluster.

For the ansatz (III.2), the evolution equations for the collective coordinates obtained from (15)–(17) become

[TABLE]

where we have dropped the tilde on $f_{m}$, and $\Sigma_{m}^{2}=\frac{1}{N_{m}}\sum_{i\in\mathcal{C}_{m}}\left(\omega^{(m)}_{i}-\Omega_{m}\right)^{2}$ is the variance of the frequencies in cluster $m$. In the following, we consider all-to-all networks, unless stated otherwise, and consider (19)–(20). Therefore, for $M$ peaks in the frequency distribution, there are $2M$ equations of motion. By introducing phase difference variables, $F_{m}=f_{m+1}-f_{m}$, we reduce the dimension of the system to $2M-1$ degrees of freedom. This suggests that chaos may be possible as long as $M\geq 2$. However, as we will show, chaos is only possible if $M\geq 3$.

In the thermodynamic limit, $N\to\infty$, with a multimodal natural frequency distribution of the form (2), $N_{m}/N\to\gamma_{m}$ (the weight of cluster $m$), and the evolution equations for the collective coordinates (19)–(20) become

[TABLE]

Note that for $M=1$, i.e., unimodal, normally distributed frequencies with $\gamma_{1}=1$ and $\Omega_{1}=0$, (21) recovers the single cluster evolution equation (9), and (22) is identically zero.

III.3 Slow-fast splitting of the shape and phase coordinates

Each synchronized cluster viewed in isolation contains oscillators with normally distributed natural frequencies. Therefore, the instantaneous order parameter for each cluster is given in the thermodynamic limit (cf. (11)) by

[TABLE]

Expressing the evolution equations (21)–(22) for the collective coordinates $\beta_{m}$ and $f_{m}$ in terms of $r_{m}$, we obtain

[TABLE]

In the case that each cluster remains tightly clustered for all time, we have $r_{m}(t)=1-\epsilon_{m}(t)$, with $0<\epsilon_{m}(t)\ll 1$ for all $t$. This is ensured provided the $\sigma_{m}$ are sufficiently small, $K$ is sufficiently large (i.e., the condition (10) is satisfied for each $\sigma_{m}$), and the means $\Omega_{m}$ are sufficiently far apart relative to the coupling strength (i.e., condition (12) fails and global synchronization does not occur). Expanding (23)–(24) in powers of $\epsilon$ yields

[TABLE]

We can view the order parameters $r_{m}$ as describing the intracluster dynamics and the phase coordinates $f_{m}$ as describing the intercluster dynamics. Since $\epsilon_{m}\ll 1$, the evolution equations (25)–(26) for $r_{m}$ and $f_{m}$ reveal a time-scale splitting of the dynamics, whereby the order parameters $r_{m}$ evolve slowly, whereas the phase variables $f_{m}$ evolve on a fast time scale. The intercluster dynamics is, to first-order, decoupled from the intracluster dynamics (cf. (26)). Hence, the intercluster dynamics obeys a reduced, renormalized Kuramoto model. Since the reduced intercluster dynamics has $M-1$ degrees of freedom (taking into account a change to phase difference variables), chaos is only possible if $M\geq 4$. We label this type of chaos where clusters remain localized, with only small changes in their order parameter, as phase chaos. However, it is possible that one or more of the clusters intermittently break-up, such that $r_{m}\ll 1$ and $\epsilon_{m}$ is not small anymore. In such a case, there is significant interplay between the intracluster and intercluster dynamics. This will be discussed in more detail in §V.

IV Four clusters: collective phase chaos

Phase chaos is typically observed in systems with multimodal natural frequency distributions with at least four peaks [cf. Fig. 3(b)]. The simplest case is to take four oscillators with natural frequencies equally spaced between $-1$ and $1$ and let them interact to produce chaotic dynamics. Popovych, Maistrenko, and Tass (2005); Maistrenko, Popovych, and Tass (2005); Miritello, Pluchino, and Rapisarda (2009); Carlu, Ginelli, and Politi (2018) One may then consider $N$ oscillators distributed over these four distinct natural frequencies, or, more generally, consider the natural frequency distribution of $M$ distinct mean frequencies $\Omega_{m}$

[TABLE]

where $\delta(x)$ denotes the Dirac delta-function. If $N$ frequencies, $\omega_{1},\dots,\omega_{N}$, are distributed equiprobably onto the $M$ mean frequencies $\Omega_{m}$, with $N$ divisible by $M$, then each mean frequency $\Omega_{m}$ is populated by $N/m$ oscillators with $\omega_{i}=\Omega_{m}$. That is, we can relabel such that $\omega_{1},\dots,\omega_{N/M}=\Omega_{1}$, $\omega_{N/M+1},\dots,\omega_{2N/M}=\Omega_{2}$, and so on. The Kuramoto model (1) for oscillators with natural frequency $\Omega_{m}$ and all-to-all coupling in this case becomes

[TABLE]

Since the coupling is all-to-all, oscillators with the same natural frequency will synchronize, such that $\bm{\phi}^{(m)}(t)=f_{m}(t)\bm{1}_{N/M}$ and (28) becomes

[TABLE]

which is of the exact form as the Kuramoto model for $M$ oscillators. Hence, chaos is expected for arbitrarily many oscillators if their natural frequencies are distributed according to (27) with $M\geq 4$ with equally spaced $\Omega_{m}$. Note that the evolution equation for the phases $f_{m}$ (29) is equivalent to the collective coordinate equations for $M$ clusters (25)–(26) in the limit $\epsilon_{m}\to 0$, which is the limit of perfectly synchronized clusters, with identical phases within each cluster.

Considering the Dirac $\delta$-function as the limit of normal distributions, i.e., $\delta(x)=\lim_{\sigma\to 0}\mathcal{N}(0,\sigma)$, we expect multimodal distributions of the form (2) with $M\geq 4$ to yield phase chaos for sufficiently small $\sigma$ and sufficiently large spacings between peaks in the natural frequency distribution, $|\Omega_{m+1}-\Omega_{m}|$. Our focus in this section is to explore the collective dynamics of the Kuramoto model for natural frequency distributions $g(\omega)$ of the form (2) with identical weights $\gamma_{m}=1/4$, identical standard deviations $\sigma_{m}=\sigma$, and equally spaced means $\Omega_{m}=-1+2(m-1)/3$ for $m=1,\dots,4$, that is

[TABLE]

where $P$ is the standard normal distribution.

We now numerically explore these cases for the full Kuramoto model (1) with $N=100$ oscillators; and shall compare our results with the reduced collective coordinate description (19)–(20) for $N=100$ oscillators as well as with the reduced collective coordinate description (21)–(22) in the thermodynamic limit of infinitely many oscillators. The collective coordinate systems involve 7 degrees of freedom: four shape parameters $\beta_{m}$ and three phase-difference variables $f_{m+1}-f_{m}$ (the evolution equations, however, are written for $f_{m}$ and hence are 8-dimensional). We shall see that the collective coordinate equations provide a reduced model that allows for a quantitative description of the chaotic dynamics of the Kuramoto model, and, in particular, for the estimation of the Lyapunov exponents of the full Kuramoto model. We compute and compare the time-averaged order parameter $\bar{r}$ and the leading Lyapunov exponent $\lambda$ across a multitude of cases for different coupling strengths $K$ and for different standard deviations of the natural frequency distribution $\sigma$.

Before discussing the results on the leading Lyapunov exponent we shall describe the dependence of $\bar{r}$ on $K$ and $\sigma$, shown in the left column of Fig. 4. Shown is the order parameter $\bar{r}$ for the full Kuramoto model (1) with $N=100$ oscillators [Fig. 4(a)], the 8D collective coordinate model (19)–(20) with $N=100$ oscillators [Fig. 4(c)], and the 8D collective coordinate model with infinitely many oscillators (21)–(22) [Fig. 4(e)]. We see good quantitative agreement between all three models throughout most of the parameter space. All three models show transitions from $\bar{r}\approx 0.45$ to $\bar{r}\approx 1$ near $K\approx 1.58$, which is the transition from four synchronized clusters to global synchronization with one synchronized cluster. This transition can be predicted by the collective coordinate ansatz, using the single cluster ansatz (6) applied to the full distribution $g(\omega)$. The transition curve is given by the condition (12) for global synchronization, and is shown by the dashed, approximately vertical, curves in Fig. 4. The transition from the incoherent state ($\bar{r}\approx 0$) to the synchronized cluster state ($\bar{r}\approx 0.45$) is predicted by the line $K=\sigma\sqrt{2e}$ (dot-dashed in Fig. 4), which derives from condition (10) for the collective coordinate ansatz. However, this line does not accurately capture the transition from incoherence to synchronized clusters in the full Kuramoto model with $N=100$ oscillators [cf. Fig. 4(a)] for which the transition occurs at lower values of $K$. This discrepancy is due to the fact that the collective coordinate models (19)–(20) and (21)–(22) do not account for partial synchronization of the clusters. In the full Kuramoto model (1) the transition from the incoherent state to a partially synchronized state is a soft second-order phase transition whereby, upon increasing the coupling strength, more and more oscillators with natural frequencies close to the mean frequency mutually synchronize until at a critical coupling strength all oscillators in a cluster have synchronized. Although this can be quantitatively described by the collective coordinate ansatz Gottwald (2015); Hancock and Gottwald (2018) we knowingly do not account for this in our simulations here to limit the computational cost of the parametric sweep.

It is remarkable that the collective coordinate models — (19)–(20) for $N=100$ oscillators [Fig. 4(d)] and (21)–(22) for $N\to\infty$ [Fig. 4(f)] — reproduce the leading Lyapunov exponent $\lambda$ of the full Kuramoto model (1) [Fig. 4(b)] with good quantitative agreement. In particular, there is a chaotic “bubble” within the region with four synchronized clusters (between the dot-dashed and dashed curves) whose width shrinks as $\sigma$ increases. The occurrence of partial synchronization of clusters in the full Kuramoto model with $N=100$ results in a positive Lyapunov exponent above and near to the dot-dashed line in Fig. 4(b), which is not captured by the collective coordinate models [Fig. 4(d) and Fig. 4(f)]. This difference is due to complex interactions between the synchronized clusters and the small number of oscillators that do not synchronize, which are not accounted for by the collective coordinate models.

In the limit as $\sigma\to 0$, the dynamics of four interacting clusters becomes equivalent to the dynamics of four interacting oscillators (cf. (29)), which has been studied extensively by Maistrenko et al. Maistrenko, Popovych, and Tass (2005) and Popovych et al. Popovych, Maistrenko, and Tass (2005) Following the approach of previous studies, we consider the first four Lyapunov exponents of the collective coordinate model (21)–(22). For small values of $\sigma<10^{-2}$ we obtain Lyapunov exponents that are qualitatively the same as those observed for four individual oscillators (compare Fig. 5(a) with Fig. 1(a) in Ref. (15)). Therefore, for these small values of $\sigma$ the bifurcation sequence is essentially the same as for four individual oscillators. At $K=K_{\text{sn}}\approx 0.91$ there is a saddle-node bifurcation, which transitions from quasiperiodic to periodic dynamics. At $K=K_{\text{td}}\approx 0.94$ there is a transition to chaos via the Afraimovich-Shilnikov torus destruction scenario. Afraimovich and Shilnikov (1991) At $K=K_{\text{cr}}\approx 1.22$ the chaotic attractor is destroyed in a boundary crisis, yielding a chaotic saddle. Lastly, at $K=K_{c}\approx 1.58$ the transition to global synchronization occurs. There are many periodic regions observed within the chaotic region $K_{\text{td}}<K<K_{\text{cr}}$, and also near $K=1.5$, which correspond to the resonances discussed by Maistrenko et al. Maistrenko, Popovych, and Tass (2005) The resonances within the chaotic region can also be observed within the chaotic bubble shown in the right plots of Fig. 4, evident as white bands ($\lambda_{1}=0$) that extend approximately vertically from the horizontal axis $\sigma=0$ (most clearly seen in Fig. 4(f) which has the highest resolution). The resonances near $K=1.5$ can be seen in the right plots of Fig. 4 as thin bands of positive largest Lyapunov exponent. For larger values of $\sigma$, such as $\sigma=5\times 10^{-2}$ shown in Fig. 5(b), we see similar dynamics, but there are some key differences. First, the chaotic window is smaller, and is punctuated by a large periodic region near $K=1.05$. In addition, there appears to be only one resonance near $K=1.5$.

The complex bifurcation structure shown in Fig. 5 also explains the discontinuous transition curves between different shades of gray in the plots for $\bar{r}$ (left plots of Fig. 4). These transitions are due to bifurcations between different stable chaotic, periodic, and quasiperiodic states.

V Three clusters: chaos via intermittent cluster desynchronization

For three clusters, as discussed previously, if the reduction and time-scale splitting shown in (25)–(26) is valid, the dynamics is essentially phase dynamics of three oscillators, excluding chaotic dynamics because there are only two degrees of freedom (recall that due to the phase-gauge invariance of the Kuramoto model, we may assume without loss of generality that $\sum_{i}f_{i}=0$). However, the time-scale splitting requires $\epsilon_{m}=1-r_{m}\ll 1$ for all time. If this is not true, e.g. one cluster intermittently desynchronizes, then chaos is possible.

As an example, consider the trimodal natural frequency distribution shown in Fig. 6. Simulating the 6D collective coordinate model, (23)–(24), for $K=1.205$ we find a positive largest Lyapunov exponent, $\lambda=0.036$, as well as time-averaged cluster-wise order parameters $\bar{r}_{1}=0.989$, $\bar{r}_{2}=0.981$, $\bar{r}_{3}=0.918$. Hence the system is both chaotic and collectively organized. While $\bar{r}_{3}$ is close to one, and, hence, the cluster would be considered synchronized, the time series for $r_{3}(t)$, shown in Fig. 7(c), intermittently dips to values around $0.5$, showing that the cluster intermittently desynchronizes, with oscillators spreading over the entire circle. Therefore, we cannot say that $\epsilon_{3}$ is close to zero for all time, meaning the time-scale splitting is invalid for $r_{3}$, and, hence, chaos is possible.

This intermittent desynchronization phenomenon predicted by our collective coordinate reduction is confirmed in the full Kuramoto model (1). For $N=200$ oscillators, with natural frequencies drawn equiprobably from the distribution $g(\omega)$ shown in Fig. 6, we compute the leading Lyapunov exponent as $\lambda=0.039$, which is within $10\%$ of the Lyapunov exponent computed using collective coordinates in the thermodynamic limit, (23)–(24), which has $\lambda=0.036$. Furthermore, the time-series of $r_{1},r_{2},r_{3}$, shown in Fig. 7(d–f), are qualitatively similar to those shown in Fig. 7(a–c). In particular, $r_{1}$ and $r_{2}$ remain close to $1$ for all time, whereas $r_{3}$ experiences intermittent dips. The dips occur in the collective coordinate model (23)–(24) with an average period of $48.9$, compared to an average period of $52.6$ in the full Kuramoto model. For the full Kuramoto model (1) with $N=1000$ oscillators, which is closer to the thermodynamic limit, the dips occur at the same frequency as with $N=200$, i.e., with a period of $52.6$, and the time series of $r_{1,2,3}$, shown in Fig. 7(g–i), are even more similar to the collective coordinate model in the thermodynamic limit (23)–(24), shown in Fig. 7(a–c), in that the dynamics between the dips becomes more regular, with high frequency oscillations and a slow negative trend. The collective coordinate model is representative of the full Kuramoto model, and has the advantage of being more analytically tractable.

We now investigate more closely the nature of this type of chaotic dynamics and how it is generated. We first describe qualitatively the dynamics of a single desynchronization event in the full Kuramoto model (1). We then show that these desynchronization events can be resolved by considering further reductions of the collective coordinate equations (23)–(24). This collective coordinate reduction is then used to show that chaos via intermittent desynchronization is a robust phenomenon.

We describe the dynamics of a desynchronization event qualitatively using the snapshots of the phases of oscillators shown in Fig. 8(b–g) (Multimedia view), which correspond to the red points marked on the time series of $r_{3}$ shown in Fig. 8(a). In the lead-up to a dip in $r_{3}$, the second and third clusters are phase-locked, with an approximately constant phase difference $F_{2}=f_{3}-f_{2}$. However, each time the first cluster passes by the second cluster, the second cluster slows down, which causes a small increase in the phase separation between the second and the third clusters, implying a small increase in $F_{2}$, as shown in Fig. 8(b) and Fig. 8(c) as a small increase in separation between the second and third clusters. Eventually, a critical point is reached, such that the oscillators in the third cluster that are furthest from the second cluster [those with the highest natural frequencies, closest to $i=1000$ in Fig. 8(b–i)] begin to desynchronize with the rest of the oscillators in the cluster, as shown in the transition from Fig. 8(b) to Fig. 8(d). This desynchronization results in the oscillators in the third cluster wrapping around and covering the entire circle, and corresponds to a sharp dip in $r_{3}$, as shown in Fig. 8(a). The desynchronization of the third cluster occurs as a traveling front, starting first with the oscillators with highest natural frequency, traveling down to the oscillators with the lowest frequency. The oscillators in the third cluster eventually cover the entire circle, and those with the highest natural frequencies (furthest to the right in Fig. 8) overtake those with the lowest natural frequencies (furthest to the left in Fig. 8), meaning they experience additional revolutions during each “dip” event. Once the oscillators in the third cluster with lowest natural frequencies catch up with the second cluster, the third cluster re-synchronizes, as shown in Fig. 8(f) and Fig. 8(g), once again becoming phase-locked with the second cluster, and the process repeats.

We now use collective coordinate reductions to analyze the dynamical scenario described above. As seen in Fig. 7, $r_{1}(t)$ and $r_{2}(t)$ are close to one for all time, demonstrating that the time-scale splitting remains valid for those variables. This suggests that we may set $r_{1}(t)=\bar{r}_{1}$ and $r_{2}(t)=\bar{r}_{2}$ as constant in the collective coordinate equations (23)–(24). Then the collective coordinate model reduces to a system of four fast variables, $r_{3},f_{1,2,3}$, with three degrees of freedom (again since, without loss of generality, $\sum_{i}f_{i}=0$). The evolution equations become

[TABLE]

where $F_{m}=f_{m+1}-f_{m}$ and $\Delta\Omega_{m}=\Omega_{m+1}-\Omega_{m}$. We see good agreement in the time series plots of $r_{3}$ for the 3D system (30)–(32), shown as dashed black in Fig. 9) and for the full 6D collective coordinate system (23)–(24), shown as solid gray in Fig. 9. In both models, $r_{3}$ experiences the same oscillations at the start and end, and both have a significant dip to $r_{3}\approx 0.4$ between $t=10$ and $t=30$. Furthermore, the Poincaré sections through the plane $F_{1}=0$, shown in the $(r_{3},F_{2})$-plane in Fig. 10, are similar for both models.

To explain the pronounced dips in $r_{3}$ in more detail, observe that for the time series of $r_{3}$, shown in Fig. 7(c), in between the sharp dips, $r_{3}$ exhibits small oscillations and a small negative trend. To explain this, let us assume that $r_{3}$ is constant, so the dynamics (30)–(32) reduces to a 2D system for $F_{1}$ and $F_{2}$, given by (31)–(32), with $r_{3}$ being a parameter. For $r_{3}>r_{c}\approx 0.981$, this 2D system (31)–(32) has one stable and one unstable limit cycle, as demonstrated in Fig. 11(a) for $r_{3}=1$ by the thick solid and dashed red curves, respectively. The gray arrows in Fig. 11(a) are the 2D velocity field. As $r_{3}$ decreases, the stable and unstable limit cycles move toward each other, as demonstrated in Fig. 11(b) for $r_{3}=0.981$. At $r_{3}=r_{c}$, the stable and unstable limit cycles annihilate via a saddle-node bifurcation, and the dynamics is topologically equivalent to quasiperiodic rotation on the torus. We observe in Fig. 11(c) that trajectories of the full 6D collective coordinate system (23)–(24), projected onto the $F_{1},F_{2}$ plane, follow curves that closely match the limit cycles corresponding to constant $r_{3}$. The tracer (whose trajectory is shown in thin black) slowly advances upward in between the lower limiting stable limit cycle corresponding to $r_{3}=1$ (lower thick red curve), and the upper limiting stable limit cycle corresponding to $r_{3}=r_{c}$ (upper thick red curve). This slow advance upward corresponds to the slow decay of $r_{3}$ in between the sharp dips.

Expanding further, starting at a time $t_{0}$ when $r_{3}\approx 1$, a tracer in the full 6D collective coordinate model (23)–(24) will have a trajectory in the $F_{1},F_{2}$ plane that is very similar to the limit cycle obtained from the assumption that $r_{3}$ is constant (equal to $r_{3}(t_{0})$). However, while $r_{3}$ is approximately constant, it decreases slightly over one period of the limit cycle. We can approximate the decrease in $r_{3}$ by computing $\Delta r_{3}=r_{3}(t_{0}+T)-r_{3}(t_{0})$, where $T=T(r_{3}(t_{0}))$ is the period of the stable limit cycle, denoted by $\mathcal{C}_{r_{3}}$, of the 2D system (31)–(32) with $r_{3}=r_{3}(t_{0})$ held constant. Here $r_{3}(t_{0}+T)$ is found by integrating (30) along the stable limit cycle $\mathcal{C}_{r_{3}}$. This is valid under the assumption that $r_{1,2,3}$ are all constant between $t=t_{0}$ and $t=t_{0}+T$. Note that the values $r_{3}(t_{0}+T)$ and $\Delta r_{3}$ are independent of the initial locations of $F_{1},F_{2}$ on the limit cycle $\mathcal{C}_{r_{3}}$. We find that $\Delta r_{3}<0$ for all $r_{3}>r_{c}$, and so it is inevitable that $r_{3}$ will eventually reach the critical value, $r_{c}$, where the stable limit cycle bifurcates.

The scenario of chaotic dynamics through intermittent desynchronization events is a robust phenomenon, occurring for a range of parameters of the natural frequency distribution (2). We show this by investigating the effect of varying $\sigma_{3}$. As $\sigma_{3}$ decreases, we observe that the average time interval between dips in $r_{3}$ increases, and at a critical value of $\sigma_{3}=\sigma_{c}\approx 0.035$ the dips no longer occur. For $\sigma_{3}<\sigma_{c}$, $r_{3}$ remains close to 1 for all time, and so the slow-fast splitting found in §III.3 is valid, and the dynamics is non-chaotic. The value of $\sigma_{c}$ can be estimated using the collective coordinate system (30)–(32). Consider $\Delta r_{3}$, the change in $r_{3}$ over the stable limit cycle that exists under the assumption that $r_{3}$ is constant. The distribution of $\Delta r_{3}$ across a range of $r_{3}$ and $\sigma_{3}$ values is shown in Fig. 12. The turning point of the curve $\Delta r_{3}=0$ (dashed black in Fig. 12) yields the approximation $\sigma_{c}^{*}=0.038$ (solid black line in Fig. 12) for $\sigma_{c}$. For $\sigma_{3}>\sigma_{c}^{*}$, $\Delta r_{3}$ is negative for all values of $r_{3}$ that have a stable limit cycle. Hence $r_{3}$ decreases after each period of the limit cycle until reaching the saddle-node bifurcation (solid gray curve in Fig. 12). For $\sigma_{3}<\sigma_{c}^{*}$, the curve $\Delta r_{3}=0$ (dashed black in Fig. 12) indicates the locations of fixed points of the map $r_{3}\mapsto r_{3}+\Delta r_{3}$, with the right-most fixed point being stable. The presence of these stable fixed points indicates a periodic solution of the three-dimensional system (30)–(32). Therefore, $\sigma_{c}^{*}$ represents a bifurcation between periodic dynamics and intermittent desynchronization dynamics, i.e., it is an approximation for $\sigma_{c}$. Note that $\sigma_{c}^{*}=0.038$ slightly over-predicts $\sigma_{c}=0.035$, which is due to the inaccuracies that occur from making the assumption that $r_{3}$ is constant over the period of the limit cycle $\mathcal{C}_{r_{3}}$, when, as we have seen, it is both oscillating and slowly decreasing [cf. Fig. 8(a)]. A similar approach can be used to determine critical values at which chaos ceases to occur when other parameters in the natural frequency distribution are varied, such as the distance between peaks and the proportion of oscillators in each cluster.

We now explain why the transition into desynchronization occurs on a fast time-scale, as observed in Fig. 8(a), using the collective coordinate equations (30)–(32) for the three-cluster interactions. The intercluster interaction term between the second and third cluster in $\dot{r}_{3}$ is $G(r_{3})\cos F_{2}$, where $G(r_{3})$ scales as $r_{3}^{2}\sqrt{-2\log{r_{3}}}$ [see (30)], which is positive for $0<r_{3}<1$, and equal to zero at $r_{3}=0$ and $r_{3}=1$. On the stable limit cycles, $F_{2}$ oscillates around $\pi/2$. Therefore, the interaction term is small while on the limit cycles, but when $r_{3}<r_{c}$, and the saddle-node bifurcation occurs, $F_{2}$ increases away from $\pi/2$, and so $\dot{r}_{3}$ becomes strongly negative, explaining the sharp decline of $r_{3}$. At the point where $F_{2}$ crosses $3\pi/2$, the sign of $\cos(F_{2})$ changes, and so $\dot{r}_{3}$ becomes strongly positive, until $F_{2}$ once again approaches $\pi/2$, at which point $r_{3}$ once again becomes slow, and the system relaxes to a limit cycle corresponding to $r_{3}\approx 1$. This restarts the cycle of slow decay followed by a sharp decline and recovery.

We have established how chaos is generated through the delicate interaction of three clusters using the collective coordinate framework. As a summary, chaos occurs as a sensitivity between the entry and exit locations to the regular limit cycle zone. This sensitivity is shown by the infinite, fractal accumulation of folds in the zoomed in Poincaré section through the plane $F_{1}=0$ (cf. Fig. 10). The folds accumulate in the small region with $r_{3}\approx 1$ and $F_{2}\approx\pi/2$, corresponding to the regular limit cycle dynamics and slow, predictable decay of $r_{3}$. While we have shown that chaos is possible for trimodal natural frequency distributions, it is a rare phenomenon. In the process of finding the natural frequency distribution shown in Fig. 6, we computed the maximal Lyapunov exponent for $5\times 10^{4}$ randomly drawn sets of natural frequency distribution parameters $(\bm{\Omega},\bm{\sigma},\bm{\gamma})$ and coupling strengths $K$ that produce synchronized clusters, and found that only $90$ cases were chaotic (with a positive Lyapunov exponent), i.e., only $0.18\%$.

VI Two clusters: no chaos

For two clusters, $M=2$, the thermodynamic limit collective coordinate equations (23)–(24) become

[TABLE]

where $F=f_{2}-f_{1}$ is the phase difference of the two clusters, and $\Delta\Omega=\Omega_{2}-\Omega_{1}$. Hence, it appears that there are three degrees of freedom, and chaos is theoretically possible.

We now show, using the collective coordinate approach, that chaos is not possible. In particular, phase chaos is not possible as it would reduce the dimension of the system to 1D with both $r_{1}$ and $r_{2}$ being constant. The case of intermittent desynchronization leads to decoupled 1D slow and 2D fast dynamics, excluding the possibility of chaos.

Let us begin with excluding the possibility of phase chaos. If the time scale splitting between $r_{1},r_{2}$ (slow) and $F$ (fast) is valid, i.e., if the two clusters remain synchronized for all time, we can average the slow $r_{1},r_{2}$ dynamics over the fast dynamics $F$. Assuming $r_{1}$ and $r_{2}$ are constant, the dynamics of $F$ can be solved analytically, with solution

[TABLE]

where $\kappa=Kr_{1}r_{2}$, $B=\sqrt{\Delta\Omega^{2}-\kappa^{2}}$, $C=\arctan\left(\frac{-\kappa+\Delta\Omega\tan\frac{F_{0}}{2}}{B}\right)$, and $F(0)=F_{0}$. Note that $F$ is a periodic function, with period $T=2\pi/B$. For the bimodal frequency distribution shown in Fig. 13(a), where the peaks have very little overlap, the approximate solution (36) [dashed black in Fig. 13(b)] for the phase difference closely matches the time series of $F$ of the collective coordinate model (33)–(35) [solid red in Fig. 13(b)].

Furthermore, since $F(t)$ ranges from [math] to $\pi$, we can choose, without loss of generality, our starting time such that $F_{0}=\pi/2$, and so $C=\arctan\frac{\Delta\Omega-\kappa}{B}$. It can be shown that

[TABLE]

which implies that

[TABLE]

Therefore, $\cos F(t)$ is an odd periodic function, and so its average over one period, $\langle\cos F(t)\rangle$, is zero. This means the dynamics of the time-averaged variables $\bar{r}_{1},\bar{r}_{2}$ becomes decoupled from one another. The dynamics for each cluster is equivalent to the single cluster ansatz equation (9), with $K$ replaced by $K\gamma_{i}$ for $i=1,2$. Hence, $r_{1}(t)$ [solid blue curve in Fig. 13(c)] and $r_{2}(t)$ [solid red curve in Fig. 13(c)] oscillate around the stable equilibria, $r_{1}^{*}$ and $r_{2}^{*}$ (dashed blue and red resp.), obtained from the respective single cluster ansatz equations, and phase chaos cannot occur.

Now we go on to exclude the case that one cluster intermittently desynchronizes, like in the three cluster case discussed in the previous section. This occurs for the natural frequency distribution shown in Fig. 13(d), where the second cluster intermittently desynchronizes, approaching $r_{2}\approx 0.4$, as shown by the solid red curve in Fig. 13(f). In this case, the dynamics of $r_{1}$, which remains close to 1 for all time, is slower than $r_{2}$ and $F$. This is confirmed in Fig. 13(e), where it is shown that the time evolution of $F$ given by numerical simulation of (33)–(35) (solid red) is not well approximated by the function (36) (dashed black), which assumes perfect time-scale splitting. Therefore, we may not assume time-scale separation between $r_{2}$ and $F$. We have an effective 2D fast system for $r_{2}$ and $F$. This 2D system has a stable limit cycle in cases with two clusters that do not globally synchronize. In turn, the dynamics of $r_{1}$ cannot be chaotic, since the time-averaged dynamics is a 1D system with time-periodic forcing.

The only other possibility is that both clusters intermittently desynchronize. However, since it is the intercluster terms in $\dot{r}_{1}$ and $\dot{r}_{2}$ that drive the push away from the single cluster ansatz equilibrium, and both intercluster terms are multiples of $r_{1}r_{2}\cos F$, it follows that $r_{1}$ cannot rapidly decay without $r_{2}$ also rapidly decaying. If one, say $r_{1}$, decays faster than the second, $r_{2}$, then it will asymptote toward $r_{1}=0$, and so it has no effect on the second cluster. The second cluster is then governed by the single cluster ansatz, and will either approach the stable synchronized state, or approach $r_{2}=0$, depending on whether $r_{2}$ crosses the unstable fixed point of the single cluster ansatz equation while the first cluster is desynchronizing. In either case, the dynamics is regular, and stationary in the long run. If both $r_{1}$ and $r_{2}$ decay at the same rate, then the system possesses a symmetry, which further reduces the effective dimension, excluding the possibility of chaos.

The above discussion used the thermodynamic limit. In finite size networks, however, chaos can occur for bimodal natural frequency distributions. This occurs due to sampling effects. In our numerical simulations of finite size networks, we found that it is typical that when chaos occurs, a small group of oscillators, with natural frequencies at one or the other extreme of the distribution (i.e., very high or very low), do not synchronize with the other oscillators corresponding to the same peak in the natural frequency distribution. This group of “rogues” may either constitute a set of incoherent oscillators or another small cluster. In either case, the system must be considered as having more than two clusters, which agrees with our results obtained in §IV and §V. We find fewer chaotic cases as the number of oscillators increases, which confirms that the issue is a finite-size effect. It is important to note that we have found no bimodal cases with finite $N$ that are chaotic and do not have unsynchronized rogue oscillators.

VII Summary and outlook

VII.1 Summary

Employing detailed numerical simulations guided by analytical results from a collective coordinate reduction we have established necessary conditions for collective chaos in the Kuramoto model with multimodal natural frequency distributions. We have shown that phase chaos can occur provided there are at least four peaks in the natural frequency distribution. This is due to a time-scale splitting between slow intracluster collective coordinates and fast intercluster collective coordinates, which reduces the Kuramoto model to $M-1$ active degrees of freedom, where $M$ is the number of peaks in the natural frequency distribution.

For three peaks in the natural frequency distribution, we have shown that chaos can occur via intermittent desynchronization of clusters. When a cluster desynchronizes, its intracluster collective coordinate becomes fast, resulting in an additional active degree of freedom. Through the slow-fast splitting, the collective coordinate description has allowed us to study the intricate dynamics of intermittent desynchronization, and show that it is a robust phenomenon.

For two peaks in the natural frequency distribution, the collective coordinate description has allowed us to rule out the possibility of chaos.

We have shown that for both phase chaos and chaos via intermittent desynchronization, the reduced collective coordinate description can be used to quantitatively predict the leading Lyapunov exponent, and, hence, regions of the parameter space where chaos occurs.

However, it is important to note that these results are primarily for the thermodynamic limit. For finite size networks, even bimodal natural frequency distributions can be chaotic. In those cases, there are rogue oscillators that do not synchronize with the rest of their cluster. These rogues can be treated as separate clusters, each of which requiring its own additional collective coordinate, increasing the number of active degrees of freedom, and opening up the possibility of chaos.

VII.2 Outlook

In our numerical simulations, we have observed regions in the parameter space of multimodal natural frequency distributions with four peaks that exhibit multistability, including natural frequency distributions that yield both strange attractors and limit cycles, depending on the initial condition. For example, Fig. 14 shows that for $K=0.95$ and multimodal distributions like Fig. 1, a second stable branch exists for $\sigma>0.022$ for the full Kuramoto model (1) with $N=100$ oscillators (green squares). This multistability is well reproduced by the collective coordinate model (19)–(20) with $N=100$ (red +’s) and by the collective coordinate model in the thermodynamic limit (21)–(22) (not shown). On the lower branch, the dynamics is periodic, and has the property that $r_{1}(t)=r_{4}(t+T/2)$ and $r_{2}(t)=r_{3}(t+T/2)$, where $T$ is the period of the system. On the upper stable branch there is no such relation between the cluster order parameters. Further study is required to understand this phenomenon, and the bifurcations that control it. Since the reduced collective coordinate models are more analytically tractable than the full Kuramoto model and accurately predict the existence of multistability, they may be used to provide deeper insight into this phenomenon.

Here we have considered all-to-all networks with synchronized clusters that result from distinct peaks in the natural frequency distribution. However, synchronized clusters can also occur due to the network topology. Future studies should consider whether topological clusters can yield chaos. Furthermore, chaos could result from a combination of frequency clustering and topological clustering. For example, a bimodal natural frequency distribution and a network with two clusters, which would result in four synchronized clusters of oscillators and the three degrees of freedom required for chaos.

Acknowledgements.

We wish to acknowledge support from the Australian Research Council, Grant No. DP180101991.