Repulsively Coupled Kuramoto-Sakaguchi Phase Oscillators Ensemble Subject to Common Noise

TL;DR

This paper investigates the effects of common noise and repulsive coupling on identical Kuramoto-Sakaguchi oscillators, revealing that observed clustering in simulations is a numerical artifact and not a true dynamical state.

Contribution

The study demonstrates that cluster formation in this model is a numerical artifact, supported by analytical and numerical analysis, contrasting with previous claims of stable multiclusters.

Findings

Clustering observed in simulations is a numerical artifact.

Two-cluster states are shown to be non-attractive in the model.

Multiclusters can occur naturally in more general oscillators like Van der Pol.

Abstract

We consider the Kuramoto-Sakaguchi model of identical coupled phase oscillators with a common noisy forcing. While common noise always tends to synchronize the oscillators, a strong repulsive coupling prevents the fully synchronous state and leads to a nontrivial distribution of oscillator phases. In previous numerical simulations, a formation of stable multicluster states has been observed in this regime. However we argue here, that because identical phase oscillators in the Kuramoto-Sakaguchi model form a partially integrable system according to the Watanabe-Strogatz theory, the formation of clusters is impossible. Integrating with various time steps reveals that clustering is a numerical artifact, explained by the existence of higher order Fourier terms in the errors of the employed numerical integration schemes. Monitoring the induced change in certain integrals of motion we quantify these errors. We support these observations by showing, on the basis of the analysis of the corresponding Fokker-Planck equation, that two-cluster states are non-attractive. On the other hand, in ensembles of general limit cycle oscillators, such as Van der Pol oscillators, due to an anharmonic phase response function, as well as additional amplitude dynamics, multiclusters can occur naturally.

Figures

Equations

$${\dot{\phi }}_k=$$

$$+{\sigma }_1{\eta }_1(t)\sin {\phi }_k+{\sigma }_2{\eta }_2(t)\cos {\phi }_k.$$

$${\dot{\phi }}_k=\omega +\mathrm{Im}[(Z{e}^{i\gamma }+\sigma \xi )e^{-i{\phi }_k}],$$

$$\lambda =-\cos \gamma -\frac{{\sigma }^2}{2}.$$

$${\dot{\phi }}_k=\omega (t)+\text{Im}[H(t)e^{-i{\phi }_k}],k=1,\dots ,N$$

$$\mathcal{M}$$

$${\mathcal{M}}^{-1}$$

$$\dot{z}$$

$$\dot{\beta }$$

$${\dot{\Phi }}_1=$$

$$+\sigma \sin {\Phi }_1{\eta }_1(t)+\sigma \cos {\Phi }_1{\eta }_2(t)$$

$${\dot{\Phi }}_2=$$

$$+\sigma \sin {\Phi }_2{\eta }_1(t)+\sigma \cos {\Phi }_2{\eta }_2(t).$$

$$\dot{\delta }=$$

$$\displaystyle+\sigma\eta_{1}(t)\cos\Phi_{1}-\sigma\eta_{2}(t)\sin\Phi_{1}\big{]}~{}.$$

$${\lambda }_1=$$

$$+\sigma ⟨{\eta }_1(t)\cos {\Phi }_1⟩-\sigma ⟨{\eta }_2(t)\sin {\Phi }_1⟩,$$

$$\frac{\partial P(\Delta \Phi )}{\partial t}=-\frac{\partial }{\partial \Delta \Phi }\{[(p_2-p_1)\sin \gamma (1-\cos \Delta \Phi )$$

$$-\cos \gamma \sin \Delta \Phi ]P\}+{\sigma }^2\frac{{\partial }^2}{\partial \Delta {\Phi }^2}\left[(1-\cos \Delta \Phi )P\right].$$

$$P(\Delta \Phi )\sim \exp \left[\frac{\Delta \Phi (p_2-p_1)\sin \gamma }{{\sigma }^2}\right]{\left|\sin \frac{\Delta \Phi }{2}\right|}^{-2\left(\frac{\cos \gamma }{{\sigma }^2}+1\right)}.$$

$$P(\Delta \Phi )$$

$$\begin{array}{c}⟨{\eta }_2(t)\sin {\Phi }_1⟩=\int \delta (t-t^{\prime })⟨\frac{\delta \sin {\Phi }_1}{\delta {\eta }_2}⟩d{t}^{\prime }=\\ ⟨\frac{d(\sin {\Phi }_1)}{d{\Phi }_1}\frac{\delta {\Phi }_1}{\delta {\eta }_2}⟩=⟨\frac{1}{2}\sigma {\cos }^2{\Phi }_1⟩=\frac{\sigma }{4}.\end{array}$$

$${\lambda }_1=-p_1\cos \gamma -(1-p_1){\int }_0^{2\pi }\cos (\Delta \Phi +\gamma )P(\Delta \Phi )d\Delta \Phi -\frac{{\sigma }^2}{2}.$$

$$\lambda_{1}=\left\{\begin{array}[]{lr}p_{1}+p_{2}(\sigma^{2}-1)-\frac{\sigma^{2}}{2},&\sigma^{2}<2\;;\\ 1-\frac{\sigma^{2}}{2},&\sigma^{2}\geq 2\;.\end{array}\right.$$

$$\begin{array}{rl}{\lambda }_1+{\lambda }_2& =-\cos \gamma -{\sigma }^2-\cos \gamma {\int }_0^{2\pi }\cos \Delta \Phi P(\Delta \Phi )d\Delta \Phi \\ & +(p_2-p_1)\sin \gamma {\int }_0^{2\pi }\sin \Delta \Phi P(\Delta \Phi )d\Delta \Phi \end{array}$$

$$\mathcal{U}(\rho ,\Phi )=\frac{1}{N}\sum _j\log \frac{1-2\rho \cos ({\phi }_j-\Phi )+{\rho }^2}{1-{\rho }^2}.$$

$$\dot{\rho }=-{\mathcal{U}}_{\rho },\dot{\Phi }=-{\mathcal{U}}_{\Phi },$$

$$C_k=\frac{S_{k,k+2}}{S_{k,k+3}}\cdot \frac{S_{k+1,k+3}}{S_{k+1,k+2}},\text{ where }{S}_{ij}=\sin \frac{{\phi }_i-{\phi }_j}{2}.$$

$$e_k=S_{k,k+2}(t)S_{k+1,k+3}(t)S_{k,k+3}(0)S_{k+1,k+2}(0)-$$

$$S_{k,k+3}(t)S_{k+1,k+2}(t)S_{k,k+2}(0)S_{k+1,k+3}(0).$$

$${\text{Err}}_{\text{WS}}(t)=$$

$$\text{where }k=2,...N,$$

$${\text{Err}}_{\text{MMS}}(t)=$$

$$Z_k=\frac{1}{N}\sum _m{e}^{ik{\phi }_m}.$$

$$g_k(\vec{\phi })=-\frac{1}{2}\sum _m{f}_m(\vec{\phi })\frac{\partial }{\partial {\phi }_m}{f}_k(\vec{\phi }).$$

$${\dot{\phi }}_k$$

Full text

Coupled oscillators can synchronize, if the nature of coupling makes their phases attract each other (attractive coupling); or desynchronize, if there is repulsion of the phases (repulsive, or inhibitory coupling). This effect is well captured by the Kuramoto model of coupled phase oscillators. Another way to synchronize oscillators is to act on them with a common force – in particular, even common random noise will bring the phases together, an effect known as noise-induced synchronization. An interplay of repulsive coupling with common noise can be nontrivial, since the two effects act against each other. We argue in this paper that the Kuramoto model in such a case cannot form multiple clusters, due to its special integrability based on the Watanabe-Strogatz theory. The only possible end state for an evolution under Kuramoto-Sakaguchi model with repulsive coupling and weak noise is complete asynchrony of the oscillators. However, clusters can be observed in numerical simulations, because a discretization of the dynamics in general breaks the integrability. We also show that in a more realistic model like coupled Van der Pol oscillators, clusters can typically form due to the anharmonic phase response intrinsic to these oscillators.

I Introduction

The Kuramoto model, since its formulation in 1975 by Y.Kuramoto Kuramoto75 , has been vastly successful across scientific disciplines in describing naturally occurring phenomena in coupled oscillatory systems. On the one hand its simple mathematical form allows the analytical solution of a mean field theory for coupled oscillators with non-identical natural frequencies in the infinite system size limit Kuramoto75 ; kura86 . On the other hand, it is still able to generate complex behavior, e.g. chimeras chimera1 ; chimera2 , chaos kurachaos etc. The significance of the Kuramoto model to complexity science is akin to that of a simple model organism in the study of genetics, as its simple mathematical form nevertheless provides qualitative and quantitative insight into a variety of complex phenomena, especially as a paradigm for synchronization.

Since the original publications, there has been a plethora of literature dedicated to modelling systems by coupled Kuramoto oscillators. A particularly idealised and simple model is that of the globally coupled identical oscillators. For such a system a powerful Watanabe-Strogatz (WS) theory ws94 exists, according to which identical phase oscillators under identical forcing in their frequencies and first harmonics evolve under a Lie group belonging to the general class of Möbius group action 2009Chaos..19d3104M . Such systems possess a hidden low-dimensional dynamics. Moreover, it implies partial integrability, as it guarantees the existence of constants of motion, which are preserved throughout the dynamics ws94 . The system is only “partially” integrable because there are still 3 variables that are not constant and follow nontrivial dynamical equations. Closely related to the WS approach is the Ott-Antonsen theory OA2018 ; OA2019 , which allows an exact formulation of the mean field dynamics (with 2 degrees of freedom) for an infinite-sized population of non-identical oscillators.

Beyond the discovery of periodic components of a complex dynamics, the theoretical significance of the WS integrability is still not quite clear. For example, the integrability cannot be stated with certainty to exist in synchronous states, because in this case a cluster is in complete synchrony the same way as if the synchronized system had been “non-integrable”. However, in situations with partial synchrony, including the evolution towards aforementioned fully synchronized state, the WS integrability leads to a multiplicity of regimes due to the integrals of motion, and appear to have additional periodic components in the dynamics (cf. application of WS integrability to chimera states in Ref. PRzZ2, ). Recently, a perturbation theory based on the WS integrability has been suggested Vlasov-Rosenblum-Pikovsky-16 , which allows one to follow the evolution of the integrals under perturbations. One of the important consequences of the WS approach, which we explore below, is that it leads to the restriction on possible not fully synchronous states – namely, it excludes the formation of several clusters.

This study is inspired by the results of Gil et al. Mikhailov , where clustered states were observed in a set of identical oscillators under common multiplicative noise and repulsive coupling. In general, globally coupled identical oscillators can demonstrate configurations of different structural complexity: complete synchrony (one cluster state), partial synchrony (a nontrivial continuous distribution of phases, where all individual phases can be different), clusters (several groups of fully synchronized oscillators), chimeras (when one or several cluster coexist with partially synchronous oscillators), and solitary states solitarystate (when only one oscillator with a different phase exists apart from the fully synchronous cluster). Under strong repulsive coupling, a fully synchronized cluster becomes unstable, and it is not evident a priori which of these aforementioned configurations will be observed. Therefore, Gil et al. Mikhailov conducted numerical simulations and reported that common noise generally causes clustering in globally repulsively coupled interactions.

We call this claim into question based on the following reasons. First, clustering is indeed observed in some models of globally coupled identical phase oscillators cluster1 ; cluster2 ; syncbook , but always in situations with complex interaction functions, i.e, when the coupling term includes higher harmonics of the coupled phases. However, no such higher harmonics were present in the interaction term in the model proposed by Gil et al. Secondly, recent studies of the competition between common noise and repulsive coupling revealed non-trivial distributions for identical and non-identical oscillators Pimenova16 ; Mason , but no clustering has been observed. Gil et al., on the other hand, reported that for identical oscillators, clusters formed without a threshold, at any noise intensity. Finally, because the model used by Gil et al. can be fully described by the WS theory, there are restrictions due to the general properties of the Möbius transform governing the dynamics 2009Chaos..19d3104M , i.e. clusters are not allowed to appear (see Section III.1 below). Therefore, a thorough numerical and analytical study is needed to resolve the conflict between numerical findings Mikhailov and known theories.

In this paper we will thoroughly study the formation and stability of clusters in numerical experiments. After formulation of the problem (section II) we will present the WS approach and show that clustering is impossible. This conclusion will be further supported by an analytic and numerical investigation of the Lyapunov exponents of oscillators evolving in the field of two fully synchronized clusters (section III.2). In terms of numerics, the exact integrability is not preserved by the standard numerical schemes, both for deterministic and stochastic equations. In section IV we analyse errors in different numerical methods in terms of the change in the constants of motion which should be preserved by the Möbius group action, and show that this can lead to the formation of clusters as a numerical artifact. The WS approach is restricted to oscillators which couple to common external forces in their first harmonics. Furthermore, the Langevin equations must be understood in the Stratonovich interpretation in order for the usual differential calculus used in the WS approach to be applicable. Since we attribute the formation of clusters to the violation of WS integrability, we test this hypothesis in section IV.4 by including higher order terms in the phase dynamics, and in section V by studying repulsively coupled Van der Pol oscillators under common additive noise, for which higher order Fourier terms and multiplicative noise are naturally present in the phase reduced dynamics.

II Problem Formulation

We study a population of identical phase oscillators with phases $\{\varphi_{k}\}$ subject to a global coupling of Kuramoto-Sakaguchi type kura86 and common phase-dependent noise terms

[TABLE]

Here $\eta_{1,2}(t)$ are Gaussian white noise terms, with $\langle\eta_{i}(t)\rangle=0$, $\langle\eta_{i}(t)\eta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime})$ .

The parameters $\sigma_{1}$ and $\sigma_{2}$ parametrize the noise strengths for the two noise terms. The Langevin Eq.(1) is to be interpreted in the Stratonovich sense so that WS theory is applicable. Indeed, when interpreted in the sense of Itô calculus, for $\sigma_{1}\neq\sigma_{2}$ a noise-induced drift, i.e. Stratonovich shift, exists, and is proportional to the second harmonics in $\varphi_{k}$ which violates the conditions of WS theory. The phase shift parameter $\gamma$ parametrizes the degree of repulsion and attraction in the coupling term. In particular, when $\gamma=0$, the coupling is purely attractive, for $\gamma=\pi$, it is purely repulsive, and for $\gamma=\pi/2$ the coupling is neutral. Because it is always possible to rescale time, we assume, without loss of generality, that the phases couple with unit strength to the mean field in Eq.(1).

Models of type (1) have been analysed in Ref. Nagai-Kori-10, for the case of one noise term, and in Gil et al. Mikhailov for two noise terms. In the latter work it has been argued, that model (1) is the proper approximation after phase reduction, for a population of weakly coupled identical Stuart-Landau oscillators with a common additive noise term which is isotropic in the complex plane, making the phase equations invariant under rotation. In this case, $\sigma_{1}=\sigma_{2}=\sigma$ and one can rewrite (1) as

[TABLE]

where $Z=\frac{1}{N}\sum\limits_{j=1}^{N}{e^{{i}\varphi_{j}}}$ is the Kuramoto mean field and $\xi=-\eta_{1}+i\eta_{2}$ is complex Gaussian white noise. There exist well-known results on the stability of the completely synchronous cluster $\varphi_{1}=\varphi_{2}=\ldots=\varphi_{N}$ for model (2): To quantify the degree of stability of a fully synchronous cluster, which corresponds to $|Z|=1$, one calculates the transversal Lyapunov exponent (in previous literature also known as “evaporation” or “split Lyapunov exponent” Pimenova16 ), which describes the evolution of oscillator phases slightly deviated from the cluster. It is their average exponential rate of approach towards a cluster (or the rate of moving away from the cluster if the exponent is positive). The expression for this exponent is  Pimenova16

[TABLE]

For a negative Lyapunov exponent, complete synchronization is stable, i.e. it attracts nearby phases that are perturbed from it. According to this formula, for strong enough noise, the cluster of complete synchrony, with $|Z|=1$, is stable. For repulsive coupling, with $\cos\gamma<0$ and weak enough noise, the cluster is unstable.

While in Ref. Pimenova16, mainly the statistical properties of $|Z|$ have been analysed, Ref. Mikhailov, focuses on the occurrence of clusters, which are distinct attractive subgroups of oscillators with identical phases within the groups. The main goal of this paper is twofold: (i) to demonstrate that the occurrence of clusters in system (2) is impossible, and (ii) to identify numerical artefacts that may nevertheless lead to cluster formation in simulations.

III Theory

III.1 Application of the Watanabe-Strogatz theory

Our basic equation (2) in the Stratonovich interpretation belongs to a class of problems for which a theory, first developed by Watanabe and Strogatz  ws94 in 1994, which we shall call WS theory, is applicable. WS theory reduces the N-dimensional dynamics of a system of identically driven identical phase oscillators

[TABLE]

where $\omega(t)$ and $H(t)$ are arbitrary real and complex-valued functions of time, respectively, to a three-dimensional dynamical system preserving $N-3$ independent constants of motion. It is evident, that Eq. (2) belongs to class (3). One must stress here, that qualitative arguments below are applicable to any common force acting on the oscillators, not necessarily to the white Gaussian noise case mainly treated in this paper. It can be colored noise, or a chaotic/quasi-periodic/periodic force, or any combination of these functions of time.

At the heart of WS theory (see Refs. ws94, ; PR08, for a detailed presentation) is a coordinate transformation $\mathcal{M}$ formally belonging to the class of Möbius maps, which is a type of fractional linear transformation, mapping the complex unit circle one-to-one to itself. Specifically, $\mathcal{M}$ and its inverse $\mathcal{M}^{-1}$ can be written as

[TABLE]

Here $\{\varphi_{k}\}$ are the original phase coordinates, $z$ is a complex parameter of absolute value smaller than 1 and the parameter $\beta$ is a rotation angle. If the $\varphi_{k}$ evolve according to (3) and $z$ and $\beta$ evolve according to

[TABLE]

then the $\psi_{k}=\mathcal{M}^{-1}(\varphi_{k})$ are time independent constants of motion. The constants of motion are determined by the actual phases $\varphi_{k}$ and three time-dependent, real-valued parameters, the amplitude and angle of $z$ and by $\beta$. It is possible to impose conditions on the constants which make the Möbius transform unique. For instance, if no majority cluster exists one can choose $z$ and $\beta$ such that $\langle\exp(i\psi_{k})\rangle=0$ and $\textrm{arg}\left(\langle\exp(i2\psi_{k})\rangle\right)=0$ [ws94, ]. The Möbius transform (4) consists essentially of a common rotation of the angles $\psi_{n}$ by the angle $\beta$ and a subsequent contraction along the circle into the direction of the angle of $z$. Indeed, $|z|$ can be used as a measure of synchronization akin to the Kuramoto order parameter as both become equal to unity at complete synchronization. The existence of the constants of motion implies the system is integrable. However, it must be stressed that WS integrability only holds for phase oscillator models of the form Eq. (3), e.g. the Kuramoto-Sakaguchi model for globally coupled phase oscillators or the Winfree model winfree with a harmonic phase response function.

Equations (6) in principle allow for a numerical integration of the system which automatically conserves all the constants $\psi_{k}$. However, typically the forcing term $H(t)$ contains the Kuramoto order parameter $Z=\langle\exp(i\varphi_{k})\rangle_{k}$ for which the $\varphi_{k}=\mathcal{M}(\psi_{k})$ must be calculated at every time step. Only for constants of motion which are uniformly distributed on the unit circle the parameter $z$ approximates the actual mean field $Z$ to the order $1/\sqrt{N}$ where $N$ is the number of oscillators PR08 PRzZ2 .

Expressions (4) and (5) alone already rule out the formation of clustered states from non-clustered initial conditions. Two oscillators with initially different phases are mapped to a single point only if $|z|\to 1$, in which case $z_{k}=\exp({i}\varphi_{k})$ tends to $z$ for all points of the circle except a singular “solitary state” point with $z+\exp\left(i(\psi_{k}+\beta)\right)=0$ [solitarystate, ] for which the Möbius transform is not unique at $|z|\to 1$. Therefore, only a single cluster attracting different phases can exist at a time. Nevertheless, this only prohibits the formation of multiple clusters but not their existence under this model. There is not restriction for oscillators to be in one or several distinct clusters with identical phases within each cluster and stay in that configuration.

III.2 Linear stability analysis of a two-cluster state

In section III.1 above we have shown that clusters cannot appear from non-clustered initial conditions. The same arguments can be applied when the dynamics evolves from an initial multicluster state. Here, either the multicluster remains with the same partition, or the fully synchronized state with maximally one additional cluster or oscillator in a solitary state appears. Imperfect clusters, i.e. configurations with phases very close to one another, can also dissolve or contract depending on their dynamical stability. It is therefore instructive to look on the stability of the multicluster states. According to the WS theory, one expects that not more than one of the clusters can be asymptotically attracting. Otherwise multiclusters would also form from non-clustered initial conditions, a phenomenon which is forbidden by the argument in section III.1. In this section we provide a linear stability analysis of the two-cluster state under the stochastic evolution given by model (2) which confirms our expectation. We write equations (2) for a two-cluster state as

[TABLE]

Here $\Phi_{1}$ and $\Phi_{2}$ are the phases of the two clusters. $\Delta\Phi=\Phi_{2}-\Phi_{1}$ is their phase difference. Parameters $p_{1}$ and $p_{2}=1-p_{1}$ are their relative population sizes.

To evaluate the stability of the two-cluster state in terms of a small perturbation from one of the clusters, say, cluster 1, we perturb two oscillators belonging to cluster 1 by pulling them by a small amount in opposite directions away from the cluster, i.e. $\varphi_{1,2}=\Phi_{1}\pm\delta$. For small $\delta$, linearization yields

[TABLE]

This allows us to express stability of cluster 1 via the split/evaporation Lyapunov exponent splitLE_Kaneko-89 ; splitLE_Pikovsky2001b ; splitLE_pikovsky_politi_2016 as

[TABLE]

where $\langle\cdot\rangle$ indicates time average, which in this case also equals to the ensemble average, because as we will see in (11), the probability distribution of the phase $\Delta\Phi$ is stationary.

While the Stratonovich shift for the Langevin Eqs.(7) happens to be zero, that is not the case anymore when Eq.(8) is considered as well. It is important to keep this in mind and choose a correct integration scheme when Eqs.(7) and (8) are integrated numerically. To calculate the Lyapunov exponent analytically, we need to know the probability distribution of $\Delta\Phi$ and the averages $\langle\eta_{1}(t)\cos\Phi_{1}\rangle,\;\langle\eta_{2}(t)\sin\Phi_{1}\rangle$. First, we write a two-dimensional Fokker-Planck equation for $\Phi_{1}$ and $\Delta\Phi$ corresponding to the Langevin equations (7) under Stratonovich interpretation HANGGI1982207 . Then, integrating the joint density of $\Phi_{1}$ and $\Delta\Phi$ over the variable $\Phi_{1}$, using the fact that the probability distribution of $\Phi_{1}$ is rotational symmetric, i.e. uniform, we obtain a closed equation for the probability distribution $P(\Delta\Phi)$ of the phase difference

[TABLE]

Note that this probability density function is defined and restricted on the open interval $(0,2\pi)$ since the two clusters cannot cross each other. The stationary solution has the form

[TABLE]

A closed expression for the normalized probability density is possible when $\gamma=\pi$, i.e. the repulsion between the oscillators is maximal. In this case,

[TABLE]

where $B(x,y)$ is the Beta function. The shape of the probability density function in the general case (11) is shown in Fig. 1.

In the expression (11), a nonzero phase shift parameter $\gamma$ introduces a curious asymmetry in form of the exponential factor which is not a periodic function, i.e. when we consider the distribution on $(0,2\pi)$ wrapped around the circle, the derivative is not continuous at the singular absorbing point $\Delta\Phi=0=2\pi$. The critical noise strength $\sigma_{c}$ beyond which the two clusters are synchronized to become one cluster, corresponds to the case where $P(\Delta\Phi)$ becomes a delta function $\delta(\Delta\Phi)$. Formally, this corresponds to divergence of the integral of the probability density (11). This happens if the exponent of $|\sin(\Delta\Phi/2)|$ is smaller than $-1$, and from this we can calculate the critical noise strength to be $\sigma_{c}^{2}/2=-\cos\gamma$.

In addition to the distribution of the phase difference $\Delta\Phi$, one needs to calculate the averages $\langle\eta_{1}(t)\cos\Phi_{1}\rangle$ and $\langle\eta_{2}(t)\sin\Phi_{1}\rangle$ to evaluate (9). Since $\eta_{1}(t)$ and $\eta_{2}(t)$ are independent Gaussian white noise processes and $\Phi_{1}[\eta_{1},\eta_{2}]$ is a functional of both $\eta_{1}$ and $\eta_{2}$ this can be accomplished by virtue of the Furutsu-Novikov formula Furutsu63 ; Novikov64

[TABLE]

Similarly $\langle\eta_{1}(t)\cos\Phi_{1}\rangle=-\frac{\sigma}{4}$. A general expression for the Lyapunov exponent $\lambda_{1}$, which describes the stability of the cluster 1, is therefore

[TABLE]

Lyapunov exponent (14) can even be analytically represented for the case $\gamma=\pi$

[TABLE]

Exchanging $p_{1}$ and $p_{2}$, we obtain the Lyapunov exponent $\lambda_{2}$ of the other cluster. From this special case one can easily see that when $\sigma^{2}<2$, i.e. when a fully synchronized one-cluster state is unstable and the two-cluster Lyapunov exponents are well defined, they satisfy $\lambda_{1}+\lambda_{2}=0$.

Through direct numerical evaluation of the Lyapunov exponent $\lambda_{1}$ in Fig. 2 , we obtain a confirmation of the above analytical result.

For the case of Kuramoto-Sakaguchi model with general phase shift parameter $\gamma$, using (9) and the corresponding expression for $\lambda_{2}$, we obtain for the sum

[TABLE]

After applying integration by parts for the product of three functions, $\langle\sin\Delta\Phi\rangle$ can be written in terms of $\langle\cos^{2}(\Delta\Phi/2)\rangle$ and $\langle\cos\Delta\Phi\rangle$. After simple algebra, the relation $\lambda_{1}+\lambda_{2}=0$ for the generic Kuramoto-Sakaguchi model can be shown. This means that for two narrowly distributed groups of repulsively coupled oscillators with common multiplicative noise, the larger group will dissolve while the smaller group is attractive (Fig. 2). Simultaneous attraction into two clusters is not possible.

IV Numerics

As we have argued above, the WS theory of integrability prevents a multicluster state from ever occurring in model (2), and a linear stability analysis via Fokker-Planck formulation has confirmed it in the case of two clusters in the repulsive regime. But the observation of clustering in simulations by Gil et al. Mikhailov may be attributed to numerical artefacts, as one cannot expect WS integrability to be preserved by standard numerical methods. In this section we explore how different numerical integration methods for integrating deterministic and stochastic equations affect WS integrability and clustering. First we discuss methods which quantify the errors occurred from a deviation from integrability and measure the degree of clustering.

IV.1 Numerical evaluation of the constants of motion

As has been outlined in section III.1, the constants of motion of the system can be determined via the Möbius transformation (4) of the $N$ phases $\{\varphi_{j}\}$. In practice, one must first determine the complex variable $z=\rho e^{{i}\Phi}$ which characterises the transformation. Watanabe and Strogatz (see section 4.2 in Ref. ws94, ) have shown that this can be accomplished with the help of a potential function

[TABLE]

The proper value of $z$ corresponds to the minimum of this function with respect to its modulus $\rho$ and to its argument $\Phi$. The easiest way to determine the minimum is by integrating

[TABLE]

until the steady state is established. The angles $\psi_{k}+\beta$ can then be obtained with the Möbius transformation (4). To avoid determination of the value of $\beta$, it is convenient to consider only the differences $\psi_{j}-\psi_{1}$, $j=2,\ldots,N$ as constants of motion. The disadvantage of this method lies in the necessity of solving the minimization problem for the potential (17), which can be performed with finite accuracy only. There exists, however, another possibility to determine the constants of motion.

Marvel et al. (see Sec. V in Ref. 2009Chaos..19d3104M, ) have demonstrated that the cross ratio of four complex numbers on the unit circle is a preserved quantity under the Möbius transformation. For any four phases $\varphi_{k},\varphi_{k+1},\varphi_{k+2},\varphi_{k+3}$ (not necessarily ordered on the circle), the constant of motion $C_{k}$ is defined as

[TABLE]

Our method of checking the conservation of these quantities is based on (LABEL:eq:constantMMS3), but we find it appropriate to avoid calculating fractions, because as the phases synchronize, the denominators can be very small or vanish. Instead we calculate the errors of the form

[TABLE]

In summary, we test for integrability in numerical schemes by calculating the following errors containing changes in the conserved quantities under Möbius action

[TABLE]

IV.2 Numerical evaluation of clustering

In numerical simulations of model (2) (details shall be outlined below) we may observe different clustered states, illustrated in Fig. 3.

We quantify the formation of synchronized clusters with the help of the Kuramoto-Daido mean fields

[TABLE]

After long integration time, the first order mean field $Z_{1}$ for repulsive coupling is either small if noise is present, or vanishes completely in the deterministic case. The second order parameter $R_{2}=|Z_{2}|$ is maximal and equal to 1 for two fully synchronized clusters of arbitrary sizes with phase difference $\pi$ between them. Altogether, the degree of the formation of two clusters can be measured by a growth of $R_{2}$ approaching values close to one. We will henceforth use the evolution of $R_{2}$ as an indication for a two-cluster state.

IV.3 Deterministic evolution

We first explore how well the WS integrability is preserved in numerical simulation of deterministic equations. Here the original Kuramoto-Sakaguchi model is not optimal. After a short initial transient, the evolution of $R_{1}=|Z_{1}|$ effectively comes to a halt as soon as a steady state is reached, i.e. $R_{1}=|Z_{1}|$ is zero for repulsive coupling or unity for attractive coupling. Instead, we integrate a model of type (3) with a prescribed modulated time-dependent forcing $\omega(t)=0.2\sin(1.752t)$, $H(t)=0.4\cos(2.33t)\cdot Z$, and $N=10$, designed to ensure the state remains nontrivial (see Fig. 4). Integrating this deterministic equation, we use the standard Runge-Kutta method of 4th order (RK4) and the first-order Euler method.

First, comparing Fig. 4 panels (a) and (b), where the two methods (19) and (20) of determining the constants of motion are used, we can conclude that, while the errors in the constants of motion are similar for large steps, the calculation of the errors via $\text{Err}_{\text{WS}}$ (19) does not allow for a proper estimation of very small errors, due to the necessity of a minimization procedure which can be performed only with finite precision. Therefore, for the rest of the paper we calculate the errors using only $\text{Err}_{\text{MMS}}$ (20).

The second observation is that in all the cases the errors grow in time roughly linear, with prefactors depending on the integration step $h$: $\text{Err}_{\text{MMS,RK4}}\sim h^{4.94}t$ for the RK4 method, and $\text{Err}_{\text{Euler}}\sim h^{0.99}t$ for the Euler method, indicating a drift of the constants of motion. This is consistent with the fact that RK4 makes an error of $h^{5}$ in each time step and for Euler it is $h$.

In Fig. 5 we present the results for the integration of the deterministic equation (1) with $\omega=\sigma_{1}=\sigma_{2}=0$, $N=100$ and $\gamma=0.54\pi$ (slightly repulsive). Here we use rather large integration steps to make the clustering effect visible during a relatively short transient time interval, before the main order parameter becomes very small and the dynamics stops. One can see that for $h>0.6$ the order parameter $R_{2}$, which measures formation of a two-cluster state, grows to macroscopic values.

For instructive purposes, we explore here which type of perturbations are introduced by the numerical integration methods to the original dynamics (1) when noise is not present. The simplest case is to estimate the perturbations introduced by the Euler method. The Euler method models a continuous-time dynamical system $\dot{\varphi}_{k}=f_{k}(\vec{\varphi})$ up to the order $h$ as a map $\varphi_{k}(t+h)=\varphi_{k}(t)+hf_{k}(\vec{\varphi}(t))$. Then we might ask, what continuous equation is integrated by the same map correctly up to the order $h^{2}$. Looking for this equation in the form of $\dot{\varphi}_{k}=f_{k}(\vec{\varphi})+hg_{k}(\vec{\varphi})$, we find

[TABLE]

Substituting for the Kuramoto-Sakaguchi model $f_{k}=\omega+\textrm{Im}\left[Ze^{{i}\left(\gamma-\varphi_{k}\right)}\right]$, we obtain a modified equation where the error in the Euler integration is part of the dynamics

[TABLE]

Here $Z_{2}=\left\langle\exp\left(i2\varphi_{m}\right)\right\rangle_{m}$ is the second Kuramoto-Daido mean field. One can see, that in addition to the new coupling terms proportional to $\sin\varphi$ or $\cos\varphi$, which do not violate the WS integrability, terms proportional to $\sin 2\varphi,\;\cos 2\varphi$ appear, which violate the WS integrability and may result in the formation of two clusters. Clusters of higher orders can presumably be attributed to terms $\sim\sin 3\varphi$ etc. appearing in higher order errors in $h$.

IV.4 Stochastic evolution

Throughout this paper we interpret the stochastic system (1) in the Stratonovich sense. However, the additional drift term needed to transform it into Itô interpretation, vanishes for the case when the two noise terms correspond to an isotropic noise in the complex plane, i.e. $\sigma_{1}=\sigma_{2}$. Therefore, the numerical schemes for both Stratonovich and Itô interpretations can be used to integrate the phases in the case of two noise terms of equal strength. On the other hand, the two noise terms in (1) do not commute. It has been shown, e.g. in Refs. Burrage, ; BURRAGE1998161, ; MR2052268, , that the strong order of convergence of all higher order integration methods for stochastic differential equations with non-commutative noise cannot be higher than 0.5. Strong order of convergence is defined by the average error made by the time-discretized approximation of the stochastic integration scheme in approximating each individual path of the continuous-time process. Therefore we restrict ourselves to using only low-order integration schemes. Only in the case of noise in a single direction (i.e. $\sigma_{2}=0$), high-order methods like the stochastic Runge-Kutta method Burrage are used. In Fig. 6 we show the results in the case of two relatively strong noise terms $\sigma_{1}=\sigma_{2}=0.1$ . The integration is performed with the Euler-Heun scheme for different time steps $h$. Due to the rotational invariance preserved by two noise terms of equal strength, we have set $\omega=0$. One can clearly see the formation of two clusters, indicated by values of $R_{2}$ growing close to one, on a time scale $\sim h^{-1}$. For a weaker noise $\sigma=0.01$, clustering appears much slower. Only initial growth of the second order parameter can be observed at the maximal integration time of $T=4\times 10^{5}$. This dependence of the clustering time scale on the integration step size demonstrates that clustering in this system is a numerical artifact.

In fact, when we break WS integrability by including a term in the stochastic dynamics (2) proportional to the error in the deterministic integration scheme as in Eq.(IV.3), i.e.

[TABLE]

we observe robust clustering under dynamics Eq. (IV.4) at a similar time scale as in the original system (2) for an integration time step of $h=\epsilon$ (see Fig.7).

In Fig. 8, we compare different integration schemes applied to models with one or two noise terms. Here for the cases of two noise terms (like in Fig. 6) and of one noise (where we set $\omega=10$ because the rotational symmetry is broken), we present results for the Euler-Heun scheme, suitable for the Stratonovich interpretation of the stochastic differential equation. Additionally, we show the results of the stochastic Runge-Kutta scheme (SRK), which is suitable for the one-noise case only, because of the non-commutativity of the two noise terms mentioned above. One can see that all plots are qualitatively similar, with only some quantitative differences. As one would expect, the conservation of the constants of motion under the SRK scheme is the best, and here also the growth of the second order parameter is rather weak on the chosen time interval $t<2\cdot 10^{5}$. We have also performed simulations with the Euler-Maruyama scheme with Stratonovich shift for the model (1) in the Stratonovich interpretation, both with one noise term and with two noise terms (where the Stratonovich shift is zero), all of which yield quantitatively identical results to the Euler-Heun scheme and are therefore not shown.

We can conclude this section by stating that in general, numerical schemes do not conserve the integrals of motion of the system, and eventually may lead to formation of clusters. Because the methods for integrating stochastic differential equations have typically lower order than the deterministic ones, the clusters may be more easily observed in the integration of noisy equations. In the deterministic case, clustering may not be observed at all if a zero mean field steady state is reached first. The presence of noise can prolong the time within which the constants of motion continue to drift and their deviation from their initial values continues to grow until at some point fully synchronized multiclusters are formed.

As mentioned in section III, the best way to avoid the numerical artefacts of clustering is to integrate the Watanabe-Strogatz equations (6), but to accomplish this one has to perform multiple Möbius transforms at each time step for a full time series of the mean field, which may be quite computationally expensive. Furthermore, discretization errors are still present in integrating the low-dimensional dynamics of the Möbius group parameters. Only multicluster formation would be guaranteed to no longer occur.

V Oscillators with naturally occurring clusters under repulsive coupling

and common noise - the Van der Pol oscillators

Unlike the Kuramoto model, more realistic oscillator models such as the Van der Pol oscillator have limit cycles which intrinsically contain higher order Fourier terms and additional amplitude dynamics. Under common additive noise and repulsive coupling, formation of multiclusters is no longer forbidden and could now naturally occur. We consider $N$ identical repulsively all-to-all coupled Van der Pol oscillators subject to additive common Gaussian white noise in one direction

[TABLE]

Here $b>0$ is the repulsive coupling strength, $a$ parametrizes the nonlinearity of the Van der Pol oscillators, $\sigma$ is the noise strength, and $\xi(t)$ is a Gaussian white noise force. Using phase reductionwinfree the additive noise term will become multiplicative with the linear phase response function as a factor.

With a similar approach to that of section III.2, one can determine the Lyapunov exponents for the two-cluster state in (24) numerically by integrating a perturbation from one of the clusters in the linearized dynamics of the two cluster system. Contrary to the case of the Kuramoto model, presented in Fig. 2, now in Fig. 9 we see that the two-cluster state with $p_{1}\approx p_{2}$ is locally stable, which is confirmed in Fig. 10 by direct simulations of Eq.(24). Here, we defined the Kuramoto order parameters using the phases defined by virtue of Poincaré sections. One oscillator has been chosen as a reference, and the moments of time $\overline{t}_{1},\overline{t}_{2},\ldots$ at which it crosses half-line $(x>0,y=0)$ have been determined. Then the phase differences of all other oscillators to the reference one at time $\overline{t}_{m}$ were defined as $2\pi(t_{m}^{(k)}-\overline{t}_{m})/(\overline{t}_{m+1}-\overline{t}_{m})$. Here $t_{m}^{(k)}-\overline{t}_{m}$ is the time needed for an oscillator with index $k$ to reach the Poincaré section from its position at time $\overline{t}_{m}$.

In general, clustering strongly depends on the level of nonlinearity, described by $a$. For large $a$ ($a=1$, $b=0.01$), in the deterministic case three clusters can be observed. In the presence of noise, the picture is not so clear as several different cluster states may appear depending on the realization of initial conditions and of the noise, but a tendency at least to temporal formation of clusters is clearly observed. For small values of nonlinearity parameter $a$, typically non-clustered states are observed both with and without noise. This is to be expected, since the Van der Pol oscillator with a weakly anharmonic limit cycle has comparably much smaller higher order Fourier terms in its phase response function. In general, dynamic complexity of systems like (24) with non negligible amplitude dynamics can be very high, with chimera-like states becoming possible (i.e. where clusters coexist with dispersed elements), and a full characterization is beyond the scope of this paper.

From the above observation we can therefore conclude that there exists a qualitative difference between the dynamics of the phase oscillator model (e.g. Fig. 2), and that of the more general oscillator model with additional amplitude dynamics (e.g. Fig. 9), specifically under a repulsive coupling and common noise: while under Kuramoto-Sakaguchi model clusters are not allowed to form, under Van der Pol model they are naturally forming and are stablized by the common stochastic forces.

VI Conclusions

In this study, we apply the Watanabe-Strogatz (WS) theory ws94 to the Kuramoto-Sakaguchi model of repulsively coupled phase oscillators under common noise, studied previously in Ref. Mikhailov,  . Our main result is that although both the WS theory and the stability analysis of clustered states exclude appearance of clusters as observed in Mikhailov, , these observations can be generally explained as artefacts from the finite accuracy of numerical simulations. The correct long term behavior for repulsively coupled phase oscillators under common noise is either an incoherent state with no clustering (when the common noise has weaker effect compared to the repulsive coupling) or a completely coherent state (when the common noise has a stronger effect compared to repulsive coupling). We study the numerical errors of different deterministic and stochastic schemes by monitoring the evolution of the constants of motion which are conserved under the exact dynamics. It should be stressed that the conclusions of WS theory only apply to a restricted class of phase oscillators which approximate weakly coupled, weakly nonlinear limit cycle oscillators. Violation of WS integrability occurs naturally in general coupled oscillator systems. We show that in the case of repulsively coupled Van der Pol oscillators noise-induced or deterministic clustering can indeed be easily observed in regimes of larger nonlinearity. Due to the limitation of the Kuramoto-Sakaguchi system in describing real-world oscillator models or even more complicated coupled systems of differential equations, in terms of numerics, this paper presents only a cautionary tale. For most types of high dimensional coupled differential equations, a hidden low-dimensional dynamics such as present in Kuramoto-type system is not available, nor do there often exist integrals of motion. For these systems, often the only way to measure or to gauge numerical errors is by using integration steps which are as small as possible, and to compare results under various degrees of numerical accuracies.

Acknowledgments

This paper is developed within the scope of the IRTG 1740/TRP 2015/50122-0, funded by the DFG/ FAPESP. We thank Denis S. Goldobin and Michael Zaks for valuable discussions. C.Z. also acknowledges the financial support from China Scholarship Council (CSC). Work of A.P. is supported by Russian Science Foundation (Grant Nr. 17-12-01534).