Synchronization of coupled noisy oscillators: Coarse-graining from continuous to discrete phases

We consider the parameter range $J\gg K$ and $J\gg\eta$, so that the time scale of the internal dynamics of each oscillator dominates over (i.e. is shorter than) that of the interactions between the oscillators. Then, after a fast transient defined by the internal dynamics, we have that $\left|A_{s}\right|\approx 1$. After that, the phase of each oscillator varies as a function of time on a slower time scale defined by the interactions (albeit with very rapid fluctuations). On this longer time scale we can write $A_{s}\approx e^{i\phi_{s}}$. The dynamics specified by Eq. (1) can then be reduced to the phase equation Kuramoto2

Here we have defined a Kuramoto order parameter $R$ which follows directly from Eq. (3),

from which we extract the real phase variable $\psi$, and

where $r\in[0,1]$ and $\phi_{s}\in[0,2\pi]$ are also real. The noise $\xi_{s}(t)=-\sin\phi_{s}\zeta_{R}^{s}\left(t\right)+\cos\phi_{s}\zeta_{I}^{s}\left(t\right)$ is again Gaussian and white, with zero mean and correlation function that follows directly from Eqs. (5) and (6),

The order parameter $R$ can be written as

where we have introduced the density of oscillators with phase $\phi$,

In the thermodynamic limit $N\rightarrow\infty$,

where $\rho\left(\phi,t\right)d\phi$ is the probability that the phase of an oscillator lies in the interval $[\phi,\phi+d\phi]$ at time $t$.

In the thermodynamic limit, the stochastic phase equation (8) (which is an equation of the Langevin form) can be replaced by a nonlinear Fokker-Planck equation (see Ref. Acebron for a detailed derivation using the path integral formalism),

where the second derivative term on the right is the diffusion term, and where the drift contains

with

Since all of our oscillators have the same frequency ($\omega=0$ in the moving frame), synchronization in this framework corresponds to the tendency of the oscillators to have the same phase. The desynchronized state corresponds to a uniform distribution of phases,

This choice obeys the normalization condition

The dynamics described by the phase equation (8) contains two competing trends: the fluctuations, which tend to desynchronize the system and stabilize the state described in Eq. (18), and the attractive interactions that tend to synchronize the system. When the coupling strength is weaker than a critical value ($K<K_{c}$), the desynchronized state is stable, while an interaction stronger than this critical value ($K>K_{c}$) destabilizes the desynchronized state, leading to a phase transition to synchronization. Calculations of the critical point have been widely documented in the literature Acebron ; Kuramoto2 , and we see no point in repeating them here. Instead, here we consider the continuous limit of the coarse-graining procedure presented in Sec. III and detailed in Appendix A. This limit is in turn detailed in Appendix B, where we show that this critical value is given by

As would be expected, stronger fluctuations require a stronger interaction for synchronzation to occur, that is, the critical value of the coupling parameter increases with increasing noise intensity $\eta$.

In the thermodynamic limit the absence of synchronization is associated with a vanishing order parameter, $R=0$ . Away from the thermodynamic limit $R$ fluctuates to values away from zero because of the finite number of oscillators, but remains small as long as $N$ is not exceedingly small. The occurrence of self-organization is characterized by a non-vanishing order parameter. Near the critical point (20), the mean-field evolution of the order parameter can be approached via the normal form (see Refs.Acebron ; Kuramoto2 ; Bonilla and Appendix B)

with

where $f^{\prime}(0)\equiv\left.(df(r)/dr)\right|_{r=0}$. For the original Kuramoto model Kuramoto1 , $f=1$, the transition to synchronization is always supercritical because $\beta>0$. The presence of a nonlinear interaction $f(r)$ might induce a change in sign of $\beta$ and might therefore lead to a subcritical transition. The function $f(r)$ affects the critical coupling strength in Eq. (20) simply as a renormalization of this critical value. The main role of $f(r)$ is to increase or decrease the dispersion of the phases around a common value in the synchronized regime. The qualitative effect of including the function $f$ can be seen in the interaction term in the phase equation (8). In the synchronized phase a function $f(r)$ that increases with increasing $r$ increases the attractive interactions among the oscillators, leading to a decrease of the dispersion of phases around a common value. In contrast, if $f(r)$ decreases with increasing $r$, this tendency will be weakened and, since fluctuations lead to a greater dispersion of the phases, the synchronization will be coarser. A decreasing function might thus be a better candidate for a successful increasingly coarse coarse-graining procedure.

For numerical simulations we must choose a particular form for the function $f(r)$. As stated above, our goal of coarse-graining the phase space is better achieved if we choose a decreasing function of $r$. We have tested a number of different functions and have determined that the particular form and parameters in it are not as important as is simply the choice of almost any decreasing function of $r$. We use the exponentially decreasing function

with $a$ a positive parameter. Note that when $a\rightarrow\infty$, $f\rightarrow 1$ and we recover the standard Kuramoto model.

In the upper panel of Fig.1 we show, after a transient, the state of an ensemble of $5000$ interacting oscillators governed by Eq. (1). The points represent the complex amplitude of each oscillator in the complex plane. As expected for large $J$ ($J=200$ in all our simulations), after a transient the amplitudes settle around a circle of unit radius (red line). We have chosen a coupling strength $K\gg K_{c}$ (see caption) so that the array of oscillators is in the synchronized regime. The system exhibits coarse phase clustering, that is, the synchronization is not perfect. This is attributed to the fluctuations that tend to disperse the phases, and to the rapid decrease of the function (23).

The bottom panel of Fig.1 displays the structure of a phase cluster. The figure shows the histogram of the phases of the oscillators in a cluster obtained from a numerical simulation of Eq. (1). The red continuous curve is a numerical computation of the steady state of the nonlinear Fokker-Planck equation (15) . This steady state fits the simulation results very well even though the nonlinear Fokker-Planck equation relies on two approximations: the reduction of the amplitude and phase equation Eq. (1) to Eq. (8) for only the phases, and the mean field assumption that $N\to\infty$ although we have a finite number of oscillators. The dashed curve shows an analytic estimation of the steady state solution of the nonlinear Fokker-Planck equation, obtained from the normal form approach (see Appendix B):

Here the phase $\psi$ is a constant that defines the phase-cluster location in the unit circle and is fixed by the initial condition and the particular realization of the fluctuations. For finite $N$, the location of this phase in the unit circle is random. This random phase is not captured by the steady state solution of Eq. (15) because in the mean field approximation (14) the limit $N\rightarrow\infty$ is applied before the limit $t\rightarrow\infty$. In general these limits do not conmute, and finite size effects are lost in the mean field equation Eq. (15). In any case, sufficiently large $N$ allows for the formation of a phase-cluster which is associated with a non-vanishing value of the order parameter $R$. This fact is independent of the order of the limits.

Equation (26) can be written in the form

where we define the transition rates

Equation (30) is a master equation for a discrete state system, where the transitions between the $M$ states are governed by the rates (31). Only nearest-neighbor transitions ($j\rightarrow j\pm 1$) are allowed. We stress that Eq. (30) is nonlinear because the rates (31) depend on the $P_{j}$ via $\Omega_{j}$ [see Eqs. (27) and(28)].

We have arrived at the following point in this development: the stochastic dynamics of a single Kuramoto oscillator described by the continuum equation Eq. (1) may be converted to the problem of a continuous-time Markov chain of $M$ states, with periodic boundary conditions. These states model the internal structure of each stochastic oscillator. The interaction between the oscillators is contained in the dependence of the rates (31) on the global state of the system. An ensemble of $N$ of these $M$-state oscillators is characterized by the densities

where $N_{j}\left(t\right)$ is the number of oscillators in state $j\in\left\{0,...,M-1\right\}$ at time $t$. The rates (31) depend on the global state $\left\{n_{j}\left(t\right)\right\}_{j=0}^{M-1}$ via the $\Omega_{j}$, which depend on the finite size order parameter

In the thermodynamic limit,

and the order parameter (32) converges to the mean-field version (28). In this limit, the infinite size ensemble of $M$-state units is described by the nonlinear master equation (30), in the same spirit as the nonlinear Fokker-Planck equation (15) description of the continuum units. If $N$ is finite, we need to work with a set of coupled Langevin equations, as in Ref. Pinto ; Rosas . However, except for our numerical simulations, finite size effects are beyond the scope of this paper; here we focus on the mean-field theory. We note that first taking the limit $N\rightarrow\infty$, and then the limit $t\to\infty$ to study the steady states of (30) in general does not commute with taking these two limits in the opposite order Pinto .

The continuous-time Markov-chain modeling scheme rests on the assumption that the rates (31) are positive, which is not trivially satisfied. This requirement imposes constraints: $w_{j\rightarrow j\pm 1}$ is positive for any initial condition and realization of the fluctuations only if

where $\Omega_{max}$ is the maximum value of $\left|\Omega_{j}\right|$ for any $j$ in any realization of the evolution of the system. An evident bound is obtained from Eq. (27) by noting that the maximum value of the sine is unity,

Fixing the noise intensity and using the coupling strength $K$ as the control parameter, we impose the condition

Equations (35) and (34) ensure that condition (33) is satisfied. In the limit $\Delta\phi\rightarrow 0$ (continuous phase space), there is no upper limit on the coupling ($K_{max}\rightarrow\infty$), as expected from our earlier analysis.

The desynchronized state (29) becomes unstable when the coupling strength $K$ exceeds the critical value (see Appendix A)

This converges to the value (20) when $\Delta\phi\rightarrow 0$ (continuous state space limit). At the opposite extreme, when $\Delta\phi=\pi$ (which corresponds to only two states on the circle, ($M=2$), Eq. (36) gives $K_{c}=\infty$, so that no instability to a synchronized state is observed. To observe synchronization within the framework of the continuous-time Markov chain picture, we must impose the condition

which implies [see Eqs. (25) and (35)]

For the standard Kuramoto model, $f(r)=1$ and $F(r)=r$, that is, $f(0)/F_{max}=1$. Coarse-graining for this model requires that $M>\pi/\arctan(1/2)$, that is, $M\geq 7$. Discretization with smaller $M$ within the Kuramoto model can not be interpreted as a Markov chain (the resulting “transition rates” may not be positive).

Coarser graining can be obtained by moving away from the Kuramoto model and considering ratios $f(0)/F_{max}<1$. For the model function (23), maximizing the coupling $F(r)$ [cf. Eq. (10)] yields

This allows us to arbitrarily decrease the number of states in each unit. The limiting case $M=3$ can be reached by choosing $a<e/6\cong 0.453$.

Although we can arbitrarily manipulate the number of states by using the function $f(r)$, the limiting case $M=3$ exhibits anomalies in the bifurcation structure that deserve a separate analysis. In the next subsection we discuss the 3-state case in detail. For $M\geq 4$ the transition to synchronization is described by the normal form for the mean field order parameter (see Appendix A)

where

Clearly, as $\Delta\phi\to 0$ or equivalently $M\to\infty$, $\alpha_{M}\rightarrow\alpha$ and $\beta_{M}\rightarrow\beta$ as defined in Eq. (22). Moreover, the bifurcation picture is quite similar to that of the continuous phase oscillator. For decreasing $f$-functions the bifurcation is supercritical, while increasing $f$-functions may induce subcriticality. Figure 2 displays a comparison between the direct numerical simulation of Eq. (1) presented as a five-bar histogram, the steady state solution of the Fokker-Planck equation (15) obtained numerically (red line), and the direct numerical simulation of the continuous-time Markov chain model for $M=5$ (black squares). The 5-state periodic Markov chain gives a surprisingly good picture of the dynamics displayed by the model Eq. (1).

The case $M=3$ requires a separate treatment because its behavior is completely different from those of the discrete $M\geq 4$ oscillators. In the $M=3$ case the mean field order parameter near the critical point and at the lowest nonlinear order obeys the normal form (see Appendix A)

where $R^{*}$ is the complex conjugate of $R$, and

which corresponds to a transcritical bifurcation. Separating magnitude and phase, that is, $R=re^{i\psi}$, we obtain

The desynchronized state $r=0$ with an arbitrary phase $\psi$ is stable for $K<K_{c}$ and unstable for $K>K_{c}$. Also, here we have the unstable fixed points

These are hyperbolic points in a two-dimensional phase space that undergo a phase shift from one of the angles in $\psi_{-}$ to one in $\psi_{+}$ when they cross the critical point (at the critical point $K=K_{c}$, $r_{-}=r_{+}=0$ because $\alpha_{3}=0$). The unstable manifold of $\left\{r_{-},\psi_{-}\right\}$ is unstable to perturbations of the modulus $r$, but stable to phase perturbations. In contrast, $\left\{r_{+},\psi_{+}\right\}$ is stable to perturbations of the modulus $r$, but unstable to phase perturbations.

The nonlinear saturation of the instability occurs with the inclusion of higher nonlinear orders, which are not captured in the weakly nonlinear analysis used to deduce Eq. (40). We analyze it numerically and again use the model function (23) for $f(r)$. Fig. 3 displays the nullclines of the dynamical system

where the nonlinear functions $G_{0}$ and $G_{1}$ are obtained from the mean field master equation (30), with $M=3$ and the normalization condition $P_{2}=1-P_{0}-P_{1}$. The nullclines correspond to the curves $G_{0}\left(P_{0},P_{1}\right)=0$ and $G_{1}\left(P_{0},P_{1}\right)=0$ in the bidimensional phase space $\left(P_{0},P_{1}\right)$. Moreover, due to the fact that $\left(P_{0},P_{1}\right)$ are normalized probabilities, the physically accessible phase space is restricted to the right triangle defined by the region bounded by the lines $P_{0}>0$, $P_{1}>0$, and $P_{0}+P_{1}<1$. The fixed points correspond to the intersections of the two nullclines. Therefore, the bifurcation scenario is the following: For low enough coupling strength, the only stable equilibrium is the desynchronized state, $r=0$ (Fig. 3a). At some point, say $K=\widehat{K}_{c}<K_{c}$, six new fixed points appear simultaneously by saddle-node bifurcation, three of them stable and three of them unstable. In the region $\widehat{K}_{c}<K<K_{c}$, there are four stable fixed points (Fig. 3b). At the critical point $K=K_{c}$ (Fig.3c), the three unstable fixed points collide with the desynchronized state, destabilizing it. For $K>K_{c}$, there remain only three stable fixed points, related to one another by synchronization with a phase shift $\Delta\phi$ (Fig. 3d). The normal form (40) only describes the collision between the desynchronized state and the three hyperbolic points. A neater representation of the bifurcation scenario is shown in Fig. 4, where we have plotted the equilibrium order parameter modulus of $r$ as a function of the coupling strength. The nonvanishing-$r$ branches represent three steady states with different phases $\psi$. For the parameters in this section, $\widehat{K}_{c}=1.548\ldots$ and $K_{c}=1.654\ldots$.

Note that for $M=3$ (the three-state system) the phase invariance occurs in discrete steps $\Delta\phi$. This should be quite general and, from the structure of the expansion performed in Appendix A, one can see that,

where, in general, the coefficient $\Gamma$ should be computed for higher orders (perhaps thus revealing the existence of an extra time scale, slower than the one used for the expansion). In any case, it is intuitively rather clear that the distributions obtained from the discrete model will be invariant only to phase transitions occurring at discrete steps $\Delta\phi$.