text:
# Synchronization in populations of globally coupled oscillators with inertial effects

## I Introduction

The dynamical behavior of large populations of nonlinearly coupled oscillators may describe many phenomena in Physics, Biology and Medicine, . In particular synchronization of mean-field coupled phase oscillators with different natural frequencies is nicely illustrated by Kuramoto’s well-known and extesively analyzed model . To describe certain biological phenomena, inertial effects should be added to this model. In , Ermentrout revisited the special problem of self-synchronization in populations of fireflies of a certain kind (the Pteroptyx malaccae). Compared to observed behavior, the approach to oscillator synchronization as described by the Kuramoto model seems to be too fast. Thus a more appropriate adaptive frequency model has been proposed in , where the natural frequency of an oscillator is a new independent variable, which is allowed to vary with time. Thus an oscillator is described by its phase and frequency. From the mathematical standpoint, the new model is governed by a system of coupled second-order differential equations containing inertial terms, in contrast to the system of first-order differential equations governing the Kuramoto model. Indeed inertia slows down synchronization and this may result in better agreement between theory and experimental measurements. Other possible biological applications of Ermentrout-type models include after-effects in alterations of circadian cycles in mammalians, cf.

A different set of applications for oscillators with inertia include power systems described by the swing equations , or by Hamilton equations . An important technologically relevant application is the study of superconducting Josephson junctions arrays . Here inertial terms describe the effect of nonzero electrical capacitance. Such effect is often far from being negligible, and it is absent in the Kuramoto model, .

In this paper, we consider the model equations of Ref. ,

$`\dot{\theta _j}`$ $`=`$ $`\omega _j`$ (1)
$`m\dot{\omega _j}`$ $`=`$ $`\omega _j+\mathrm{\Omega }_j+Kr_N\mathrm{sin}(\psi _N\theta _j)+\xi _j(t),`$ (3)
$`j=1,\mathrm{},N,`$
where $`\theta _j`$, $`\omega _j`$ and $`\mathrm{\Omega }_j`$ denote phase, frequency and natural frequency of the $`j`$th oscillator, respectively. The natural frequencies are distributed with probability density $`g(\mathrm{\Omega })`$, which may have a single maximum (unimodal distribution), or several peaks (multimodal distribution). The positive parameters $`m`$ and $`K`$ are the “inertia” and the coupling strength, respectively. The complex order parameter defined by

$$r_Ne^{i\psi _N}=\frac{1}{N}\underset{j=1}{\overset{N}{}}e^{i\theta _j},$$
(4)
measures phase synchronization: $`r_N>0`$ if the oscillators are synchronized and $`r_N=0`$ if not. Finally, $`\xi _j`$’s are independent identically distributed Gaussian white noises, with $`\xi _j=0,\xi _i(t)\xi _j(s)=2D\delta _{ij}\delta (ts)`$. White noise terms were not included in . When the inertial terms vanish, $`m=0`$, Eqs. (3) and (4) are exactly the Kuramoto model.

Typically, N is very large, and oscillator synchronization is conveniently analyzed in the limiting case of infinitely many oscillators. In this limit, models with mean-field coupling are described by an evolution equation for the one-oscillator probability density, $`\rho (\theta ,\omega ,\mathrm{\Omega },t)`$, . For the present model this equation is

$`{\displaystyle \frac{\rho }{t}}`$ $`=`$ $`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{^2\rho }{\omega ^2}}`$ (5)
$``$ $`{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{\omega }}[\left(\omega +\mathrm{\Omega }+Kr\mathrm{sin}(\psi \theta )\right)\rho ]\omega {\displaystyle \frac{\rho }{\theta }},`$ (6)
where the order parameter is now given by

$`re^{i\psi }={\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{i\theta }\rho (\theta ,\omega ,\mathrm{\Omega },t)g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\omega 𝑑\theta .`$ (7)
Equations (6) and (7) should be supplemented with appropriate initial and boundary data ($`\rho `$ is $`2\pi `$-periodic in $`\theta `$ and has suitable decay behavior as $`\omega \pm \mathrm{}`$) plus the normalization condition,

$$_0^{2\pi }_{\mathrm{}}^+\mathrm{}\rho (\theta ,\omega ,\mathrm{\Omega },t)𝑑\omega 𝑑\theta =1.$$
(8)
Differentiating $`_0^{2\pi }_{\mathrm{}}^+\mathrm{}\rho (\theta ,\omega ,\mathrm{\Omega },t)𝑑\omega 𝑑\theta `$ with respect to time, and then using Eq. (6) itself, together with periodicity in $`\theta `$ and decay in $`\omega `$, we find that the left side of (8) is time independent. Normalization to unity of the initial probability density then implies (8) for the solution of (6).

In this paper, we study oscillator synchronization and transition from incoherence to synchronization in the model (6) - (8). The incoherent solution of (6) - (8) (or simply incoherence) is a stationary solution which is independent of $`\theta `$. This solution asigns equal probability to all angles and has $`r=0`$ (no order), so it corresponds to lack of oscillator synchronization. There are synchronized solutions which branch off from incoherence as the coupling among oscillators is increased. These bifurcations describe the synchronization transitions, which we have analyzed and compared to the corresponding ones in the Kuramoto model. Our main results are that inertia: (i) may stabilize incoherence, making it harder to synchronize oscillators, and (ii) it may harden the synchronization transition. In the Kuramoto model ($`m=0`$) or with oscillators with identical natural frequencies, the synchronization transition is soft (supercritical bifurcation), whereas it may become hard (subcritical bifurcation) if the distribution of natural frequencies has a nonzero spread (unimodal Lorentzian distribution) or several peaks (e.g., a discrete bimodal distribution). The methods we have used in our analysis are similar to those previously employed in the Kuramoto model : linear stability of incoherence, bifurcation analysis, high-frequency singular perturbations and numerical solutions. An important difference is that now we do not have an explicit functional form for stationary solutions (as it was the case for the Kuramoto model). This has led us to use mode-coupling expansions of the solution and solving the corresponding mode-coupling equations. Solutions of these equations in close form are not always accesible, so that we have introduced some closure assumptions. The results of these uncontrolled assumptions have been compared to direct simulations or to approximate amplitude equations and found reasonable in the limit of small inertia.

The rest of the paper is as follows. In Section II, we find the incoherent solution and study its linear stability for several natural frequency distributions. Results are compared with those obtained in the massless case . It is found that the critical coupling needed to destabilize incoherence increases with $`m`$ for “unimodal” frequency distributions of the Lorentzian type. The critical coupling is independent of $`m`$ when $`g(\mathrm{\Omega })=\delta (\mathrm{\Omega })`$. In this case the time needed to reach synchronization increases as $`m`$ increases. If $`g(\mathrm{\Omega })=[\delta (\mathrm{\Omega }\mathrm{\Omega }_0)+\delta (\mathrm{\Omega }+\mathrm{\Omega }_0)]/2`$ (discrete bimodal distribution), the critical coupling may grow or decrease with $`m`$ depending on the values of $`\mathrm{\Omega }_0`$. In Section III, we construct other stationary solutions by two procedures: an amplitude expansion for solutions branching off from incoherence and a general expansion in Hermite polynomials which is appropriately truncated. An exact analytical solution is obtained if $`g(\mathrm{\Omega })=\delta (\mathrm{\Omega })`$ (cf. ), while analytical approximations for small $`m`$ and $`\mathrm{\Omega }`$ are available in the general case. We have observed that inertia tends to harden the synchronization transition: in the Kuramoto model ($`m=0`$) or with oscillators with identical natural frequencies, the synchronization transition is soft (supercritical bifurcation), whereas it becomes hard (subcritical bifurcation) in the cases of unimodal Lorentzian or discrete bimodal frequency distributions. In Section IV, we obtain approximations to stable time-dependent solutions of Eq. (6) in the “high frequency limit”, $`\mathrm{\Omega }\mathrm{}`$ . There are partially synchronized nonlinearly stable solutions of standing wave type, as in the Kuramoto model ). Finally, numerical results are presented in Section V, and compared to the approximate or exact solutions of previous Sections. Two Appendices at the end are devoted to technical details.

## II Linear stability of the incoherent solution

The incoherent solution is a $`\theta `$-independent stationary solution of (6). Its order parameter is $`r=0`$ according to (7). Then Eqs. (6), decay as $`\omega \pm \mathrm{}`$ and the normalization condition (8) yield the incoherent solution:

$$\rho _0(\omega ,\mathrm{\Omega })=\frac{1}{2\pi }\sqrt{\frac{m}{2\pi D}}e^{\frac{m}{2D}(\omega \mathrm{\Omega })^2}.$$
(9)
To analyze its linear stability, let us consider a small disturbance about incoherence,

$$\rho (\theta ,\omega ,\mathrm{\Omega },t)=\rho _0(\omega ,\mathrm{\Omega })+\epsilon \eta (\theta ,\omega ,\mathrm{\Omega },t)+O(\epsilon ^2),$$
(10)
where $`\epsilon 1`$. Normalization of $`\rho (\theta ,\omega ,\mathrm{\Omega },t)`$ then implies

$$_0^{2\pi }_{\mathrm{}}^+\mathrm{}\eta (\theta ,\omega ,\mathrm{\Omega },t)𝑑\omega 𝑑\theta =0.$$
(11)
We now introduce (10) into (6) and (7) and equate like terms in $`\epsilon `$. To order $`\epsilon `$, the result is

$`{\displaystyle \frac{\eta }{t}}+\omega {\displaystyle \frac{\eta }{\theta }}{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{\omega }}[(\omega \mathrm{\Omega })\eta ]{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{^2\eta }{\omega ^2}}=`$ (12)
$`{\displaystyle \frac{K\frac{\rho _0}{\omega }}{m}}{\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle _{\mathrm{}}^+\mathrm{}}\eta (\varphi ,\omega ,\mathrm{\Omega },t)`$ (13)
$`\times \mathrm{sin}(\varphi \theta )g(\mathrm{\Omega })d\mathrm{\Omega }d\omega d\varphi .`$ (14)
We now insert a trial solution

$$\eta (\theta ,\omega ,\mathrm{\Omega },t)=e^{\lambda t}\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}b_n(\omega ;\mathrm{\Omega },\lambda )e^{in\theta }$$
(15)
(which is $`2\pi `$-periodic in $`\theta `$) into (14), thereby obtaining

$`{\displaystyle \frac{d^2b_n}{d\omega ^2}}+{\displaystyle \frac{m(\omega \mathrm{\Omega })}{D}}{\displaystyle \frac{db_n}{d\omega }}+{\displaystyle \frac{m(1m\lambda inm\omega )b_n}{D}}=`$ (16)
$`{\displaystyle \frac{\pi mK(i\delta _{n,1}i\delta _{n,1})\frac{\rho _0}{\omega }}{D}}1,b_n,`$ (17)
where we have defined the scalar product

$`\phi ,\psi ={\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle _{\mathrm{}}^+\mathrm{}}\overline{\phi (\omega ,\mathrm{\Omega })}\psi (\omega ,\mathrm{\Omega })g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\omega .`$ (18)
Notice that $`b_n=\overline{b_n}`$ and that $`b_n=0`$ because of the normalization condition (11).

Equation (17) can be transformed into a nonhomogeneous parabolic cylinder equation by the following change of variable:

$`b_n(\omega ;\mathrm{\Omega },\lambda )=\mathrm{exp}\left[{\displaystyle \frac{m(\omega \mathrm{\Omega })^2}{4D}}\right]\beta _n(w;\mathrm{\Omega },\lambda ),`$ (19)
$`w=\sqrt{{\displaystyle \frac{m}{D}}}(\omega \mathrm{\Omega }+2nDi).`$ (20)
Inserting (19) and (20) into (17), we obtain

$`{\displaystyle \frac{d^2\beta _n}{dw^2}}+\left[{\displaystyle \frac{1}{2}}{\displaystyle \frac{w^2}{4}}m(\lambda +in\mathrm{\Omega }+n^2D)\right]\beta _n=`$ (21)
$`i\pi K{\displaystyle \frac{\rho _0}{\omega }}e^{\frac{1}{4}(w2i\sqrt{mD})^2}1,e^{\frac{1}{4}(w2i\sqrt{mD})^2}\beta _1\delta _{n,1}.`$ (22)
(Recall that $`d\omega =\sqrt{D/m}dw`$ when using the definition of scalar product). Let us assume now that $`n\pm 1`$ and that $`\mathrm{\Omega }`$ is a fixed real number. Then the right hand side of (22) is zero and the resulting equation has the following eigenvalues

$`\lambda _{p,n}(\mathrm{\Omega })={\displaystyle \frac{p}{m}}n^2Din\mathrm{\Omega },p=0,1,2,\mathrm{},`$ (23)
associated to the eigenfunctions

$`\beta _{p,n}(w;\mathrm{\Omega },\lambda _{p,n})=D_p(w)=2^{\frac{p}{2}}e^{\frac{w^2}{4}}H_p\left({\displaystyle \frac{w}{\sqrt{2}}}\right),`$ (24)
which are independent of $`n`$ and $`\mathrm{\Omega }`$. In this formula, $`D_p(w)`$ and $`H_p(x)`$ are the parabolic cylinder function and the Hermite polynomial of index $`p`$, respectively . The eigenvalues $`\lambda _{p,n}(\mathrm{\Omega })`$ of (23), with $`n=\pm 1,\pm 2,\mathrm{}`$, $`p=0,1,\mathrm{}`$ and $`\mathrm{\Omega }`$ belonging to the support of $`g(\mathrm{\Omega })`$, constitute the continuous spectrum of the linear stability problem. In fact, a nonhomogeneous linear problem with a homogeneous part given by (22) cannot be solved for an arbitrary source term if $`\lambda =\lambda _{p,n}`$. Notice that the continuous spectrum lies to the left side of the imaginary axis if $`D>0`$ and $`n0`$. Then the “eigenvalues” (23) have negative real parts (and therefore correspond to stable modes). As we have already said, the neutrally stable modes with $`n=0`$ have zero amplitude due to the normalization condition (11).

If $`n=1`$, we can solve (17) by means of an expansion in eigenfunctions $`D_p(w)`$, $`p=0,1,2,\mathrm{}`$. To obtain the generalized Fourier coefficients of $`\beta _1`$, we multiply both sides of (17) by $`D_p(w)`$ and integrate over $`w`$. As $`_{\mathrm{}}^{\mathrm{}}D_p(x)D_n(x)𝑑x=\sqrt{2\pi }p!\delta _{pn}`$ (orthogonality condition, cf. §7.711.1 of ), the result is

$`\beta _1(\omega ;\mathrm{\Omega },\lambda )`$ $`=`$ $`{\displaystyle \frac{i\pi K}{m}}1,e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}\beta _1`$ (25)
$`\times `$ $`{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p\rho _0^{}𝑑w}{\sqrt{2\pi }p!\left(\frac{p}{m}+\lambda +i\mathrm{\Omega }+D\right)}}D_p(w),`$ (26)
where

$`\rho _0^{}(w)`$ $`=`$ $`{\displaystyle \frac{\rho _0}{\omega }}|_{\omega =\mathrm{\Omega }i2D+(D/m)^{\frac{1}{2}}w}`$ (27)
$`=`$ $`{\displaystyle \frac{m(wi2\sqrt{mD})}{(2\pi )^{\frac{3}{2}}D}}e^{\frac{1}{2}(wi2\sqrt{mD})^2}.`$ (28)
Once we have found $`\beta _1`$, we can calculate the scalar product $`1,e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}\beta _1`$. Since this scalar product appears as a factor in both sides of the resulting expression, we can divide by it, thereby obtaining an eigenvalue equation for $`\lambda `$:

$`1={\displaystyle \frac{i\pi K\sqrt{D}}{\sqrt{2\pi m^3}}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{p!}}{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p\rho _0^{}𝑑w`$ (29)
$`\times {\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_pdw`$ (30)
$`\times {\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{g(\mathrm{\Omega })}{\frac{p}{m}+\lambda +i\mathrm{\Omega }+D}}d\mathrm{\Omega }.`$ (31)
In Appendix A, we show that this equation may be rewritten as

$`1={\displaystyle \frac{K}{4\pi \sqrt{mD}}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{𝒜_p(\sqrt{mD})𝒜_p^{}(\sqrt{mD})}{p!}}`$ (32)
$`\times {\displaystyle _{\mathrm{}}^+\mathrm{}}{\displaystyle \frac{g(\mathrm{\Omega })}{\frac{p}{m}+\lambda +i\mathrm{\Omega }+D}}d\mathrm{\Omega },`$ (33)
where $`𝒜_p(x)`$ is defined as

$`𝒜_p(x)={\displaystyle _{\mathrm{}}^+\mathrm{}}D_p(w)e^{\left(\frac{w}{2}ix\right)^2}𝑑w.`$ (34)
The result of evaluating this integral is (cf. Appendix A):

$`𝒜_p(x)=i^p\sqrt{2\pi }e^{\frac{x^2}{2}}x^p(x>0).`$ (35)
Inserting (35) in (33), we obtain

$`1={\displaystyle \frac{Ke^{mD}}{2}}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^p\left(1+\frac{p}{mD}\right)}{p!}}`$ (36)
$`\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{g(\mathrm{\Omega })d\mathrm{\Omega }}{\lambda +D+i\mathrm{\Omega }+\frac{p}{m}}}.`$ (37)
As $`m0`$, this equation coincides with that obtained for the Kuramoto model, . Equation (37) can be rewritten in terms of incomplete gamma functions as follows :

$`{\displaystyle \frac{2D}{K}}=e^{mD}{\displaystyle _{\mathrm{}}^{\mathrm{}}}[mD\gamma (m(\lambda +D+i\mathrm{\Omega }),mD)`$ (38)
$`\gamma (1+(\lambda +D+i\mathrm{\Omega })m,mD)]{\displaystyle \frac{g(\mathrm{\Omega })d\mathrm{\Omega }}{(mD)^{m(\lambda +D+i\mathrm{\Omega })}}}`$ (39)
$`={\displaystyle _{\mathrm{}}^{\mathrm{}}}\left[1{\displaystyle \frac{(\lambda +i\mathrm{\Omega })me^{mD}\gamma (m(\lambda +D+i\mathrm{\Omega }),mD)}{(mD)^{(\lambda +D+i\mathrm{\Omega })m}}}\right]`$ (40)
$`\times g(\mathrm{\Omega })d\mathrm{\Omega },`$ (41)
where

$`\gamma (a,x)={\displaystyle _0^x}e^tt^{a1}𝑑t.`$ (42)
¿From now on, we analyze Eq. (37) for special frequency distributions:

(a) Unimodal frequency distribution, $`g(\mathrm{\Omega })=\delta (\mathrm{\Omega })`$.

In this case we show that if $`\text{Re}\lambda =0`$, then $`\text{Im}\lambda =0`$. Thus, the eigenvalues that may acquire a positive real part are real. Then the critical coupling is obtained by setting $`\lambda =0`$. By subtracting from Eq. (37) its complex conjugate, we obtain

$`0=\text{Im}(\lambda )f(\text{Im}\lambda ,m,D),`$ (43)
where

$$f(\text{Im}\lambda ,m,D)=\underset{p=0}{\overset{\mathrm{}}{}}\frac{(mD)^p\left(1+\frac{p}{mD}\right)}{p!\left[(\text{Im}\lambda )^2+\left(D+\frac{p}{m}\right)^2\right]}.$$
(44)
Notice than the even function $`f(\text{Im}\lambda ,m,D)`$ decreases monotonically with Im$`(\lambda )>0`$. On the other hand, $`f`$ tends to zero, as $`\text{Im}\lambda +\mathrm{}`$, and

$$f\frac{m}{D}(mD)^{mD}\gamma (mD,mD)>0,\text{as }\text{Im}\lambda 0.$$
(45)
Thus $`f`$ does not vanish at finite values of Im$`\lambda `$, and therefore the only solution of (43) is Im$`\lambda =0`$. Setting now $`\lambda =0`$, Equation (37) yields the critical coupling $`K=K_c`$,

$`{\displaystyle \frac{K_c}{2D}}e^{mD}\left[{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^n}{n!}}\right]=1,\text{and therefore}`$ (46)
$`K_c=2D.`$ (47)
Figures 1 and 2 show the largest eigenvalue $`\lambda `$ as a function of $`m`$ and $`K`$. To compute $`\lambda `$ numerically, we used from Eq. (10), and (15), that the amplitude of the order parameter is

$$rCe^{\lambda t},$$
(48)
close to incoherence. Then, the goal is to simulate the evolution of the system, choosing the initial condition sufficiently close to the incoherent solution, and obtain numerically the amplitude order parameter $`r(t)`$. Fig. 2 shows that different eigenvalue curves (for different $`m`$) intersect the horizontal axis, $`\lambda =0`$, at the same value of $`K`$, as expected from Eq. (47).

(b) Unimodal Lorentzian frequency distribution, $`g(\mathrm{\Omega })=\frac{\epsilon /\pi }{\epsilon ^2+\mathrm{\Omega }^2}`$.

Eq. (37) becomes

$`1={\displaystyle \frac{K}{2D}}e^{mD}[{\displaystyle \frac{D}{\lambda +D+\epsilon }}+{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^n}{n!}}`$ (49)
$`\times {\displaystyle \frac{n+mD}{m(\lambda +D+\epsilon )+n}}]={\displaystyle \frac{K}{2}}me^{mD}`$ (50)
$`\times (mD)^{m(\lambda +D+\epsilon )}[(mD)^{m(\lambda +D+\epsilon )1}e^{mD}`$ (51)
$`\gamma (m(\lambda +D+\epsilon ),mD)({\displaystyle \frac{\lambda +D+\epsilon }{D}}1)]`$ (52)
An explicit solution for $`\lambda `$ cannot in general be found. Thus, we consider several limiting cases corresponding to physically interesting parameter choices. In the small noise limit, $`D1`$, we consider the cases $`(i)`$ $`m=O(1)`$ fixed, and $`(ii)`$ $`mD=1`$. It is remarkable that the expansion

$$\gamma (a,x)=e^xx^a\underset{n=0}{\overset{\mathrm{}}{}}\frac{x^n}{(a)_{n+1}},$$
(53)
where $`(a)_k=a(a+1)\mathrm{}(a+k1)`$, $`k=1,2,\mathrm{}`$, holds in both cases. If $`x=mD0`$, $`a=mD+m(\lambda +\epsilon )>x`$, (53) holds as a convergent expansion, . If $`x=mD=1`$, and $`a=1+m(\lambda +\epsilon )\mathrm{}`$ (with fixed $`\lambda `$ and $`\epsilon `$ of order 1), (53) holds as an asymptotic expansion, . Inserting (53) in Eq. (52), we obtain

$`{\displaystyle \frac{2D}{K}}=1+{\displaystyle \frac{xa}{a}}\left[1+{\displaystyle \frac{x}{a+1}}+{\displaystyle \frac{x^2}{(a+1)(a+2)}}+O\left({\displaystyle \frac{x^3}{a^3}}\right)\right]`$ (54)
Similarly to the unimodal case, it is possible to prove that $`\lambda `$ is always real. To this purpose, notice that replacing $`\lambda +\epsilon `$ in eq. (52) with $`\lambda `$ in eq. (37) (setting $`\mathrm{\Omega }=0`$) we obtain the same equation. The critical coupling $`K=K_c`$ is then found by setting $`\lambda =0`$ in (54). In case $`(i)`$, we have

$$K_c=2\epsilon (m\epsilon +1)+\frac{2(2+3m\epsilon )}{2+m\epsilon }D+O\left(D^2\right),$$
(55)
In the limit of vanishing mass, we recover the result $`K_c=2(D+\epsilon )`$ valid for the Kuramoto model . Another important limit is $`\epsilon =0`$, which reproduces the unimodal distribution. We find $`K_c=2D`$, independently of mass. Thus the spread in frequency distribution plays an important role in synchronizing populations of oscillators affected by inertia.

In case $`(ii)`$, Eq. (54) yields

$$1=\frac{Kx}{2Da(a+1)}+O\left(a^3\right),$$
(56)
from which

$$\lambda =\left[\epsilon +\frac{1}{2m}\left(3\pm \sqrt{1+2Km}\right)\right].$$
(57)
This quantity is always real and vanishes for $`K=K_c=2\epsilon (m\epsilon +3)+\frac{4}{m}`$. Note that $`K_c`$ grows roughly linearly with $`m`$. Thus oscillator synchronization is made harder by increasing inertia in the limit of vanishing noise. This behavior is slightly different from that described in . There numerical simulations seemed to show that incoherence remains stable up to a critical coupling, which was independent of $`m`$. The singular nature of the limit $`D0`$ makes the cause of this discrepancy unclear, although we should mention that no stability analysis was conducted in . In the opposite limit $`m\mathrm{}`$, $`\lambda \epsilon `$, and incoherence is always stable.

The stability diagram in the parameter space ($`\epsilon `$,$`K`$) is shown in Figure 3 for $`m=0.2`$, and compared to that of the Kuramoto model ($`m=0`$). This diagram is obtained from Eq. (52) with $`\lambda =0`$, for fixed $`D`$ and $`m`$. In this figure, we have also plotted the evolution of the order parameter amplitude for the parameter values marked in the stability diagram by (1) to (4). In all cases, the initial condition is taken sufficiently close to the incoherent solution, $`r=0`$.

(c) Bimodal frequency distribution, $`g(\mathrm{\Omega })=\frac{1}{2}[\delta (\mathrm{\Omega }\mathrm{\Omega }_0)+\delta (\mathrm{\Omega }+\mathrm{\Omega }_0)]`$.

Eq. (37) becomes

$`{\displaystyle \frac{Ke^{mD}}{2D}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^n}{n!}}{\displaystyle \frac{(n+mD)[m(\lambda +D)+n]}{[m(\lambda +D)+n]^2+m^2\mathrm{\Omega }_0^2}}=1.`$ (58)
In the high frequency limit, $`\mathrm{\Omega }_0\mathrm{}`$, we can find an analytical formula for $`\lambda `$ by inserting the following asymptotic expansion for the incomplete gamma function in (41) ,

$$\gamma (a,x)\frac{e^xx^a}{a},a\mathrm{},$$
(59)
where $`a=m(\lambda +D+i\mathrm{\Omega }_0)`$, and $`x=mD`$. The result is

$$1=\frac{K}{4}\left(\frac{1}{\lambda +D+i\mathrm{\Omega }_0}+\frac{1}{\lambda +Di\mathrm{\Omega }_0}\right),$$
(60)
which yields

$$\lambda =D+\frac{K}{4}+\frac{K}{4}i\sqrt{16\mathrm{\Omega }_0^2K^2}.$$
(61)
Re$`\lambda =0`$ gives the same critical coupling as the Kuramoto model, $`K_c=4D`$, for the same bimodal frequency distribution .

Figure 4 shows the stability diagram, which is obtained from Eq. (58) with Re$`\lambda =0`$, for $`D=1`$ and three different mass values $`m=0.1,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}6}`$. The stability diagram of the Kuramoto model ($`m=0`$) is also depicted. Notice how the curves corresponding to differerent masses tend to $`K=4D`$ as $`\mathrm{\Omega }_0\mathrm{}`$, as expected. Fig. 5 displays the evolution of the order parameter amplitude in different regions of the stability diagram corresponding to $`m=0.8`$ and $`D=1`$.

In conclusion, increasing $`m`$, $`D`$, $`\epsilon `$ (inertia, noise, and frequency spread) makes it more difficult to synchronize the oscillator population via stationary bifurcations from incoherence.

## III Mode-coupling equations and stationary solutions

Inspired by the previous linear stability analysis, we shall expand the distribution function using a basis of parabolic cylinder functions (or, equivalently, Hermite polynomials) of unit mean square norm ,

$`\rho (\theta ,\omega ,\mathrm{\Omega },t)=\left({\displaystyle \frac{2\pi D}{m}}\right)^{\frac{1}{4}}e^{\frac{m\omega ^2}{4D}}{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}c_n(\theta ,\mathrm{\Omega },t)\psi _n(\omega ).`$ (62)
In (62), we have defined

$`\psi _n(\omega )=\left(n!\sqrt{{\displaystyle \frac{2\pi D}{m}}}\right)^{\frac{1}{2}}D_n\left(\sqrt{{\displaystyle \frac{m}{D}}}\omega \right)`$ (63)
$`=\left(n!\mathrm{\hspace{0.17em}2}^n\sqrt{{\displaystyle \frac{2\pi D}{m}}}\right)^{1/2}H_n\left(\sqrt{{\displaystyle \frac{m}{2D}}}\omega \right)e^{m\omega ^2/4D},`$ (64)
so that $`_{\mathrm{}}^{\mathrm{}}\psi _n\psi _p𝑑\omega =\delta _{np}`$. The functions $`c_n(\theta ,\mathrm{\Omega },t)`$ are $`2\pi `$-periodic in $`\theta `$, and we have

$$_0^{2\pi }c_0(\theta ,\mathrm{\Omega },t)𝑑\theta =1,$$
(65)
as it follows from the normalization of $`\rho (\theta ,\omega ,\mathrm{\Omega },t)`$.

We shall find a system of mode-coupling equations for the coefficient functions $`c_n`$. Then we shall try to find stationary solutions for different frequency distributions $`g(\mathrm{\Omega })`$. This is not so easy in the general case, so that we shall follow a standard approach: We shall recognize in the system of equations for the coefficient functions a particular stationary solution corresponding to incoherence. Then we shall try to find a bifurcation equation for other stationary solutions branching off from incoherence. We shall see that even this requires different approximation schemes in order to succeed.

### A Mode-coupling equations

Let us insert Eq. (62) into the Fokker-Planck equation (6). We then obtain the following hierarchy of coupled partial differential equations for $`c_n(\theta ,\mathrm{\Omega },t)`$,

$`{\displaystyle \frac{c_0}{t}}=\sqrt{{\displaystyle \frac{D}{m}}}{\displaystyle \frac{c_1}{\theta }},`$ (66)
$`{\displaystyle \frac{c_1}{t}}=c_0{\displaystyle \frac{1}{m}}c_1\sqrt{2}\sqrt{{\displaystyle \frac{D}{m}}}{\displaystyle \frac{c_2}{\theta }},`$ (67)
$`\mathrm{}`$ (68)
$`{\displaystyle \frac{c_n}{t}}=\sqrt{n}c_{n1}{\displaystyle \frac{n}{m}}c_n\sqrt{n+1}\sqrt{{\displaystyle \frac{D}{m}}}{\displaystyle \frac{c_{n+1}}{\theta }}.`$ (69)
Here we have defined the operator

$`f=\sqrt{{\displaystyle \frac{D}{m}}}\left[{\displaystyle \frac{}{\theta }}{\displaystyle \frac{\mathrm{\Omega }+Kr\mathrm{sin}(\psi \theta )}{D}}\right]f.`$ (70)
In terms of the functions $`c_n`$, the equation for the order parameter becomes

$`re^{i\psi }={\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{i\theta }c_0(\theta ,\mathrm{\Omega },t)g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\theta .`$ (71)
$`r\mathrm{sin}(\psi \theta )={\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{sin}(\varphi \theta )c_0(\varphi ,\mathrm{\Omega },t)g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\varphi .`$ (72)
### B Incoherence and bifurcating stationary solutions

The incoherent $`\theta `$-independent solution, $`c_n=C_n(\mathrm{\Omega })`$, depends only on $`\mathrm{\Omega }`$. It can readily be obtained from the above hierarchy of equations, by ignoring all derivatives. The result is

$$C_n(\mathrm{\Omega })=\frac{1}{2\pi }\frac{1}{\sqrt{n!}}\left(\frac{m}{D}\right)^{n/2}\mathrm{\Omega }^n.$$
(73)
Inserting this into the expansion in Eq. (62), we indeed recover Eq. (9). In particular setting $`\mathrm{\Omega }=0`$, this yields the simplest incoherent solution, $`C_n=1/(2\pi )`$, corresponding to the unimodal frequency distribution. Other interesting stationary solutions are partially synchronized distributions, which depend on $`\theta `$. Notice that eq. (66) implies that $`c_1`$ does not depend on $`\theta `$.

In the following we analyze stationary solutions bifurcating from incoherence for different frequency distributions.

(a) Discrete unimodal frequency distribution, $`g(\mathrm{\Omega })=\delta (\mathrm{\Omega })`$.

Let us look for stationary solutions with finitely many nonzero coefficients $`c_n=0`$ for $`n>N`$ and $`g(\mathrm{\Omega })=\delta (\mathrm{\Omega })`$. Eq. (69) yields $`c_N=0`$, and therefore, $`c_N=K_Ne^{\frac{K}{D}r\mathrm{cos}(\psi \theta )}`$ ($`K_N=`$constant). Inserting this result in (69) for $`n=N1`$, we find $`c_{N1}=\sqrt{N}c_N/m`$. The corresponding solution is not periodic in $`\theta `$ unless $`K_N=0`$. Repeating this argument, we obtain the stationary solution

$$c_0(\theta )=\frac{e^{\frac{K}{D}r\mathrm{cos}(\psi \theta )}}{_0^{2\pi }e^{\frac{K}{D}r\mathrm{cos}(\psi \theta )}𝑑\theta },$$
(74)
where $`c_n=0`$ for $`n>0`$ and the normalization condition has been used. Observe that the resulting distribution, $`\rho (\theta ,\omega )=\left(\frac{2\pi D}{m}\right)^{1/4}c_0(\theta )\psi _0(\omega )e^{m\omega ^2/4D}`$, is factorized with respect to its two arguments, $`\theta `$ and $`\omega `$, cf. . It is remarkable that $`c_0(\theta )`$ is independent of $`m`$, and coincides with that obtained for the Kuramoto model. Therefore, from (71) the order parameter does not depend on inertia, and the bifurcation diagram for $`r`$ is exactly the same as that for the Kuramoto model.

(b) General frequency distributions.

For general frequency distributions, we could try to find stationary solutions of the mode-coupling equations (69) which bifurcate from incoherence. It is however more direct to work with the stationary Fokker-Planck equation (6) as follows. We consider stationary solutions as functions of a fixed value of the synchronization parameter $`r`$. Then we expand these solutions in power series of $`r`$ and write a hierarchy of equations for the coefficient functions. The first coefficient function should be the incoherent solution of synchronization parameter $`r=0`$. Inserting the power series probability density function into Eq. (7), we find the amplitude equation for stationary solutions bifurcating from incoherence. This procedure is explained in Appendix B. We quote here the result:

$$r=\frac{Kr}{2D}\alpha +\frac{(Kr)^3}{6}\beta +O((Kr)^4),$$
(75)
where

$`\alpha =e^{mD}[{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{1+\frac{\mathrm{\Omega }^2}{D^2}}}g(\mathrm{\Omega })d\mathrm{\Omega }`$ (76)
$`+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^p(1+\frac{p}{mD})^2}{p!}}`$ (77)
$`\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{(1+\frac{p}{mD})^2+\frac{\mathrm{\Omega }^2}{D^2}}}g(\mathrm{\Omega })d\mathrm{\Omega }],`$ (78)
and $`\beta `$ can be calculated numerically by solving a system of differential equations. As the inertia vanishes, $`mD0+`$, $`\alpha `$ becomes

$`\alpha ={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{g(\mathrm{\Omega })d\mathrm{\Omega }}{1+\frac{\mathrm{\Omega }^2}{D^2}}}{\displaystyle \frac{m}{D}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{\mathrm{\Omega }^2}{1+\frac{\mathrm{\Omega }^2}{D^2}}}g(\mathrm{\Omega })𝑑\mathrm{\Omega }+O(m^2D^2).`$ (79)
Notice that $`\alpha `$ in Equation (78) coincides with $`2D/K`$ in Eq. (37) provided $`\lambda =0`$ and $`g(\mathrm{\Omega })=g(\mathrm{\Omega })`$. This means that the critical coupling for bifurcation towards stationary synchronized states is obtained at $`K\alpha =2D`$, no matter what the symmetric frequency distribution may be.

Before interpreting the amplitude equation (75), we shall outline a procedure to obtain its coefficients based upon an uncontrolled closure assumption which is accurate for small values of $`mD`$. Consider the expansion (62). We expect that the coefficients $`c_n`$ of stationary solutions close to incoherence do not differ much from the coefficients of the latter. In view of the functional form (73), we anticipate that the coefficients $`c_n`$ approach zero as $`n\mathrm{}`$ faster for smaller values of the mass. Thus we shall now consider stationary solutions for general frequency distributions, such that $`c_n=0`$ in (69), for all $`n3`$. By solving equation (69) for such a stationary solution with $`n=2`$, we obtain

$$c_2(\theta ,\mathrm{\Omega })=\sqrt{\frac{m}{2D}}[\mathrm{\Omega }+Kr\mathrm{sin}(\psi \theta )]c_1(\mathrm{\Omega }).$$
(80)
We now insert this expression into (67). The resulting equation is solved for a $`c_0`$, which is $`2\pi `$-periodic in $`\theta `$ and obeys the normalization condition (8). We find

$`c_0(\theta ,\mathrm{\Omega })={\displaystyle \frac{e^{\frac{Kr}{D}\mathrm{cos}(\psi \theta )}\phi (\theta ,\mathrm{\Omega })}{Z(\mathrm{\Omega })}}`$ (81)
$`\phi (\theta ,\mathrm{\Omega })={\displaystyle _0^{2\pi }}[1mKr\mathrm{cos}(\psi \theta \eta )]`$ (82)
$`\times e^{\frac{1}{D}[\mathrm{\Omega }\eta +Kr\mathrm{cos}(\psi \theta \eta )]}d\eta ,`$ (83)
$`Z(\mathrm{\Omega })={\displaystyle _0^{2\pi }}e^{\frac{Kr}{D}\mathrm{cos}(\psi \theta )}\phi (\theta ,\mathrm{\Omega })𝑑\theta ,`$ (84)
In (81), $`r`$ should be determined so that

$`r={\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}c_0(\theta ,\mathrm{\Omega })\mathrm{cos}(\psi \theta )g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\theta ,`$ (86)
holds. The function $`c_1(\mathrm{\Omega })`$ can be obtained by integrating (67) with respect to $`\theta `$ and using the normalization condition for $`c_0`$ together with the $`2\pi `$-periodicity in $`\theta `$ of $`c_0`$ and $`c_2`$:

$$c_1=\frac{1}{2\pi }\sqrt{\frac{m}{D}}\left[\mathrm{\Omega }+Kr_0^{2\pi }\mathrm{sin}(\psi \theta )c_0(\theta ,\mathrm{\Omega })𝑑\theta \right].$$
(87)
Then, the stationary distribution can be approximated by

$`\rho (\theta ,\omega ,\mathrm{\Omega })\left({\displaystyle \frac{2\pi D}{m}}\right)^{1/4}[c_0(\theta ,\mathrm{\Omega })\psi _0(\omega )+c_1(\mathrm{\Omega })\psi _1(\omega )`$ (88)
$`+c_2(\theta ,\mathrm{\Omega })\psi _2(\omega )]e^{m\omega ^2/4D}.`$ (89)
Notice that a nonvanishing $`c_1(\mathrm{\Omega })`$ in (62) implies that the probability density is no longer even in $`\omega `$ ($`\psi _1`$ is an odd function of $`\omega `$). By (87), this occurs in the synchronized phase for the case of nonidentical oscillators. For the bimodal frequency distribution, Fig. 6 shows $`c_1`$ as function of $`\mathrm{\Omega }_0`$ for two different values of $`m`$. Both the approximate expression (87) and results of direct numerical simulations are depicted. Notice that the agreement between our approximation and the numerical result improves as $`m`$ decreases. Thus we observe that for each fixed nonzero $`\mathrm{\Omega }`$ and each fixed $`\theta `$, the distribution function is no longer peaked at $`\omega =0`$. However, the instantaneous frequency distribution, defined by

$$_0^{2\pi }_{\mathrm{}}^{\mathrm{}}\rho (\theta ,\omega ,\mathrm{\Omega })g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\theta ,$$

may turn out to be even in $`\omega `$. This is certainly true for the approximate stationary distribution (89), for

$$_{\mathrm{}}^{\mathrm{}}c_1(\mathrm{\Omega })g(\mathrm{\Omega })𝑑\mathrm{\Omega }=\sqrt{\frac{m}{D}}\frac{Kr}{2\pi }_0^{2\pi }_{\mathrm{}}^{\mathrm{}}c_0(\theta ,\mathrm{\Omega })\mathrm{sin}(\psi \theta )g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\theta =0,$$

as it follows from the definition of the order parameter and the expansion (62).

Another indication that the exact instantaneous frequency distribution may be even in $`\omega `$ is that the average frequency tends to zero as $`t+\mathrm{}`$. This would occur if the stable stationary distribution is even in $`\omega `$ (although, admittedly, distributions which are not even in $`\omega `$ may have zero mean). The result can be shown directly from the Fokker-Planck equation (6). We multiply that equation by $`\omega g(\mathrm{\Omega })`$ and integrate with respect to all variables. Then we obtain the following equation for the mean value $`\omega `$:

$$\frac{d}{dt}\omega +\frac{\omega }{m}=_{\mathrm{}}^{\mathrm{}}\mathrm{\Omega }g(\mathrm{\Omega })𝑑\mathrm{\Omega }.$$

The right hand side of this expression is zero if $`g(\mathrm{\Omega })=g(\mathrm{\Omega })`$. Then $`\omega `$ tends to zero exponentially fast as $`t+\mathrm{}`$.

Let us now go back to the problem of obtaining the coefficients in the amplitude equation (75). Let us expand (81) in powers of $`r`$ and insert the result in (86). We obtain an approximation to the coefficients of the amplitude equation (86). $`\alpha `$ is again given by (75), and the expression for $`\beta `$ is

$`\beta ={\displaystyle _{\mathrm{}}^{\mathrm{}}}[{\displaystyle \frac{3\frac{\mathrm{\Omega }^2}{D^2}\frac{3}{2}+\frac{39}{4}m\frac{\mathrm{\Omega }^2}{D}}{(1+\frac{\mathrm{\Omega }^2}{D^2})^2(4+\frac{\mathrm{\Omega }^2}{D^2})D^3}}`$ (90)
$`+{\displaystyle \frac{3m^2\mathrm{\Omega }^2\left(1+\frac{\mathrm{\Omega }^2}{2D^2}\right)}{2D^3\left(1+\frac{\mathrm{\Omega }^2}{D^2}\right)^2}}]g(\mathrm{\Omega })d\mathrm{\Omega }.`$ (91)
We now interpret the amplitude equation (75) in the usual way. Notice that the nonzero solution is approximately given by

$$Kr\sqrt{\frac{3(2DK\alpha )}{K\beta D}}.$$

Then the critical value of $`K`$ is $`K^{}=2D/\alpha `$, and the sign of $`\beta `$ determines the direction of the bifurcating branch of stationary solutions. Assume $`\alpha >0`$. Then the bifurcating stationary solution exists for $`K>K^{}`$ if $`\beta <0`$ (supercritical bifurcation). If $`\beta >0`$, we have a subcritical bifurcation and the partially synchronized stationary solution exists for $`K<K^{}`$. The critical coupling tends to infinity as $`\alpha 0+`$, and the bifurcation does not occur at positive couplings if $`\alpha <0`$. Given the formulas (78) or (79), $`\alpha `$ could become negative, depending on the value of $`m`$. Let us now analyze the bifurcation diagram corresponding to two different frequency distributions.

* Unimodal Lorentzian frequency distribution.

For an unimodal Lorentzian frequency distribution, the coefficients $`\alpha `$ and $`\beta `$ of (79) and (91) are approximately given by

$`\alpha ={\displaystyle \frac{D}{D+\gamma }}(1m\epsilon ),`$ (92)
$`\beta ={\displaystyle \frac{3}{4}}{\displaystyle \frac{1}{(D+\epsilon )^2(2D+\epsilon )}}[1+{\displaystyle \frac{m}{D}}\epsilon (3D+\epsilon )].`$ (93)
Apparently, the sign of $`\alpha `$ could again be negative provided $`m\epsilon `$ is sufficiently large. However, evaluation of the exact expression (78) shows that $`\alpha >0`$. Figure 7 shows $`\alpha `$ as a function of $`\epsilon `$ for two different masses. Notice that the approximate expression for $`\alpha `$ fits better the numerical result as $`m`$ and $`\epsilon `$ decrease.

In contrast with $`\alpha `$, $`\beta `$ can really change sign for a Lorentzian frequency distribution. Fig. 8 shows the good qualitative agreement between the approximate expression for $`\beta `$ and its numerical evaluation from the exact equations (B22).

There is a critical mass $`m_c`$ ($`\epsilon `$ and $`D`$ are kept fixed) for which $`\beta =0`$. This mass separates the supercritical and subcritical bifurcation regions. Setting $`\beta =0`$ in eq. (93) yields

$$m_c=\frac{D}{\epsilon (3D+\epsilon )}.$$
(94)
Figure 9 shows the qualitative agreement between the approximate expression (94), and the numerical result obtained from the solution of (B22).

Note that in the massless case, the stationary solution always branches off supercritically, independently of $`\epsilon `$. In the presence of inertia and for appropriate values of $`\epsilon `$, we have found a subcritical bifurcation. This prediction is illustrated by Figures 10 (parameters corresponding to a supercritical bifurcation) and 11 (parameters corresponding to a subcritical bifurcation). Figure 10 is a typical supercritical bifurcation diagram (order parameter amplitude versus $`K`$). In the subcritical case, Figure 11 shows the different evolution of the order parameter amplitude for two different initial conditions. Notice the bistability between incoherence and the stable synchronized solution typical of a subcritical bifurcation.
* Bimodal frequency distribution

For the bimodal frequency distribution, $`\alpha `$ may be zero. This means that there is a critical frequency $`\mathrm{\Omega }_0^c(m)`$ above which $`K^{}=\mathrm{}`$. In fact, we have shown in Fig. 4, looking to the branch corresponding to $`\lambda =0`$, that in case of a bimodal frequency distribution, it is possible to find some values of $`\mathrm{\Omega }_0`$ without its corresponding $`K`$ in the stability diagram of incoherence, due to the existence of an horizontal asymptote in the space parameter $`(K,\mathrm{\Omega }_0)`$. This means that there is not a finite coupling $`K^{}`$ where the stationary solution branches off. On the other hand, this never occurs in the bimodal Kuramoto model \[$`m=0`$, $`K_c/D=2(1+\mathrm{\Omega }_0^2/D^2)`$\], because in this case there is not such an horizontal asymptote.

Eq. (79) shows that the critical frequency at which $`\alpha =0`$ is

$$\frac{\mathrm{\Omega }_0^{\mathrm{}}}{D}=\sqrt{\frac{1}{mD}}.$$
(95)
Figure 12 compares the previous approximate expression to the exact critical frequency obtained from (78). Note that the approximation improves as $`m`$ decreases. Fig. 13 shows that $`K^{}`$ increases as $`\mathrm{\Omega }_0`$ does. $`K^{}`$ becomes infinity for $`\mathrm{\Omega }_0\mathrm{\Omega }_0^{\mathrm{}}`$.

For $`\alpha >0`$, there is another important critical frequency. In this case, the sign of $`\beta `$ decides whether the bifurcation is subcritical ($`\beta >0`$) or supercritical ($`\beta <0`$). The sign of $`\beta `$ depends on $`m`$ and $`\mathrm{\Omega }_0`$. For small mass, we can use the approximations (79) and (91) (ignoring terms which are quadratic in the mass). Then we find that the critical frequency at which $`\beta =0`$ is

$$\frac{\mathrm{\Omega }_0^c}{D}=\sqrt{\frac{1}{2+\frac{13}{8}mD}}.$$
(96)
We have solved numerically the system of equations (B22), and compared with the approximation (91) for $`\beta `$. In Figure 14 the coefficient $`\beta `$ is plotted as a function of $`\mathrm{\Omega }_0`$ for different masses. Similarly, Figure 15 shows how the critical frequency $`\mathrm{\Omega }_0^c`$ varies as a function of $`m`$. Note that the analytical approximation improves as $`m`$ goes to zero (as expected).

Notice that $`\mathrm{\Omega }_0^c<\mathrm{\Omega }_0^{\mathrm{}}`$. Then the branch of synchronized stationary states bifurcates subcritically from incoherence at $`K^{}=\mathrm{}`$, provided $`\mathrm{\Omega }_0>\mathrm{\Omega }_0^{\mathrm{}}`$.

Summarizing, inertia favors the subcritical character of bifurcations describing the transition from incoherence to the partial synchronized state. In fact in the bimodal case, the transition of incoherence to synchronization will most likely occur as a subcritical bifurcation as inertia increases. For a continuous unimodal Lorentzian frequency distribution, inertia may turn subcritical the supercritical bifurcation to stationary synchronized states, which is always found in the masslesss Kuramoto model.

## IV High frequency limit

The high frequency limit of Eq. (6) can be analyzed by means of a multiscale method. For the Kuramoto model , this method lead to the result that the probability density is (to leading order) a superposition of different components corresponding to the different peaks of the oscillator frequency distribution. To apply this method to the present model, we shall assume that the frequency distribution $`g(\mathrm{\Omega })`$ has $`m`$ maxima located at $`\mathrm{\Omega }_0\nu _l,l=1,\mathrm{}.,m`$, so that $`g(\mathrm{\Omega })d\mathrm{\Omega }`$ tends to the distribution

$$\mathrm{\Gamma }(\nu )\underset{l=1}{\overset{m}{}}\alpha _l\delta (\nu \nu _l)d\nu ,$$
(97)
as $`\mathrm{\Omega }_0\mathrm{}`$. In order to simplify the calculations, we change variables to a comoving frame:

$$\beta =\theta \mathrm{\Omega }t\theta \frac{\nu }{ϵ}t,$$
(98)
and let $`\omega ^{}=\omega \mathrm{\Omega }`$, where

$$ϵ=\frac{1}{\mathrm{\Omega }_0}1$$
(99)
Then we obtain the following equations:

$`{\displaystyle \frac{\rho _j}{t}}=D{\displaystyle \frac{^2\rho _j}{\omega ^2}}{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{\omega ^{}}}(\rho _jU_j)\omega ^{}{\displaystyle \frac{\rho _j}{\beta }},`$ (100)
$`U_j=\omega ^{}+\text{Im}K\{{\displaystyle \underset{l=1}{\overset{m}{}}}\alpha _le^{i(\nu _l\nu _j)t/ϵ}`$ (101)
$`\times {\displaystyle }e^{i(\beta ^{}\beta )}\rho (\beta ^{},\omega ^{},t,\nu _l;ϵ)d\beta ^{}d\omega ^{}\},`$ (102)
$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _0^{2\pi }}\rho _j(\beta ^{},t,\omega ^{};ϵ)𝑑\beta ^{}𝑑\omega ^{}=1,`$ (103)
where $`\rho _j=\rho (\beta ,t,\omega ^{},\nu _j;ϵ)`$, and $`\rho \alpha _j\rho _j\delta (\nu \nu _j)`$, . We now define fast and slow time scales, $`\tau =t/ϵ`$ and $`t`$, and make the Ansatz:

$$\rho =\underset{n=0}{\overset{2}{}}\rho ^{(n)}(\beta ,\omega ^{},\tau ,t,\nu )ϵ^n+O(ϵ^3).$$
(104)
Inserting this into the governing equations, we obtain the following hierarchy:

$`{\displaystyle \frac{\rho _j^{(0)}}{\tau }}=0,`$ (105)
$`{\displaystyle \frac{\rho _j^{(1)}}{\tau }}={\displaystyle \frac{\rho _j^{(0)}}{t}}+D{\displaystyle \frac{^2\rho _j^{(0)}}{\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{\omega ^{}}}(\omega ^{}\rho _j^{(0)})`$ (106)
$`\omega ^{}{\displaystyle \frac{\rho _j^{(0)}}{\beta }}{\displaystyle \frac{K}{m}}{\displaystyle \frac{}{\omega ^{}}}\{\rho _j^{(0)}[\text{Im}(\alpha _je^{i\beta }Z_j^{(0)}`$ (107)
$`+{\displaystyle \underset{lj}{}}\alpha _je^{i(\nu _l\nu _j)\tau }e^{i\beta }Z_l^{(0)})]\},`$ (108)
where

$$Z_j^{(0)}(t)=_{\mathrm{}}^{\mathrm{}}_0^{2\pi }e^{i\eta }\rho _j^{(0)}(\eta ,\omega ,t,\nu _j)𝑑\eta 𝑑\omega .$$
(109)
Elimination of secular terms yields the following condition:

$``$ $`{\displaystyle \frac{\rho _j^{(0)}}{t}}+D{\displaystyle \frac{^2\rho _j^{(0)}}{\omega ^2}}\omega ^{}{\displaystyle \frac{\rho _j^{(0)}}{\beta }}`$ (110)
$``$ $`{\displaystyle \frac{1}{m}}{\displaystyle \frac{}{\omega ^{}}}\left\{\rho _j^{(0)}\left[\omega ^{}+K\alpha _j\text{Im}(e^{i\beta }Z_j^{(0)})\right]\right\}=0.`$ (111)
Notice that this equation corresponds to the equation of a unimodal frequency distribution, already studied in . Its solution evolves towards the following stationary state as time evolves:

$`\rho _j^{(0)}(\beta ,\omega ^{})`$ $`=`$ $`\sqrt{{\displaystyle \frac{m}{2\pi D}}}e^{(m/2D)\omega ^2}`$ (112)
$`\times `$ $`{\displaystyle \frac{e^{(K\alpha _j/D)R_j\mathrm{cos}(\mathrm{\Psi }_j\beta )}}{_0^{2\pi }e^{(K\alpha _j/D)R_j\mathrm{cos}(\mathrm{\Psi }_j\beta ^{})}𝑑\beta ^{}}},`$ (113)
where

$`R_je^{i\mathrm{\Psi }_j}={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _0^{2\pi }}e^{i\eta }\rho _j^{(0)}(\eta ,\omega )𝑑\eta 𝑑\omega \underset{t\mathrm{}}{lim}Z_j^{(0)}(t).`$ (114)
Incoherence of a component corresponds to $`R_j=0`$. We know that as $`K\alpha _j`$ surpasses $`2D`$, a stable synchronized solution bifurcates from incoherence. In the particular case of a symmetric bimodal frequency distribution, the bifurcation value is $`K=4D`$, independently of the inertia $`m`$. This agrees with the results of the linear stability analysis for $`\mathrm{\Omega }_0\mathrm{}`$. All the results previously obtained for the Kuramoto model can be applied to the present model without any modification . Thus, both a stable standing wave solution (SW) and an unstable travelling wave solution (TW) bifurcate supercritically from incoherence at $`K=4D`$ . In Figure 16 illustrates the comparison between the asymptotic solution (113) and the result of direct numerical simulation for a large enough value of the frequency $`\mathrm{\Omega }_0`$. Finally, Fig. 17 shows the global bifurcation diagram of the bimodal case for the case of positive $`\alpha `$ and $`\beta `$. As explained in the previous Section, the branch of synchronized stationary states bifurcates subcritically from incoherence at $`K^{}=\mathrm{}`$, provided $`\mathrm{\Omega }_0>\mathrm{\Omega }_0^{\mathrm{}}`$. Thus $`K^{}`$ and the end of the SW and TW branches in Fig. 17 should extend to infinity as $`\mathrm{\Omega }_0+\mathrm{}`$.

## V Numerical results

Four different numerical methods have been used. Numerical simulations of the system of Langevin equations (3) were carried out for a large number of oscillators (N=20 000), using an Euler method. A standard finite difference method was used to solve numerically the Fokker-Planck equation (6), or the system of partial differential equations (69). In addition, we have used a simple spectral method, which generalizes the one proposed in . The idea is to solve a set of ordinary differential equations for moments of $`\rho `$:

$`(x_i^j)_k:={\displaystyle _0^{2\pi }}c_i(\theta ,\mathrm{\Omega }_k,t)\mathrm{cos}[j(\psi \theta )]𝑑\theta ,`$ (115)
$`(y_i^j)_k:={\displaystyle _0^{2\pi }}c_i(\theta ,\mathrm{\Omega }_k,t)\mathrm{sin}[j(\psi \theta )]𝑑\theta ,`$ (116)
$`r={\displaystyle _0^{2\pi }}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}\rho (\theta ,\omega ,\mathrm{\Omega },t)\mathrm{cos}(\psi \theta )`$ (117)
$`\times g(\mathrm{\Omega })d\mathrm{\Omega }d\omega d\theta .`$ (118)
The coefficient functions $`c_i(\theta ,\mathrm{\Omega },t)`$ are the same as in the parabolic cylinder expansion (62). The integral in (118) will be approximated by a suitable quadrature formula, as

$`r{\displaystyle \underset{q=1}{\overset{Q}{}}}\alpha _q{\displaystyle _0^{2\pi }}c_0(\mathrm{\Omega }_q,\theta ,t)\mathrm{cos}(\psi \theta )𝑑\theta ={\displaystyle \underset{q=1}{\overset{Q}{}}}\alpha _q(x_0^1)_q.`$ (119)
For instance, the Gauss-Laguerre quadrature has been chosen case for the case of a Lorentzian frequency distribution. The system of ordinary differential equations is given by:

$`(\dot{x}_i^j)_k=j\sqrt{i}\sqrt{{\displaystyle \frac{D}{m}}}(y_{i1}^j)_k+\sqrt{i}{\displaystyle \frac{1}{\sqrt{mD}}}\mathrm{\Omega }_k(x_{i1}^j)_k`$ (120)
$`+\sqrt{{\displaystyle \frac{i}{mD}}}{\displaystyle \frac{Kr}{2}}\left[(y_{i1}^{j+1})_k(y_{i1}^{j1})_k\right]{\displaystyle \frac{i}{m}}(x_i^j)_k`$ (121)
$`+\sqrt{i+1}\sqrt{D}mj(y_{i+1}^j)_kj\dot{\psi }(y_i^j)_k,`$ (122)
$`(\dot{y}_i^j)_k=j\sqrt{i}\sqrt{{\displaystyle \frac{D}{m}}}(x_{i1}^j)_k+\sqrt{i}{\displaystyle \frac{1}{\sqrt{mD}}}\mathrm{\Omega }_k(y_{i1}^j)_k`$ (123)
$`+\sqrt{{\displaystyle \frac{i}{mD}}}{\displaystyle \frac{Kr}{2}}\left[(x_{i1}^{j1})_k(x_{i1}^{j+1})_k\right]{\displaystyle \frac{i}{m}}(y_i^j)_k`$ (124)
$`\sqrt{i+1}\sqrt{D}mj(x_{i+1}^j)_k+j\dot{\psi }(x_i^j)_k,`$ (125)
$`i=1,\mathrm{},N,j=1,\mathrm{},M,`$ (126)
$`(\dot{x}_0^j)_k=j\sqrt{{\displaystyle \frac{D}{m}}}(y_1^j)_kj\dot{\psi }(y_0^j)_k,`$ (127)
$`(\dot{y}_0^j)_k=j\sqrt{{\displaystyle \frac{D}{m}}}(x_1^j)_k+j\dot{\psi }(x_0^j)_k,`$ (128)
$`i=0,j=1,\mathrm{},M.`$ (129)
In order to numerically simulate this system, it is necessary to truncate the hierarchy after a reasonable number of modes $`N`$, and $`M`$. These numbers will depend on the inertia $`m`$, the spread, and the coupling strength. They should be chosen large enough, so that the numerical results do not depend on $`N`$ and $`M`$.

Fig. 18 shows the ratio $`c_n^{max}/c_0^{max}`$ as a function of $`n`$, for the stationary solution and for various mass values. $`c_n^{max}`$ is the maximum value of the stationary coefficient $`c_n(\theta ,\mathrm{\Omega }_0=1)`$. It has been obtained from a finite-difference solution of Eq. (6). Note that this ratio increases as inertia does. Thus the truncation approximation (89) ceases to make sense for larger values of inertia. Whether truncating (at higher order) the moment equations (115) and (116) yields good numerical results, is tested in Fig. 19. This figure shows how closely the previous system (with moments of order 4 or 10) resembles the direct solution of the Langevin equations. We notice that the system containing moments of order 10 is rather close to the solution of the Langevin equations, but it does not contain the fluctuations unavoidable in stochastic methods. If we are interested in solving the nonlinear Fokker-Planck equation, the system of moment equations seems a good alternative to solving Langevin equations for a large number of oscillators.

## VI conclusions

We have investigated synchronization properties of a model of globally and nonlinearly coupled phase oscillators, where the effects of white noise, inertia and spread in the natural frequency distribution are all considered. The linear stability of the incoherent solution is rigorously analyzed. Stationary and time-dependent solutions of the standing and travelling wave type are obtained by a variety of perturbative methods. These include finding and approximately solving mode-coupling equations for expansions of the probability distribution in parabolic cylinder functions, finding amplitude equations near bifurcation points, and hierarchy-closure assumptions. Numerical simulations of the different equations and the original model favorably agree with the different perturbative results.

Inertia changes the stability boundaries of the incoherent solution in a non trivial way. In the case of an unimodal Lorentzian frequency distribution, incoherence is stabilized, but the effect of mass completely drops out if the frequency spread vanishes. For a discrete bimodal frequency distribution, both stability boundaries and the character of the transition from incoherence to synchronized states depend on the values of the natural frequency $`\mathrm{\Omega }_0`$ and on inertia. The effect of inertia on the stationary solution is dramatic in some cases. In general, inertia hardens the synchronization transition: it may render subcritical (hard) originally supercritical (soft) transitions (in unimodal Lorentzian frequency distributions), or it increases the region in parameter space where the transition is subcritical (discrete bimodal frequency distributions). Analytical as well as numerical calculations confirm these findings.

## VII Acknowledgments

This work was supported, in part, by the European Union TMR contract ERB FMBX-CT97-0157, by UNESCO under contract UVO–ROSTE 875.629.9, by the Spanish DGES grant PB98-0142-C04-C01, and by the GNFM of the Italian C.N.R. All numerical simulations were conducted at CASPUR–Rome.

## A Integrals of special functions and eigenvalue equation

In (31) there appear two integrals over $`w`$:

$`A={\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p\rho _0^{}𝑑w,`$ (A1)
$`B={\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p𝑑w`$ (A2)
The second integral is directly $`𝒜_p(\sqrt{mD})`$. The first integral equals

$`A={\displaystyle \frac{m}{(2\pi )^{\frac{3}{2}}D}}{\displaystyle _{\mathrm{}}^+\mathrm{}}(wi2\sqrt{mD})e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p(w)𝑑w,`$ (A3)
because of (28). An equivalent expression is

$`A={\displaystyle \frac{2m}{(2\pi )^{\frac{3}{2}}D}}{\displaystyle _{\mathrm{}}^+\mathrm{}}D_p(w){\displaystyle \frac{}{w}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}𝑑w.`$ (A4)
This can be written as

$`A={\displaystyle \frac{2m}{(2\pi )^{\frac{3}{2}}D}}{\displaystyle _{\mathrm{}}^+\mathrm{}}D_p(w){\displaystyle \frac{}{w}}e^{\frac{1}{4}\left(wi2x\right)^2}𝑑w`$ (A5)
$`={\displaystyle \frac{m}{i(2\pi )^{\frac{3}{2}}D}}{\displaystyle _{\mathrm{}}^+\mathrm{}}D_p(w){\displaystyle \frac{}{x}}e^{\frac{1}{4}\left(wi2x\right)^2}𝑑w`$ (A6)
provided we set $`x=\sqrt{mD}`$ after performing the calculations. Thus we see that

$`A={\displaystyle \frac{m}{i(2\pi )^{\frac{3}{2}}D}}𝒜_p^{}(\sqrt{mD}).`$ (A7)
This expression and the previous one for $`B`$ yield (33).

Let us now obtain $`𝒜_p(x)`$ for integer $`p`$ and $`x>0`$. By using that $`D_p(x)`$ is even for even $`p`$ and odd for odd $`p`$, we can write

$`𝒜_p(x)=e^{x^2}{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\frac{w^2}{4}+iwx}D_p(w)𝑑w`$ (A8)
$`=e^{x^2}{\displaystyle _{\mathrm{}}^+\mathrm{}}e^{\frac{w^2}{4}}(\mathrm{cos}xw+i\mathrm{sin}xw)D_p(w)𝑑w`$ (A9)
$`=2e^{x^2}{\displaystyle _0^+\mathrm{}}D_p(w)e^{\frac{w^2}{4}}f(xw)𝑑w,`$ (A10)
where $`f(xw)=\mathrm{cos}xw`$ for $`p`$ even and $`f(xw)=i\mathrm{sin}xw`$ for $`p`$ odd. We can find (cf. §7.741 of ):

$`{\displaystyle _0^+\mathrm{}}D_p(w)e^{\frac{w^2}{4}}\mathrm{cos}xwdw=(1)^n\sqrt{{\displaystyle \frac{\pi }{2}}}e^{\frac{x^2}{2}}x^{2n},`$ (A11)
$`{\displaystyle _0^+\mathrm{}}D_p(w)e^{\frac{w^2}{4}}\mathrm{sin}xwdw=(1)^n\sqrt{{\displaystyle \frac{\pi }{2}}}e^{\frac{x^2}{2}}x^{2n+1},`$ (A12)
provided $`x>0`$. From these formulas, we obtain

$`𝒜_{2n}(x)=(1)^n\sqrt{2\pi }e^{\frac{x^2}{2}}x^{2n},`$ (A13)
$`𝒜_{2n+1}(x)=i(1)^n\sqrt{2\pi }e^{\frac{x^2}{2}}x^{2n+1}.`$ (A14)
These two expressions are equivalent to (35).

## B Calculation of the coefficients $`\alpha `$ and $`\beta `$

We will now calculate coefficients $`\alpha `$ and $`\beta `$ of Eq. (75). The idea is simply to expand the right hand side in the definition of the order parameter,

$$r=_0^{2\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\mathrm{cos}(\psi \theta )\rho (\theta ,\omega ,\mathrm{\Omega },t)g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\omega 𝑑\theta .$$
(B1)
as a power series in $`r`$. To this end, we fix $`r`$ and $`\psi `$ in Eq. (6), expand its stationary solution in powers of $`r`$ and insert the result in (B1). As we do not have a closed formula for the stationary solution of (6), we shall first derive a hierarchy of equations for its coefficients in the series in powers of $`r`$. Thus we have

$`\rho (\theta ,\omega ,\mathrm{\Omega };r)={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{G_n(\theta \psi ,\omega ,\mathrm{\Omega })}{n!}}r^n,`$ (B2)
$`r={\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{r_n}{n!}}r^n.`$ (B3)
Then

$$r_n=_0^{2\pi }_{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}\mathrm{cos}\theta G_n(\theta ,\omega ,\mathrm{\Omega })g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\omega 𝑑\theta ,$$
(B4)
where we have absorbed $`\psi `$ in the definition of $`G_n`$. Comparison with eq. (75), establishes that $`r_0=r_2=0`$, $`r_1=K\alpha /2D`$, and $`r_3=K^3\beta `$. From (B2) and (6), we obtain

$$\frac{D}{m^2}\frac{^2G_n}{\omega ^2}+\frac{1}{m}\frac{}{\omega }[(\omega \mathrm{\Omega })G_n]\omega \frac{G_n}{\theta }=\frac{n}{m}K\mathrm{sin}\theta \frac{G_{n1}}{\omega },$$
(B5)
with $`G_10`$. Since we are trying to find solutions bifurcating from incoherence, we should have $`G_0=\rho _0(\omega ,\mathrm{\Omega })`$, i.e. the $`\theta `$-independent incoherent solution (9). This directly confirms that $`r_0=0`$.

The unknowns $`G_n(\theta ,\omega ,\mathrm{\Omega })`$ are $`2\pi `$-periodic in $`\theta `$, so that we may Fourier expand them as

$$G_n(\theta ,\omega ,\mathrm{\Omega })=\underset{l=\mathrm{}}{\overset{\mathrm{}}{}}Z_n^l(\omega ,\mathrm{\Omega })e^{il\theta },$$
(B6)
where

$$Z_n^l(\omega ,\mathrm{\Omega })=\frac{1}{2\pi }_0^{2\pi }e^{il\theta }G_n(\theta ,\omega ,\mathrm{\Omega })𝑑\theta .$$
(B7)
Here $`Z_0^l0`$ for $`l0`$ because $`G_0`$ is the $`\theta `$-independent incoherent solution. The condition that the probability density function be real yields $`\overline{G_n}=G_n`$, which in turn implies $`Z_n^l=\overline{Z_n^l}`$. Inserting these expressions into (B4), we find

$$r_n=2\pi _{\mathrm{}}^{\mathrm{}}_{\mathrm{}}^{\mathrm{}}Re(Z_n^1)g(\mathrm{\Omega })𝑑\mathrm{\Omega }𝑑\omega .$$
(B8)
The unknowns $`Z_n^l`$ satisfy the following hierarchy of equations:

$$\frac{D}{m^2}\frac{d^2Z_n^l}{d\omega ^2}+\frac{1}{m}\frac{d}{d\omega }[(\omega \mathrm{\Omega })Z_n^l]il\omega Z_n^l=\frac{n}{m}\frac{K}{2i}\frac{d}{d\omega }[Z_{n1}^{l1}Z_{n1}^{l+1}].$$
(B9)
The normalization condition for $`\rho `$ and the incoherent solution together with (B7) imply

$$_{\mathrm{}}^{\mathrm{}}Z_n^0g(\omega )𝑑\omega =\delta _{n0}.$$
(B10)
In order to obtain $`\alpha `$, and $`\beta `$, we should show that $`r_2=0`$ and calculate $`r_1`$ and $`r_3`$. We do this by means of (B9), with $`n=1,2,3`$.

(a) Case n=1. Eq. (B9) becomes

$`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{d^2Z_1^1}{d\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{d}{d\omega }}[(\omega \mathrm{\Omega })Z_1^1]i\omega Z_1^1`$ (B11)
$`={\displaystyle \frac{1}{m}}{\displaystyle \frac{K}{2i}}{\displaystyle \frac{dZ_0^0}{d\omega }},`$ (B12)
where $`Z_0^0=G_0`$ is the incoherent solution (9), and we have used $`Z_0^2=0`$. We can solve Eq. (B12) by means of the solution of Eq. (17), in which we set $`\lambda =0`$ and replace $`b_1(\omega ,\mathrm{\Omega })`$ with $`Z_1^1(\omega ,\mathrm{\Omega })`$, and the right-hand side with $`i\frac{mK}{2D}dZ_0^0/d\omega `$. Then we obtain

$`Z_1^1(\omega ;\mathrm{\Omega })`$ $`=`$ $`{\displaystyle \frac{iK}{2m}}e^{\frac{m}{4D}(\omega \mathrm{\Omega })^2}`$ (B13)
$`\times `$ $`{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{_{\mathrm{}}^{\mathrm{}}e^{\left(\frac{w}{2}i\sqrt{mD}\right)^2}D_p\frac{dZ_0^0}{d\omega }𝑑w}{\sqrt{2\pi }p!\left(\frac{p}{m}+i\mathrm{\Omega }+D\right)}}D_p(w),`$ (B14)
The previous analogy allows us to calculate $`r_1`$ from (B8) and Equations (31), (34) and (35) (with $`\lambda =0`$). The result is

$`r_1={\displaystyle \frac{KD}{2}}e^{mD}[{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{D^2+\mathrm{\Omega }^2}}g(\mathrm{\Omega })d\mathrm{\Omega }`$ (B15)
$`+{\displaystyle \underset{p=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(mD)^p(1+\frac{p}{mD})^2}{p!}}`$ (B16)
$`\times {\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{1}{(D+\frac{p}{m})^2+\mathrm{\Omega }^2}}g(\mathrm{\Omega })d\mathrm{\Omega }].`$ (B17)
(b) Case n=2. We now show that $`r_2=0`$.
In order to analyze $`r_2`$, it is necessary to consider $`Z_2^1`$, which is the solution to the following equation:

$$\frac{D}{m^2}\frac{d^2Z_2^1}{d\omega ^2}+\frac{1}{m}\frac{d}{d\omega }[(\omega \mathrm{\Omega })Z_2^1]i\omega Z_2^1=\frac{2}{m}\frac{K}{2i}\frac{d}{d\omega }[Z_1^0Z_1^2].$$
(B18)
We shall show that $`Z_1^0=Z_1^2=0`$. If this is so, the resulting homogeneous equation can be transformed to the parabolic cylinder equation with a quadratic potential having complex coefficients. Its solution cannot decay to zero both as $`\omega \pm \mathrm{}`$ unless it is identically zero. The reason $`Z_1^0`$ and $`Z_1^2`$ are zero is similar. The equation for $`Z_1^0`$ is homogeneous ($`Z_0^l=0`$ for $`l0`$), and it has a solution $`C(\mathrm{\Omega })e^{m(\omega \mathrm{\Omega })^2/(2D)}`$. The normalization condition (B10) then implies $`C0`$. $`Z_1^2`$ again obeys a homogeneous equation which can be transformed into the parabolic cylinder equation with a complex potential. The only solution which decays to zero as $`\omega \pm \mathrm{}`$ is again $`Z_1^2=0`$.

(c) Case n=3.
The equations for $`Z_3^1`$ are

$`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{d^2Z_3^1}{d\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{d}{d\omega }}[(\omega \mathrm{\Omega })Z_3^1]i\omega Z_3^1={\displaystyle \frac{3}{m}}{\displaystyle \frac{K}{2i}}{\displaystyle \frac{d}{d\omega }}[Z_2^0Z_2^2],`$ (B19)
$`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{d^2Z_2^2}{d\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{d}{d\omega }}[(\omega \mathrm{\Omega })Z_2^2]2i\omega Z_2^2={\displaystyle \frac{2}{m}}{\displaystyle \frac{K}{2i}}{\displaystyle \frac{d}{d\omega }}Z_1^1,`$ (B20)
$`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{d^2Z_2^0}{d\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{d}{d\omega }}[(\omega \mathrm{\Omega })Z_2^0]={\displaystyle \frac{2}{m}}{\displaystyle \frac{K}{2i}}{\displaystyle \frac{d}{d\omega }}[Z_1^1Z_1^1]={\displaystyle \frac{2K}{m}}Im{\displaystyle \frac{d}{d\omega }}Z_1^1,`$ (B21)
$`{\displaystyle \frac{D}{m^2}}{\displaystyle \frac{d^2Z_1^1}{d\omega ^2}}+{\displaystyle \frac{1}{m}}{\displaystyle \frac{d}{d\omega }}[(\omega \mathrm{\Omega })Z_1^1]i\omega Z_1^1={\displaystyle \frac{1}{m}}{\displaystyle \frac{K}{2i}}{\displaystyle \frac{d}{d\omega }}Z_0^0.`$ (B22)
This system of equations has to be solved numerically in order to obtain the coefficient $`\beta `$.