Analytical approach to synchronous states of globally coupled noisy rotators

We study populations of globally coupled noisy rotators (oscillators with inertia) allowing a nonequilibrium transition from a desynchronized state to a synchronous one (with the nonvanishing order parameter). The newly developed analytical approaches resulted in solutions describing the synchronous state with constant order parameter for weakly inertial rotators, including the case of zero inertia, when the model is reduced to the Kuramoto model of coupled noise oscillators. These approaches provide also analytical criteria distinguishing supercritical and subcritical transitions to the desynchronized state and indicate the universality of such transitions in rotator ensembles. All the obtained analytical results are confirmed by the numerical ones, both by direct simulations of the large ensembles and by solution of the associated Fokker–Planck equation. We also propose generalizations of the developed approaches for setups where different rotators parameters (natural frequencies, masses, noise intensities, strengths and phase shifts in coupling) are dispersed.


Introduction
The Kuramoto model of globally coupled phase oscillators [1] is a paradigmatic model to study synchronization phenomena [2][3][4][5]. In the thermodynamic limit, stationary (in a proper rotating reference frame) synchronized states can be found analytically for an arbitrary distribution of natural frequencies [6]. The idea is that for any subgroup of oscillators having a certain natural frequency, a stationary distribution under a constant mean field and the corresponding subgroup order parameter can be found analytically. Then the solution of the general self-consistent problem can be expressed via integrals of these subgroup order parameters. Other analytical approaches for the Kuramoto model may be not restricted to stationary states, but usually have other restrictions since they apply only for identical oscillators like the Watanabe-Strogatz theory [7], or for a population with a Lorentzian distribution of natural frequencies like the Ott-Antonsen ansatz [8].
An important generalization of the Kuramoto model considers an ensemble of globally coupled rotators (that is, oscillators with inertia). It is in particular relevant for modeling power grid networks [9,10]. Synchronization features of globally coupled rotators, both deterministic and noisy, have been widely studied [11-14, 16, 15]; in particular, an approach similar to the analytical description [6] was developed [3,17,18]. However, for noisy coupled rotators, so far the stationary distributions were found only numerically. Similarly, there are no analytical formulae expressing the order parameter as a function of other parameters for the Kuramoto model with noise.
In this paper, we present a fully analytical approach to the problem of globally coupled noisy rotators for small intertia (small 'masses'). We analyze a stationary solution of the Fokker-Planck (Kramers) equation describing identical noisy rotators driven by the mean field, and obtain a closed expression for the order parameter for this subgroup. Then, for an arbitrary distribution of natural frequencies, a stationary solution for the global order parameter is expressed in an analytical parametric form. Based on this analysis, we derive the asymptotic form with respect to the small order parameter and use it to classify the transition to synchrony as a supercritical or a subcritical one.
The paper is organized as follows. We introduce the model of coupled noisy rotators in section 2, where we also formulate stationary equations for the distribution density. A solution of these equation in the limit of small inertia terms is presented in section 3 (another derivation of this solution is given in appendix). Important particular cases of noise-free rotators and of noisy coupled oscillators (i.e. the case of vanishing inertia term) are discussed in section 4. There, we also give general expressions valid for small order parameters, i.e. close to the transition point. Solutions for several popular distributions of the natural frequencies (e.g. Gaussian, Lorentzian, etc) with numerical examples are presented in section 5. In section 6 we show how the approach can be used for generic ensembles where not only natural frequencies, but other relevant parameters of rotators (masses, noise strengths, etc) are distributed. We conclude with section 7.

Noisy coupled rotators
In this paper we consider an ensemble of N globally coupled rotators characterized by their angles j n and velocities j n  (n=1,2, K, N). The rotators are coupled via the complex mean field The unit of time is chosen so that the coefficient at the friction term ( j n  ) is one and all the parameters are normalized by the friction coefficient. Parameter μ describes the mass of rotators (more precisely, μ should be called the moment of inertia of a rotator); below, we focus on the overdamped (small-masses) case of small μ. Parameter ε is the coupling strength. Parameters ω n describe torques acting on rotators; we assume them to be distributed with a density g(ω). Note, that our approach can be straightforwardly generalized to the case where other parameters are distributed as discussed in section 6. Because uncoupled rotators have mean angular velocities ω n , we speak of 'natural frequencies' instead of 'torques'. We do so also to make transparent a relation to the particular case μ=0, where system(2) is nothing else but the standard Kuramoto model for globally coupled phase oscillators with a distribution g(ω) of natural frequencies ω. The rotators are acted upon by the independent white Gaussian noise forcings sz t n ( ) with amplitude σ (which is assumed to be the same for all rotators), zero means z á ñ = t 0 n ( ) , and auto-correlations nn denotes the Kronecker delta, and d t ( ) is the Dirac δ-function). Whereas equations (2) are used for numerical simulations below, the analytical approach is developed in the thermodynamics limit  ¥ N . In the mean-field coupling models, synchronization is one of the fundamental effects which is relevant to many systems. A transition to synchrony can be fully characterized in terms of the Kuramoto order parameter (1). Physically, the amplitude r of this complex parameter describes the synchrony level of the elements belonging to the population considered. Nonvanishing r indicates the collective synchronization. So, it is important to develop an universal analytical approach allowing one to calculate r value and predict the global evolution of an ensemble consisting of macroscopically large number of interacting units.
Note that model (2) is not the most general one-it might include a phase shift in coupling (which corresponds to the Kuramoto-Sakaguchi overdamped system). We postpone the discussion of this case to section 6. Furthermore, for simplicity of presentation we consider in sections 2-5 only symmetric unimodal frequency distributions; the discussion of asymmetric distributions is also postponed to section 6.
We now consider the limit of infinitely large number of elements, i.e.  ¥ N . Furthermore, we look for stationary synchronous states, i.e. those with constant modulus of the order parameter and with a uniformly rotating angle: . Such states are possible only in the thermodynamics limit where the finitesize fluctuations vanish. It is convenient to introduce a new angle variable related to the angle of the mean field q j y q j = -= = -W u , Because all the deterministic forces in this reference frame are constant, this density evolves toward the stationary, time-independent, one. As the only parameters governing this distribution are detuning n w = -W and forcing term A=ε r, we seek this stationary distribution as a function depending on these parameters explicitly, q n P u A , , 0 ( | ). Then, the system takes the form The system of (4) and (5) is in fact a self-consistent system for determining the unknown order parameter r (frequency Ω is obtained from the condition that r is real). Below, we restrict ourselves to symmetric distributions only, w n w n ). Then, one may chose Ω=ω 0 , which corresponds to the symmetric stationary density q n q n = --- . In this case, the integral in (5) is real, which proofs that ω 0 is a correct value of the frequency of the mean field Ω. Note that there could also be other appropriate values of Ω besides ω 0 . These different values correspond to different synchronous branches. Now, (5) provides a parametric solution of the self-consistent problem since both r and ε are represented as functions of a free parameter A, The only remaining act here is finding solution P 0 of(4). In [17,19], a method of matrix continuous fractions was employed to solve(4) numerically, another numerical approach is described in [18]. In the next section 3, we present an analytical solution for the overdamped (small-mass) case.

Stationary distribution of the phases in the limit of small masses
Here, we present one of the analytical methods for obtaining the stationary density P 0 from(4) based on series representation. The second method is given in appendix.

Matrix representation of the Fokker-Planck equation
According to [12,17], we represent a stationary solution q n P u , 0 ( | ) of(4) as a double series in the parabolic cylinder functions F u p ( ) in u and the Fourier modes q e q i in θ, Here functions F u p ( ) are defined as Thus, the main challenge is to find the quantity n a A , 0,1 ( )via solving (7); this coefficient is nothing else but the order parameter for the group of oscillators having natural frequency w n + ; 0 therefore, we call it subgroup order parameter. It is then used in (8) to find the total order parameter r via averaging over the subgroup frequency ν. Below, we also use the symmetry property following from (7), For the fixed parameters A and ν, the set of equations for a pq can be formally solved by virtue of the matrix continuous fractions method [17,19]. According to this method, the expression for the subgroup order parameter n a A , 0,1 ( )is as follows: )is the element of the infinite matrix S, which is given by the recurrence formula where I is the identical matrix, and the infinite diagonal matrix D and the infinite tri-diagonal matrix D are defined as The main steps of the numerical procedure based on relations (9)-(11) employed for finding the order parameter with desired accuracy are described in detail in [17]. Although this approach for the computation of the steady-state value of r is efficient and has many potential applications, there are several limitations in its use. In particular, this method has high time cost for the case of slowly descending distributions w g ( ), e.g. for the Lorentz distribution. So, it is important to develop an analytical approximation for the calculation of the value r in a nonequilibrium stationary state.

Small mass approximation
Expression (10) allows for a perturbative expansion in the small parameter μ. Below we omit o(μ) assuming mass to be small enough. One can consider the resulting equalities approximately valid if for sufficiently small μ. In the first order in μ, we obtain In order to evaluate the inverse matrix -D 1 , we denote the principal minors as and introduce a cutoff (cyclic) harmonic number  +¥ d for Fourier modes. Then, the elements of the inverse matrix -D 1 take the form According to the Laplace theorem, for all jk (for convenience, we take = ), the following representation holds (14), we obtain for  j k that The principal minor M jk (13) is found from a recurrent equation with fixed j, i.e. with fixed bottom right corner (17)

Recurrent equation
Here I z and K z denote the principal branches of the modified Bessel functions of the first and second kind, respectively, of the order z. Using properties of the modified Bessel functions, we evaluate the limit Employing both (18) and (19), we obtain Similarly, we find a representation of the principal minor M kj , using its expansion at fixed k, i.e. at fixed upper left corner Using expressions (21) and (19), we find  (11)-(13), (16), (20), and (22) (23) is based on the method of moments for the density and on elimination of the velocity; it is presented in appendix.

Limiting cases
Above, we have derived a general expression for the subgroup order parameter (23). Here, we discuss its form in several important particular cases.

Noise-free case
For purely deterministic rotators, we have to set σ → 0. To find the noise-free limit of general expression (23) So, it is neccesary to calculate limits at s  0 of the two fractions of the modified Bessel functions in (24). The limits can be evaluated by virtue of the known uniform asymptotic expansions of the modified Bessel functions (see formula (9.7.7) in [20]) or using recurrence relations for these special functions (see formula (10.29.1) in [21]). In particular, it follows that

Massless rotators: Kuramoto model
In the massless case (μ=0), the problem (2) is in fact the Kuramoto model with noise. The analytical expression of the local order parameter is exact in the thermodynamics limit (we can now omit the first index) as follows (this formula has been independently derived in [22]): A more elaborated application of the theory to the Kuramoto model with noise will be presented elsewhere.

Synchronization transition
Here, we employ the general parametric representation of the order parameter as a function of the coupling constant, with n a A , 0,1 ( )given by (23) to characterize the synchronization transition, i.e. to characterize states with the order parameter close to zero. All formulas below are valid in the first order in small mass μ like (23) and are therefore approximate, but for the sake of simplicity we write them as exact relations below.
At the transition point the amplitude A vanishes. Thus, to unfold the nature of the transition, we expand r(A) in the Taylor series for small values of A. This expansion contains only odd powers of A since (23) is an odd To obtain the three initial coefficients in expansion(30), we rewrite (24) In the massless case μ=0, the expressions for C 0 , C 1 have been obtained in [23].
The nontrivial branch of solutions e r ( ) starts at e = C 1 .  This expression coinsides, in the first order in μ, with the result for deterministic rotators obtained in [16].
It is instructive to compare this critical value with the expression for the stability loss of the asynchronous state r=0 (equation (24) in [12]). In our notations, this general formula for the imaginary part of the eigenvalue x, valid also for large masses μ, reads Here the synchronization transition is discontinuous (a first-order transition). We see, that the type of the transition is determined by the sign of C 1 . Because this coefficient depends on the mass μ, there is a critical value of the mass at which the type of the transition changes, This formula is obtained from (31) under condition C 1 =0. The expression is valid only if the value of m* is small enough. Concluding this section we would like to mention, that the full bifurcation analysis of the synchronization transition should also include stability analysis of all branches. This at the moment is missing, as we have analytic expressions from the stability of the trivial asynchronous state only. Therefore we use above the term 'transitions' understood at the physical level of rigor, and avoided term 'bifurcations' which would suggest mathematically rogortous statements.

Examples
Here we present explicit calculations of the synchronization transition for several commonly explored distributions of the frequencies. Each of these distributions is characterized by a width γ, and in all cases the result depends on dimensionless parameters x g s = 2 , ε/γ and μγ.

Gaussian distribution
The Gaussian distribution of the rotators' frequencies is where erfc is the complementary error function. We illustrate several cases with supercritical and subcritical transition in figure 1. Here parameter μ is not too small, nevertheless the correspondence between the analytical and numerical results is very good. One can see that for a narrow distribution (γ=0.1), the synchronization transition is supercritical (solid curve), whereas for the broader distributions, it is subcritical.  figure 1 for the Gaussian distribution. The transition is supercritical for a narrow distributions of natural frequencies and subcritical for larger values of γ.

Bimodal distribution
For a bimodal distribution  we have m < 0 * , which means that in this range of parameters the stationary branch bifurcates subcritically for any masses.
Noteworthy, for the bimodal distribution, nonstationary (time-periodic) solutions can dominate the transition, as have been demonstrated in [12,23], and the above-presented analysis of stationary states solves only a part of the problem.

Uniform distribution
Here we consider a uniform distribution In the noise-free case, i.e. at σ=0, it follows from (4.1) that The critical value of the coupling is, in the first order in μ, e g mg p p = + 4 1 2 . A remarkable feature of this distribution is that it demonstrates a discontinuous transition without hysteresis for μ=σ=0 [24]: at ε=4γ/π the order parameter r jumps from zero to π/4. Both small noise and small inertia destroy this degeneracy, although in different ways. For small values of μ and vanishing noise σ→0, the transition is subcritical, see figure 3(a), where we compare the analytical result(35) with direct numerical simulations. For small noise, in the massless case, the transition is supercritical, but it turns to be subcritical beyond the critical mass m* given by(34). This situation is illustrated in figures 3(b), (c). Direct numerical simulations clearly show a hysteresis and regions of bistability, where both the asynchronous state and the stationary synchronous branch are stable. Here, one can also see different effects of finite-size fluctuations on the transitions to synchrony and back: the asynchronous state is much more sensitive to these fluctuations, which results in a shift of the transition point to smaller values of coupling.

Generalizations
Above we focused on the case where the rotators differ solely by their natural frequencies ω n , and the distribution of these frequencies is symmetric. Often, it is desirable to consider a more general situation, where all the parameters governing the dynamics of rotators are different (see a similar generalization of the Kuramoto model in [25]): The goal is to find uniformly rotating solutions y = =W r const,  in the thermodynamic limit. For the given distribution of individual parameters m w e a s , , , , n n n n n , and for values of global parameters E and B, we look for a solution r, Ω (multiple solutions are also possible). Similar to the consideration above, it is more convenient to fix Ω and A=rE, and to find the solution in a parametric form where w m e s a G , , , , ( ) is a joint distribution density over the parameters of the problem. Expression (38) yields the representation of the solution in the explicit parametric form Clearly, if only some of the parameters are distributed, general expressions(37), (38) can be simplified.

Conclusion
In conclusion, we have developed an analytical description of stationary synchronous regimes in a population of noisy globally coupled rotators ('oscillators with inertia'). The main analytical formula is expression (23) for the subgroup order parameter of oscillators having detuning ν to the frequency of the mean field. This expression contains only normalized detuning ν/σ 2 and forcing A/σ 2 , and is valid for small masses μ. We also discussed different limiting cases which can be straightforwardly derived from this formula. For example, for massless rotators, i.e. for the standard Kuramoto oscillators, we obtain an exact expression (25). This provides an analytical formula for stationary solutions for the Kuramoto model (or, more generally, for the Kuramoto-Sakaguchi model) with noise. Furthermore, we obtained an analytical expression for the critical mass of rotators, at which a change of the synchronization transition type (supercritical versus subcritical) occurs. Moreover, we applied general analytical formulations containing integrals over a general distribution of natural frequencies, to several commonly used distributions. There, the expressions for the critical mass reduce to algebraic analytical formulae (see section 5). Our approach is restricted to stationary solutions only, it does not capture possible regimes with a nonconstant modulus of the order parameter. The latter are essential for multimodal distributions of natural frequencies, in particular for the bimodal distribution considered in section 5.3, where periodic regimes dominate the transition to synchrony. Another restriction of our approach is that it does not provide stability of the nontrivial solutions: the corresponding linearized equations have to be explored numerically (even stability analysis of the trivial asynchronous state is rather involved, see [12]).
Finally, we have shown that the approach can be directly generalized to the case where not only the natural frequencies of rotators are different, but also their masses (see [17]), coupling strengths and phase shifts in coupling (see [25]), or even noise intensities. Such a generalization can be useful for applications of mean field theory to random network couplings, e.g. the diversity of incoming degrees can be modeled as a distribution of effecting coupling constants.

Appendix. Derivation of the subgroup order parameter by virtue of moment expansion
We start with a general reduction of a second-order stochastic equation mj j j sz We employ the moment method described in [26]. First, we introduce the moments in the velocity:     Finally, in the first order in μ, we obtain the following Fokker-Planck equation for the distribution density of the phases r j j º t W t , , Next, similarly to the procedure described in section 2, we transform to the rotating reference frame by introducing q j y = -, where y = W  . In this reference frame we set n q = -F A sin and look for a stationary solution r q 0 ( ), which satisfies the following equation