for convenience assume that $\hat d_p$ is small, in which case $(\theta
^+_i-\theta ^-_i)$ is small, and we can use the approximation, $$\theta _i^+\cong \theta ^-_i+\hat d_ph(\theta ^-_i); \ p=0,1 \ .$$

Following Kuramoto we introduce the complex valued order parameter $R(t)$, $$R(t)=\frac{1}{N} \sum ^N_{j=1} e^{i\theta _j(t)} \ ,$$ in terms of which Eq. (1) can be rewritten as $$d\theta _i/dt=\omega _i+(K/N)Im[e^{-i\theta _i}R(t)]-h(\theta
_i)\Delta (t)+n(t) \ .$$

In our analysis in Sec. III we will take the limit $N\rightarrow
\infty$ useful for approximating the situation where $N\gg 1$. In that limit it is appropriate to describe the system state by a continuous distribution function $f(\theta ,\omega ,t)$, where $$\int ^{2\pi }_0f(\theta ,\omega ,t)\frac{d\theta}{2\pi }=1 \ ,$$ and the fraction of oscillators with angles and natural frequencies in the ranges $(\theta , \theta +d\theta )$ and $(\omega ,\omega +d\omega )$ is $f(\theta ,\omega ,t)g(\omega
)d\omega d\theta /2\pi $. The conservation of the number of oscillators then gives the time evolution equation for $f(\theta
,\omega ,t)$, $$\frac{\partial f}{\partial t}+\frac{\partial}{\partial \theta
}\left\{ f\left[ \omega +KIm(R(t)e^{-i\theta })-h(\theta )\Delta
(t)\right]\right\}=\xi \frac{\partial ^2f}{\partial \theta ^2} \ ,$$ $$R^*(t)=\int d\omega f_1(\omega ,t)g(\omega ) \ ,$$ where $R^*$ denotes the complex conjugate of $R$, $f(\omega
,\theta ,t)\equiv 1$ 