Dynamical equivalence between Kuramoto models with first- and higher-order coupling

TL;DR

This paper demonstrates that the Kuramoto model with simple higher-order coupling is dynamically equivalent to the standard model, allowing properties of the former to be derived from the latter, simplifying analysis of complex interactions.

Contribution

It establishes a dynamical equivalence between Kuramoto models with simple higher-order coupling and the standard model, providing a new perspective for analyzing complex oscillator interactions.

Findings

Higher-order coupling can be reduced to the standard Kuramoto model

Properties of higher-order models can be derived from the standard model

Simplifies analysis of complex coupled oscillator systems

Abstract

The Kuramoto model with high-order coupling has recently attracted some attention in the field of coupled oscillators in order, for instance, to describe clustering phenomena in sets of coupled agents. Instead of considering interactions given directly by the sine of oscillators' angle differences, the interaction is given by the sum of sines of integer multiples of these angle differences. This can be interpreted as a Fourier decomposition of a general $2{\pi}$-periodic interaction function. We show that in the case where only one multiple of the angle differences is considered, which we refer to as the "Kuramoto model with simple $q$th-order coupling," the system is dynamically equivalent to the original Kuramoto model. In other words, any property of the Kuramoto model with simple higher-order coupling can be recovered from the standard Kuramoto model.

Figures

Equations

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

$$\begin{array}{rl}{\pi }_q:{\mathbb{T}}^n& ⟶{\mathbb{T}}^n\\ ({\theta }_1,...,{\theta }_n)& ⟼(q{\theta }_1,...,q{\theta }_n)\mathrm{mod}2\pi ,\end{array}$$

$$\frac{\mathrm{d}}{\mathrm{d}t}{\pi }_q{\Theta }_{{\mathbf{\psi }}^{\ast }}$$

$$=q{\omega }_i-\frac{qK}{n}\sum _{j=1}^n\sin \left({\pi }_q{\Theta }_{{\mathbf{\psi }}^{\ast },i}-{\pi }_q{\Theta }_{{\mathbf{\psi }}^{\ast },j}\right),$$

$$q^{-1}{\mathbf{\eta }}^{\ast }+\frac{2\pi }{q}\mathbf{\rho },$$

$$\frac{\mathrm{d}}{\mathrm{d}t}$$

$$=q^{-1}\left[q{\omega }_i-\frac{qK}{n}\sum _{j=1}^n\sin \left({\Phi }_{{\mathbf{\eta }}^{\ast },i}-{\Phi }_{{\mathbf{\eta }}^{\ast },j}\right)\right]$$

$$={\omega }_i-\frac{K}{n}\sum _{j=1}^n\sin \left[q\left({\pi }_q^{-1}{\Phi }_{{\mathbf{\eta }}^{\ast },i}-{\pi }_q^{-1}{\Phi }_{{\mathbf{\eta }}^{\ast },j}\right)\right].$$

$${\pi }_q\left({\pi }_q^{-1}{\Phi }_{{\mathbf{\eta }}^{\ast }}\right)$$

$${\pi }_q^{-1}\left({\pi }_q{\Theta }_{{\mathbf{\psi }}^{\ast }}\right)$$

$$r_q(\mathbf{\theta })$$

$$n_1$$

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

Full text

Along the last decades, the Kuramoto model has attracted a lot of interest in the field of synchronization of coupled oscillators. Its simple formulation and the variety of synchronization phenomena that it describes make it a good candidate to investigate such phenomena both numerically and analytically. The sinusoidal interaction in the Kuramoto model can be seen as the first order of the Fourier decomposition of a general coupling function. It is then natural to extend the coupling to higher orders, such that the dynamics depend on the sines of multiples of the angle differences. In this manuscript, we show that considering only the angle difference or a unique integer multiple of it in the sinusoidal coupling is qualitatively equivalent.

I Introduction

In the context of synchronization of coupled dynamical systems, the Kuramoto model Kuramoto (1975, 1984) has drawn a lot of attention within the last decades. Strogatz (2000); Acebrón et al. (2005); Dörfler and Bullo (2014) Synchony is observed in many real systems, ranging from the brain’s oscillatory pacemaker cells establishing the circadian rythm, Lu et al. (2016) to synchronous machines connected to the high-voltage AC electrical grid. Bergen and Vittal (2000); Dörfler et al. (2013) This popularity brought the topic to a point where the remaining open questions are both hard and poorly rewarding to answer. To describe more realistic systems, some generalized versions of the Kuramoto model have been considered, as, for instance, meshed interaction graphs, higher order dynamics, Dörfler et al. (2013) or directed interactions. Restrepo et al. (2006); Delabays et al. (2019) One of these generalizations is to consider higher-order couplings. Hansel et al. (1993); Daido (1996); Skardal et al. (2011); Komarov and Pikovsky (2013); Li et al. (2014); Wang et al. (2017); Eydam and Wolfrum (2017); Li et al. (2019a) While the Kuramoto model is defined as

[TABLE]

for $i\in\{1,...,n\}$, where $\theta_{i}\in\mathbb{R}$ is the $i^{\rm th}$ oscillator’s angle, $\omega_{i}\in\mathbb{R}$ is its natural frequency, and $K>0$ is the coupling strength, the Kuramoto model with $q^{\rm th}$-order coupling, for $q\in\mathbb{Z}_{>0}$, is defined as

[TABLE]

for $i\in\{1,...,n\}$. The sum over $\ell$ can be seen as a truncated Fourier decomposition of a general $2\pi$-periodic coupling function.

To this day, most of the works about this version of the model has been limited to second-order couplings ($K_{1,2}\neq 0$), which already exhibits behaviors significantly different from the original Kuramoto model. Komarov and Pikovsky (2013); Li et al. (2014); Wang et al. (2017); Eydam and Wolfrum (2017) While most descriptions of the synchronous states of this model, for $q=2$ and large $n$, have been performed numerically, Li et al. (2014); Wang et al. (2017) an analytical approach, based on a self-consistency equation for the order parameters of the system, is given in Ref. Komarov and Pikovsky, 2013. The Kuramoto model with $q^{\rm th}$-order coupling has been used to describe clustering phenomena in systems of coupled synchronized oscillators. Clustered synchronous states of some particular versions of the Kuramoto model with second-order coupling and their dynamics are described in Refs. Hansel et al., 1993; Daido, 1996.

When one of the coupling orders largely dominates the others, a simplifying assumption is to consider the case where $K_{\ell}=0$ for all $\ell\neq q$, which we refer to as the Kuramoto model with simple $q^{\rm th}$-order coupling. Skardal et al. (2011); Li et al. (2019a) In this case, Eq. (2) reduces to

[TABLE]

for $i\in\{1,...,n\}$, which is the dynamical system considered in this manuscript. Note that the Kuramoto model, Eq. (1), can be seen as the Kuramoto model with simple first-order coupling. An analytical description of the transient dynamics of Eq. (3) and its synchronous clustered states are given in Ref. Skardal et al., 2011.

In this manuscript, we show that the dynamical systems described by Eqs. (1) and (3) are qualitatively equivalent, in the sense that both have the same dynamical properties (fixed points, linear stability, basin stability, order parameter) up to a projection from the state space to itself. Doing so, we make rigorous the claim below Eq. (3) in Ref. Li et al., 2014 that “In the case with $K_{1}=0$ (or $K_{2}=0$), the model is reduced to the original [Kuramoto model]…”. It also explains the striking similarity between Figure 2 in Ref. Li et al., 2019a (Kuramoto model with simple second-order coupling) and Figure 2 in Ref. Li et al., 2019b (Kuramoto model with simple first-order coupling). In Sec. II, we see that a direct relation can be drawn between the two models, allowing to translate any property of one model to the other. Therefore, a thorough investigation of the Kuramoto model with simple $q^{\rm th}$-order coupling is not needed, as any of its properties can be derived from properties of the original Kuramoto model. A selection of such properties and some implication for potential applications are detailed in Sec. III.

Remark. In Eq. (3), we consider all-to-all coupling, and we will do so for the whole manuscript for sake of readability. Nevertheless, our results can be straightforwardly extended to any coupling graph. We comment on this at the end of Sec. III.

II Dynamical equivalence

The Kuramoto model is usually considered as a dynamical system on the torus $\mathbb{T}^{n}$. We consider the variables $\theta_{i}$ as elements of $\mathbb{S}^{1}$, which we parametrize as $[0,2\pi)$ with periodic boundary conditions. We show now that the two following dynamical systems

[TABLE]

describe the same dynamics, up to rescaling the variables $\theta_{i}$ by a factor $q$. The rescaling is performed by the covering map from $\mathbb{T}^{n}$ to itself, which we define parametrically as

[TABLE]

where the modulo is applied elementwise.

More formally, we show below that $\pi_{q}$ sends any solution of Eq. (5) to a solution of Eq. (4), and that any lift of a solution of Eq. (4) through $\pi_{q}$ is a solution of Eq. (5). As $\pi_{q}$ is a smooth covering map, all properties of a solution of Eq. (5) are preserved in its projection [and similarly for the lift of a solution of Eq. (4)]. This is what we mean by dynamical equivalence.

II.1 Projecting

Let $\Theta_{\bm{\psi}^{*}}\colon\mathbb{R}\to\mathbb{T}^{n}$ be the solution of Eq. (5) with initial conditions $\bm{\psi}^{*}$. We verify that the projection $\pi_{q}\Theta_{\bm{\psi}^{*}}(t)$ solves Eq. (4),

[TABLE]

with initial conditions $\pi_{q}\Theta_{\bm{\psi}^{*}}(0)=q\bm{\psi}^{*}$.

II.2 Lifting up

The other way is a bit more intricate, because there are multiple preimages through $\pi_{q}^{-1}$ for each element of $\mathbb{T}^{n}$. Let $\Phi_{\bm{\eta}^{*}}\colon\mathbb{R}\to\mathbb{T}^{n}$ be the solution of Eq. (4) with initial conditions $\bm{\eta}^{*}\in\mathbb{T}^{n}$. The preimage $\pi_{q}^{-1}\bm{\eta}^{*}$ is a set of $q^{n}$ points, one of them being $q^{-1}\bm{\eta}^{*}\in[0,2\pi/q)^{n}\subset\mathbb{T}^{n}$, whose $i^{\rm th}$ component is $\eta_{i}/q$. The other $q^{n}-1$ can be constructed as

[TABLE]

with $\bm{\rho}\in\{0,1,...,q-1\}^{n}$. Each one of these points can be chosen as a starting point for the lifting.

Now, whichever element $\bm{\psi}^{*}\in\pi_{q}^{-1}\bm{\eta}^{*}$ we choose as starting point for the lifting, there is a unique smooth lifting $\pi_{q}^{-1}\Phi_{\bm{\eta}^{*}}(t)$ of the trajectory $\Phi_{\bm{\eta}^{*}}(t)$ satisfying the two following properties:

1. (i)

$\pi_{q}^{-1}\Phi_{\bm{\eta}^{*}}(0)=\bm{\psi}^{*}$; 2. (ii)

$\pi_{q}\left[\pi_{q}^{-1}\Phi_{\bm{\eta}^{*}}(t)\right]=\Phi_{\bm{\eta}^{*}}(t)$, for all $t\in\mathbb{R}$.

The time derivative of the $i^{\rm th}$ component of this lifting is

[TABLE]

The lifting then solves Eq. (5) with initial conditions $\bm{\psi}^{*}$, and this is true independently of the choice of representative $\bm{\psi}^{*}\in\pi_{q}^{-1}\bm{\eta}^{*}$. The preimage, by $\pi_{q}^{-1}$, of a solution $\Phi_{\bm{\eta}^{*}}$ of Eq. (4) is then a set of solutions of Eq. (5) differing from one another by a shift $2\pi\bm{\rho}/q$, with $\bm{\rho}\in\{0,1,...,q-1\}^{n}$.

II.3 Dynamical equivalence

It is now clear that

[TABLE]

and the unique smooth lifting of $\pi_{q}\Theta_{\bm{\psi}^{*}}$ such that $\pi_{q}^{-1}\left(\pi_{q}\Theta_{\bm{\psi}^{*}}\right)(0)=\bm{\psi}^{*}$ is exactly

[TABLE]

Properties of the solutions are then preserved by the projection $\pi_{q}$ as well as by its correponding lifting. We verify this for some dynamical properties in Sec. III.

III Consequences on clustering

Now that we established the dynamical equivalence between the Kuramoto model with simple first- and higher-order coupling, we review some results known for the original Kuramoto model and translate them in the Kuramoto model with $q^{\rm th}$-order coupling, in order to derive some results about clustering in the latter.

III.1 Fixed points

A given fixed point of the original Kuramoto model corresponds to $q^{n}$ fixed points of the Kuramoto model with $q^{\rm th}$-order coupling. Each of these fixed points differ by a shift $2\pi\bm{\rho}/q$, with $\bm{\rho}\in\{0,1,...,q-1\}^{n}$. When the natural frequencies are small ($\omega_{i}\ll K/n$), the synchronous state $\bm{\theta}^{*}\in[0,2\pi)^{n}$ of Eq. (4) is such that all angles are close to each other. For any $\bm{\rho}\in\{0,1,...,q-1\}^{n}$, the point $q^{-1}\bm{\theta}^{*}+2\pi\bm{\rho}/q\in[0,2\pi)^{n}$ is a synchronous state for Eq. (5). The integer vector $\bm{\rho}$ describes the clustering pattern of the corresponding synchronous state. If $\rho_{i}=\rho_{j}$, oscillators $i$ and $j$ are close to each other (at least for rather small natural frequencies), and approximately $2\pi/q$ appart from an oscillator $k$ such that $\rho_{k}=\rho_{i}\pm 1$. Then each oscillators with the same value in $\bm{\rho}$ form a cluster. The number of different values in $\bm{\rho}$ gives the total number of clusters in the synchronous state under consideration. From the point of view of the dynamics however, the clustering pattern has no effect. The $2\pi/q$ shifts introduced by the vector $\bm{\rho}$ leave Eq. (5) unchanged. The only information that is lost between the $q^{\rm th}$- and first-order coupling Kuramoto models is the clustering pattern. But this is only a combinatorial problem that can be addressed independently of dynamical considerations. In Fig. 1, for instance, each color corresponds to a different clustering pattern. At the end of the trajectory, the orange line has five oscillators in one cluster ($\{1,2,3,4,6\}$), with angles in $[0,\pi)$ and oscillator $5$ forming a cluster by itself, with its angle in $[\pi,2\pi)$. Similarly, for the green line, oscillators $1$ and $3$ are in one cluster with angles in $[0,\pi)$ and the others are in the cluster with angles in $[\pi,2\pi)$.

III.2 Linear stability

One can verify that the linear stability of Eq. (5) at a fixed point $\bm{\theta}^{*}$ is identical to the linear stability of Eq. (4) at $\pi_{q}\bm{\theta}^{*}$. More precisely, the Jacobian matrices ${\cal J}_{q}(\bm{\theta}^{*})$ of Eq. (5) at $\bm{\theta}^{*}$, and ${\cal J}_{1}(\pi_{q}\bm{\theta}^{*})$ of Eq. (4) at $\pi_{q}\bm{\theta}^{*}$ are equal. Thus, for a given distribution of natural frequencies, all clustered states have the same linear stability properties.

III.3 Order parameter

The order parameter $r_{1}$ is a quantity describing the level of coherence between oscillators’ angles in the Kuramoto model. It has been a major ingredient to analyze synchronization in this model. Kuramoto (1975); Aeyels and Rogge (2004); Mirollo and Strogatz (2005); Verwoerd and Mason (2008) To take clustering into account, in the Kuramoto model with higher-order coupling, it has been generalized Skardal et al. (2011) to the $q^{\rm th}$ order parameter

[TABLE]

It directly translates from the Kuramoto model with simple $q^{\rm th}$-order coupling to its equivalent Kuramoto model with first-order coupling. Namely, the order parameter $r_{q}(\bm{\theta}^{*})$ is equal to $r_{1}(\pi_{q}\bm{\theta}^{*})$. A large $q^{\rm th}$ order parameter $r_{q}$ indicates that the current state of Eq. (5) is clustered, but does not give any information about the clustering pattern, because it is blind to any angle shift of $2\pi/q$. The order parameter $r_{q}$ takes the same value $r=r_{1}(\bm{\psi}^{*})$ for each element of the preimage $\pi_{q}^{-1}\bm{\psi}^{*}$.

More generally, $r_{p}\approx 1$ for $1\leq p\leq q$ indicates that the system is clustered in $p$ equidistance clusters. Furthermore, if $r_{p}\approx 1$, it implies that $r_{kp}\approx 1$ for any $k\in\mathbb{N}$. Looking at all the order parameters with $1\leq p\leq q$ can give more information about the clustering pattern. For instance in the special case where $q=2$, if $r_{2}=1$, then one can verify that the number of oscillators in each cluster $n_{1}$ and $n_{2}$ respectively are given by

[TABLE]

In more general cases, however, it is not possible, as far as we can tell, to determine the number of oscillators in each cluster only based on the order parameters.

III.4 Synchronization

It is known Dörfler and Bullo (2011) for the Kuramoto model, Eq. (4), that if the coupling is sufficiently large to grant the existence of a synchronous state [$K>\max_{i,j}(\omega_{i}-\omega_{j})$], then there exists a value $\gamma_{\rm max}\in(\pi/2,\pi]$ such that the system synchronizes if all initial angles are in an arc of length at most $\gamma_{\rm max}$. In the Kuramoto model with simple $q^{\rm th}$-order coupling, this translates as follows. First, for the same value of $K$, the system synchronizes to the single cluster fixed point if all initial angles are in an arc of length $\gamma_{\rm max}/q$. Second, if the initial conditions $\bm{\psi}^{*}$ are such that all angles of $\bm{\psi}^{*}-2\pi\bm{\rho}^{*}/q$ are contained in an arc of length $\gamma_{\rm max}/q$, then the system synchronizes to the state with clusters given by $\bm{\rho}^{*}\in\{0,1,...,q-1\}^{n}$.

III.5 Basins of attraction

By equivalence of the dynamics, the basin of attraction of fixed point $\bm{\theta}^{*}$ of Eq. (5) is a copy of the basin of attraction of the fixed point $\pi_{q}\bm{\theta}^{*}$ of Eq. (4), rescaled by a factor $q^{-1}$. Namely, its volume is $q^{-n}$ times the volume of the basin of attraction of $\pi_{q}\bm{\theta}^{*}$. We illustrate this in Fig. 2, where we show the basins of attraction for the systems Eq. (4) (left panel) and Eq. (5) (right panel).

III.6 Basin escape

Suppose we introduce an additive noise in Eq. (5) to account for unpredictable perturbation of the environment. This will eventually lead our system to jump from a synchronous state to another. DeVille (2012); Tyloo et al. (2019) For the Kuramoto model, Eq. (4), such jumps bring the system from a synchronous state to a translate of itself, where some angles slip and accumulate integer multiples of $2\pi$. Lifting up such a trajectory to the Kuramoto model with simple $q^{\rm th}$-order coupling, the jumps then occur between the basins of attraction of different clustered states, where some angles accumulate an integer multiple of $2\pi/q$. Additive noise in Eq. (5) is then a possible mechanism for cluster formation.

If the noise has sufficiently small amplitude, the system remains for a long time in a neighborhood of a synchronous state, until the noise generates a sequence of perturbations that make it jump to another synchronous state. As pointed out by previous research on the Kuramoto model, DeVille (2012); Hindes and Schwartz (2018) the expected time between two jumps is exponential in (i) the inverse of the natural frequencies’ distribution width and (ii) the potential difference between the initial synchronous states and the closest $1$-saddle. It can also be related to (iii) the distance (in the state space) between the stable synchronous state and the closest $1$-saddle. Tyloo et al. (2019) By the equivalence derived in Sec. II, the expected time between jumps from a clustering pattern to another then follows the same exponential dependence (i)-(iii).

The equivalence is clearly seen in Fig. 3, where the angle trajectories in the two systems seem to be rescaled copies of each other, while they are obtained by two different simulations. The difference being that in the original Kuramoto model, the system jumps from the unique synchronous state to itself, by periodicity of the phase space, while in the Kuramoto model with $q^{\rm th}$-order coupling, the jumps occur between different clustered states.

III.7 Information storage

It has been acknowledged Kiss et al. (2005); Niyogi and English (2009) that the dynamics of some real-world systems can be described rather accurately by Eq. (5). In particular, it was proposed that such systems could be used as information storage, Ashwin and Borresen (2004); Skardal et al. (2011) where the value of each lower order parameter is a piece of information. Ref. Skardal et al., 2011 defines a forcing that allows to determine the clusterd state to which the system synchronizes and thus to control the value of $r_{p}$ for $1\leq p<q$. In this scope, the dynamical equivalence shown in Sec. II indicates that the stability properties of the system does not depend on its states. This is of major importance for such an application as it guarantees that the reliability of the storage system does not depend on the information it contains, which would render such an application much more complicated.

III.8 Generalization

More generally, instead of non-oriented, homogeneous, all-to-all couplings, our argument can be straightforwardly extended to the Kuramoto model with interactions given by any graph, weighted or not, directed or not. Our argument also shows dynamical equivalence between the two following, more general versions of the Kuramoto model with higher-order coupling:

[TABLE]

for $i\in\{1,...,n\}$ and $p\in\mathbb{Z}_{>0}$.

IV Conclusion

The main consequence of the dynamical equivalence presented in this manuscript is that any property of the original Kuramoto model, Eq. (4), can be lifted to the Kuramoto model with simple $q^{\rm th}$-order coupling, Eq. (5). The only discrepancy being the multiplicity of the elements of the lifting from Eq. (5) to Eq. (4).

To summarize, we showed that the Kuramoto model with first-order coupling is dynamically equivalent to the Kuramoto model with simple $q^{\rm th}$-order coupling. As a matter of fact, any dynamical property of the latter can be derived from the corresponding property of the original Kuramoto model. The behavior of the Kuramoto model with higher-order coupling qualitatively changes only if at least two different coupling orders are considered. Clustering occurs in the Kuramoto model with simple $q^{\rm th}$-order couping because of the choice of coupling function. But the dynamics are blind to the clustering pattern as each synchronous states is dynamically equivalent.

To take into account clustered states whose characteristics (linear stability, basin shape and size,…) differ, other models should be used. Some promising examples are, for instance, the more general Kuramoto model with $q^{\rm th}$-order coupling [Eq. (2)] or some dynamical systems with bounded confidence. Lorenz (2007)

Acknowledgments

This work has been supported by the Swiss National Science Foundation under grant 200020_182050.