I am trying to derive the critical coupling strength for synchronisation in a network of phase oscillators with noisy input.

I am following the steps outlined in Sakaguchi, Hidetsugu. "Cooperative phenomena in coupled oscillator systems under external fields." Progress of theoretical physics 79.1 (1988): 39-46.

However I am stuck in an expansion and I would be grateful if someone had an insight to help.

I have a problem to derive the expansion in Eq. 25 but i will also copy here the equations for convenience.

We consider phase oscillators with all-to-all coupling in the thermodynamic limit, and thereby have expressed the mean-field Fokker-Planck for the phase distribution

$$ \frac{\partial n(\psi ; \omega)}{\partial t}=-\frac{\partial}{\partial \psi}\{(\omega-K \sigma \sin \psi) n(\psi ; \omega)\}+D \frac{\partial^{2}}{\partial \psi^{2}} n(\psi ; \omega), $$

where $\psi$ denotes the phase and $\omega$ the natural frequency, $\sigma$ the synchronisation order parameter, and $K$ the coupling strength. We compute the stationary solution: $$ \begin{aligned} n(\psi ; \omega)=& \exp \left(\frac{-K \sigma+\omega \psi+K \sigma \cos \psi}{D}\right) n(0 ; \omega) \\ & \times\left\{1+\frac{\left(e^{-2 \pi \omega / D}-1\right) \int_{0}^{\phi} e^{(-\omega \phi-K \sigma \cos \phi) / D} d \phi}{\int_{0}^{2 \pi} e^{(-\omega \phi-K \sigma \cos \phi) / D} d \phi} \right\}, \end{aligned} $$

and we are about to write a self-consistent equation for the synchronisation order parameter $\sigma$. We have from previously (Eq.9) $$ \begin{aligned} \sigma \exp \left(i \phi_{0}\right) &=\int_{-\infty}^{\infty} d \omega g(\omega) \int_{0}^{2 \pi} d \psi n(\psi ; \omega) \exp (i \psi), \end{aligned} $$ where $g(\omega)$ denotes the natural frequency distribution, and $\phi_0$ stands for the mean phase.

We write the self consistent equation for the order parameter $$ \sigma=\int_{-\infty}^{\infty} d \omega g\left(\omega+\omega_{0}\right) \int_{0}^{2 \pi} d \psi n(\psi ; \omega) \exp (i \psi), $$ and to find the critical coupling strength we proceed to expand the expression for the order parameter in powers of $K \sigma/D$

$$ \begin{aligned} \sigma=& K \sigma\left[\frac{1}{2} \int_{-\infty}^{\infty} g\left(D \omega+\omega_{0}\right) \frac{d \omega}{\omega^{2}+1}-\frac{K^{2} \sigma^{2}}{4 D^{2}} \int_{-\infty}^{\infty} g\left(D \omega+\omega_{0}\right)\left\{\frac{1}{\omega^{2}+4}-\frac{\omega}{\left(\omega^{2}+1\right)^{2}}\right\} d \omega\right.\\ &\left.+O\left(\left(\frac{K \sigma}{D}\right)^{4}\right)\right]. \end{aligned} $$

So my question is how can I derive this expansion. The author mentions that they consider that $g(\omega+\omega_0)$ is symmetric about 0, and that the imaginary part of the self-consistent equation for the order parameter is zero. So the expansion is only for the real part.

I would be grateful for any insight!