Noise-induced stabilization of collective dynamics

TL;DR

This paper demonstrates that small additive noise can unexpectedly stabilize unstable collective oscillations in coupled oscillator systems, revealing noise-induced bifurcations and partial synchrony.

Contribution

It introduces a semi-analytical approach to explain how small noise stabilizes unstable regimes in globally coupled oscillators, supported by numerical and stability analyses.

Findings

Small white noise stabilizes unstable collective oscillations.

Two distinct noise-induced bifurcations lead to partial synchrony.

The macroscopic Fokker-Planck approach confirms the numerical results.

Abstract

We illustrate a counter-intuitive effect of an additive stochastic force, which acts independently on each element of an ensemble of globally coupled oscillators. We show numerically and semi-analytically that a very small white noise is able to stabilize an otherwise linearly unstable collective periodic regime. Microscopic simulations reveal two noise-induced bifurcations of different nature towards self-consistent partial synchrony. We develop a macroscopic treatment solving the corresponding nonlinear Fokker-Planck equation by means of a perturbative approach. The associated linear stability analysis confirms the results anticipated by the numerics. We also argue about the generality of the phenomenon.

Figures

Equations

$${\dot{\varphi }}_j=\frac{1}{N}\sum _{m=1}^{N}G({\varphi }_m-{\varphi }_j)+{\xi }_j(t),$$

$$G(\varphi )=\sin (\varphi +{\gamma }_1)+a\sin (2\varphi +{\gamma }_2).$$

$${\tilde{P}}_D(k)=\frac{1}{N}\sum _{j=1}^N{e}^{ik{\varphi }_j}.$$

$${\dot{\varphi }}_j=-\frac{i}{2}\left[{\tilde{P}}_D(1)e^{-i({\varphi }_j-{\gamma }_1)}+{\tilde{P}}_D(2)e^{-i(2{\varphi }_j-{\gamma }_2)}-\text{c.c.}\right].$$

$${\nu }_D=⟨{\dot{\varphi }}_j⟩,$$

$${\Omega }_D=⟨\dot{\Phi }⟩\text{where}\Phi =\arg \left[{\tilde{P}}_D(1)\right].$$

$$\sigma =\sqrt{⟨{\left(|{\tilde{P}}_D(1)|-⟨|{\tilde{P}}_D(1)|⟩\right)}^2⟩},$$

$$\frac{\partial P_D}{\partial t}$$

$$+D\frac{{\partial }^2{P}_D}{\partial {\theta }^2},$$

$${\tilde{P}}_D(k,t)={\int }_0^{2\pi }\mathrm{d}\psi P_D(\psi ,t)e^{ik\psi }.$$

$$\displaystyle B\bigr{[}\theta,\tilde{P}_{D}(1,t),\tilde{P}_{D}(2,t)\bigr{]}:=\int{\rm d}\psi G(\psi-\theta)P_{D}(\psi,t)$$

$$=-\frac{i}{2}\left({\tilde{P}}_D(1,t)e^{-i(\theta -{\gamma }_1)}+a{\tilde{P}}_D(2,t)e^{-i(2\theta -{\gamma }_2)}-\text{c.c.}\right).$$

$$-H_{D}=\left\{\Omega_{D}-B\bigr{[}\theta,\tilde{P}_{D}(1),\tilde{P}_{D}(2)\bigr{]}\right\}P_{D}(\theta)+D\frac{\partial P_{D}}{\partial\theta}(\theta)$$

$$\frac{\partial u}{\partial t}=-\frac{1}{2\pi }\frac{\partial }{\partial \theta }\int \mathrm{d}\psi G(\psi -\theta )u(\psi ,t)+D\frac{{\partial }^{2}u}{\partial {\theta }^2}$$

$$u(\theta ,t)=\frac{1}{2\pi }\sum _{k=-\infty }^{+\infty }\tilde{u}(k,t)e^{-i\theta k}$$

$$\tilde{u}(k,t)={\int }_0^{2\pi }u(\theta ,t)e^{ik\theta }\mathrm{d}\theta ,$$

$$\frac{\partial u}{\partial t}(\theta ,t)=\frac{1}{2\pi }\frac{\partial }{\partial t}\sum _{k=-\infty }^{+\infty }\tilde{u}(k,t)e^{-i\theta k}$$

$$=\frac{1}{4\pi }\left(\tilde{u}(1,t)e^{-i(\theta -{\gamma }_1)}+2a\tilde{u}(2,t)e^{-i(2\theta -{\gamma }_2)}+\text{c.c.}\right)$$

$$-\frac{D}{2\pi }\sum _{k=-\infty }^{+\infty }{k}^2\tilde{u}(k,t)e^{-i\theta k}.$$

$$\frac{\partial \tilde{u}(1,t)}{\partial t}$$

$$\frac{\partial \tilde{u}(2,t)}{\partial t}$$

$$\frac{\partial \tilde{u}(k,t)}{\partial t}$$

$$\mathrm{Re}({\lambda }_1)=\frac{1}{2}\cos {\gamma }_1-D<0.$$

$$\frac{\partial V(\theta)}{\partial\theta}=\Omega_{D}-B\bigr{[}\theta,\tilde{P}_{D}(1,t),\tilde{P}_{D}(2,t)\bigr{]}\;,$$

$$V(\theta )={\Omega }_D\theta -\frac{1}{2}{\tilde{P}}_D(1)e^{-i(\theta -{\gamma }_1)}-\frac{a}{4}{\tilde{P}}_D(2)e^{-i(2\theta -{\gamma }_2)}+\text{c.c.}$$

$$P_{0}(\theta)=\frac{-H_{0}}{\Omega_{0}-B\bigr{[}\theta,\tilde{P}_{0}(1),\tilde{P}_{0}(2)\bigr{]}}\;.$$

$$\eta=\frac{H_{0}}{P_{0}(\theta)}p(\theta)-P_{0}^{\prime}(\theta)-\Bigl{\{}\omega-B\bigr{[}\theta,\tilde{p}(1),\tilde{p}(2)\bigr{]}\Bigr{\}}P_{0}(\theta)$$

$$p(\theta)\!=\!\frac{P_{0}(\theta)}{H_{0}}\Biggl{(}\eta+P_{0}^{\prime}(\theta)+\Bigl{\{}\omega-B\bigr{[}\theta,\tilde{p}(1),\tilde{p}(2)\bigr{]}\Bigr{\}}P_{0}(\theta)\Biggr{)}\;.$$

$$\frac{\partial }{\partial t}u(\theta ,t)$$

$$\displaystyle-P_{D}(\theta)B\bigr{[}\theta,\tilde{u}(1,t),\tilde{u}(2,t)\bigr{]}\Biggr{)}+D\frac{\partial^{2}}{\partial\theta^{2}}u(\theta,t).$$

$$u(\theta ,t)=\frac{1}{2\pi }\sum _{k=-\infty }^{\infty }\tilde{u}(k,t)e^{-ik\theta }.$$

$$\dot{\mathbf{\tilde{\mathbf{u}}}}=\bigl{[}\mathcal{S}(\mathbf{\tilde{P}_{D}},\Omega)-D\mathcal{C}\bigr{]}\mathbf{\tilde{u}}$$

$${\mathcal{S}}_{mn}$$

$$+\frac{1}{2}{e}^{\text{i}{\gamma }_1}[{\tilde{P}}_D(1){\delta }_{m-1,n}+{\tilde{P}}_D(m-1){\delta }_{1,n}]$$

Full text

I Introduction

Typically, noise decreases the coherence of a dynamical system, by blurring, for instance, a perfect periodicity, or smoothing the fractal structure of low-dimensional chaos. Furthermore, noise can destabilize an attractor, when sufficiently large fluctuations allow overcoming an (effective) energy barrier. Such effects naturally occur whenever the unavoidable presence of stochastic forces is included into an otherwise deterministic evolution. Sometimes, however, noise may unexpectedly have the opposite effect of either increasing the overall coherence or stabilizing a given dynamical regime. Well known examples are the stabilization of the inverted pendulum Simons and Meerson (2009), stochastic resonance Benzi et al. (1981); McDonnell and Abbott (2009) and coherence resonance Lindner et al. (2004); Pikovsky and Kurths (1997).

In this paper we discuss another such instance, where a finite but small amount of white noise, acting independently on an ensemble of identical oscillators, stabilizes self-consistent partial synchrony (SCPS), an ubiquitous, collective regime Clusella et al. (2016) observed in ensembles of identical oscillators. This phenomenon differs from standard noise-induced bifurcations (see, e.g.  Aumaître et al. (2007); Lucke (1989)), where the Lyapunov exponent measuring the local stability of a given regime changes sign because of the fluctuations induced by a (multiplicative) stochastic force. Here, the noise acts on the microscopic level, while the regime we are interested in is a collective state. Macroscopically, SCPS corresponds to a rotation of the probability density of oscillator phases. Its dynamics, controlled by a nonlinear continuity (Liouville-type) equation, takes place within an infinite-dimensional functional space. Depending on the control parameters, one or more eigenvalues of the linearized equations (the so-called Floquet exponents) may have a positive real part, implying that SCPS is unstable and cannot be thereby maintained indefinitely.

On the macroscopic level, the effect of noise is described by a diffusion operator, which adds up to the continuity equation, transforming it into a nonlinear Fokker-Planck equation. The impact on the overall dynamics may be relatively trivial, as the regularization of the switching dynamics (see the next section for its definition), but also much less so, when it stabilizes SCPS as shown in this paper with the help of both numerical simulations and semi-analytical calculations.

We mostly focus on an ensemble of Kuramoto-Daido phase oscillators Daido (1993a); *Daido-93a; *Daido-96, whose evolution, in the absence of noise, is entirely controlled by the coupling function $G(\phi)$. The simplest such example is the Kuramoto-Sakaguchi model Sakaguchi and Kuramoto (1986) where $G(\phi)$ is a sinusoidal function. However, in the last years it has been understood that a purely harmonic coupling is rather special: a few macroscopic variables suffice to describe the collective dynamics Watanabe and Strogatz (1993); *Watanabe-Strogatz-94; Ott and Antonsen (2008). At the same time, it has emerged that the addition of a second harmonic suffices to enrich the resulting phenomenology. For instance, two- (and three-) cluster states Hansel et al. (1993); Kori and Kuramoto (2001) have been found in ensembles of identical oscillators, while a high degree of multistability has been observed in the presence of diversity Komarov and Pikovsky (2014). Furthermore, this setup is the minimal one where SCPS can spontaneously emerge Clusella et al. (2016).

The addition of noise to the bi-harmonic setup has been already studied in Ref. Vlasov et al. (2015) for parameter values where, however, no “shear” phenomena such as SCPS can arise. Here, we indeed explore a region where neither the fully synchronous, nor the splay state are stable in the deterministic limit. Numerical simulations of the microscopic equations performed for different noise levels show that a noise amplitude $1.2\cdot 10^{-4}$ (see the next section for a proper definition) can stabilize SCPS.

For smaller noise, a second regime is present: it is characterized by a pulsating bimodal distribution of phases, which can be traced back to the existence of two-cluster states. The transition between the two regimes is controlled by a bifurcation that can be either super or subcritical, this meaning that there exists a bistability region where, depending on the initial conditions, either of the two regimes can be attained.

The details of the bifurcation diagram have been unraveled with the help of a perturbative approach, the noise amplitude being the smallness parameter. As a result, we have been able to estimate the deviations induced by the noise on the actual shape of the probability density of phases, and then to solve the corresponding linearized equations, to determine the stability properties of SCPS. These detailed studies of the biharmonic model are accompanied by a similar (but purely numerical) investigation of an ensemble of weakly coupled Rayleigh oscillators, where an analogous stabilization of SCPS emerges, when a small noise is added to the microscopic evolution equations.

Altogether, the paper is organized as follows. In Sec. II we introduce the basic model, the biharmonic coupling function, and recall the phase diagram observed in the purely deterministic limit. In Sec. III we reconstruct the scenario in the presence of noise by performing direct simulations of the oscillators. There we give pictorial representations of the relevant regimes and confirm the typical signature of SCPS: a difference between the mean frequency of the distribution and that of the single oscillators. Sec. IV is devoted to a thorough perturbative analysis of the macroscopic equations. This includes a general remark on the fact that, in the presence of noise, the fully synchronous state becomes conceptually indistinguishable from SCPS. Finally, a brief analysis of Rayleigh oscillators is discussed in Sec. V together with general remarks about future perspectives.

II The model

In this section we introduce the main reference model: an ensemble of identical, globally-coupled Kuramoto-Daido oscillators,

[TABLE]

where $\xi_{j}(t)$ is a zero-average white noise such that $\langle\xi_{j}(t)\xi_{m}(t^{\prime})\rangle=2D\delta_{jm}\delta(t-t^{\prime})$. The coupling function $G(\phi)$ is assumed to have a biharmonic shape

[TABLE]

The dynamics of the phase oscillators can be characterized with the help of the order parameters

[TABLE]

The first two parameters suffice to express the evolution equation (1) in the more compact form

[TABLE]

A first important observable used throughout the paper to identify the different regimes is the microscopic frequency,

[TABLE]

where $\langle\cdot\rangle$ stands for temporal average, which is the same for all oscillators (no symmetry breaking). Furthermore it is convenient to introduce the mean-field frequency of the ensemble, defined as

[TABLE]

The analysis carried out in Ref. Clusella et al. (2016) has revealed two symmetric parameter regions where, in the absence of noise, neither the fully synchronous regime nor the splay state are stable. Therein, two dynamical regimes have been identified: self-consistent partial synchrony (SCPS) and a switching dynamics (SD). In the SCPS regime, the probability distribution of oscillator phases rotates with a constant velocity without changing shape. In the SD regime (first discussed in Ref. Hansel et al. (1993); Kori and Kuramoto (2001)), a two-cluster regime is attained, characterized by oscillations of the cluster widths and of their separation. This regime originates from a nontrivial form of instability. The behavior of infinitesimal perturbations is controlled by the inter-cluster exponent, which measures the response to perturbations of the mutual distance between the two clusters; and by two intra-cluster exponents, which quantify the growth rate of the two cluster-widths. In the biharmonic model the former exponent is negative, while the two latter ones have opposite sign. This means that while the width of one cluster increases, the other decreases exponentially. When the width of the wider cluster becomes of order 1, nonlinearities induce an exchange in the order of the two clusters, so that it starts decreasing. Therefore, the overall stability is controlled by the sum of the two intracluster exponents, which is negative in this biharmonic model. This means that both cluster-widths oscillate between a value of order 1 and a minimal value, which becomes progressively smaller: this is nothing but the convergence towards an heteroclinic cycle. When the minimal width becomes smaller than the computer accuracy, a spurious convergence to the cluster state occurs. A small amount of disorder among the oscillators (of order $10^{-12}$) eliminates this artificial effect, giving rise to periodic oscillations which, however, depend weakly on the amount of disorder.

For $a=0.2$ and $\gamma_{2}=\pi$, the parameter regions characterized by these nontrivial regimes are the intervals $1.159\simeq\operatorname{acos}(0.4)<\gamma_{1}<\pi/2$ and $3\pi/2<\gamma_{1}<2\pi-\operatorname{acos}(0.4)\simeq 5.124$. Here, we will focus on the first one since the second is analogous under the transformations $\phi\to-\phi$ and $\gamma_{1}\to 2\pi-\gamma_{1}$. In Fig. 1 we show the dependence of various observables on $\gamma_{1}$ in the purely deterministic case (see the red circles). In panel (a) we see that upon decreasing $\gamma_{1}$, SCPS emerges from the splay state at $\pi/2$ and becomes unstable at $\gamma_{1}\simeq 1.401$ through a subcritical Hopf bifurcation. For yet smaller $\gamma_{1}$-values, a switching dynamics is observed, until the perfectly synchronous state becomes stable below $\gamma_{1}=1.159$. In panel (b), we can appreciate the typical signature of SCPS: a difference between the microscopic $\nu_{D}$ and macroscopic $\Omega_{D}$ frequency. The SD is also characterized by a frequency difference, although it is so small it cannot be well appreciated in Fig. 1b. Finally, the behavior of the microscopic frequency is plotted in panel (c), where one can again recognize the loss of stability of SCPS upon decreasing $\gamma_{1}$.

III Microscopic approach

The very existence of SD reveals that the deterministic dynamics is extremely sensitive to the presence of disorder. It is therefore crucial to construct a phase diagram which includes the noise strength. The results of detailed simulations performed for different values of $\gamma_{1}$ and $D$ (and different system sizes) are reported in Fig. 2. Five regions can be recognized: on the left synchronous dynamics is observed, while the splay state is found on the right of the diagram; in between, region A corresponds to stable SCPS, while region B to stable SD; finally, the two latter regimes coexist within region C.

On a more quantitative level, we have monitored $\tilde{P}_{D}(1)$, $\Omega_{D}$, and $\nu_{D}$, for $D=1.2\cdot 10^{-4}$, upon varying $\gamma_{1}$. The results are reported in Fig. 1 (see the black crosses), where the noise is so small that whenever SCPS is deterministically stable, no appreciable changes are observed. An important difference with the deterministic case is that SCPS seems to extend down to the region where full synchrony is stable. Additionally, we see that the transition from SCPS to full synchrony, expected around $\gamma_{1}\approx 1.159$, is smoothed out. We anticipate that this is because it is no longer a true bifurcation. As shown in the next section, the onset of SCPS out of full synchrony can be seen as the tilting of a washboard potential beyond the point where a minimum is present. In the presence of noise, an otherwise $\delta$-distribution is broadened, making it possible to have phase jumps even in the synchronous regime. In other words, the synchronous state is not exactly synchronous and the microscopic and macroscopic frequencies differ from one another. Additionally, in Fig. 1(b) we see a curious finite-size effect induced by noise. In the asynchronous (splay-state) regime, $\tilde{P}_{D}(1)$ is not strictly zero for finite $N$. One can thereby determine its phase and compute the corresponding growth rate, which, in the presence of noise, coincides with the frequency $\Omega_{D}$ of the relaxation oscillations (no such oscillations are generated when $D=0$).

In order to illustrate the difference between SCPS and SD, in Fig. 3 we have plotted a snapshot of the probability density $P_{D}(\theta)$ for two different points in parameter space, both falling in the region where SCPS is unstable for $D=0$ (see the green circles in Fig. 2). Panel (a) corresponds to the point inside region A: here, the probability density shifts rigidly with a constant velocity, as expected for SCPS. At variance with the deterministic case, the microscopic quasiperiodicity is obviously lost, due to the presence of noise: it still holds true that the microscopic average frequency differs from the macroscopic one. Panel (b) corresponds to the green point inside region B: here, the distribution is bimodal and the two peaks breath – a reminiscence of the different stability of the two clusters.

A more accurate characterization of SCPS and SD is obtained by looking at the time evolution of the order parameter $|\tilde{P}_{D}(1)|$. In Fig. 4 we present the trace for three points, all corresponding to the same $\gamma_{1}$-value and different noise amplitudes (see the purple squares in Fig. 2). The black dotted curve corresponds to zero noise (a very small quenched randomness has been added to avoid the spurious collapse onto a two-cluster state). Strong, periodic fluctuations are observed, associated to the alternation of contraction and expansion of the cluster widths. Upon increasing the noise amplitude, periodic oscillations are still observed, the amplitude of which decreases while their period shrinks (see the red dashed curve in Fig. 4). A logarithmic reduction of the period was already noticed in Hansel et al. (1993): it is due to the fact that the noise prevents a cluster to become too thin. The additional fluctuations are finite-size effects which decrease upon increasing the system size.

For a still small but larger noise ($D=1.2\cdot 10^{-4}$), the oscillations practically disappear: the fluctuations exhibited by the blue continuous curve are just manifestations of finite-size effects which decrease upon increasing the number of oscillators. In fact, for this noise amplitude any reminiscence of SD is lost, as confirmed by the series of crosses plotted in Fig. 1.

Altogether, the direct simulations suggest that the transition from SCPS to SD corresponds to a Hopf bifurcation, beyond which a constant order parameter starts exhibiting periodic oscillations, which are reminiscent of the presence of the unstable two-cluster state. One can interpret SD as a sort of more structured SCPS, since also in this case there is a difference between the microscopic and macroscopic frequency. In order to validate this interpretation, we have determined the time averaged Kuramoto order parameter $\langle|\tilde{P}_{D}(1)|\rangle$ and the mean-field frequency $\Omega_{D}$ for $\gamma_{1}=1.34$ and an adiabatic increase of $D$ from region B to A (along the left dashed line in Fig. 2). In the top panel of figure 5 we see that $\langle|\tilde{P}_{D}(1)|\rangle$ progressively decreases upon increasing the noise: this is because noise tends to glue together the two clusters. A signature of a true transition can be seen in panel (b) where the (temporal) standard deviation of $|\tilde{P}_{D}(1)|$,

[TABLE]

is plotted for different numbers of oscillators. Upon increasing $N$, $\sigma$ clearly approaches zero above a critical noise strength. Finally, in panel (c) we see that $\Omega_{D}-\nu_{D}$ is different from zero both above and below the bifurcation, confirming that the qualitative difference between SD and SCPS is just the periodic modulation of the Kuramoto order parameter.

The transition scenario from SD to SCPS remains essentially unchanged so long as $\gamma_{1}\lesssim 1.4$. For larger $\gamma_{1}$-values, in region C, both SCPS and SD are stable for arbitrarily small noise so that one expects the Hopf bifurcation separating B from C to be subcritical. This scenario is confirmed by the discontinuous jumps observed in Fig. 6, where the average order parameter, the standard deviation $\sigma$, and the frequency difference are plotted while adiabatically increasing $D$ starting from the SD regime (along the right vertical dashed line in Fig. 2). In practice SD suddenly disappears through a saddle-node bifurcation, where it collides with a similar unstable regime.

Upon increasing $D$, these two bifurcation lines bounding region C merge together in a tricritical point, a Bautin bifurcation. The identification of this point requires some extra effort since, close to it, the amplitude of the discontinuity progressively vanishes. By comparing simulations performed by adiabatically increasing and decreasing $\gamma_{1}$, we estimate the tricritical point to be located at $\gamma_{1}\simeq 1.3885$, $D\simeq 7\times 10^{-5}$.

IV Macroscopic description

In the previous section we have seen that a small noise stabilizes SCPS with the exception of the tiny parameter region B, where periodic oscillations are observed. Here, we approach the problem from the macroscopic point of view, extending the method introduced in Ref. Clusella et al. (2016) to account for the presence of microscopic noise.

In the thermodynamic limit, $N\to\infty$, the evolution of the probability density $P_{D}(\theta,t)$ of oscillators with phase $\theta$ at time $t$ is controlled by the nonlinear Fokker-Planck equation

[TABLE]

where $D$ is the diffusion term. This equation refers to a rotating frame $\theta=\phi-\Omega_{D}t$. For a properly selected $\Omega_{D}$-value this equation admits a stationary solution $P_{D}(\phi)$, which corresponds to the SCPS regime. It is convenient to introduce the Fourier representation $P_{D}(\phi)=\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}\tilde{P}_{D}(k,t)\operatorname{\rm e}^{-ik\phi}$ where

[TABLE]

In fact, by using this notation, the integral in Eq. (IV) can be simplified, making use of Eq. (2),

[TABLE]

Notice that the $\tilde{P}_{D}(k,t)$ coefficients coincide with the order parameters introduced in Eq. (3).

A stationary solution of Eq. (IV) can be obtained by setting the time derivative of $P_{D}$ equal to zero and thereby integrating the r.h.s. to obtain

[TABLE]

where $H_{D}$ is the probability flux. The mean frequency $\nu_{D}$ of the single oscillators can be expressed in terms of the flux as $\nu_{D}=\Omega_{D}+2\pi H_{D}$.

Two simple solutions are characterized by a zero flux $H_{D}=0$: the splay state and the full synchrony. The solution corresponding with $H_{D}\neq 0$ corresponds to SCPS. We now discuss in detail these three cases.

IV.1 Stability of the splay state

The splay state is characterized by $P_{D}(\theta,t)=1/(2\pi)$ and $\Omega_{D}=0$. Let us consider an infinitesimal perturbation $u(\theta,t)$ of $P_{D}$. The corresponding linearized equation is

[TABLE]

where we have used that B\bigr{[}\theta,\tilde{P}(1),\tilde{P}(2)\bigr{]}=0. Using the Fourier expansion

[TABLE]

where

[TABLE]

and solving the integral term as in Eq. (5), the time evolution of $u$ reads

[TABLE]

The evolution equations of the Fourier modes $\tilde{u}(k,t)$ are then given by

[TABLE]

complemented by the complex conjugate equations for the negative $k$-modes. Manifestly, the equations are diagonal. Since the eigenvalues corresponding to the eigenfunctions with $|k|>2$ are negative real numbers, the stability is determined only by the eigenvalues corresponding to the two first Fourier modes, $\lambda_{1}=\frac{1}{2}\operatorname{\rm e}^{i\gamma_{1}}-D$ and $\lambda_{2}=a\operatorname{\rm e}^{i\gamma_{2}}-4D$. In the particular case $\gamma_{2}=\pi$, the second eigenvalue has $\operatorname{Re}(\lambda_{2})<0$ so that the stability of the splay state is controlled only by the first mode,

[TABLE]

The critical curve where $\operatorname{Re}(\lambda_{1})=0$ corresponds to the right almost vertical line reported in Fig. 2.

IV.2 The synchronous state and self-consistent partial synchrony

In the deterministic case, the fully synchronous state is characterized by a $\delta$-like distribution and there is a well defined stability boundary for this solution. In order to understand what happens once noise is added, it is convenient to look at Eq. (IV) as if the velocity field were given a priori. It corresponds to a standard Fokker-Planck equation in a washboard potential defined as

[TABLE]

so that

[TABLE]

apart from an arbitrary additional constant. Because of the tilting, as soon as $D>0$, any stationary solution is characterized by a non-zero flux $H_{D}$, even if the potential has well defined minima. In other words, the fully synchronous state becomes formally equivalent to SCPS and no transition can be any longer found. Nevertheless, one can still identify a sort of critical line separating the regime where the potential has local minima and the flux is thereby driven by the noise, from the one where no minima exist and the flux is the result of a deterministic current. In figure 7 we show how the potential $V(\theta)$ changes shape upon varying $\gamma_{1}$ for a fixed noise strength (the unknown parameters $\Omega_{D}$, $\tilde{P}_{D}(1,t)$ and $\tilde{P}_{D}(2,t)$ have been determined with the help of direct simulations). Using this technique we have reconstructed the “transition” line reported in Fig. 2 (see the left quasi-vertical line).

Having understood that, once noise is added, SCPS and full synchrony are one and the same regime, now we focus on the procedure to determine the shape of the stationary distribution. For $D=0$, i.e. in the deterministic case discussed in Ref. Clusella et al. (2016), the probability density can be determined from Eq. (6) without the need of performing any integration,

[TABLE]

This solution is valid if the denominator has no zeros, i.e. if there are no minima in the corresponding potential. The above expression depends on two complex variables $\tilde{P}_{0}(1)$, $\tilde{P}_{0}(2)$ and two scalars $\Omega_{0}$, $H_{0}$. Since we are free to choose the phase of the distribution $P_{0}$, we can assume that $\tilde{P}_{0}(1)$ is real. Moreover, $H_{0}$ can be obtained by imposing the normalization condition. Therefore, the determination of $P_{0}$ requires finding a fixed point in a four-dimensional space: this problem was tackled and solved numerically in Ref. Clusella et al. (2016).

In the general case $D>0$, the stationary solution $P_{D}$ must be obtained by integrating the ODE (6). One could obtain an explicit expression for $P_{D}$ by introducing a Jacobi-Anger expansion. Such expression would involve series with terms depending on Bessel functions, making both the analytic and numerical treatment highly complex and, thus, inappropriate to effectively study $P_{D}$ and its stability properties.

It is more convenient to develop a perturbative formalism to investigate the problem in a semi-analytic way in the limit of small noise, i.e. $D\ll 1$. At first order, the probability density can be written as $P_{D}(\phi)=P_{0}(\phi)+Dp(\phi)$. Moreover, since we expect variations of the macroscopic frequency as well as of the flux, we assume $\Omega_{D}=\Omega_{0}+D\omega$ and $H_{D}=H_{0}+D\eta$. Upon replacing these assumptions into Eq. (6) and retaining first order corrections in $D$, it is found that

[TABLE]

where $\tilde{p}(k)$ stands for the $k$-th Fourier mode of $p$, while the prime denotes a derivative with respect to $\theta$. This equation can be easily solved for $p(\theta)$, obtaining

[TABLE]

In order to complete the identification of the solution, it is necessary to determine $\tilde{p}(1),\tilde{p}(2)$ and the two scalars $\omega$, $\eta$. Since $\tilde{P}(1)$ in equation (9) is real, so has to be $\tilde{p}(1)$, while $\eta$ can be obtained by imposing the “normalization” condition $\int p(\psi){\rm d}\psi=0$. Therefore we are facing a problem of the same complexity as in Eq. (9); it can be solved using similar procedures. An example of $P_{0}$ and of the corresponding variation $p$ obtained by solving Eqs. (9,10), is shown in Fig. 8. From the shape of $p$, we see that the effect of noise is to deplete the left shoulder of the density and to raise a bit its tails.

IV.3 Stability analysis of SCPS

An accurate estimate of $P_{D}$ is a necessary requisite for a reliable stability analysis. In fact, only the stationary solution is marginally stable against a rigid translation.

Upon linearizing Eq. (IV) around $P_{D}(\theta)$, we find that an infinitesimal perturbation $u(\theta,t)$ satisfies the equation

[TABLE]

The spectrum of this linear operator is purely point-like and can be determined by approximating the infinite-dimensional operator with finite matrices of increasing size. The most effective method consists in expanding $u(\theta,t)$ into Fourier modes

[TABLE]

By making use of Eq. (5), we can rewrite the evolution equation as

[TABLE]

where $\mathbf{\tilde{u}}=\{\tilde{u}(k,t)\}_{k=-\infty}^{\infty}$, $\mathbf{\tilde{P}_{D}}=\{\tilde{P}_{D}(k)\}_{k=-\infty}^{\infty}$, while $\mathcal{S}$ and $\mathcal{C}$ are two matrices defined as follows

[TABLE]

[TABLE]

By recalling that $\tilde{P}_{D}(m)=\tilde{P}_{0}(m)+D\tilde{p}(m)$, we can insert this perturbative expression into Eq. (11), obtaining

[TABLE]

The stability of SCPS is determined by the eigenvalues $\{\Lambda^{(k)}_{D}\}_{k}$ of the linear operator defined by this equation. Observe that $\dot{\tilde{u}}(0,t)=0$, while $\tilde{u}(0,t)$ is decoupled from all other Fourier modes. This is a consequence of norm conservation, which implies that one zero eigenvalue is always present. Another obvious property of the linear operator is that exchanging positive with negative components is equivalent to taking the complex conjugate: this is due to the probability density being a real quantity.

In the limit $D\to 0$ the evolution operator reduces to $\mathcal{S}(\mathbf{\tilde{P}_{0}},\Omega_{0})$. Its eigenvalues $\Lambda^{(k)}_{0}$ are responsible for the instability of SCPS in the zero-noise limit, reported in the previous section.

The perturbative corrections to the evolution operator defined in Eq. (11) are composed of two contributions: the first one is due to the change of shape of the probability density induced by noise; the second one is the direct consequence of the diffusion term, which strongly damps short-wavelength Fourier modes (see the negative diagonal terms of $\mathcal{C}$ proportional to $m^{2}$). The eigenvalues $\{\Lambda^{(k)}_{D}\}_{k}$ can be determined perturbatively by following standard procedures. Upon expanding the eigenvalues up to first order in $D$, $\Lambda^{(k)}_{D}=\Lambda^{(k)}_{0}+D\lambda^{(k)}$, one can indeed write Wilkinson (1965)

[TABLE]

where $\mathbf{W}_{0}$ and $\mathbf{V}_{0}$ are the left-, resp. right-hand eigenvectors associated to $\Lambda^{(k)}_{0}$.

In principle, the matrix indices range from $-\infty$ to $+\infty$. One can however obtain sufficiently accurate estimates of the largest eigenvalue by restricting the analysis to a finite range $[-L,L]$, provided that $L$ is large enough. We have verified that $L=100$ suffices to reconstruct the relevant part of the stability spectrum. In practice, one first needs to solve numerically the Eqs. (9,10) to determine the stationary solution of the Fokker-Planck equation (up to first order in $D$). Then, left and right eigenvectors of the unperturbed evolution operator are obtained, so that we are finally able to determine the corrections to the eigenvalues.

The results are illustrated in Fig. 9, where the eigenvalues are plotted for $D=0$ (red dots) and $D=1.7\cdot 10^{-5}$ (the arrows denote the shift of each eigenvalue). The spectrum is composed of pairs of complex conjugate eigenvalues; two of them are characterized by a strictly zero imaginary part: the first one corresponds to the zero exponent, which follows from the invariance of the solution under a rigid phase shift and is present in the noisy regime as well; the second one is the most negative eigenvalue (see the leftmost red dot) - very weakly affected by noise.

In the deterministic case the imaginary parts of the eigenvalues are almost equispaced. They can be used to parametrize the eigenvectors, as $\mathrm{Im}(\Lambda_{D})$ is proportional to the wavenumber $k$ (see Ref. Clusella et al. (2016)). Moreover, for $D=0$ the real parts decrease exponentially for increasing $k$ Clusella et al. (2016). This is a crucial difference with mono-harmonic systems such as the Kuramoto-Sakaguchi model, where the conservation laws enforced by the Watanabe-Strogatz theorem Watanabe and Strogatz (1993); *Watanabe-Strogatz-94 imply the existence of infinitely many strictly-imaginary eigenvalues.

The effect of noise is to basically shift all eigenvalues towards more negative values: in fact, all of the arrows reported in Fig. 9 are practically horizontal and point to the left. The stabilizing effect depends strongly on $k$ (because of the diffusive term in Eq. IV); it is approximately proportional to $k^{2}$. In particular, it is very weak for the two pairs of unstable eigenvalues (see also the inset). Nevertheless, even for a maximally unstable $\gamma_{1}$ (the value chosen in Fig. 9), a tiny amount of noise is sufficient to stabilize one of the two pairs of unstable directions (a noise amplitude about 7 times larger would fully stabilize SCPS).

Additional information can be obtained by looking at the eigenvectors. Unsurprisingly they are localized in the region where the probability density is concentrated. Moreover, the higher the imaginary part of an eigenvalue, the larger the number of oscillations of the corresponding eigenvector: this is a manifestation of the above mentioned (approximate) relationship beteen imaginary parts and wavenumbers. The eigenvectors corresponding to the unstable directions, $\mathbf{V}_{0}^{(1)}$ and $\mathbf{V}_{0}^{(2)}$ are plotted in Fig. 10 for $D=0$, separating the real from the imaginary component (see the solid lines). There we also see that $\mathbf{V}_{0}^{(1)}$ is also reminiscent of the first and second derivatives of $P_{0}$ respectively (compare red continuous and black dashed lines in panels (a) and (b) in Fig. 10). The analogy extends to $\mathbf{V}_{0}^{(2)}$, if compared with the third and fourth derivatives of $P_{0}$, although it is much more qualitative (see panels (c) and (d) in Fig. 10).

Upon tracking the eigenvalues with positive real part for increasing $D$ we can identify where the stabilization of SCPS takes place. The transition points are plotted as blue triangles in figure 2 and show a good agreement with the direct numerical simulations.

V Discussion and conclusions

In the previous section we have performed a detailed stability analysis of SCPS in the presence of noise. This collective regime corresponds to a stationary state in a suitably moving frame, where its evolution is described by a nonlinear Fokker-Planck equation. The nonlinearities play a double role: on the one hand they contribute to a self-consistent determination of the underlying potential; on the other hand they contribute to the stability of the state itself. In fact, the dynamics is not purely drift-diffusion driven; as shown in Clusella et al. (2016) there may be unstable directions. In the previous section we have, however, seen that even a very small noise is sufficient to stabilize the collective dynamics.

It is therefore natural to ask whether this scenario is peculiar of the biharmonic setup. We check this point by studying another model, an ensemble of mean-field coupled Rayleigh oscillators, where SCPS has been observed and found to lose stability in a purely deterministic setup Clusella et al. (2016). We show that a small amount of noise is again able to stabilize SCPS. The Rayleigh oscillator model reads as

[TABLE]

where $X=N^{-1}\sum_{m}x_{m}$ and $Y=N^{-1}\sum_{m}\dot{x}_{m}$ are the mean field contributions to the coupling, and $\varepsilon=0.05$ is the coupling strength. We assume again white noise $\eta_{j}(t)$ with $\langle\eta_{j}(t)\rangle=0$ and $\langle\eta_{j}(t)\eta_{m}(t^{\prime})\rangle=2D\delta_{jm}\delta(t-t^{\prime})$. The parameter $\zeta$ determines the stability of the limit-cycle. In this work we discuss the case $\zeta=5$ for which there is a strong attraction. Therefore $\gamma$ is the main control parameter that is going to be tuned. An appropriate order parameter is

[TABLE]

where

[TABLE]

is the root mean square of the time evolution. Therefore $\rho=1$ when there is fully synchrony and $\rho=0$ when the oscillators are distributed uniformly. In the deterministic case, in the range $\gamma\in[-0.7,0.2]$ a wide number of dynamical regimes are observed (see red circles in figure 11). In particular, SCPS is observed for $-0.18\lesssim\gamma\lesssim 0.05$. Above $\gamma\simeq 0.05$ the system converges to a homogeneous nine-cluster state. Similarly to the biharmonic model, SCPS loses stability towards a two-cluster state for $\gamma\simeq-0.18$. This cluster state finally converges to full synchrony at $\gamma\simeq-0.57$. On the one hand, a small noise ($D=1\cdot 10^{-6}$) does not substantially affect the regions where full synchrony and the nine-cluster states are stable. On the other hand, it is once again able to stabilize SCPS in the entire interval up to full synchrony (see the black crosses in Fig. 11).

The study of two models of phase oscillators has shown that a small amount of microscopic noise stabilizes self-consistent partial synchrony. It is natural to ask whether this effect extends to other types of collective dynamics. In globally coupled identical maps, collective chaos can be observed Shibata et al. (1999). In such a setup, it was found that an additive noise of the type considered in this paper can reduce the dimensionality of the collective dynamics Shibata et al. (1999); Takeuchi and Chaté (2013). Considering that a chaotic evolution can be seen as a sort of wandering process among different unstable periodic orbits, it is tempting to interpret this reduction of dimensionality as a progressive stabilization of the dynamics along various directions. It will be instructive to further investigate this interpretation.

Acknowledgement

This work has been financially supported by the EU project COSMOS (642563).