Analytical approach to synchronous states of globally coupled noisy rotators

The Kuramoto model of globally coupled phase oscillators [1] is a paradigmatic model to study synchronization phenomena [2–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.

In this paper we consider an ensemble of N globally coupled rotators characterized by their angles φ_n and velocities ${\dot{\varphi }}_{n}$ (n = 1,2, ..., N). The rotators are coupled via the complex mean field

$$\begin{eqnarray}&&R\equiv {{r}{\rm{e}}}^{{\rm{i}}\psi }=\displaystyle \frac{1}{N}\sum _{n=1}^{N}{{\rm{e}}}^{{\rm{i}}{\varphi }_{n}}\end{eqnarray} \tag{ 1 }$$

and obey equations of motion

$$\begin{eqnarray}\begin{array}{rcl}\mu {\ddot{\varphi }}_{n}+{\dot{\varphi }}_{n} & = & {\omega }_{n}+\displaystyle \frac{\varepsilon }{N}\sum _{n^{\prime} =1}^{N}\sin \left({\varphi }_{n^{\prime} }-{\varphi }_{n}\right)+\sigma {\zeta }_{n}\left(t\right)\\ & = & {\omega }_{n}+\varepsilon r\sin \left(\psi -{\varphi }_{n}\right)+\sigma {\zeta }_{n}\left(t\right).\end{array}\end{eqnarray} \tag{ 2 }$$

The unit of time is chosen so that the coefficient at the friction term ($\sim {\dot{\varphi }}_{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 $\sigma {\zeta }_{n}\left(t\right)$ with amplitude σ (which is assumed to be the same for all rotators), zero means $\langle {\zeta }_{n}\left(t\right)\rangle =0$, and auto-correlations $\langle {\zeta }_{n}\left({t}_{1}\right){\zeta }_{n^{\prime} }\left({t}_{2}\right)\rangle =2{\delta }_{{nn}^{\prime} }\delta \left({t}_{1}-{t}_{2}\right)$ (where ${\delta }_{{nn}^{\prime} }$ denotes the Kronecker delta, and $\delta \left(t\right)$ is the Dirac δ-function). Whereas equations (2) are used for numerical simulations below, the analytical approach is developed in the thermodynamics limit $N\to \infty$.

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\to \infty$. Furthermore, we look for stationary synchronous states, i.e. those with constant modulus of the order parameter and with a uniformly rotating angle: $r=\mathrm{const},\dot{\psi }={\rm{\Omega }}$. Such states are possible only in the thermodynamics limit where the finite-size fluctuations vanish. It is convenient to introduce a new angle variable related to the angle of the mean field ${\theta }_{n}={\varphi }_{n}-\psi ,{\dot{\theta }}_{n}={u}_{n}={\dot{\varphi }}_{n}-{\rm{\Omega }}$. Then, model (2) can be formally rewritten as an infinite-dimensional system

$$\begin{eqnarray}&&\mu {\ddot{\theta }}_{n}+{\dot{\theta }}_{n}={\omega }_{n}-{\rm{\Omega }}-\varepsilon r\sin \left({\theta }_{n}\right)+\sigma {\zeta }_{n}\left(t\right),\quad r=\langle {{\rm{e}}}^{{\rm{i}}{\theta }_{n}}\rangle .\end{eqnarray} \tag{ 3 }$$

Stochastic equation (3) allows one to use the Fokker–Planck (Kramers) equation for the probability density $P(\theta ,u,t| \omega )$ for a subset of rotators having natural frequency ω and to express the mean field as the integral over this density,

$$\begin{eqnarray*}&&\begin{array}{l}{\partial }_{t}P+{\partial }_{\theta }\left({uP}\right)+{\mu }^{-1}{\partial }_{u}\left[\left(-u+\omega -{\rm{\Omega }}-\varepsilon r\sin \theta \right)P\right]={\mu }^{-2}{\sigma }^{2}{\partial }_{u}^{2}P,\\ r=\displaystyle \int {\rm{d}}\omega \,g(\omega )\displaystyle \int {\rm{d}}\theta \,\displaystyle \int {\rm{d}}u\,{{\rm{e}}}^{{\rm{i}}\theta }P(\theta ,u,t\,| \,\omega ).\end{array}\end{eqnarray*}$$

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 $\nu =\omega -{\rm{\Omega }}$ and forcing term A = ε r, we seek this stationary distribution as a function depending on these parameters explicitly, ${P}_{0}(\theta ,u\,| \,\nu ,A)$. Then, the system takes the form

$$\begin{eqnarray}&&u{\partial }_{\theta }{P}_{0}={\mu }^{-1}{\partial }_{u}\left(\left(u-\nu +A\sin \theta \right){P}_{0}\right)+{\mu }^{-2}{\sigma }^{2}{\partial }_{u}^{2}{P}_{0},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&r=\int {\rm{d}}\nu \,g\left({\rm{\Omega }}+\nu \right)\int {\rm{d}}\theta \,\int {\rm{d}}u\,{e}^{i\theta }{P}_{0}\left(\theta ,u\,| \,\nu ,A\right),\quad A=\varepsilon r.\end{eqnarray} \tag{ 5 }$$

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, $g({\omega }_{0}+\nu )=g({\omega }_{0}-\nu )$. Then, one may chose Ω = ω₀, which corresponds to the symmetric stationary density ${P}_{0},{P}_{0}(\theta ,u\,| \,\nu ,A)={P}_{0}(-\theta ,-u\,| -\nu ,A)$. In this case, the integral in (5) is real, which proofs that ω₀ is a correct value of the frequency of the mean field Ω. Note that there could also be other appropriate values of Ω besides ω₀. 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,

$$\begin{eqnarray*}&&r(A)=\int {\rm{d}}\nu \,{\rm{d}}\theta \,{\rm{d}}u\,g\left({\omega }_{0}+\nu \right){{\rm{e}}}^{{\rm{i}}\theta }{P}_{0}\left(\theta ,u\,| \,\nu ,A\right),\quad \varepsilon (A)=\displaystyle \frac{A}{r(A)}.\end{eqnarray*}$$

The only remaining act here is finding solution P₀ 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.

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 in appendix.

3.1. Matrix representation of the Fokker–Planck equation

According to [12, 17], we represent a stationary solution ${P}_{0}\left(\theta ,u\,| \,\nu \right)$ of (4) as a double series in the parabolic cylinder functions ${{\rm{\Phi }}}_{p}\left(u\right)$ in u and the Fourier modes ${{\rm{e}}}^{{\rm{i}}{q}\theta }$ in θ,

$$\begin{eqnarray}&&{P}_{0}\left(\theta ,u\,| \,\nu ,A\right)={\left(2\pi \right)}^{-1/2}{{\rm{\Phi }}}_{0}\left(u\right)\sum _{p=0}^{+\infty }\sum _{q=-\infty }^{+\infty }{a}_{{pq}}\left(\nu ,A\right){{\rm{\Phi }}}_{p}\left(u\right){{\rm{e}}}^{{\rm{i}}{q}\theta }.\end{eqnarray} \tag{ 6 }$$

Here functions ${{\rm{\Phi }}}_{p}\left(u\right)$ are defined as

$$\begin{eqnarray*}&&{{\rm{\Phi }}}_{p}\left(u\right)=\sqrt{\displaystyle \frac{\varkappa }{{2}^{p}p!\sqrt{\pi }}}\exp \left[-{\varkappa }^{2}{u}^{2}/2\right]{H}_{p}\left(\varkappa u\right),\end{eqnarray*}$$

where $\varkappa =\sqrt{\mu /2{\sigma }^{2}}$ and ${H}_{p}\left(\varkappa u\right)$ are the Hermite polynomials. Substituting expansion (6) in (4), (5) and using orthogonality conditions for the basis functions, we obtain the infinite system for unknown coefficients ${a}_{{pq}}\left(\nu ,A\right)$ where the order parameter r is just proportional to an integral over one of the expansion coefficients,

$$\begin{eqnarray}&&\sigma \sqrt{\displaystyle \frac{p+1}{\mu }}{{\rm{i}}{q}{a}}_{p+1,q}+\displaystyle \frac{p}{\mu }{a}_{{pq}}-\sigma \sqrt{\displaystyle \frac{p}{\mu }}\left[\left(\displaystyle \frac{\nu }{{\sigma }^{2}}-{\rm{i}}{q}\right){a}_{p-1,q}-{\rm{i}}\displaystyle \frac{A}{2{\sigma }^{2}}\left({a}_{p-1,q+1}-{a}_{p-1,q-1}\right)\right)=0,\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&r=\sqrt{2\pi }\int {\rm{d}}\nu \,g\left({\omega }_{0}+\nu \right){a}_{0,1}^{* }\left(\nu ,A\right).\end{eqnarray} \tag{ 8 }$$

Thus, the main challenge is to find the quantity ${a}_{\mathrm{0,1}}(\nu ,A)$ via solving (7); this coefficient is nothing else but the order parameter for the group of oscillators having natural frequency ${\omega }_{0}+\nu ;$ 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), ${a}_{{pq}}\left(-\nu ,A\right)={\left(-1\right)}^{p}{a}_{{pq}}^{* }\left(\nu ,A\right)$.

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 ${a}_{\mathrm{0,1}}\left(\nu ,A\right)$ is as follows:

$$\begin{eqnarray}&&{a}_{\mathrm{0,1}}\left(\nu ,A\right)=\displaystyle \frac{{S}_{\mathrm{1,0}}\left(\nu ,A\right)}{\sqrt{2\pi }{S}_{\mathrm{0,0}}\left(\nu ,A\right)},\end{eqnarray} \tag{ 9 }$$

where ${S}_{{jk}}\left(\nu ,A\right)$ is the element of the infinite matrix S, which is given by the recurrence formula

$$\begin{eqnarray}&&{\bf{S}}={\tilde{{\bf{D}}}}^{-1}\left({\bf{I}}-\mu {\bf{D}}{\left({\bf{I}}-\displaystyle \frac{\mu }{2}{\bf{D}}{\left({\bf{I}}-\displaystyle \frac{\mu }{3}{\bf{D}}{\left({\bf{I}}-\ldots \right)}^{-1}\tilde{{\bf{D}}}\right)}^{-1}\tilde{{\bf{D}}}\right)}^{-1}\tilde{{\bf{D}}}\right),\end{eqnarray} \tag{ 10 }$$

where I is the identical matrix, and the infinite diagonal matrix D and the infinite tri-diagonal matrix $\tilde{{\bf{D}}}$ are defined as

$$\begin{eqnarray}&&{D}_{{jk}}={\rm{i}}{k}{\delta }_{{jk}},\quad {\tilde{D}}_{{jk}}=\left({\rm{i}}{\sigma }^{2}k-\nu \right){\delta }_{{jk}}+{\rm{i}}\displaystyle \frac{A}{2}\left({\delta }_{j,k-1}-{\delta }_{j,k+1}\right).\end{eqnarray} \tag{ 11 }$$

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 $g\left(\omega \right)$, 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.

3.2. 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

$$\begin{eqnarray}&&{\bf{S}}={\tilde{{\bf{D}}}}^{-1}\left({\bf{I}}-\mu {\bf{D}}\tilde{{\bf{D}}}\right).\end{eqnarray} \tag{ 12 }$$

In order to evaluate the inverse matrix ${\tilde{{\bf{D}}}}^{-1}$, we denote the principal minors as

$$\begin{eqnarray}{M}_{{jk}}=\tilde{{\bf{D}}}\left[\begin{array}{cccc}j & j+1 & \ldots & k\\ j & j+1 & \ldots & k\end{array}\right]=\left|\begin{array}{cccc}{\tilde{D}}_{{jj}} & {\tilde{D}}_{j,j+1} & \ldots & {\tilde{D}}_{{jk}}\\ {\tilde{D}}_{j+1,j} & {\tilde{D}}_{j+1,j+1} & \ldots & {\tilde{D}}_{j+1,k}\\ \ldots & \ldots & \ldots & \ldots \\ {\tilde{D}}_{{kj}} & {\tilde{D}}_{k,j+1} & \ldots & {\tilde{D}}_{{kk}}\end{array}\right|,\end{eqnarray} \tag{ 13 }$$

and introduce a cutoff (cyclic) harmonic number $d\to +\infty$ for Fourier modes. Then, the elements of the inverse matrix ${\tilde{{\bf{D}}}}^{-1}$ take the form

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{D}}_{{jk}}^{-1} & = & \mathop{\mathrm{lim}}\limits_{d\to +\infty }{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{j-k}\displaystyle \frac{{M}_{-d,k-1}{M}_{j+1,d}}{\det \tilde{{\bf{D}}}}{\rm{for}}j\geqslant k,\\ {\tilde{D}}_{{jk}}^{-1} & = & {\left(-1\right)}^{j+k}{\tilde{D}}_{{kj}}^{-1}\qquad \qquad \qquad {\rm{for}}j\lt k.\end{array}\end{eqnarray} \tag{ 14 }$$

According to the Laplace theorem, for all j ≥ k (for convenience, we take ${M}_{k+1,k}=1$ and ${M}_{k+2,k}=0$), the following representation holds

$$\begin{eqnarray}\begin{array}{rcl}\det \tilde{{\bf{D}}} & = & {M}_{{kj}}{M}_{-d,k-1}{M}_{j+1,d}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k,j-1}{M}_{-d,k-1}{M}_{j+2,d}\\ & & +{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k+1,j}{M}_{-d,k-2}{M}_{j+1,d}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{4}{M}_{k+1,j-1}{M}_{-d,k-2}{M}_{j+2,d}.\end{array}\end{eqnarray} \tag{ 15 }$$

Substituting (15) in (14), we obtain for $j\geqslant k$ that

$$\begin{eqnarray}\begin{array}{rcl}{\tilde{D}}_{{jk}}^{-1} & = & \mathop{\mathrm{lim}}\limits_{d\to +\infty }{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{j-k}\left({M}_{{kj}}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k+1,j}\displaystyle \frac{{M}_{-d,k-2}}{{M}_{-d,k-1}}\right.\\ & & {\left.+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k,j-1}\displaystyle \frac{{M}_{j+2,d}}{{M}_{j+1,d}}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{4}{M}_{k+1,j-1}\displaystyle \frac{{M}_{-d,k-2}}{{M}_{-d,k-1}}\displaystyle \frac{{M}_{j+2,d}}{{M}_{j+1,d}}\right)}^{-1}.\end{array}\end{eqnarray} \tag{ 16 }$$

The principal minor M_jk (13) is found from a recurrent equation with fixed j, i.e. with fixed bottom right corner

$$\begin{eqnarray}&&{M}_{{kj}}=\left({\rm{i}}{\sigma }^{2}k-\nu \right){M}_{k+1,j}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k+2,j}.\end{eqnarray} \tag{ 17 }$$

Recurrent equation (17) possesses a solution

$$\begin{eqnarray}&&{M}_{{kj}}={\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{j-k+1}\displaystyle \frac{{I}_{-j-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right){K}_{-k+1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)-{K}_{-j-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right){I}_{-k+1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{-j-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right){K}_{-j-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)-{K}_{-j-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right){I}_{-j-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 18 }$$

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

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{d\to +\infty }\displaystyle \frac{{I}_{d-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{K}_{d-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)}=0.\end{eqnarray} \tag{ 19 }$$

Employing both (18) and (19), we obtain

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{d\to +\infty }\displaystyle \frac{{M}_{-d,k-2}}{{M}_{-d,k-1}}=-{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{-1}\displaystyle \frac{{I}_{-k+1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{-k-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 20 }$$

Similarly, we find a representation of the principal minor M_kj, using its expansion at fixed k, i.e. at fixed upper left corner

$$\begin{eqnarray*}&&{M}_{{kj}}=\left({\rm{i}}{\sigma }^{2}j-\nu \right){M}_{k,j-1}+{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{2}{M}_{k,j-2},\end{eqnarray*}$$

which yields

$$\begin{eqnarray}&&{M}_{{kj}}={\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{j-k+1}\displaystyle \frac{{I}_{k-1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right){K}_{j+1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)-{K}_{k-1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right){I}_{j+1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{k-1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right){K}_{k+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)-{K}_{k-1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right){I}_{k+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 21 }$$

Using expressions (21) and (19), we find

$$\begin{eqnarray}&&\mathop{\mathrm{lim}}\limits_{d\to +\infty }\displaystyle \frac{{M}_{j+2,d}}{{M}_{j+1,d}}=-{\left({\rm{i}}\displaystyle \frac{A}{2}\right)}^{-1}\displaystyle \frac{{I}_{j+1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{j+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 22 }$$

Combination of expressions (11)–(13), (16), (20), and (22) with the general formula (9) results in the closed-form formula to the first order in the mass parameter μ

$$\begin{eqnarray}&&{a}_{\mathrm{0,1}}\left(\nu ,A\right)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\left(1-\mu \displaystyle \frac{{\sigma }^{2}}{\pi }\displaystyle \frac{\sin \left({\rm{i}}\pi \tfrac{\nu }{{\sigma }^{2}}\right)}{{I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right){I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right).\end{eqnarray} \tag{ 23 }$$

An alternative derivation of our main result (23) is based on the method of moments for the density and on elimination of the velocity; it is presented in appendix.

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

4.1. Noise-free case

For purely deterministic rotators, we have to set σ → 0. To find the noise-free limit of general expression (23), it is convenient to rewrite it in an equivalent form

$$\begin{eqnarray}&&{a}_{\mathrm{0,1}}\left(\nu ,A\right)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\left(1-\mu \displaystyle \frac{A}{2}\left(\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}-\displaystyle \frac{{I}_{-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)\right).\end{eqnarray} \tag{ 24 }$$

So, it is neccesary to calculate limits at $\sigma \to 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

$$\begin{eqnarray*}&&{\left(\mathop{\mathrm{lim}}\limits_{\sigma \to 0}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)}^{2}+2{\rm{i}}\displaystyle \frac{\nu }{A}\left(\mathop{\mathrm{lim}}\limits_{\sigma \to 0}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)-1=0.\end{eqnarray*}$$

From this quadratic equation we arrive at

$$\begin{eqnarray*}&&\mathop{\mathrm{lim}}\limits_{\sigma \to 0}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}=-{\rm{i}}\displaystyle \frac{\nu }{A}\pm \sqrt{1-\displaystyle \frac{{\nu }^{2}}{{A}^{2}}}.\end{eqnarray*}$$

The sign is chosen for reasons of continuity with respect to ν finiteness at ν → ±∞ for domains ν > A and ν < −A, and equality of expression to unity at ν = 0 for domain $\left|\nu \right|\leqslant A$. After evaluating the limits, one obtains

$$\begin{eqnarray*}\begin{array}{rcl}\mathop{\mathrm{lim}}\limits_{\sigma \to 0}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)} & = & \left\{\begin{array}{ll}-{\rm{i}}\tfrac{\nu }{A}+\sqrt{1-\tfrac{{\nu }^{2}}{{A}^{2}}}, & \nu \gt -A,\\ -{\rm{i}}\tfrac{\nu }{A}-\sqrt{1-\tfrac{{\nu }^{2}}{{A}^{2}}}, & \nu \leqslant -A,\end{array}\right.\\ \mathop{\mathrm{lim}}\limits_{\sigma \to 0}\displaystyle \frac{{I}_{-1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)} & = & \left\{\begin{array}{ll}-{\rm{i}}\tfrac{\nu }{A}+\sqrt{1-\tfrac{{\nu }^{2}}{{A}^{2}}}, & \nu \leqslant A,\\ -{\rm{i}}\tfrac{\nu }{A}-\sqrt{1-\tfrac{{\nu }^{2}}{{A}^{2}}}, & \nu \gt A.\end{array}\right.\end{array}\end{eqnarray*}$$

The resulting expression for the subgroup order parameter (in the first order in μ) is

$$\begin{eqnarray*}{a}_{\mathrm{0,1}}\left(\nu ,A\right)=\displaystyle \frac{1}{\sqrt{2\pi }A}\left\{\begin{array}{ll}{\rm{i}}\left(-\nu -\sqrt{{\nu }^{2}-{A}^{2}}\right)\left(1+{\rm{i}}\mu \sqrt{{\nu }^{2}-{A}^{2}}\right), & \nu \lt -A,\\ \left(-{\rm{i}}\nu +\sqrt{{A}^{2}-{\nu }^{2}}\right), & \left|\nu \right|\leqslant A,\\ {\rm{i}}\left(-\nu +\sqrt{{\nu }^{2}-{A}^{2}}\right)\left(1-{\rm{i}}\mu \sqrt{{\nu }^{2}-{A}^{2}}\right), & \nu \gt A.\end{array}\right.\end{eqnarray*}$$

4.2. 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]):

$$\begin{eqnarray}&&{a}_{1}\left(\nu ,A\right)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 25 }$$

Correspondingly, the expression for the full order parameter is also exact

$$\begin{eqnarray}&&r=\int {\rm{d}}\nu \,g\left({\omega }_{0}+\nu \right)\displaystyle \frac{{I}_{1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 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 the Lorentz distribution

$$\begin{eqnarray}&&g\left({\omega }_{0}+\nu \right)=\displaystyle \frac{\gamma }{\pi \left({\nu }^{2}+{\gamma }^{2}\right)},\end{eqnarray} \tag{ 27 }$$

the resulting expression for r reads

$$\begin{eqnarray}&&r=\displaystyle \frac{{I}_{1+\tfrac{\gamma }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{\tfrac{\gamma }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}.\end{eqnarray} \tag{ 28 }$$

A more elaborated application of the theory to the Kuramoto model with noise will be presented elsewhere.

4.3. Synchronization transition

Here, we employ the general parametric representation of the order parameter as a function of the coupling constant,

$$\begin{eqnarray}&&r=\sqrt{2\pi }\int {\rm{d}}\nu \,g\left({\omega }_{0}+\nu \right){a}_{0,1}^{* }\left(\nu ,A\right),\quad \varepsilon =\displaystyle \frac{A}{r}\end{eqnarray} \tag{ 29 }$$

with ${a}_{\mathrm{0,1}}\left(\nu ,A\right)$ 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 function of A

$$\begin{eqnarray}&&r\left(A\right)={C}_{0}A+{C}_{1}{A}^{3}+{C}_{2}{A}^{5}+\ldots .\end{eqnarray} \tag{ 30 }$$

To obtain the three initial coefficients in expansion (30), we rewrite (24) using recurrence relations for the modified Bessel functions (see formula (10.29.1) in [21])

$$\begin{eqnarray*}&&{a}_{\mathrm{0,1}}\left(\nu ,A\right)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\left(1-{\rm{i}}\mu \nu -{\rm{i}}\mu A{\mathfrak{I}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)\right).\end{eqnarray*}$$

Then, introducing an auxiliary variable ${\sigma }^{2}y=\nu$ in expression (29) and taking into account that the value r is real, we get

$$\begin{eqnarray*}&&r={\sigma }^{2}\int {\rm{d}}y\,g\left({\omega }_{0}+{\sigma }^{2}y\right)\left({\mathfrak{R}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)+\mu {\mathfrak{I}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)\left({\sigma }^{2}y+A{\mathfrak{I}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)\right)\right).\end{eqnarray*}$$

Next, we use the expansion of the modified Bessel function (see formula (10.25.2) in [21])

$$\begin{eqnarray*}&&{I}_{\kappa }(z)={\left(z/2\right)}^{\kappa }\sum _{\kappa ^{\prime} =0}^{\infty }\displaystyle \frac{{\left(z/2\right)}^{2\kappa ^{\prime} }}{\kappa ^{\prime} !{\rm{\Gamma }}(\kappa +\kappa ^{\prime} +1)},\end{eqnarray*}$$

where Γ(κ) denotes the gamma function. As a result, we arrive at the following relations

$$\begin{eqnarray*}\begin{array}{rcl}{\mathfrak{R}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right) & = & \displaystyle \frac{A}{2{\sigma }^{2}\left(1+{y}^{2}\right)}-\displaystyle \frac{\left(1-2{y}^{2}\right){A}^{3}}{4{\sigma }^{6}{\left(1+{y}^{2}\right)}^{2}\left(4+{y}^{2}\right)}+\displaystyle \frac{\left(3-17{y}^{2}+4{y}^{4}\right){A}^{5}}{8{\sigma }^{10}{\left(1+{y}^{2}\right)}^{3}\left(4+{y}^{2}\right)\left(9+{y}^{2}\right)}+\ldots ,\\ {\mathfrak{I}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}{y}}\left(\tfrac{A}{{\sigma }^{2}}\right)}\right) & = & -\displaystyle \frac{{yA}}{2{\sigma }^{2}\left(1+{y}^{2}\right)}+\displaystyle \frac{\left(5-{y}^{2}\right){{yA}}^{3}}{8{\sigma }^{6}{\left(1+{y}^{2}\right)}^{2}\left(4+{y}^{2}\right)}-\displaystyle \frac{\left(23-24{y}^{2}+{y}^{4}\right){{yA}}^{5}}{16{\sigma }^{10}{\left(1+{y}^{2}\right)}^{3}\left(4+{y}^{2}\right)\left(9+{y}^{2}\right)}+\ldots ,\end{array}\end{eqnarray*}$$

which yield:

$$\begin{eqnarray}\begin{array}{rcl}{C}_{0} & = & \displaystyle \frac{1}{2}\displaystyle \int {\rm{d}}y\,\displaystyle \frac{g\left({\omega }_{0}+{\sigma }^{2}y\right)}{1+{y}^{2}}\left(1-\mu {\sigma }^{2}{y}^{2}\right),\\ {C}_{1} & = & -\displaystyle \frac{1}{8{\sigma }^{4}}\displaystyle \int {\rm{d}}y\,\displaystyle \frac{g\left({\omega }_{0}+{\sigma }^{2}y\right)}{{\left(1+{y}^{2}\right)}^{2}\left(4+{y}^{2}\right)}\left(2\left(1-2{y}^{2}\right)-\mu {\sigma }^{2}{y}^{2}\left(13+{y}^{2}\right)\right),\\ {C}_{2} & = & \displaystyle \frac{1}{16{\sigma }^{8}}\displaystyle \int {\rm{d}}y\,\displaystyle \frac{g\left({\omega }_{0}+{\sigma }^{2}y\right)\left(2\left(3-17{y}^{2}+4{y}^{4}\right)+\mu {\sigma }^{2}{y}^{2}\left(-113+32{y}^{2}+{y}^{4}\right)\right)}{{\left(1+{y}^{2}\right)}^{3}\left(4+{y}^{2}\right)\left(9+{y}^{2}\right)}.\end{array}\end{eqnarray} \tag{ 31 }$$

In the massless case μ = 0, the expressions for C₀, C₁ have been obtained in [23].

The nontrivial branch of solutions $r(\varepsilon )$ starts at

$$\begin{eqnarray*}&&{\varepsilon }_{c}^{\left(1\right)}=1/{C}_{0}.\end{eqnarray*}$$

In the noiseless limit, $\sigma \to 0$, the critical value of the coupling parameter can be expressed as

$$\begin{eqnarray*}&&{\left.{\varepsilon }_{c}^{\left(1\right)}\right|}_{\sigma =0}=\displaystyle \frac{2}{\pi g\left({\omega }_{0}\right)-\mu }.\end{eqnarray*}$$

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

$$\begin{eqnarray}&&1=\displaystyle \frac{\mu \hat{\varepsilon }{{\rm{e}}}^{\mu {\sigma }^{2}}}{2}\sum _{p=0}^{\infty }\displaystyle \frac{{\left(-\mu {\sigma }^{2}\right)}^{2}\left(1+\tfrac{p}{\mu {\sigma }^{2}}\right)}{p!}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\displaystyle \frac{g(\omega )}{\mu {\sigma }^{2}+p+{\rm{i}}(\mu \omega +x)}.\end{eqnarray} \tag{ 32 }$$

For small μ, this expression can be simplified. Assuming $x={x}_{0}+\mu {x}_{1}+{\mu }^{2}{x}_{2}$ and substituting this in (32), we find x₀ = x₂ = 0. For a unimodal frequency distribution symmetric around ω₀, a solution ${x}_{1}=-{\omega }_{0}$ exists, which yields

$$\begin{eqnarray*}&&{\hat{\varepsilon }}^{-1}=\displaystyle \frac{{\sigma }^{2}}{2}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\displaystyle \frac{g(\omega )}{{\sigma }^{4}+{\left(\omega -{\omega }_{0}\right)}^{2}}-\displaystyle \frac{\mu }{2}{\int }_{-\infty }^{\infty }{\rm{d}}\omega \,\displaystyle \frac{g(\omega ){\left(\omega -{\omega }_{0}\right)}^{2}}{{\sigma }^{4}+{\left(\omega -{\omega }_{0}\right)}^{2}}.\end{eqnarray*}$$

One can easily see that $\hat{\varepsilon }={\varepsilon }_{c}^{\left(1\right)}$, which means that the branch of stationary solutions joins the axis r = 0 in the $\left(r,\varepsilon \right)$ plane exactly where the instability of the asynchronous state 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

$$\begin{eqnarray*}&&r={C}_{0}^{2}\sqrt{\left({\varepsilon }_{c}^{\left(1\right)}-\varepsilon \right)/{C}_{1}}\end{eqnarray*}$$

existing for $\varepsilon \gt {\varepsilon }_{c}^{\left(1\right)}$.

Subcritical transition occurs for C₁ > 0. Here, the branch of solutions exists for $\varepsilon \lt {\varepsilon }_{c}^{\left(1\right)}$. If ${C}_{2}\lt 0$, one can estimate that this branch spreads till the minimal value at

$$\begin{eqnarray*}&&{\varepsilon }_{c}^{\left(2\right)}\approx {\left({C}_{0}-\displaystyle \frac{{C}_{1}^{2}}{4{C}_{2}}\right)}^{-1}.\end{eqnarray*}$$

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₁. Because this coefficient depends on the mass μ, there is a critical value of the mass at which the type of the transition changes,

$$\begin{eqnarray}&&{\mu }^{* }=\displaystyle \frac{2}{{\sigma }^{2}}\left(\int {\rm{d}}y\,\displaystyle \frac{g\left({\omega }_{0}+{\sigma }^{2}y\right)\left(1-2{y}^{2}\right)}{{\left(1+{y}^{2}\right)}^{2}\left(4+{y}^{2}\right)}\right){\left(\int {\rm{d}}y\displaystyle \frac{g\left({\omega }_{0}+{\sigma }^{2}y\right){y}^{2}\left(13+{y}^{2}\right)}{{\left(1+{y}^{2}\right)}^{2}\left(4+{y}^{2}\right)}\right)}^{-1}.\end{eqnarray} \tag{ 33 }$$

This formula is obtained from (31) under condition C₁ = 0. The expression is valid only if the value of ${\mu }^{* }$ 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.

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 $\xi =\gamma /{\sigma }^{2}$, ε/γ and μγ.

5.1. Gaussian distribution

The Gaussian distribution of the rotators' frequencies is

$$\begin{eqnarray*}&&g\left(\omega \right)=\displaystyle \frac{1}{\sqrt{2\pi {\gamma }^{2}}}{{\rm{e}}}^{-\tfrac{{\left(\omega -{\omega }_{0}\right)}^{2}}{2{\gamma }^{2}}}.\end{eqnarray*}$$

For this distribution, the coefficients are

$$\begin{eqnarray*}\begin{array}{rcl}{C}_{0} & = & \displaystyle \frac{1}{{\sigma }^{2}}\left(\sqrt{\displaystyle \frac{\pi }{8{\xi }^{2}}}{{\rm{e}}}^{\tfrac{1}{2{\xi }^{2}}}\mathrm{erfc}\left(\displaystyle \frac{1}{\sqrt{2}\xi }\right)\left(1+\mu {\sigma }^{2}\right)-\displaystyle \frac{\mu {\sigma }^{2}}{2}\right),\\ {\mu }^{* } & = & \displaystyle \frac{1}{{\sigma }^{2}}\displaystyle \frac{2\pi {{\rm{e}}}^{\tfrac{1}{2{\xi }^{2}}}\left({\xi }^{2}-1\right)\mathrm{erfc}\left(\tfrac{1}{\sqrt{2}\xi }\right)+2\xi \left(\sqrt{2\pi }\left(1-2{\xi }^{2}\right)+3\pi {{\rm{e}}}^{\tfrac{2}{{\xi }^{2}}}\xi \mathrm{erfc}\left(\tfrac{\sqrt{2}}{\xi }\right)\right)}{{\xi }^{2}\left(\sqrt{2\pi }\xi \left(10+{\xi }^{2}\right)-4\pi \left({{\rm{e}}}^{\tfrac{1}{2{\xi }^{2}}}\mathrm{erfc}\left(\tfrac{1}{\sqrt{2}\xi }\right)+3{{\rm{e}}}^{\tfrac{2}{{\xi }^{2}}}\mathrm{erfc}\left(\tfrac{\sqrt{2}}{\xi }\right)\right)\right)},\end{array}\end{eqnarray*}$$

where $\mathrm{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.

Zoom In Zoom Out Reset image size

5.2. Lorentz distribution

In the case of the Lorentz distribution (27), the expressions for characterization of the synchronization transition are as follows:

$$\begin{eqnarray*}&&{C}_{0}=\displaystyle \frac{1}{2{\sigma }^{2}}\displaystyle \frac{1-\mu {\sigma }^{2}\xi }{1+\xi },\quad {C}_{1}=\displaystyle \frac{1}{8{\sigma }^{6}}\displaystyle \frac{-1+\mu {\sigma }^{2}\xi \left(3+\xi \right)}{{\left(1+\xi \right)}^{2}\left(2+\xi \right)},\quad {\mu }^{* }=\displaystyle \frac{1}{{\sigma }^{2}}\displaystyle \frac{1}{\xi \left(\xi +3\right)}.\end{eqnarray*}$$

Different cases of supercritical and subcritical transitions are illustrated in figure 2. The results are very close to the results shown in figure 1 for the Gaussian distribution. The transition is supercritical for a narrow distributions of natural frequencies and subcritical for larger values of γ.

Zoom In Zoom Out Reset image size

5.3. Bimodal distribution

For a bimodal distribution

$$\begin{eqnarray*}&&g\left(\omega \right)=\displaystyle \frac{1}{2}\delta \left(\omega -{\omega }_{0}-\gamma \right)+\displaystyle \frac{1}{2}\delta \left(\omega -{\omega }_{0}+\gamma \right),\end{eqnarray*}$$

we find

$$\begin{eqnarray*}&&{C}_{0}=\displaystyle \frac{1}{2{\sigma }^{2}}\displaystyle \frac{1-\mu {\sigma }^{2}{\xi }^{2}}{1+{\xi }^{2}},\quad {C}_{1}=\displaystyle \frac{1}{8{\sigma }^{6}}\displaystyle \frac{\left(4+13z\right){\xi }^{2}-2+\mu {\sigma }^{2}{\xi }^{4}}{{\left(1+{\xi }^{2}\right)}^{2}\left(4+{\xi }^{2}\right)},\quad {\mu }^{* }=\displaystyle \frac{2}{{\sigma }^{2}}\displaystyle \frac{1-2{\xi }^{2}}{{\xi }^{2}\left(13+{\xi }^{2}\right)}.\end{eqnarray*}$$

Furthermore, here the integral in (29) can be evaluated in terms of the modified Bessel functions,

$$\begin{eqnarray*}&&r={\mathfrak{R}}\left(\displaystyle \frac{{I}_{1+{\rm{i}}\xi }\left(\tfrac{A}{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\xi }\left(\tfrac{A}{{\sigma }^{2}}\right)}\left(1-\mu \displaystyle \frac{{\sigma }^{2}}{\pi }\displaystyle \frac{\sin \left({\rm{i}}\pi \xi \right)}{{I}_{-{\rm{i}}\xi }\left(\tfrac{A}{{\sigma }^{2}}\right){I}_{{\rm{i}}\xi }\left(\tfrac{A}{{\sigma }^{2}}\right)}\right)\right),\end{eqnarray*}$$

Remarkably, for $\gamma \gt {\sigma }^{2}/\sqrt{2}$ we have ${\mu }^{* }\lt 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.

5.4. Uniform distribution

Here we consider a uniform distribution

$$\begin{eqnarray*}g\left(\omega \right)=\left\{\begin{array}{ll}\tfrac{1}{2\gamma }, & \left|\omega -{\omega }_{0}\right|\leqslant \gamma ,\\ 0, & \left|\omega -{\omega }_{0}\right|\gt \gamma .\end{array}\right.\end{eqnarray*}$$

Using (31) and (33), we obtain

$$\begin{eqnarray}\begin{array}{rcl}{C}_{0} & = & \displaystyle \frac{1}{{\sigma }^{2}}\displaystyle \frac{\left(1+\mu {\sigma }^{2}\right)\arctan \xi -\mu {\sigma }^{2}\xi }{2\xi },\\ {C}_{1} & = & -\displaystyle \frac{1}{8{\sigma }^{6}}\displaystyle \frac{\left(1+2\mu {\sigma }^{2}\right)\xi +\left(1+{\xi }^{2}\right)\left(\left(1+2\mu {\sigma }^{2}\right)\arctan \tfrac{\xi }{2}-\left(1+3\mu {\sigma }^{2}\right)\arctan \xi \right)}{\xi \left(1+{\xi }^{2}\right)},\\ {\mu }^{* } & = & -\displaystyle \frac{1}{{\sigma }^{2}}\displaystyle \frac{\xi +\left({\xi }^{2}+1\right)\left(\arctan \tfrac{\xi }{2}-\arctan \xi \right)}{2\xi +\left({\xi }^{2}+1\right)\left(2\arctan \tfrac{\xi }{2}-3\arctan \xi \right)}.\end{array}\end{eqnarray} \tag{ 34 }$$

In the noise-free case, i.e. at σ = 0, it follows from (4.1) that

$$\begin{eqnarray}r\left(A\right)=\left\{\begin{array}{ll}A\left(\tfrac{\pi }{4\gamma }-\tfrac{\mu }{3}\left(1-\tfrac{A}{\gamma }\right)\tfrac{\gamma -A+2\sqrt{{\gamma }^{2}-{A}^{2}}}{\gamma +\sqrt{{\gamma }^{2}-{A}^{2}}}\right), & A\leqslant \gamma ,\\ \tfrac{1}{2}\left(\sqrt{1-\tfrac{{\gamma }^{2}}{{A}^{2}}}+\tfrac{A}{\gamma }\arcsin \left(\tfrac{\gamma }{A}\right)\right), & A\gt \gamma .\end{array}\right.\end{eqnarray} \tag{ 35 }$$

The critical value of the coupling is, in the first order in μ, ${\varepsilon }_{c}^{\left(1\right)}=4\gamma \left(1+2\mu \gamma /\pi \right)/\pi$.

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 ${\mu }^{* }$ 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.

Zoom In Zoom Out Reset image size

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]):

$$\begin{eqnarray}\begin{array}{rcl}{\mu }_{n}{\ddot{\varphi }}_{n}+{\dot{\varphi }}_{n} & = & {\omega }_{n}+\displaystyle \frac{{\epsilon }_{n}}{N}\sum _{n^{\prime} =1}^{N}\sin \left({\varphi }_{n^{\prime} }-{\varphi }_{n}-{\beta }_{n}\right)+{\sigma }_{n}{\zeta }_{n}\left(t\right)\\ & = & {\omega }_{n}+{\epsilon }_{n}r\sin \left(\psi -{\varphi }_{n}-{\beta }_{n}\right)+{\sigma }_{n}{\zeta }_{n}\left(t\right).\end{array}\end{eqnarray} \tag{ 36 }$$

Here, we introduce two global parameters: E is the average strength of coupling, and B is the characteristic phase shift in coupling

$$\begin{eqnarray*}&&{\epsilon }_{n}=E+{\varepsilon }_{n},\quad {\beta }_{n}=B+{\alpha }_{n}.\end{eqnarray*}$$

The goal is to find uniformly rotating solutions $r=\mathrm{const},\dot{\psi }={\rm{\Omega }}$ in the thermodynamic limit. For the given distribution of individual parameters ${\mu }_{n},{\omega }_{n},{\varepsilon }_{n},{\alpha }_{n},{\sigma }_{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 $E=E(A,{\rm{\Omega }}),B=B(A,{\rm{\Omega }}),r=A/E(A,{\rm{\Omega }})$.

To accomplish this, we introduce the rotating variable ${\theta }_{n}={\varphi }_{n}-\psi -B-{\alpha }_{n}$. Then (36) takes form

$$\begin{eqnarray*}&&\mu {\ddot{\theta }}_{n}+{\dot{\theta }}_{n}={\omega }_{n}-{\rm{\Omega }}-A\varepsilon \sin \left({\theta }_{n}\right)+\sigma {\zeta }_{n}\left(t\right).\end{eqnarray*}$$

The stationary distribution for this Langevin equation was analyzed in section 3 above, it yields the following subgroup order parameter

$$\begin{eqnarray}&&{a}_{\mathrm{0,1}}\left(\omega ,\mu ,\varepsilon ,\sigma \,| \,A,{\rm{\Omega }}\right)=\displaystyle \frac{1}{\sqrt{2\pi }}\displaystyle \frac{{I}_{1+{\rm{i}}\tfrac{\omega -{\rm{\Omega }}}{{\sigma }^{2}}}\left(\tfrac{A\varepsilon }{{\sigma }^{2}}\right)}{{I}_{{\rm{i}}\tfrac{\omega -{\rm{\Omega }}}{{\sigma }^{2}}}\left(\tfrac{A\varepsilon }{{\sigma }^{2}}\right)}\left(1-\mu \displaystyle \frac{{\sigma }^{2}}{\pi }\displaystyle \frac{\sin \left({\rm{i}}\pi \tfrac{\omega -{\rm{\Omega }}}{{\sigma }^{2}}\right)}{{I}_{-{\rm{i}}\tfrac{\omega -{\rm{\Omega }}}{{\sigma }^{2}}}\left(\tfrac{A\varepsilon }{{\sigma }^{2}}\right){I}_{{\rm{i}}\tfrac{\omega -{\rm{\Omega }}}{{\sigma }^{2}}}\left(\tfrac{A\varepsilon }{{\sigma }^{2}}\right)}\right).\end{eqnarray} \tag{ 37 }$$

Substituting this solution in the expression for the global order parameter r, we obtain

$$\begin{eqnarray}&&{{re}}^{-{\rm{i}}{B}}=\int {\rm{d}}\omega \,{\rm{d}}\mu \,{\rm{d}}\varepsilon \,{\rm{d}}\sigma \,{\rm{d}}\alpha \,{a}_{0,1}^{* }\left(\omega ,\mu ,\varepsilon ,\sigma \,| \,A,{\rm{\Omega }}\right){e}^{{\rm{i}}\alpha }G(\omega ,\mu ,\varepsilon ,\sigma ,\alpha )\end{eqnarray} \tag{ 38 }$$

where $G(\omega ,\mu ,\varepsilon ,\sigma ,\alpha )$ is a joint distribution density over the parameters of the problem. Expression (38) yields the representation of the solution in the explicit parametric form

$$\begin{eqnarray*}&&r=r(A,{\rm{\Omega }}),\quad B=B(A,{\rm{\Omega }}),\quad E=\displaystyle \frac{A}{r(A,{\rm{\Omega }})}.\end{eqnarray*}$$

Clearly, if only some of the parameters are distributed, general expressions (37), (38) can be simplified.

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 ν/σ² and forcing A/σ², 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.

We thank D Goldobin and R Toenjes for useful discussions. Results presented in sections 2 and 3 were supported by the RSF grant No. 17-12-01534. Results presented in sections 4–6 were supported by the RSF grant No. 19-12-00367. Results presented in appendix and numerical simulations presented in section 5 were supported by the RFBR grant No. 19-52-12053.

We start with a general reduction of a second-order stochastic equation

$$\begin{eqnarray*}&&\mu \ddot{\varphi }+\dot{\varphi }=F(\varphi ,t)+\sigma \zeta (t)\end{eqnarray*}$$

to a first-order equation, valid for small μ. The probability density $P(\varphi ,\dot{\varphi },t)$ obeys the corresponding Fokker–Planck equation

$$\begin{eqnarray}&&{\partial }_{t}P+{\partial }_{\varphi }\left(\dot{\varphi }P\right)+{\partial }_{\dot{\varphi }}\left({\mu }^{-1}(-\dot{\varphi }+F)P\right)={\mu }^{-2}{\sigma }^{2}{\partial }_{\dot{\varphi }}^{2}P.\end{eqnarray} \tag{ A.1 }$$

We employ the moment method described in [26].

First, we introduce the moments in the velocity:

$$\begin{eqnarray*}&&{w}_{m}\left(\varphi ,t\right)={\int }_{-\infty }^{\infty }{\rm{d}}\dot{\varphi }{\dot{\varphi }}^{m}P\left(\varphi ,\dot{\varphi },t\right).\end{eqnarray*}$$

According to (A.1), these moments obey following equations:

$$\begin{eqnarray}&&{\partial }_{t}{w}_{0}+{\partial }_{\varphi }{w}_{1}=0,\end{eqnarray} \tag{ A.2 }$$

$$\begin{eqnarray}&&{w}_{1}+\mu {\partial }_{t}{w}_{1}={{Fw}}_{0}-\mu {\partial }_{\varphi }{w}_{2},\end{eqnarray} \tag{ A.3 }$$

$$\begin{eqnarray}&&{w}_{m}+\displaystyle \frac{\mu }{m}{\partial }_{t}{w}_{m}={{Fw}}_{m-1}-\displaystyle \frac{\mu }{m}{\partial }_{\varphi }{w}_{m+1}+\left(m-1\right)\displaystyle \frac{{\sigma }^{2}}{\mu }{w}_{m-2}\ {\rm{for}}\ m\geqslant 2.\end{eqnarray} \tag{ A.4 }$$

In order to employ the smallness of parameter μ, it is convenient to introduce rescaled moments

$$\begin{eqnarray*}{w}_{m}=\left\{\begin{array}{ll}\tfrac{1}{{\mu }^{m/2}}{W}_{m} & \ {\rm{for}}\,{\rm{even}}\ m,\\ \tfrac{1}{{\mu }^{\left(m-1\right)/2}}{W}_{m} & {\rm{for}}\,{\rm{odd}}\ m.\end{array}\right.\end{eqnarray*}$$

In terms of moments W_m, equations (A.2)–(A.4) can be rewritten in a form free of singularities

$$\begin{eqnarray*}&&\begin{array}{l}{\partial }_{t}{W}_{0}+{\partial }_{\varphi }{W}_{1}=0,\\ {W}_{1}={{FW}}_{0}-{\partial }_{\varphi }{W}_{2}-\mu {\partial }_{t}{W}_{1},\\ {W}_{m}=\left(m-1\right){\sigma }^{2}{W}_{m-2}+\mu \left({{FW}}_{m-1}-\displaystyle \frac{1}{m}{\partial }_{\varphi }{W}_{m+1}-\displaystyle \frac{1}{m}{\partial }_{t}{W}_{m}\right)\ \ {\rm{for}}\,{\rm{even}}\ m\gt 0,\\ {W}_{m}=\left(m-1\right){\sigma }^{2}{W}_{m-2}+{{FW}}_{m-1}-\displaystyle \frac{1}{m}{\partial }_{\varphi }{W}_{m+1}-\displaystyle \frac{\mu }{m}{\partial }_{t}{W}_{m}\qquad \quad {\rm{for}}\,{\rm{odd}}\ m\gt 1.\end{array}\end{eqnarray*}$$

As we are looking for the first order corrections in μ, we rewrite this system keeping the relevant terms only

$$\begin{eqnarray*}&&\begin{array}{l}{\partial }_{t}{W}_{0}+{\partial }_{\varphi }{W}_{1}=0,\\ {W}_{1}={{FW}}_{0}-{\partial }_{\varphi }{W}_{2}-\mu {\partial }_{t}{W}_{1},\\ {W}_{2}={\sigma }^{2}{W}_{0}+\mu \left({{FW}}_{1}-\displaystyle \frac{1}{2}{\partial }_{\varphi }{W}_{3}-\displaystyle \frac{1}{2}{\partial }_{t}{W}_{2}\right),\\ {W}_{3}\approx 2{\sigma }^{2}{W}_{1}+{{FW}}_{2}-\displaystyle \frac{1}{3}{\partial }_{\varphi }{W}_{4},\\ {W}_{4}\approx 3{\sigma }^{2}{W}_{2}.\end{array}\end{eqnarray*}$$

Now, starting with substituting the expression for W₄ in the equation for W₃, we find

$$\begin{eqnarray*}\begin{array}{rcl}{W}_{3} & \approx & 2{\sigma }^{2}{W}_{1}+{{FW}}_{2}-{\sigma }^{2}{\partial }_{\varphi }{W}_{2}\\ {W}_{2} & \approx & {\sigma }^{2}{W}_{0}+\mu \left[-\displaystyle \frac{{\sigma }^{2}}{2}{\partial }_{t}{W}_{0}+F\left({{FW}}_{0}-{\sigma }^{2}{\partial }_{\varphi }{W}_{0}\right)-\displaystyle \frac{{\sigma }^{2}}{2}{\partial }_{\varphi }\left({{FW}}_{0}\right)\right.\\ & & \left.+\displaystyle \frac{{\sigma }^{4}}{2}{\partial }_{\varphi }^{2}{W}_{0}-{\sigma }^{2}{\partial }_{\varphi }\left({{FW}}_{0}-{\sigma }^{2}{\partial }_{\varphi }{W}_{0}\right)\right],\\ {W}_{1} & \approx & {{FW}}_{0}-{\sigma }^{2}{\partial }_{\varphi }{W}_{0}+\mu \left[-\left({\partial }_{t}F+F{\partial }_{\varphi }F\right){W}_{0}+{\sigma }^{2}\left({\partial }_{\varphi }F\right){\partial }_{\varphi }{W}_{0}\right].\end{array}\end{eqnarray*}$$

Finally, in the first order in μ, we obtain the following Fokker–Planck equation for the distribution density of the phases $\rho (\varphi ,t)\equiv {W}_{0}(\varphi ,t)$:

$$\begin{eqnarray*}&&{\partial }_{t}\rho +{\partial }_{\varphi }\left[\left(F\left(1-\mu {\partial }_{\varphi }F\right)-\mu {\partial }_{t}F\right)\rho \right]={\sigma }^{2}{\partial }_{\varphi }\left[\left(1-\mu {\partial }_{\varphi }F\right){\partial }_{\varphi }\rho \right].\end{eqnarray*}$$

The corresponding Langevin equation reads

$$\begin{eqnarray*}&&\dot{\varphi }=F-\mu \left({\partial }_{t}+F{\partial }_{\varphi }+\displaystyle \frac{{\sigma }^{2}}{2}{\partial }_{\varphi }^{2}\right)F+\sigma \sqrt{1-\mu {\partial }_{\varphi }F}\,\zeta \left(t\right).\end{eqnarray*}$$

Next, similarly to the procedure described in section 2, we transform to the rotating reference frame by introducing $\theta =\varphi -\psi$, where $\dot{\psi }={\rm{\Omega }}$. In this reference frame we set $F=\nu -A\sin \theta$ and look for a stationary solution ${\rho }_{0}(\theta )$, which satisfies the following equation

$$\begin{eqnarray}&&0=-{\partial }_{\theta }\left[\left(1+\mu A\cos \theta \right)\left(\nu -A\sin \theta \right){\rho }_{0}\right]+{\sigma }^{2}{\partial }_{\theta }\left[\left(1+\mu A\cos \theta \right){\partial }_{\theta }{\rho }_{0}\right].\end{eqnarray} \tag{ A.5 }$$

Solution of (A.5) satisfying periodicity ${\rho }_{0}\left(\theta \right)={\rho }_{0}\left(\theta +2\pi \right)$, reads (see [27])

$$\begin{eqnarray}&&{\rho }_{0}\left(\theta \,| \,\omega \right)=\displaystyle \frac{1}{Z}{{\rm{e}}}^{V\left(\theta \right)/{\sigma }^{2}}{\int }_{\theta }^{\theta +2\pi }{\rm{d}}\theta ^{\prime} \left(1-\mu A\cos \theta ^{\prime} \right){{\rm{e}}}^{-V\left(\theta ^{\prime} \right)/{\sigma }^{2}},\end{eqnarray} \tag{ A.6 }$$

where $V\left(\theta \right)=\nu \theta +A\cos \theta$, and Z is a normalization constant,

$$\begin{eqnarray*}\begin{array}{rcl}Z & = & 2{\pi }^{2}{{\rm{e}}}^{-\pi \nu /{\sigma }^{2}}\left\{2{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right){I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right)\right.\\ & & \left.-\mu A\left[{I}_{1+{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right){I}_{-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right)+{I}_{{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right){I}_{1-{\rm{i}}\tfrac{\nu }{{\sigma }^{2}}}\left(-\displaystyle \frac{A}{{\sigma }^{2}}\right)\right]\right\}\end{array}\end{eqnarray*}$$

(see formulas (6.681.3) and (8.511.4) in [28]). Similar integrals appear for the order parameter $\langle {{\rm{e}}}^{{\rm{i}}\theta }\rangle$ by virtue of (A.6); the resulting expression coincides with that for ${a}_{\mathrm{0,1}}\left(\omega \right)$ (equation (23)) up to normalization.