Analytical approach to synchronous states of globally coupled noisy rotators
V.O.Munyaev, L.A.Smirnov, V.A.Kostin, G.V.Osipov, A.Pikovsky
aa r X i v : . [ n li n . AO ] D ec Analytical approach to synchronous states ofglobally coupled noisy rotators
V O Munyaev , L A Smirnov , , V A Kostin , , G V Osipov and A Pikovsky , ‡ Department of Control Theory, Nizhny Novgorod State University, Gagarin Av. 23,Nizhny Novgorod, 603950, Russia. Institute of Applied Physics, Russian Academy of Sciences, Ul’yanova Str. 46,Nizhny Novgorod, 603950, Russia. Institute for Physics and Astronomy, University of Potsdam, Karl-Liebknecht-Str.24/25, 14476 Potsdam-Golm, Germany.E-mail: [email protected]
Abstract.
We study populations of globally coupled noisy rotators (oscillatorswith inertia) allowing a nonequilibrium transition from a desynchronized state to asynchronous one (with the non-vanishing order parameter). The newly developedanalytical approaches resulted in solutions describing the synchronous state withconstant order parameter for weakly inertial rotators, including the case of zeroinertia, when the model is reduced to the Kuramoto model of coupled noise oscillators.These approaches provide also analytical criteria distinguishing supercritical andsubcritical transitions to the desynchronized state. All the obtained analyticalresults are confirmed by the numerical ones, both by direct simulations of the largeensembles and by the solution of the associated Fokker-Planck equation. We alsopropose generalizations of the developed approaches for setups where different rotatorsparameters (natural frequencies, masses, noise intensities, strengths and phase shiftsin coupling) are dispersed.
Keywords : coupled rotators, synchronization transition, hysteresis, Kuramoto model,noisy systems
Submitted to:
New J. Phys. ‡ Author to whom any correspondence should be addressed. nalytical approach to synchronous states of globally coupled noisy rotators
1. Introduction
The Kuramoto model of globally coupled phase oscillators [1] is a paradigmatic modelto study synchronization phenomena [2–5]. In the thermodynamic limit, stationary (ina proper rotating reference frame) synchronized states can be found analytically for anarbitrary distribution of natural frequencies [6]. The idea is that for any subgroup ofoscillators having a certain natural frequency, a stationary distribution in a constantmean 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 ofthese subgroup order parameters. Other analytical approaches for the Kuramoto modelmay be not restricted to stationary states, but usually have other restrictions sincethey apply only for identical oscillators like the Watanabe–Strogatz theory [7], or for apopulation with a Lorentzian distribution of natural frequencies like the Ott–Antonsenansatz [8].An important generalization of the Kuramoto model considers an ensemble ofglobally coupled rotators (that is, oscillators with inertia). It is in particular relevantfor modelling power grid networks [9, 10]. Synchronization features of globally coupledrotators, both dterministic and noisy, have been widely studied [11–16]; in particular, anapproach similar to the analytical description [6] was developed [3, 17, 18]. However, fornoisy coupled rotators, so far the stationary distributions were found only numerically.Similarly, there are no analytical formulae for the Kuramoto model with noise.In this paper, we present a fully analytical approach to the problem of globallycoupled noisy rotators for small intertia (small “masses”). We analyse a stationarysolution of the Fokker–Planck (Kramers) equation describing identical noisy rotatorsdriven by the mean field, and obtain a closed expression for the order parameter forthis subgroup. Then, for an arbitrary distribution of natural frequencies, a stationarysolution 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 orderparameter and use it to classify the transition to synchrony as a supercritical or asubcritical one.The paper is organized as follows. We introduce the model of coupled noisyrotators in Section 2, where we also formulate stationary equations for the distributiondensity. A solution of these equation in the limit of small inertia terms is presentedin Section 3 (another derivation of this solution is given in Appendix Appendix A).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 givegeneral 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 6we show how the approach can be used for generic ensembles where not only naturalfrequencies, but other relevant parameters of rotators (masses, noise strengths, etc.) aredistributed. We conclude with Section 7. nalytical approach to synchronous states of globally coupled noisy rotators
2. Noisy coupled rotators
In this paper we consider an ensemble of N globally coupled rotators characterized bytheir angles ϕ n and velocities ˙ ϕ n ( n = 1 , , . . . , N ). The rotators are coupled via thecomplex mean field R ≡ re iψ = 1 N N X n =1 e iϕ n (1)and obey equations of motion µ ¨ ϕ n + ˙ ϕ n = ω n + εN N X n ′ =1 sin( ϕ n ′ − ϕ n ) + σξ n ( t )= ω n + εr sin( ψ − ϕ n ) + σξ n ( t ) . (2)The unit of time is chosen so that the coefficient at the friction term ( ∼ ˙ ϕ n ) is one andall the parameters are dimensionless (normalized by the friction coefficient). Parameter µ describes the mass of rotators (more precisely, µ should be called the moment ofinertia of a rotator); below, we focus on the overdamped case µ ≪
1. Parameter ε is the coupling strength. Parameters ω n describe torques acting on rotators; weassume them to be distributed with a density g ( ω ). Note, that our approach can bestraightforwardly generalized to the case where other parameters are distributed asdiscussed in Section 6. Because uncoupled rotators have mean angular velocities ω n , wespeak of “natural frequencies” instead of “torques”. We do so also to make transparenta relation to the particular case µ = 0, where system (2) is nothing else but the standardKuramoto model for globally coupled phase oscillators with a distribution g ( ω ) of naturalfrequencies ω . The rotators are acted upon by the independent white Gaussian noiseforcings σξ n ( t ) with equal amplitudes σ , zero means h ξ n ( t ) i = 0, and auto-correlations h ξ n ( t ) ξ n ′ ( t ) i = 2 δ nn ′ δ ( t − t ) (where δ nn ′ denotes the Kronecker delta, and δ ( t ) is theDirac δ -function). Whereas equations (2) are used for numerical simulations below, theanalytical approach is developed in the thermodynamics limit N → ∞ .In the mean-field coupling models, synchronization is one of the fundamental effectswhich is relevant to many systems. A transition to synchrony can be fully characterizedin terms of the Kuramoto order parameter (1). Physically, the amplitude r of thiscomplex parameter describes the synchrony level of the elements belonging to thepopulation considered. Nonvanishing r indicates the collective synchronization. So, it isimportant to develop an universal analytical approach allowing to calculate r value andpredict the global evolution of an ensemble consisting of macroscopically large numberof interacting units.Note that model (2) is not the most general one – it might include a phaseshift in coupling (which corresponds to the Kuramoto–Sakaguchi overdamped system).We postpone the discussion of this case to Section 6. Furthermore, for simplicityof presentation we consider in Sections 2–5 only symmetric unimodal frequencydistributions; the discussion of asymmetric distributions is also postponed to Section 6. nalytical approach to synchronous states of globally coupled noisy rotators N → ∞ .Furthermore, we look for stationary synchronous states, i.e., those with constantmodulus of the order parameter and with a uniformly rotating angle: r = const,˙ ψ = Ω. Such states are possible in the thermodynamics limit only where the finite-size fluctuations vanish. It is convenient to introduce a new angle variable related tothe angle of the mean field θ = ϕ − ψ , ˙ θ = u = ˙ ϕ − Ω. Then, model (2) can be formallyrewritten as an infinite-dimensional system µ ¨ θ n + ˙ θ n = ω n − Ω − εr sin( θ n ) + σξ n ( t ) , r = h e iθ n i . (3)Stochastic equation (3) allows one to use the Fokker–Planck (Kramers) equation for theprobability density P ( θ, u, t | ω ) for a subset of rotators having natural frequency ω andto express the mean field as an integral over this density, ∂ t P + ∂ θ ( uP ) + µ − ∂ u [( − u + ω − Ω − εr sin θ ) P ] = µ − σ ∂ u P,r = Z d ω g ( ω ) Z d θ Z d u e iθ P ( θ, u, t | ω ) . Because all the deterministic forces in this reference frame are constant, this densityevolves toward the stationary, time-independent, one. As the only parameters governingthis distribution are detuning ν = ω − Ω and forcing term A = εr , we seek this stationarydistribution as a function depending on these parameters explicitly, P ( θ, u | ν, A ).Then, the system takes the form u∂ θ P = µ − ∂ u (( u − ν + A sin θ ) P ) + µ − σ ∂ u P , (4) r = Z d ν g (Ω + ν ) Z d θ Z d u e iθ P ( θ, u | ν, A ) , A = εr. (5)The system (4) and (5) is in fact a self-consistent system for determining the unknownorder parameter r (frequency Ω is obtained from the condition that r is real).Below, we restrict ourselves to symmetric distributions only, g ( ω + ν ) = g ( ω − ν ).Then, one may chose Ω = ω , which corresponds to the symmetric stationary density P , P ( θ, u | ν, A ) = P ( − θ, − u | − ν, A ). In this case, the integral in (5) is real, whichproofs that ω is a correct value of the frequency of the mean field Ω. Note that therecould also be other appropriate values of Ω besides ω . These different values correspondto 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 , r ( A ) = Z d ν d θ d u g ( ω + ν ) e iθ P ( θ, u | ν, A ) , ε ( A ) = Ar ( A ) . The only remaining act here is finding solution P of (4). In Refs. [17, 19], a method ofmatrix continuous fractions was employed to solve (4) numerically, another numericalapproach is described in [18]. In the next Section 3, we present an analytical solutionfor the overdamped case µ ≪ nalytical approach to synchronous states of globally coupled noisy rotators
3. 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 from (4) based on series representation. The second method is given inAppendix Appendix A. According to [12, 17], we represent a stationary solution P ( θ, u | ν ) of (4) as a doubleseries in the parabolic cylinder functions Φ p ( u ) in u and the Fourier modes e iqθ in θ , P ( θ, u | ν, A ) = (2 π ) − / Φ ( u ) + ∞ X p =0 + ∞ X q = −∞ a pq ( ν, A ) Φ p ( u ) e iqθ . (6)Here functions Φ p ( u ) are defined asΦ p ( u ) = r κ p p ! √ π exp (cid:2) − κ u / (cid:3) H p ( κ u ) , where κ = p µ/ σ and H p ( κ u ) are the Hermite polynomials. Substituting expansion(6) in (4), (5) and using orthogonality conditions for the basis functions, we obtain theinfinite system for unknown coefficients a pq ( ν, A ) where the order parameter r is justproportional to one of the expansion coefficients, σ r p + 1 µ iqa p +1 ,q + pµ a pq − σ r pµ (cid:20)(cid:16) νσ − iq (cid:17) a p − ,q − i A σ ( a p − ,q +1 − a p − ,q − ) (cid:19) = 0 , (7) r = √ π Z d ν g ( ω + ν ) a ∗ , ( ν, A ) . (8)Thus, the main challenge is to find the quantity a , ( ν, A ) via solving (7); this coefficientis nothing else but the order parameter for the group of oscillators having naturalfrequency ω + ν ; therefore, we call it subgroup order parameter . It is then used in (8) tofind the total order parameter r via averaging over ν . Below, we also use the symmetryproperty following from (7), a pq ( − ν, A ) = ( − p a ∗ pq ( ν, A ).For the fixed parameters A and ν , the set of equations for a pq can be formally solvedby virtue of the matrix continuous fractions method [17, 19]. According to this method,the expression for the subgroup order parameter a , ( ν, A ) is as follows: a , ( ν, A ) = S , ( ν, A ) √ πS , ( ν, A ) , (9)where S jk ( ν, A ) is the element of the infinite matrix S , which is given by the recurrenceformula S = ˜ D − I − µ D (cid:18) I − µ D (cid:16) I − µ D ( I − . . . ) − ˜ D (cid:17) − ˜ D (cid:19) − ˜ D ! . (10) nalytical approach to synchronous states of globally coupled noisy rotators I is the identical matrix, and the infinite diagonal matrix D and the infinitetri-diagonal matrix ˜ D are defined as D jk = ikδ jk , ˜ D jk = (cid:0) iσ k − ν (cid:1) δ jk + i A δ j,k − − δ j,k +1 ) . (11)The main steps of the numerical procedure based on relations (9)–(11) employedfor 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 andhas 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 g ( ω ), e.g.,for the Lorentz distribution. So, it is important to develop an analytical approximationfor the calculation of the value r in a nonequilibrium stationary state. Expression (10) allows for a perturbative expansion in the small parameter µ . In thefirst order in µ , we obtain S = ˜ D − (cid:16) I − µ D ˜ D (cid:17) + o ( µ ) . (12)In order to evaluate the inverse matrix ˜ D − , we denote the principal minors as M jk = ˜ D (cid:20) j j + 1 . . . kj j + 1 . . . k (cid:21) , (13)and introduce a cutoff (cyclic) frequency d → + ∞ for Fourier modes. Then, the elementsof the inverse matrix ˜ D − take the form˜ D − jk = lim d → + ∞ (cid:18) i A (cid:19) j − k M − d,k − M j +1 ,d det ˜ D for j ≥ k, (14)˜ D − jk = ( − j + k ˜ D − kj for j < k. According to the Laplace theorem, for all j ≥ k (for convenience, we take M k +1 ,k = 1and M k +2 ,k = 0), the following representation holdsdet ˜ D = M kj M − d,k − M j +1 ,d + (cid:18) i A (cid:19) M k,j − M − d,k − M j +2 ,d + (cid:18) i A (cid:19) M k +1 ,j M − d,k − M j +1 ,d + (cid:18) i A (cid:19) M k +1 ,j − M − d,k − M j +2 ,d . (15)Substituting (15) in (14), we obtain for j ≥ k that˜ D − jk = lim d → + ∞ (cid:18) i A (cid:19) j − k M kj + (cid:18) i A (cid:19) M k +1 ,j M − d,k − M − d,k − + (cid:18) i A (cid:19) M k,j − M j +2 ,d M j +1 ,d + (cid:18) i A (cid:19) M k +1 ,j − M − d,k − M − d,k − M j +2 ,d M j +1 ,d ! − . (16)The principal minor M jk (13) is found from a reccurent equation with fixed j , i.e.,with fixed bottom right corner, M kj = (cid:0) iσ k − ν (cid:1) M k +1 ,j + (cid:18) i A (cid:19) M k +2 ,j . (17) nalytical approach to synchronous states of globally coupled noisy rotators M kj = (cid:18) i A (cid:19) j − k +1 I − j − − i νσ (cid:0) Aσ (cid:1) K − k +1 − i νσ (cid:0) − Aσ (cid:1) − K − j − − i νσ (cid:0) − Aσ (cid:1) I − k +1 − i νσ (cid:0) Aσ (cid:1) I − j − − i νσ (cid:0) Aσ (cid:1) K − j − i νσ (cid:0) − Aσ (cid:1) − K − j − − i νσ (cid:0) − Aσ (cid:1) I − j − i νσ (cid:0) Aσ (cid:1) . (18)Here I z and K z denote the modified Bessel functions of the first and second kind,respectively, of the order z . Using properties of the modified Bessel functions, weevaluate the limitlim d → + ∞ I d − i νσ (cid:0) Aσ (cid:1) K d − i νσ (cid:0) − Aσ (cid:1) = 0 . (19)Employing both (18) and (19), we obtainlim d → + ∞ M − d,k − M − d,k − = − (cid:18) i A (cid:19) − I − k +1 − i νσ (cid:0) Aσ (cid:1) I − k − i νσ (cid:0) Aσ (cid:1) . (20)Similarly, we find a representation of the principal minor M kj , using its expansionat fixed k , i.e., at fixed upper left corner, M kj = (cid:0) iσ j − ν (cid:1) M k,j − + (cid:18) i A (cid:19) M k,j − , which yields M kj = (cid:18) i A (cid:19) j − k +1 I k − i νσ (cid:0) − Aσ (cid:1) K j +1+ i νσ (cid:0) Aσ (cid:1) − K k − i νσ (cid:0) Aσ (cid:1) I j +1+ i νσ (cid:0) − Aσ (cid:1) I k − i νσ (cid:0) − Aσ (cid:1) K k + i νσ (cid:0) Aσ (cid:1) − K k − i νσ (cid:0) Aσ (cid:1) I k + i νσ (cid:0) − Aσ (cid:1) . (21)Using expressions (21) and (19), we findlim d → + ∞ M j +2 ,d M j +1 ,d = − (cid:18) i A (cid:19) − I j +1+ i νσ (cid:0) − Aσ (cid:1) I j + i νσ (cid:0) − Aσ (cid:1) . (22)Combination of expressions (11), (12), (13), (16), (20), and (22) with the generalformula (9) results in the closed-form formula to the first order in the mass parameter µ , a , ( ν, A ) = 1 √ π I i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) − µ σ π sin (cid:0) iπ νσ (cid:1) I − i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) ! + o ( µ ) . (23)An alternative derivation of our main result (23) is based on the method of momentsfor the density and on elimination of the velocity; it is presented in Appendix A.
4. 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. nalytical approach to synchronous states of globally coupled noisy rotators For purely deterministic rotators, we have to set σ →
0. To find the noise-free limit ofgeneral expression (23), it is convenient to rewrite it in an equivalent form a , ( ν, A ) = 1 √ π I i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) − µ A I i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) − I − − i νσ (cid:0) Aσ (cid:1) I − i νσ (cid:0) Aσ (cid:1) !! + o ( µ ) . (24)So, it is neccesary to calculate limits at σ → σ → I i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) = − i νA + r − ν A , ν > − A − i νA − r − ν A , ν ≤ − A ,lim σ → I − − i νσ (cid:0) Aσ (cid:1) I − i νσ (cid:0) Aσ (cid:1) = − i νA + r − ν A , ν ≤ A − i νA − r − ν A , ν > A ,The resulting expression for the subgroup order parameter (in the first order in µ ) is a , ( ν, A ) = o ( µ ) + 1 √ πA i (cid:16) − ν − √ ν − A (cid:17) (cid:16) iµ √ ν − A (cid:17) , ν < − A , (cid:16) − iν + √ A − ν (cid:17) , | ν | ≤ A , i (cid:16) − ν + √ ν − A (cid:17) (cid:16) − iµ √ ν − A (cid:17) , ν > A . 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 thermodynamicslimit (we can now omit the first index) as follows: a ( ν, A ) = 1 √ π I i νσ (cid:0) Aσ (cid:1) I i νσ (cid:0) Aσ (cid:1) . (25)Correspondingly, the expression for the full order parameter is also exact, r = Z d ν g ( ω + ν ) I − i νσ (cid:0) Aσ (cid:1) I − i νσ (cid:0) Aσ (cid:1) . (26)The subgroup order function (25) has no poles in lower half-plane, and the integral in(26) can be evaluated for suitable distributions g ( ω ) via residues. For example, for theLorentz distribution g ( ω + ν ) = γπ ( ν + γ ) , (27)the resulting expression for r reads r = I γσ (cid:0) Aσ (cid:1) I γσ (cid:0) Aσ (cid:1) . (28) nalytical approach to synchronous states of globally coupled noisy rotators
9A more elaborated application of the theory to the Kuramoto model with noise will bepresented elsewhere.
Here, we employ the general parametric representation of the order parameter as afunction of the coupling constant, r = √ π Z d ν g ( ω + ν ) a ∗ , ( ν, A ) , ε = Ar (29)with a , ( ν, A ) given by (23) to characterize the synchronization transition, i.e., tocharacterize states with the order parameter close to zero. All formulas below arevalid in the first order in small mass µ like (23) and are therefore approximate, but forthe sake of simplicity we write them as exact relations below.We expand r ( A ) in the Taylor series for small values of A . This expansion containsonly odd powers of A since (23) is an odd function of A , r ( A ) = C A + C A + C A + . . . . (30)The three initial coefficients in expansion (30) are as follows: C = 12 Z d y g ( ω + σ y )1 + y (cid:0) − µσ y (cid:1) ,C = − σ Z d y g ( ω + σ y )(1 + y ) (4 + y ) (cid:0) (cid:0) − y (cid:1) − µσ y (cid:0)
13 + y (cid:1)(cid:1) , (31) C = 116 σ × Z d y g ( ω + σ y )(1 + y ) (4 + y ) (9 + y ) (cid:0) (cid:0) − y + 4 y (cid:1) + µσ y (cid:0) −
113 + 32 y + y (cid:1)(cid:1) . In the massless case µ = 0, the expressions for C , C have been obtained in [20].The nontrivial branch of solutions r ( ε ) starts at ε (1) c = 1 /C . In the noiseless limit, σ →
0, the critical value of the coupling parameter can beexpressed as ε (1) c (cid:12)(cid:12) σ =0 = 2 πg ( ω ) − µ . This expression coinsides, in the first order in µ , with the result for deterministic rotatorsobtained in [16].It is instructive to compare this critical value with the expression for the stabilityloss of the asynchronous state r = 0 (equation (24) in [12]). In our notations, thisgeneral formula for the imaginary part of the eigenvalue x , valid also for large masses µ , reads 1 = µ ˆ εe µσ ∞ X p =0 ( − µσ ) (cid:16) pµσ (cid:17) p ! Z ∞−∞ d ω g ( ω ) µσ + p + i ( µω + x ) . (32) nalytical approach to synchronous states of globally coupled noisy rotators µ , this experssion can be simplified. Assuming x = x + µx + µ x andsubstituting this in (32), we find x = x = 0. For a unimodal frequency distributionsymmetric around ω , a solution x = − ω exists, which yieldsˆ ε − = σ Z ∞−∞ d ω g ( ω ) σ + ( ω − ω ) − µ Z ∞−∞ d ω g ( ω )( ω − ω ) σ + ( ω − ω ) + o ( µ ) . One can easily see that ˆ ε = ε (1) c , which means that the branch of stationary solutionsjoins the axis r = 0 in the ( r, ε ) plane exactly where the instability of the asynchronousstate first occurs.Depending on the sign of the coefficient C , there are two possibilities: Supercritical transition occurs for C <
0. Here one observes a continuous (second-order) transition with the solution branch r = C r(cid:16) ε (1) c − ε (cid:17) /C existing for ε > ε (1) c . Subcritical transition occurs for C >
0. Here, the branch of solutions exists for ε < ε (1) c . If C <
0, one can estimate that this branch spreads till the minimal value at ε (2) c ≈ (cid:18) C − C C (cid:19) − . 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 . Becausethis coeffictient depends on the mass µ , there is a critical value of the mass at whichthe type of the transition changes, µ ∗ = 2 σ (cid:18)Z d y g ( ω + σ y ) (1 − y )(1 + y ) (4 + y ) (cid:19) (cid:18)Z d y g ( ω + σ y ) y (13 + y )(1 + y ) (4 + y ) (cid:19) − . (33)We stress that this expression is valid only if the value of µ ∗ is numerically small: µ ∗ ≪
5. Examples
Here we present explicit calculations of the synchronization transition for severalcommonly explored distributions of the frequencies. Each of these distributions ischaracterized by a width γ , and in all cases the result depends on the dimensionlessparameter ξ = γ/σ , which measures relative influence of the frequency dispersion andnoise on the transition. nalytical approach to synchronous states of globally coupled noisy rotators The Gaussian distribution of the rotators’ frequencies is g ( ω ) = 1 p πγ e − ( ω − ω γ . For this distribution, the coefficients are C = 1 σ (cid:18)r π ξ e ξ erfc (cid:18) √ ξ (cid:19) (cid:0) µσ (cid:1) − µσ (cid:19) ,µ ∗ = 1 σ πe ξ ( ξ −
1) erfc (cid:16) √ ξ (cid:17) + 2 ξ (cid:16) √ π (1 − ξ ) + 3 πe ξ ξ erfc (cid:16) √ ξ (cid:17)(cid:17) ξ (cid:16) √ πξ (10 + ξ ) − π (cid:16) e ξ erfc (cid:16) √ ξ (cid:17) + 3 e ξ erfc (cid:16) √ ξ (cid:17)(cid:17)(cid:17) , where erfc is the complementary error function. We illustrate several cases withsupercritical and subcritical transition in Figure 1. Here parameter µ is not too small,nevertheless the correspondence between the analytical and numerical results is verygood. One can see that for a narrow distribution ( γ = 0 . r . . . . ε/γ .
62 1 .
65 1 .
68 1 .
71 1 .
74 1 . γ = 0 . γ = 0 . γ = 0 . γ = 0 . Figure 1.
Branches of the desynchronized stationary solutions r ( ε ) for the Gaussiandistribution of frequencies with different widths γ : the stationary order parameter r vs normalized coupling ε/γ is shown. The other parameters are µ = 0 . σ = 0 . µ ; markers areobtained through numerical solutions of (9). Small deviations of the analytical solutionfrom the numerical one are noticeable for large values of γ only. In the case of the Lorentz distribution (27), the expressions for characterization of thesynchronization transition are as follows: C = 12 σ − µσ ξ ξ , C = 18 σ − µσ ξ (3 + ξ )(1 + ξ ) (2 + ξ ) , µ ∗ = 1 σ ξ ( ξ + 3) . nalytical approach to synchronous states of globally coupled noisy rotators γ . r . . ε/γ .
08 2 .
12 2 .
16 2 . . γ = 0 . γ = 0 . γ = 0 . γ = 0 . Figure 2.
Branches of stationary solutions r ( ε ) for µ = 0 . σ = 0 .
05, and different γ for the Lorentz distribution of frequencies. Only the analytical solutions obtainedfrom (23) are presented, because precise numerical solution of the full problem is harddue to the broad tails of the distribution. For a bimodal distribution g ( ω ) = 12 δ ( ω − ω − γ ) + 12 δ ( ω − ω + γ ) , we find C = 12 σ − µσ ξ ξ , C = 18 σ (4 + 13 z ) ξ − µσ ξ (1 + ξ ) (4 + ξ ) , µ ∗ = 2 σ − ξ ξ (13 + ξ ) . Furthermore, here the integral in (29) can be evaluated in terms of the modified Besselfunctions, r = ℜ I iξ (cid:0) Aσ (cid:1) I iξ (cid:0) Aσ (cid:1) − µ σ π sin( iπξ ) I − iξ (cid:0) Aσ (cid:1) I iξ (cid:0) Aσ (cid:1) !! + o ( µ ) , Remarkably, for γ > σ / √ µ ∗ <
0, which means that in this range ofparameters the stationary branch bifurcates subcritically for any masses.Noteworthy, for the bimodal distribution, non-stationary (time-periodic) solutionscan dominate the transition, as have been demonstrated in [12, 20], and the above-presented analysis of stationary states solves only a part of the problem. nalytical approach to synchronous states of globally coupled noisy rotators Here we consider a uniform distribution g ( ω ) = γ , | ω − ω | ≤ γ ,0 , | ω − ω | > γ .Using (31) and (33), we obtain C = 1 σ (1 + µσ ) arctan ξ − µσ ξ ξ ,C = − σ (1 + 2 µσ ) ξ + (1 + ξ ) (cid:0) (1 + 2 µσ ) arctan ξ − (1 + 3 µσ ) arctan ξ (cid:1) ξ (1 + ξ ) , (34) µ ∗ = − σ ξ + ( ξ + 1) (cid:0) arctan ξ − arctan ξ (cid:1) ξ + ( ξ + 1) (cid:0) ξ − ξ (cid:1) . In the noise-free case, i.e., at σ = 0, it follows from (4.1) that r ( A ) = o ( µ ) + A π γ − µ (cid:18) − Aγ (cid:19) γ − A + 2 p γ − A γ + p γ − A ! , A ≤ γ ,12 r − γ A + Aγ arcsin (cid:16) γA (cid:17)! , A > γ . (35)The critical value of the coupling is ε (1) c = 4 γ (1 + 2 µγ/π ) /π + o ( µ ).A remarkable feature of this distribution is that it demonstrates a discontinuoustransition without hysteresis for µ = σ = 0 [21]: 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, thetransition 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 transitionis supercritical, but it turns to be subcritical beyond the critical mass µ ∗ given by (34).This situation is illustrated in Figure 3(b,c). Direct numerical simulations clearlyshow a histeresis and regions of bistability, where both the asynchronous state andthe stationary synchronous branch are stable. Here, one can also see different effects offinite-size fluctuations on the transitions to synchrony and back: the asynchronous stateis much more sensitive to these fluctuations, which results in a shift of the transitionpoint to smaller values of coupling.
6. 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 toconsider a more general situation, where all the parameters governing the dynamics of nalytical approach to synchronous states of globally coupled noisy rotators (c) r . . . . ε .
32 1 .
34 1 .
36 1 . N = 2 . · , → N = 2 . · , ← N = 5 . · , → N = 5 . · , ← Figure 3. (a) Curves: branches of stationary solutions r ( ε ) for different µ in thenoise-free case obtained analytically from equation (35). Connected markers: branchesobtained from direct numerical simulations of the ensemble of N = 10000 rotatorswith µ = 0 . ε gradually increases (’ → ’ label,diamond and square markers) or decreases (’ ← ’ label, round and triangle markers). Inpanels (b) and (c), results of similar numerical experiments are shown for rotators withmasses µ = 0 .
1, noise amplitude σ = 0 .
2, and two different widths of the distribution, γ = 2 in (b) and γ = 1 in (c). Comparison with theoretical predictions (dashedlines) shows particularly strong effect of finite-size fluctuations on the discontinuoustransition “asynchrony → synchrony”; the reverse transition is nearly at the pointpredicted by analytical theory even for not too large populations. rotators are different (cf. a similar generalization of the Kuramoto model in [22]): µ n ¨ ϕ n + ˙ ϕ n = ω n + ǫ n N N X n ′ =1 sin( ϕ n ′ − ϕ n − β n ) + σ n ξ n ( t )= ω n + ǫ n r sin( ψ − ϕ n − β n ) + σ n ξ n ( t ) . (36)Here, we introduce two global parameters: E is the average strength of coupling, and B is the characteristic phase shift in coupling, ǫ n = E + ε n , β n = B + α n . The goal is to find uniformly rotating solutions r = const, ˙ ψ = Ω in the thermodynamic nalytical approach to synchronous states of globally coupled noisy rotators µ n , ω n , ε n , α n , σ n , and forvalues of global parameters E and B , we look for a solution r, Ω (multiple solutions arealso possible). Similar to the consideration above, it is more convenient to fix Ω and A = rE , and to find the solution in a parametric form E = E ( A, Ω), B = B ( A, Ω), r = A/E ( A, Ω).To accomplish this, we introduce the rotating variable θ = ϕ − ψ − B − α . Then(36) takes form µ ¨ θ + ˙ θ = ω − Ω − Aε sin( θ ) + σξ ( t ) . The stationary distribution for this Langevin equation was analyzed in Section 3 above,it yields the following subgroup order parameter a , ( ω, µ, ε, σ | A, Ω) = 1 √ π I i ω − Ω σ (cid:0) Aεσ (cid:1) I i ω − Ω σ (cid:0) Aεσ (cid:1) − µ σ π sin (cid:0) iπ ω − Ω σ (cid:1) I − i ω − Ω σ (cid:0) Aεσ (cid:1) I i ω − Ω σ (cid:0) Aεσ (cid:1) ! + o ( µ ) . (37)Substituting this solution in the expression for the global order parameter r , we obtain re − iB = Z d ω d µ d ε d σ d α a ∗ , ( ω, µ, ε, σ | A, Ω) e iα G ( ω, µ, ε, σ, α ) (38)where 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 r = r ( A, Ω) , B = B ( A, Ω) , E = Ar ( A, Ω) . Clearly, if only some of the parameters are distributed, general expressions (37), (38)can be simplified.
7. Conclusion
In conclusion, we have developed an analytical description of stationary synchronousregimes in a population of noisy globally coupled rotators (“oscillators with inertia”).The main analytical formula is expression (23) for the subgroup order parameter ofoscillators having detuning ν to the frequency of the mean field. This expression containsonly normalized detuning ν/σ and forcing A/σ , and is valid for small masses µ . Wealso discussed different limiting cases which can be straightforwardly derived from thisformula. For example, for massless rotators, i.e., for the standard Kuramoto oscillators,we obtain an exact expression (25). This provides an analytic expression for stationarysolutions for the Kuramoto model (or, more generally, for the Kuramoto–Sakaguchimodel) with noise.Our approach is restricted to stationary solutions only, it does not capture possibleregimes with a non-constant modulus of the order parameter. The latter are essentialfor multimodal distributions of natural frequencies, in particular for the bimodaldistribution considered in Section 5.3, where periodic regimes dominate the transitionto synchrony. Another restriction of our approach is that it does not provide stabilityof the nontrivial solutions: the corresponding linearized equations have to be explored nalytical approach to synchronous states of globally coupled noisy rotators Acknowledgments
We thank D. Goldobin for useful discussions. Results presented in Sections 2 and 3were supported by the RSF grant No. 17-12-01534. Results presented in Sections 4,5 and 6 were supported by the RSF grant No. 19-12-00367. Results presented inAppendix Appendix A and numerical simulations presented in Section 5 were supportedby the RFBR grant No. 19-52-12053.
Appendix A. Derivation of the subgroup order parameter by virtue ofmoment expansion
We start with a general reduction of a second-order stochastic equation µ ¨ ϕ + ˙ ϕ = F ( ϕ, t ) + σξ ( t )to a first-order equation, valid for small µ . The probability density P ( ϕ, ˙ ϕ, t ) obeys thecorresponding Fokker–Planck equation ∂ t P + ∂ ϕ ( ˙ ϕP ) + ∂ ˙ ϕ (cid:0) µ − ( − ˙ ϕ + F ) P (cid:1) = µ − σ ∂ ϕ P. (A.1)We employ the moment method described in [23].First, we introduce the moments in the velocity: w m ( ϕ, t ) = Z ∞−∞ d ˙ ϕ ˙ ϕ m P ( ϕ, ˙ ϕ, t ) . According to (A.1), these moments obey following equations: ∂ t w + ∂ ϕ w = 0 , (A.2) w + µ∂ t w = F w − µ∂ ϕ w , (A.3) w m + µm ∂ t w m = F w m − − µm ∂ ϕ w m +1 + ( m − σ µ w m − for m ≥ . (A.4)In order to employ the smallness of parameter µ , it is convenient to introduce rescaledmoments w m = µ m/ W m for even m ,1 µ ( m − / W m for odd m . nalytical approach to synchronous states of globally coupled noisy rotators W m , equations (A.2)–(A.4) can be rewritten in a form free ofsingularities ∂ t W + ∂ ϕ W = 0 ,W = F W − ∂ ϕ W − µ∂ t W ,W m = ( m − σ W m − + µ (cid:18) F W m − − m ∂ ϕ W m +1 − m ∂ t W m (cid:19) for even m,W m = ( m − σ W m − + F W m − − m ∂ ϕ W m +1 − µm ∂ t W m for odd m. As we are looking for the first order corrections in µ , we rewrite this system keeping therelevant terms only ∂ t W + ∂ ϕ W = 0 ,W = F W − ∂ ϕ W − µ∂ t W ,W = σ W + µ (cid:18) F W − ∂ ϕ W − ∂ t W (cid:19) ,W = 2 σ W + F W − ∂ ϕ W + O ( µ ) ,W = 3 σ W + O ( µ ) . Now, starting with substituting the expression for W in the equation for W , we find W = 2 σ W + F W − σ ∂ ϕ W + O ( µ ) ,W = σ W + µ (cid:20) − σ ∂ t W + F (cid:0) F W − σ ∂ ϕ W (cid:1) − σ ∂ ϕ ( F W )+ σ ∂ ϕ W − σ ∂ ϕ (cid:0) F W − σ ∂ ϕ W (cid:1)(cid:21) + O (cid:0) µ (cid:1) ,W = F W − σ ∂ ϕ W + µ (cid:2) − ( ∂ t F + F ∂ ϕ F ) W + σ ( ∂ ϕ F ) ∂ ϕ W (cid:3) + O (cid:0) µ (cid:1) . Finally, in the first order in µ , we obtain the following Fokker–Planck equation for thedistribution density of the phases ρ ( ϕ, t ) ≡ W ( ϕ, t ): ∂ t ρ + ∂ ϕ [( F (1 − µ∂ ϕ F ) − µ∂ t F ) ρ ] = σ ∂ ϕ [(1 − µ∂ ϕ F ) ∂ ϕ ρ ] . The corresponding Langevin equation reads˙ ϕ = F − µ (cid:18) ∂ t + F ∂ ϕ + σ ∂ ϕ (cid:19) F + σ p − µ∂ ϕ F ξ ( t ) . Next, similarly to the procedure described in Section 2, we transform to the rotatingreference frame by introducing θ = ϕ − ψ , where ˙ ψ = Ω. In this reference frame we set F = ν − A sin θ and look for a stationary solution ρ ( θ ), which satisfies the followingequation 0 = − ∂ θ [(1 + µA cos θ ) ( ν − A sin θ ) ρ ] + σ ∂ θ [(1 + µA cos θ ) ∂ θ ρ ] . (A.5)Solution of (A.5) satisfying periodicity ρ ( θ ) = ρ ( θ + 2 π ), reads (cf [24]) ρ ( θ | ω ) = 1 Z e V ( θ ) /σ Z θ +2 πθ d θ ′ (1 − µA cos θ ′ ) e − V ( θ ′ ) /σ + o ( µ ) , (A.6) nalytical approach to synchronous states of globally coupled noisy rotators V ( θ ) = νθ + A cos θ , and Z is a normalization constant, Z = 2 π e − πν/σ (cid:26) I i νσ (cid:18) − Aσ (cid:19) I − i νσ (cid:18) − Aσ (cid:19) − µA (cid:20) I i νσ (cid:18) − Aσ (cid:19) I − i νσ (cid:18) − Aσ (cid:19) + I i νσ (cid:18) − Aσ (cid:19) I − i νσ (cid:18) − Aσ (cid:19)(cid:21)(cid:27) + o ( µ ) . (see formulas 6.681.3 and 8.511.4 in [25]). Similar integrals appear for the orderparameter h e iθ i by virtue of (A.6); the resulting expression coincides with that for a , ( ω )(equation (23)) up to normalization. [1] Kuramoto Y 1975 Self-entrainment of a population of coupled nonlinear oscillators InternationalSymposium on Mathematical Problems in Theoretical Physics ed Araki H (New York: SpringerLecture Notes Phys., v. 39) p 420[2] Acebrn J A, Bonilla L L, Vicente C J P, Ritort F and Spigler R 2005
Rev. Mod. Phys. J. Stat. Mech. - Theor. Exp. R08001[4] Pikovsky A and Rosenblum M 2015
Chaos Contemporary Physics Phys. Rev. Lett. (16) 164101[7] Watanabe S and Strogatz S H 1993
Phys. Rev. Lett. Chaos The European Physical Journal B Chaos Phys. Rev. Lett. Phys. Rev. E (3) 3437–3454[13] Olmi S, Navas A, Boccaletti S and Torcini A 2014 Phys. Rev. E (4) 042905[14] Barr J and Mtivier D 2016 Phys. Rev. Lett.
Statistical Physics of Synchronization (Cham: Springer)[16] Gao J and Efstathiou K 2018
Phys. Rev. E (4) 042201[17] Komarov M, Gupta S and Pikovsky A 2014 EPL
J. Stat. Mech. P05011[19] Risken H Z 1989
The Fokker–Planck Equation (Berlin: Springer)[20] Bonilla L L, Neu J C and Spigler R 1992
J. Stat. Phys. Phys. Rev. E Chaos J. Stat. Phys. Topics in the Theory of Random Noise, Volume 2 (New York: Gordon andBreach)[25] Gradshteyn I and Ryzhik I 2007