# Penrose method for Kuramoto model with inertia and noise

Artem Alexandrov<EMAIL_ADDRESS>Moscow Institute of Physics and
Technology, Dolgoprudny 141700, Russia Institute for Information Transmission
Problems, Moscow, 127994, Russia  Alexander Gorsky  Institute for Information
Transmission Problems, Moscow, 127994, Russia  Laboratory of Complex Networks,
Center for Neurophysics and Neuromorphic Technologies, Moscow, Russia

###### Abstract

Using the Penrose method of instability analysis, we consider the
synchronization transition in the Kuramoto model with inertia and noise with
all-to-all couplings. Analyzing the Penrose curves, we identify the appearance
of cluster and chimera states in the presence of noise. We observe that noise
can destroy chimera and biclusters states. The critical coupling describing
bifurcation from incoherent to coherent state is found analytically. To
confirm our propositions based on the Penrose method, we perform numerical
simulations.

##  1 Introduction

Synchronization phenomenon is widely observed in nature [1], which provoking a
natural desire to investigate in details the essence of it. Kuramoto model,
proposed in [2, 3], is one of the many models of synchronization and attracts
a great attention due to both simplicity and existence of non-trivial physics.
Despite a huge progress in the investigation of Kuramoto model properties,
some blank spaces still exist. These white spots both have physical motivation
and deserve for a formulation in rigorous mathematical form. Since a lot of
research devoted to Kuramoto model uses numerical simulations, we believe that
new analytic results are interesting and eligible.

In this paper we use the quite new results concerning the Kuramoto model on
weighted graphs with large number of nodes, obtained by G. Medvedev, H. Chiba
& M. Mizuhara in the series of papers [4, 5, 6, 7, 8]. The cornerstone of
their research is so-called Penrose method [9], which easily allows to capture
bifurcations and transitions between different states, appeared in the model.
We apply the proposed method for the Kuramoto model on weighted graph with
inertia and noise and discuss how the presence of noise affects
synchronization, appearance of phase structures (patterns), and bifurcations
related to transitions between different structures. Our goal is twofold.
First, we want go generalize results obtained with help of Penrose method.
Second, we want to understand how the noise affects synchronization
transition, appearance of chimeras & clusters states, and hysteresis, caused
by inertia term. The main passage of this work is that dissipative, inertial,
and noise terms are naturally interconnected to each other and should be
considered together. Armed with the Penrose method, we can treat all the
actors (inertia, noise, dissipation) simultaneously, easily switch between
different limiting cases (like zero noise limit or vanishing inertia term) and
incorporate graph structure as well.

The rest of paper is organized as follows. We start from section 2 with
discussion of previous research, emphasizing the development of methods and
establish the role of our findings. In section 3 we apply Penrose method for
model with inertia and noise and then provide the explicit simple form of the
expression for critical coupling constant. Developing these ideas, we analyze
the noise impact on different structures appearing in model in section 4.
Next, we reinforce our results by numerical simulation in section 5. We
finalize the narration in section 6.

##  2 Set-up and model properties

We consider the Kuramoto model with inertia and noise. The equations of motion
are given by

$m\ddot{\varphi}_{i}+\gamma\dot{\varphi}_{i}=\omega_{i}+\frac{2K}{N}\sum_{j=1}^{N}W_{ij}\sin\left(\varphi_{j}-\varphi_{i}\right)+\alpha_{i},$
(1)

where $m$ is the mass of $i$-th oscillator, $\gamma$ is dissipative factor,
$\omega_{i}$ is the oscillator eigenfrequency, $K$ is the coupling constant,
$N$ denotes the number of oscillators, $\alpha_{i}$ represents the noise,
which will be specified further. We start from model defined on weighted
symmetric graph with adjacency matrix $W_{ij}$ but mainly we focus on the case
of complete graph. The synchronization is measured by the order parameter,

$h(t)=\frac{1}{N}\sum_{i,j=1}^{N}W_{ij}\exp\left\\{i\varphi_{j}(t)\right\\}.$
(2)

Virtually always the parameters of model (1) are rescaled. Let us carefully
consider the rescaling procedure. First, we note that $[m]=T^{1}$, i.e. mass
term has dimensionality of time. For other parameters, we see that
$[\gamma]=T^{0}$, $[\omega_{i}]=T^{-1}$, $[K]=T^{-1}$, $[D]=T^{-1}$. The
structure of the equations tells us that parameters $m$ and $\gamma$ linked to
each other and moreover if we want to rescale equations, the damping factor
$\gamma$ intertwines to eigenfrequencies and coupling constant. The crucial
fact is that damping factor is dimensionless and we set $\gamma=1$. Let us for
a moment consider model without noise. In this case, the equations of motion
(1) can be represented as system of $2N$ first order ODEs,

$\displaystyle\dot{\varphi}_{i}$ $\displaystyle=v_{i}+\omega_{i},$ (3)
$\displaystyle m\dot{v}_{i}$
$\displaystyle=-v_{i}+\frac{2K}{N}\sum_{j=1}^{N}W_{ij}\sin\left(\varphi_{j}-\varphi_{i}\right).$

If we add noise term, the set of $2N$ ODEs eq. (3) turns into set of $2N$
SDEs,

$\displaystyle d\varphi_{i}$
$\displaystyle=\omega_{i}\,dt+v_{i}\,dt+\alpha_{i}\,dt,$ (4) $\displaystyle
dv_{i}$
$\displaystyle=-\frac{v_{i}}{m}\,dt+\frac{1}{m}\frac{2K}{N}\sum_{j=1}^{N}W_{ij}\sin\left(\varphi_{j}-\varphi_{i}\right)\,dt$

Note that original set-up of Kuramoto model with noise, proposed by Sakaguchi
in [10] contains (as should be) stochastic term only on phases, whereas in
case with inertia the noise usually is related to velocities, i.e.
$\alpha_{i}(t)$ appears in the second line of (4). The limiting cases of the
proposed model was considered in many papers. Let us briefly discuss them for
the sake of completeness. The overdamped model, i.e. with vanishing inertia
term, was considered rigorously by Lancelotti [11], who has proven that
dynamics of such model converges in $N\rightarrow\infty$ limit to Vlasov
equation (in case of zero noise) or to Fokker-Planck equation (in the presence
of noise). First results in stability analysis of the corresponding equations
were obtained by Mirollo and Strogatz [12]. Next, for our best knowledge, the
consideration of model with inertia was started by Tanaka and coauthors in the
paper [13]. In this research the authors have applied ideas developed by Levi
with coauthors in the work [14]. The crucial point of the model with inertia
is the existence of hysteresis, i.e. synchronization and desynchronization
occur at different coupling constants. The underlying mechanism of such
phenomenon is coexistence of stable fixed point and running periodic solution.
Discussing the model with inertia and noise, we should keep in mind the
existence of hysteresis.

Moving onto continuum limit with inertia, the first discussion of the limit
$N\rightarrow\infty$ for the model with noise and inertia was done by Acebron
in [15]. The further development was done by S. Olmi et. al [16] (model with
inertia) and S. Gupta et. al [17] (model with inertia and noise). In parallel
to this research, the rigorous results were obtained by H. Chiba in [18] for
model without inertia and by Dietert in the paper [19], where the Penrose
method for the Kuramoto model was introduced.

As was shown in paper [17], there are several well-known limits of the
proposed model. First, the overdamped limit corresponds to Kuramoto model with
noise (Kuramoto-Sakaguchi model) [10]. Second, overdamped limit with identical
frequencies coincides with so-called Brownian mean-field (BMF) model. Third,
limit with vanishing dissipation and identical eigenfrequencies corresponds to
the Hamiltonian mean field (HMF) model [20]. So, in case of complete graph
(all-to-all couplings) we have three parameters: distribution of
eigenfrequencies, amplitude of noise, and ratio between particle mass and
damping factor.

In this paper we focus on the continuum limit of the proposed model (4). For
the conventional Kuramoto model, the consideration of continuum limit allows
to easily observe low-dimensional dynamics. This dynamics is related to the
invariance under Möbius group transformations and in case of finite $N$ such
invariance is the cornerstone of Watanabe-Strogatz (W-S) ansatz [21, 22, 23].
In case of $N\rightarrow\infty$, the Ott-Antonsen (O-A) ansatz, proposed in
[24], captures this invariance via the reduction of infinite hierarchy of
differential equations for moments of distribution function into the one
equation. Such reduction corresponds to the existence of attracting invariant
manifold in the space of probability measures, so-called O-A manifold. It is
worth mentioning that in case of finite $N$ the O-A manifold can be not
attracting [25]. The O-A ansatz fails in the presence of noise, but the
expansion in terms of cumulants provides the insight for systematic treatment
of non-zero noise case [26]. The cumulant expansion also is one of candidates
to develop the O-A ansatz in the presence of inertia term. Another way to
investigate the case with inertia is so-called adiabatic elimination, which
works for small inertia. The summary of attempts to establish O-A ansatz for
the model with inertia and noise is given in [27].

Discussed generalizations of Kuramoto model alongside with approaches (W-S
ansatz, O-A ansatz, cumulant expansion, etc) are mostly focused on the case of
all-to-all coupling between degrees of freedom. A general theory for more
complicated architectures of couplings is the quite hot topic. Among different
perspectives, the fresh view on this problem was done by G. Medvedev, H. Chiba
and M. Mizuhara in papers [4, 5] and then devleloped in [6, 7, 8]. In these
papers, the authors use the concept of graphons in graph theory followed by
application of Penrose method. Based on these results, we state that the model
(1) dynamics converges in $N\rightarrow\infty$ limit to the following Fokker-
Planck (F-P) equation,

$\frac{\partial f}{\partial
t}+\frac{\partial}{\partial\varphi}\left\\{\left(v+\omega\right)f\right\\}+\frac{\partial}{\partial
v}\left\\{\left(-v+\mathcal{N}[f]\right)f\right\\}-D\frac{\partial^{2}f}{\partial\varphi^{2}}=0,$
(5)

where we set $m=1$ and the quantity $\mathcal{N}[f]$ is given by

$\mathcal{N}[f]=\frac{K}{i}\left\\{h(t,x)e^{-i\varphi}-\overline{h(t,x)}e^{+i\varphi}\right\\},$
(6)

and $h=h(t,x)$ is local continuum version of the order parameter, which is
sensitive to the structure of graph via its graphon limit,

$h(t,x)=\int_{-\pi}^{+\pi}d\varphi\int_{-\infty}^{+\infty}dv\int_{-\infty}^{+\infty}d\omega\int_{0}^{1}dy\,W(x,y)f(t,\varphi,v,\omega,y)e^{i\varphi}.$
(7)

Here $W=W(x,y)$ is the graphon of graph with adjacency matrix $W_{ij}$. The
function $W(x,y)$ is a symmetric function, defined on the unit square
$[0,1]\times[0,1]$ and this function raises the integral operator,

$\bm{W}[f](\bullet)=\int_{0}^{1}dy\,W(\bullet,y)f(y)$ (8)

which acts as $\bm{W}:L^{2}(I)\rightarrow L^{2}(I)$, $I\equiv[0,1]$. This is a
compact self-adjoint operator with a countable sequence of eigenvalues and
with a single accumulation point at zero. The equation (5) plays the key role
in our investigation. The first question that we ask is following: how the
presence of noise affects stability of the incoherent state?

##  3 Stability of mixing state

###  3.1 Linearized Fokker-Planck equation

We call the function $f_{0}=\delta(v)g(\omega)/2\pi$ _mixing state_. It solves
the Fokker-Planck equation (5) and corresponds to incoherent state, i.e. state
with order parameter equals to zero. We analyze stability of this incoherent
state against of small perturbations. To do it, we perform standard
derivations. First, we introduce Fourier expansion of distribution function
$f$,

$u_{j}(t,\zeta,\eta,x)=\int_{-\pi}^{+\pi}d\theta\,\int_{-\infty}^{+\infty}d\omega\,\int_{-\infty}^{+\infty}dv\,e^{i(j\theta+\zeta
v+\eta\omega)}f(t,\theta,v,\omega,x),\quad
j\in\mathbb{Z},\,\zeta\in\mathbb{R},\,\eta\in\mathbb{R}.$ (9)

In terms of Fourier modes, the Fokker-Plank equation (5) becomes

$\frac{\partial u_{j}}{\partial t}=\left(j-\zeta\right)\frac{\partial
u_{j}}{\partial\zeta}+j\frac{\partial
u_{j}}{\partial\eta}+K\zeta\left\\{hu_{j-1}-\overline{h}u_{j+1}\right\\}+Dj^{2}u_{j},$
(10)

Following the paper [7], we introduce the change of variables,

$\zeta-j=\begin{cases}-e^{-\xi_{j}},\quad\zeta-j<0,\\\
+e^{-\xi_{j}},\quad\zeta-j\geq 0.\end{cases}$ (11)

which gives us,

$\frac{\partial u_{j}}{\partial t}=\frac{\partial
u_{j}}{\partial\xi_{j}}+j\frac{\partial
u_{j}}{\partial\eta}+K\left(j-e^{-\xi_{j}}\right)\left\\{hu_{j-1}-\overline{h}u_{j+1}\right\\}+Dj^{2}u_{j},$
(12)

where the functions $u_{j}$ now depend on $\xi_{j}$ (we do not introduce new
variable in order to be clearer). Now we are ready to investigate stability of
the mixing state. To do it, we linearize the eq. (12) near the mixing state.
The linearized version of the Fokker-Planck equation becomes (see [7] for the
detailed derivation)

$\frac{\partial u_{j}}{\partial t}=\frac{\partial
u_{j}}{\partial\xi_{j}}+j\frac{\partial
u_{j}}{\partial\eta}+K\left(j-e^{-\xi_{j}}\right)\left\\{hu_{j-1}-\overline{h}u_{j+1}\right\\}+Dj^{2}u_{j}.$
(13)

Next, we represent the linearized equation as

$\frac{\partial w_{1}}{\partial t}=\frac{\partial
w_{1}}{\partial\xi_{1}}+\frac{\partial
w_{1}}{\partial\eta}+K\left(1-e^{-\xi_{1}}\right)\left\\{h\hat{g}(\eta)+hw_{0}-\overline{h}w_{2}\right\\}+Dw_{1},$
(14)

where $\hat{g}(\eta)$ is the Fourier image of $g(\omega)$ and
$w_{0}=u_{0}-\hat{g}(\eta)$ and $w_{j}=u_{j}$ for $j\neq 0$. In equivalent
form we can write

$\frac{\partial w_{1}}{\partial
t}=\bm{L}_{1}[w_{1}]+K\bm{B}[w_{1}]=\bm{S}[w_{1}],\quad\frac{\partial
w_{j}}{\partial t}=\bm{L}_{j}[w_{j}],\quad j\geq 0,\,j\neq 1,$ (15)

where we have introduced the following operators,

$\begin{gathered}\bm{L}_{j}[w](\xi,\eta,x)=\left(\frac{\partial}{\partial\xi}+j\frac{\partial}{\partial\eta}+Dj^{2}\right)w(\xi,\eta,x),\\\
\bm{B}[w](\xi,\eta,x)=\left(1-e^{-\xi}\right)\hat{g}(\eta)\bm{W}[w(0,0,\bullet)](x).\end{gathered}$
(16)

To proceed further we need to compute the resolvent of the operator
$\bm{L}_{j}$, which is done the Appendix. Having found the resolvent, we can
reformulate linear stability problem as the eigenvalue problem for the linear
operator.

###  3.2 Eigenvalue problem

We consider the eigenvalue problem,

$\lambda w=\bm{S}w=\left(\bm{L}_{1}+K\bm{B}\right)w,$ (17)

where the operators $\bm{L}_{1}$ and $\bm{B}$ are

$\bm{L}_{1}=\left(\partial_{\xi}+\partial_{\eta}-D\right),\quad\bm{B}[w]=\left(1-e^{-\xi}\right)\hat{g}(\eta)\bm{W}[w(0,0,\bullet)](x)$
(18)

and $\hat{g}(\eta)$ is the Fourier image of $g(\omega)$. The eigenvalue
problem can be rewritten as

$\left(\lambda-\bm{L}_{1}\right)w=K\bm{B}w\rightarrow
w=K\left(\lambda-\bm{L}\right)^{-1}\left(1-e^{-\xi}\right)\hat{g}(\eta)\bm{W}[w(0,0,\bullet)].$
(19)

Using the expression for the resolvent, we find

$\left(\lambda-\bm{L}_{1}\right)^{-1}\left(1-e^{-\xi}\right)\hat{g}(\eta)=\int
d\omega\,g(\omega)e^{i\eta\omega}\left(\frac{1}{\lambda-i\omega+D}-\frac{e^{-\xi}}{\lambda+1-i\omega+D}\right).$
(20)

Next, we define the function $\mathcal{D}(\lambda,\xi,\eta)$,

$\mathcal{D}(\lambda,\xi,\eta)=\int
d\omega\,g(\omega)e^{i\eta\omega}\left(\frac{1}{\lambda-i\omega+D}-\frac{e^{-\xi}}{\lambda+1-i\omega+D}\right).$
(21)

With this definition, the eigenvalue problem is written down in more compact
way,

$w=K\mathcal{D}(\lambda,\xi,\eta)\bm{W}[w(0,0,\bullet)].$ (22)

Setting $\xi=\eta=0$, we obtain

$w(0,0,x)=K\mathcal{D}(\lambda)\bm{W}[w(0,0,\bullet)].$ (23)

Let $\mu>0$ is the eigenvalue of the operator $\bm{W}$ and $V$ is the
corresponding eigenfunction. Let $G(\lambda)=\mathcal{D}(\lambda,0,0)$. Then a
root of the equation $G(\lambda)=(K\mu)^{-1}$ determines the eigenvalue
$\lambda=\lambda(\mu)$, i.e. as function of $\mu$. If this holds, then
$w(0,0,x)=V(x)$, which means that the set of roots of equation
$G(\lambda)=(K\mu)^{-1}$ coincides with the eigenvalues of operator $\bm{S}$.
Next, the equation $G(\lambda)=(K\mu)^{-1}$ has a root if the curve
$\mathcal{C}=G(it)$, $t\in\mathbb{R}$ has a positive winding number about
$(K\mu)^{-1}$. The curve $\mathcal{C}$ is called critical curve or Penrose
curve.

We see that the effect of network structure is encoded in the spectrum of
operator $\bm{W}$. Its spectrum is directly related to the spectrum of the
linear operator $\bm{S}$, that is responsible for stability of the mixing
state.

###  3.3 Penrose method in presence of noise

Let us consider the case of complete graph. In such case, the operator
$\bm{W}$ has only one eigenvalue $\mu=1$. We can represent the Penrose curve
as

$G(x+iy)=\int_{-\infty}^{+\infty}d\omega\,g(\omega)\left(\frac{1}{x+iy-i\omega+D}-\frac{1}{x+iy+1-i\omega+D}\right).$
(24)

At this stage, we can easily restore the quantity $m$ by dimensionality
analysis. Applying Sokhotski–Plemelj formula, we arrive to the final equations
for Penrose curve,

$\begin{gathered}x(t)=\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)D}{D^{2}+(\omega-t)^{2}}-m\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)\left(1+mD\right)}{\left(1+mD\right)^{2}+m^{2}\left(\omega-t\right)^{2}},\\\
y(t)=\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)(\omega-t)}{D^{2}+(\omega-t)^{2}}-m^{2}\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)(\omega-t)}{\left(1+mD\right)^{2}+m^{2}\left(\omega-t\right)^{2}}.\end{gathered}$
(25)

We plot Penrose curves for different parameters at fig. 1a.

Figure 1: (a) Penrose curves for Gaussian distribution with $\sigma=0.3$ for
different white noise amplitudes $D$, (b) critical coupling $K_{c}$ as
function of $m$, (c) $K_{c}$ as function of $D$ for different $m$

Analysis of Penrose curve allows to easily find the critical value of coupling
constant $K$. As we have stated above, the spectrum of operator $\bm{S}$
depends on the winding number of Penrose curve about $K^{-1}$. In case of
symmetric distribution $g(\omega)$, the curve $\mathcal{C}$ is symmetric about
the real axis. It intersects the real axis at a unique point $x_{0}>0$.
Furthermore, $G(0)=x_{0}$. It means that if $K<K_{c}$, the point $K^{-1}$ lies
outside the Penrose curve. For $K>K_{c}$, the point lies inside the curve and
the winding number is positive. Therefore, $K_{c}^{-1}=x_{0}$ and
$x_{0}=x(0)$.

Taking into account this fact, we obtain the explicit expression for the
critical coupling constant,

$K_{c}^{-1}=D\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)}{D^{2}+\omega^{2}}-m\left(1+mD\right)\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)}{(1+mD)^{2}+m^{2}\omega^{2}}.$
(26)

This expression is our main result. It can be also presented in the sligthly
different form

$K_{c}^{-1}=\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega
D)}{1+\omega^{2}}-\int_{-\infty}^{+\infty}\frac{d\omega}{1+\omega^{2}}g\left(\omega
D\left(1+\frac{1}{mD}\right)\right)$ (27)

The first term corresponds to the limit of $m\rightarrow 0$ and it is the
well-known critical coupling for the Kuramoto model with noise. The second
term contains both inertia and noise contributions, but its structure is the
same. In the limit $D\rightarrow 0$ the critical coupling for the model with
inertia is reproduced. Slightly more general limit $D\rightarrow 0$,
$m\rightarrow\infty$ with $mD=\mathrm{const}$ can be considered as well. The
dependence of $K_{c}$ on parameter is monotonic both for $m$ and $D$, see
fig.1b and fig.1c

Note that this expression is dissimilar from the critical coupling obtained in
[15]. The reason lies in that the authors consider the different set up with
noise coupled to velocities of particles, instead of phases. In addition, they
deal with another incoherent state, which can depend on velocities.

##  4 Noise impact on model patterns

It is known that the Kuramoto model has a quite rich set of states, distinct
from the incoherent & coherent states. Among them, chimeras and synchronized
clusters attract lot of attention. One of the possible mechanisms to realize
the formation of such structures is to consider the bimodal (two-humped)
frequency distribution $g(\omega)$. In this set-up chimeras and clusters
appear via the unfolding of mixing state bifurcation, which was discussed in
details in [6]. Noise is an additional ingredient, which can affect the
mentioned patterns. In this section we use the formalism of Penrose curves in
order to understand the impact of noise on different structures observed in
the model with inertia.

###  4.1 Bicluster states

To investigate the bicluster state, we consider bimodal Gaussian distribution,

$g(\omega;\mu,\sigma_{1},\sigma_{2})=\frac{1}{2\sqrt{2\pi}}\left[\frac{e^{-(x+\mu)^{2}/(2\sigma_{1}^{2})}}{\sigma_{1}}+\frac{e^{-(x-\mu)^{2}/(2\sigma_{2}^{2})}}{\sigma_{2}}\right]$
(28)

Our reasoning is that for Gaussian distribution we can immediately compare our
results with previously obtained in [8] for zero noise. Here we consider
symmetric bimodal distribution, i.e. $\sigma_{1}=\sigma_{2}=\sigma$.

We fix $\sigma=0.3$ and $m=1.0$. Guided by the previous research, we know that
for large enough value of $\mu$ the bicluster state arises via the Andronov-
Hopf (A-H) bifurcation. In our case the appearance of such state depends on
noise magnitude $D$. To find the critical value $\mu_{c}$, starting from that
A-H bifurcation takes place, we compute the function
$J(\mu)=\left.dy/dt\right|_{t=0}$. The critical value $\mu_{c}$ corresponds to
a root of equation $J(\mu_{c})=0$. In terms of Penrose curve, at large enough
value $\mu_{c}$, the curve develops a cusp.

Figure 2: (a) function $J=J(\mu)$ with $m=1$ for bimodal Gaussian distribution
with $\sigma=0.3$, (b) bifurcation domains in $(\mu,D)$-space

We observe that the noise gradually increases the critical value $\mu_{c}$
(see fig. 2a) and it means that large noise prevents the formation of
bicluster structure. There is a line in $(\mu,D)$-plane that separates the
domain with A-H bifurcation from the domain with first pitchfork bifurcation
(1PF) only, which is shown at fig. 2b.

Figure 3: Disapperance of the second pitchfork bifurcation (parameters are
$m=1.0$, bimodal Guassian distribution with $\mu=1.0$, $\sigma=0.3$ is
considered). Blue dot represents the point $x_{0}$ of the Andronov-Hopf, green
dot represents the second pitchfork bifurcation

In addition to A-H bifurcation, the second pitchfork bifurcation (2PF) occurs
at certain value $\mu=\mu_{0}$. It represents transition from the mixing state
to anti-phase coherent clusters. Noise affects such transition in the expected
way: a large enough value of noise magnitude can shift the corresponding point
$x_{0}$ of Penrose curve from negative semiaxis to the positive one. For the
parameters from fig. 3 this occurs at $D\approx 0.55$. At $D\approx 0.55$ the
point $x_{0}$ crosses the line $\text{Im}\,z=0$, which results disappearance
of the second pitchfork bifurcation. So, we conclude that in the domain
colored by red on fig. 2 the Penrose curves are diffeomorphic to curves shown
at fig.1. In the domain colored by blue the Penrose curves are diffeomorphic
to curve shown at fig. 3 for $D=0.7$, whereas in the green filled domain
curves are diffeomorphic to curve with $D=0.1$ from the same figure.

###  4.2 Chimeras

Figure 4: Destruction of chimera states by increasing of noise magnitude $D$

Now we break the symmetry of bimodal Gaussian distribution by setting
different values of $\sigma_{1}$ and $\sigma_{2}$. This occurs an appearance
of chimera states, i.e. the presence of synchronized and non-synchronized
oscillators in the whole population. Again, noise has a considerable influence
on the chimera states. To illustrate this statement, we choose $m=1$,
$\sigma_{1}=0.4$, $\sigma_{2}=0.2$, and $\mu=1$. Having fixed these
parameters, we vary the noise magnitude $D$ from zero to non-zero values. The
critical curves are shown at fig. 4. We see that the first bifurcation (red
dot) occurs at quite small value of coupling constant and this corresponds to
the partial synchronization of population but a sizeable fraction of
oscillators remains non-synchronized. The second bifurcation (blue dot) occurs
at a larger value of coupling constant and this corresponds to the
disappearance of chimera state: now the whole population is synchronized.

Noise makes the synchronization scenario more interesting. For small values of
noise the qualitative picture does not change, whereas from a certain value of
noise _one more_ transition occurs (middle picture on panel, see fig. 4).
Further increase of noise magnitude destroys the chimera state. The intuition
based on Penrose curves treatment is confirmed by direct numerical
simulations, which are given at fig. 5, fig. 6, and fig 7.

Figure 5: Destruction of chimera state by increasing noise magnitude $D$
($N=2000$ oscillators, $K=1.0$, bimodal Gaussian distribution with parameters
$\mu=1.0$, $\sigma_{1}=0.4$, $\sigma_{2}=0.2$, oscillators are grouped by
intrinsic frequency sign)

The fig. 5 represents how noise destroys the chimera state. For quite small
values of noise magnitude, $D=0.1$, the chimera state survives, but the slight
increase ($D=0.2$) in noise causes the decay of synchronized fraction of
population (middle picture). Further increase of noise magnitude (to $D=0.3$)
completely destructs synchronized remnants. The reverse story is depicted on
fig. 6. Here we fix the moderate value of noise magnitude, $D=0.4$, and
observe the stability of mixing state for different values of coupling
constant $K$. From the analysis of Penrose curves, we find two critical
values, $K_{c}^{1}\approx 1.60$ (bifurcation from mixing to chimera) and
$K_{c}^{2}\approx 2.48$ (from chimera to bicluster). Noise suppresses
synchronization, so clusters seem friable in compare with noiseless case.

Figure 6: Bifurcations of mixed state in presence of noise for increasing
coupling constant $K$ ($N=2000$ oscillators, $D=0.4$, bimodal Gaussian
distribution with parameters $\mu=1.0$, $\sigma_{1}=0.4$, $\sigma_{2}=0.2$,
oscillators are grouped by intrinsic frequency sign) Figure 7: Bifurcations of
mixed state in presence of noise for increasing coupling constant $K$
($N=2000$ oscillators, $D=0.6$, bimodal Gaussian distribution with parameters
$\mu=1.0$, $\sigma_{1}=0.4$, $\sigma_{2}=0.2$, oscillators are grouped by
intrinsic frequency sign)

As we have mentioned above, for the moderate value of noise, there exist three
critical points, which are determined as intersections of Penrose curve with
real axis. The corresponding transitions are shown at fig. 7. The smallest
critical coupling (for our choice of parameters it is $K_{c}^{1}\approx 2.57$)
corresponds to the bifurcation of mixed state to chimera. However, the
moderate noise significantly suppresses synchronization (the synchronized
population has order parameter $\approx 0.30$, whereas non-synchronized has
$\approx 0.05$). The next critical coupling ($K_{c}^{2}\approx 3.62$)
corresponds again to the bifurcation from chimera state to bicluster state.
The largest critical coupling corresponds to a merging of two clusters into a
fully coherent state.

##  5 Numerical simulations

We have emphasized that the model with inertia term exhibits hysteresis. This
fact should be taken into account in numerical simulations. The details of
simulations can be found in the Appendix. The phase diagrams were obtained
with unimodal Cauchy distribution,

$g(\omega;\mu,\sigma)=\frac{\sigma}{\pi(\sigma^{2}+(\omega-\mu)^{2})},$ (29)

and other parameters are specified in figure captions. First of all, we
reproduce previous results, obtained for the model with inertia only. The key
point is to see hysteresis from the phase diagram. We would like also to draw
attention to the convergence of finite $N$ simulations to analytical results,
derived in continuum limit, $N\rightarrow\infty$. Having adopted the numerical
scheme, described in [28], we observe the same issues, pointed out in [16]:
the backward critical coupling (which describes the destruction of coherent
state under the decrease of coupling constant) is almost insensitive to finite
size effects, whereas the forward critical coupling (which corresponds to the
appearance of the synchronized state by the increase of coupling constant) is
very sensitive to finite size effects. Second, we simulate the model with
inertia and noise. Again, we use the mentioned above scheme. For our surprise,
the finite size effects are not so significant in the presence of noise. We
observe that noise strongly suppresses the hysteresis, which is similar to the
effect, observed in [17].

Figure 8: Left: phase diagram ($N=40000$) for unimodal Cauchy with
$\sigma=1.0$, $m=1.0$ and zero noise, right: phase diagram ($N=3500$) with
$\sigma=0.5$, $m=1.0$ and $D=0.4$. On both plots red line corresponds to
critical coupling from (26), $\langle h\rangle$ is time-averaged order
parameter.

##  6 Discussion and conclusion

In this study we consider the Kuramoto model with noise and inertia for the
case of complete graph (all-to-all couplings). We analyze the stability of
mixing state against small fluctuations via linearization of the Fokker-Planck
equation. The linearization procedure yields the eigenvalue problem for the
particular differential operator. With help of recently obtained results, we
have adopted the Penrose method to investigate the spectrum of such operator
and have classified possible bifurcations of the mixing state. The formulation
in terms of Penrose (critical) curves allows to determine explicitly the
critical coupling for transition to the synchronized state. The presence of
noise affects on critical points corresponding to these bifurcations. It can
also affect the existence of bifurcations, suppressing the transitions from
mixing state to clusters (Andronov-Hopf bifurcation) and mixing state to anti-
phase coherent clusters (second pitchfork bifurcation). Also, we observe that
large noise destroys chimera states.

Numerical experiments clearly show the hysteresis in the absence of noise. We
have faced the known challenge with the convergence of numerics to the
analytically computed critical coupling constant and have observed that noise
strongly affects the hysteresis. We emphasize that developed in [8] method
allows to capture only the so-called ‘‘forward’’ critical coupling, which
corresponds to transition from non-synchronized state to the synchronized one.
Nevertheless, we assume that the Penrose method can provide a profound insight
into the hysteresis phenomenon in Kuramoto model with inertia.

There are several issues for the further study. First, we have not covered the
analysis of graph architecture, but it can be done in straightforward way via
the graphons formalism. Second, it is known that hysteresis behavior can be
related to the graph architecture via frequency-degree correlated
eigenfrequencies. The simplest possible example is the case of star graph,
which was investigated analytically in [28]. The analysis of synchronization
on star graph with noise [29] shows that there is a tricritcal point, where
the synchronization transition becomes continuous. For our best understanding,
it is not so clear how to realize the frequency-degree correlations in terms
of graphon-based treatment. Finally, there is one more method that allows to
indicate critical point, based on the Fisher information metric in parameter
space [30]. The critical points appear as singularities of this metric,
defined through the stationary solutions of Fokker-Planck equations. In case
of noise and inertia, the corresponding parameter space is three dimensional
and it would be interesting to probe the hysteresis via singularities.

The usage of graphons in accompany with the Penrose method seems a very
powerful tool to analyze the effects of long-range interactions in the
Hamiltonian mean-field (HMF) model with the additional dissipative term. The
inclusion of dissipation looks natural and in the limit of vanishing
dissipation all the effects known for original HMF model are reproduced. We
expect that formulation in terms of graphons will allow to tune the strength
of long-range interactions, starting from the complete graph to more sparsed
graphs. Introducing the dissipative term into HMF model, we immediately relate
this problem to the results for the networks of Josephson contacts and systems
of granulated superconductors.

The Kuramoto model is the popular simplified model for the synchronization of
neurons, see [31] for the review. The account of noise in the neuroscience
problem is standard however the inertia term has not been discussed in that
context so far. Our study suggests that combination of inertia and noise
strongly influence at least two phenomena certainly relevant for the brain
synchronization – formation and fate of chimeras [32] and interplay between
the continuous and abrupt synchronization corresponding to the second and
first order phase transitions. The latter is important for the transition from
of consciousness and unconsciousness [33].

Last point that we would like to mention is the development of neuromorphic
architectures based on spintronic devices. The spin-torque oscillators (STOs)
are prominent building blocks for neuromorphic computing. For instance, quite
recently the set of synchronized STOs was used for vovel recognition [34]. The
dynamics of STO magnetization can be treated in terms of
Landau–Lifshitz–Gilbert (LLG) equation but typically this dynamics corresponds
to relatively slow magnetization change. In case of fast dynamics, memory
effects should be taken into account. It can be done by introducing the
inertial term into LLG [35]. The LLG equation is directly related to the
dynamics of Landau-Stuart oscillators and under some assumptions the
synchronization of STOs can be treated in terms of Kuramoto model.

## Acknowledgements

AA would like to express his huge gratitude to Georgi Medvedev for fruitful
discussions and advice. Also, AA is grateful to Denis Goldobin for sharing the
relevant notes and discussion of O-A ansatz for the Kuramoto model with
inertia. This work is supported by the IDEAS Research Center. A.A. was
supported by the Foundation for the Advancement of Theoretical Physics and
Mathematics ‘‘BASIS’’ (grant №23-1-5-41-1).

## Data availability

The code snippets and the datasets generated and analysed during the current
study are available from the corresponding author on reasonable request.

##  Appendix A Numerical simulations details

All the numerical simulations were realized on ȷulıa programming language. For
numerical solution of differential equations, we use the package
`DifferentialEquations`. For the simulations with non-zero noise, the Euler-
Maryama was used. In case of zero noise, we have compared two implementations,
based on the system of 2-nd order ODEs (via `SecondOrderODEProblem`) and based
on the system of 1-st order ODEs (via `ODEProblem`). We observed that there is
no significant difference in performance. For the computation of phase
diagrams with non-zero inertia and all-to-all couplings, the following
simulation scheme was used:

  1. 1.

Generate initial conditions, i.e. set of $\varphi_{i}(t=0)$. For mixing state,
initial phases are drawn from uniform distribution on $(-\pi,+\pi]$

  2. 2.

Numerically solve set of $N$ differential equations for a fixed value of
coupling constant $K_{0}$ on the interval $[0,T]$ with predefined $T$

  3. 3.

Compute time-averaged value of order parameter:

    1. (a)

Compute sum $h(t)=\sum_{i=j}^{N}\exp\\{i\varphi_{j}(t)\\}$

    2. (b)

Compute numerically the integral $\displaystyle\langle
h\rangle=\frac{1}{T}\int_{0}^{T}dt\,|h(t)|$

    3. (c)

Store $\langle h\rangle$ in memory

  4. 4.

Shift coupling constant by a predefined small step $\delta K$,
$K_{0}\rightarrow K_{0}+\delta K$

  5. 5.

Numerically solve set of $N$ differential equations for $K_{0}+\delta K$ with
initial conditions $\varphi_{i}(T)$ (i.e. on this step we use the stored in
memory values $\varphi_{i}(T)$ from step 2) on the interval $[0,T]$ and
refresh values $\varphi_{i}(t=T)$ in memory

  6. 6.

Repeat steps 3, 4, 5 until the current value of coupling constant will reach
predefined value $K_{\max}$ (i.e. we split interval $[0,K_{\max}]$ with
predefined step $\delta K$)

We observe the quite good convergence between analytical predictions in
thermodynamic $N\rightarrow\infty$ limit and results of numerical simulations
starting from $N\sim 2000$. Surprisingly, the presence of noise improves
convergence. If not specified explicitly, we use $T=1000$ and time interval is
discretized with step $\delta T=0.05$. The choice of $K_{\max}$ is dictated by
analytical expressions for critical couplings and we set $\delta K=0.01$.

For the simulation of patterns, we find (analytically, if possible or
numerically) corresponding critical values of coupling constant and choose the
value $K^{*}$, which lies in a interval between critical values (for instance,
$K_{c}^{1}<K^{*}<K_{c}^{2}$ in case of two critical coupling constants,
predicted by Penrose curve). We generate initial conditions once and then use
them. Next, we numerically solve set of $2N$ SDE and extract values of phases
$\varphi_{i}(t=T)$, $i\in\\{1,N\\}$. Finally, we group oscillators in
accordance with the sign of their intrinsic frequencies (this is done to
improve visualization, see [8]). The number of oscillators in simulations
related to patterns is $N=2000$.

##  Appendix B Operators and their spectra

###  B.1 Resolvent

For the operator $\bm{L}_{j}$, we have the eigenvalue problem,

$\left(\lambda-\bm{L}_{j}\right)u=u.$ (30)

Now apply Fourier transform $\mathcal{F}[\bullet]$ for both sides,

$\left(\lambda-i\zeta-
ij\omega+Dj^{2}\right)\mathcal{F}[u]=\mathcal{F}[v],\,F[u](\zeta,\omega,x)=\int
d\xi\,d\eta\,e^{-i\xi\zeta-i\eta\omega}u(\xi,\eta,x).$ (31)

It gives us

$u=\left(\lambda-\bm{L}_{j}\right)^{-1}u=\mathcal{F}^{-1}\mathcal{F}[u]$

. The explicit computation of the integral is

$\int d\zeta\,d\omega\,\frac{e^{+i\xi\zeta+i\omega\eta}}{\lambda-i\zeta-
ij\omega+Dj^{2}}\int
d\xi^{\prime}\,d\eta^{\prime}\,e^{-i\zeta\xi^{\prime}-i\eta^{\prime}\omega}u(\xi^{\prime},\eta^{\prime},x)=\\\
=\int_{0}^{\infty}dt\,e^{-\lambda t-Dj^{2}t}\int
d\zeta\,d\omega\,e^{+i\xi\zeta+i\omega\eta}e^{+i\zeta t+ij\omega t}\int
d\xi^{\prime}\,d\eta^{\prime}\,e^{-i\zeta\xi^{\prime}-i\eta^{\prime}\omega}u(\xi^{\prime},\eta^{\prime},x)=\\\
=\int_{0}^{\infty}dt\,e^{-\lambda t-Dj^{2}t}\int d\xi^{\prime}\int
d\eta^{\prime}\left\\{\int
d\omega\,e^{i\omega(jt+\eta-\eta^{\prime})}\right\\}\left\\{\int
d\zeta\,e^{i\zeta(\xi+t-\xi^{\prime})}\right\\}u(\xi^{\prime},\eta^{\prime},x)=\\\
=\int_{0}^{\infty}dt\,e^{-\lambda t-Dj^{2}t}u(\xi+t,\eta+jt).$ (32)

It means that the resolvent of the operator $\bm{L}_{j}$ is given by

$\left(\lambda-\bm{L}_{j}\right)^{-1}u=\begin{cases}\displaystyle\int_{0}^{\infty}dt\,e^{-\lambda
t-Dj^{2}t}u(\xi+t,\eta+jt),\quad\operatorname{\mathrm{Re}\,}\lambda>0,\\\
\displaystyle-\int_{-\infty}^{0}dt\,e^{-\lambda
t-Dj^{2}t}u(\xi+t,\eta+jt),\quad\operatorname{\mathrm{Re}\,}\lambda<0.\end{cases}$
(33)

###  B.2 Eigenvalue problem

The eigenvalue problem is given by

$\left(\lambda-\bm{L}_{1}\right)u=K\bm{B}u\rightarrow
u=K\left(\lambda-\bm{L}\right)^{-1}\left(1-e^{-\xi}\right)\hat{g}(\eta)\bm{W}[u(0,0,\bullet)].$
(34)

Using the expression for resolvent, we find

$\left(\lambda-\bm{L}_{1}\right)^{-1}\left(1-e^{-\xi}\right)\hat{g}(\eta)=\int
d\zeta\,d\omega\frac{e^{i\xi\zeta+i\eta\omega}}{\lambda-i(\zeta+\omega)+D}\mathcal{F}\left[\left(1-e^{-\xi}\right)\hat{g}(\eta)\right]=\\\
=\int
d\zeta\,d\omega\,\frac{e^{i\zeta\xi+i\eta\omega}}{\lambda-i(\zeta+\omega)+D}\int
d\xi^{\prime}d\eta^{\prime}e^{-i\eta^{\prime}\omega-i\xi^{\prime}\zeta}\left(1-e^{-\xi^{\prime}}\right)\int
d\omega^{\prime}e^{+i\eta^{\prime}\omega^{\prime}}g(\omega^{\prime})=\\\
=\int\frac{d\zeta\,d\omega\,e^{i\zeta\xi+i\eta\omega}}{\lambda-i(\zeta+\omega)+D}\int
d\omega^{\prime}g(\omega^{\prime})\int d\xi^{\prime}\left\\{\int
d\eta^{\prime}\,e^{i\eta^{\prime}(\omega^{\prime}-\omega)}\right\\}e^{-i\xi^{\prime}\zeta}\left(1-e^{-\xi^{\prime}}\right)=\\\
=\int
d\omega\,g(\omega)e^{i\eta\omega}\left(\frac{1}{\lambda-i\omega+D}-\frac{e^{-\xi}}{\lambda+1-i\omega+D}\right).$
(35)

##  Appendix C Derivation of Penrose curve equations

Representing $\lambda=x+iy$, let us consider limit $x\rightarrow 0$. We denote
the first integral in the definition of $G$ as $G_{1}$ and the second as
$G_{2}$. So, we write $G=G_{1}-G_{2}$, where

$G_{1}=\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)}{x+iy-i\omega+D},\quad
G_{2}=\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)}{x+iy+1-i\omega+D}.$
(36)

We start from $G_{1}$ and denote $G_{1}=iI_{1}$, where

$I_{1}=\int\frac{d\omega\,g(\omega)}{ix-y+\omega+iD}=-i\int\frac{d\omega\,g(\omega)(x+D)}{(x+D)^{2}+(\omega-y)^{2}}+\int\frac{d\omega\,g(\omega)(\omega-y)}{(x+D)^{2}+(\omega-y)^{2}},$
(37)

and take the limit of $x\rightarrow 0$,

$\lim\limits_{x\rightarrow 0}G_{1}=i\lim\limits_{x\rightarrow
0}I_{1}=+D\int\frac{d\omega\,g(\omega)}{D^{2}+(\omega-y)^{2}}+i\int\frac{d\omega\,g(\omega)(\omega-y)}{D^{2}+(\omega-y)^{2}}.$
(38)

For $G_{2}$, we denote $G_{2}=iI_{2}$, where

$I_{2}=\int\frac{d\omega\,g(\omega)}{ix-y+i+\omega+iD}=-i\int\frac{d\omega\,g(\omega)(x+D+1)}{(x+D+1)^{2}+(\omega-y)^{2}}+\int\frac{d\omega\,g(\omega)(\omega-y)}{(x+D+1)^{2}+(\omega-y)^{2}},$
(39)

which gives us in the limit of $x\rightarrow 0$,

$\lim\limits_{x\rightarrow 0}G_{2}=i\lim\limits_{x\rightarrow
0}I_{2}=\int\frac{d\omega\,g(\omega)(D+1)}{(D+1)^{2}+(\omega-y)^{2}}+i\int\frac{d\omega\,g(\omega)(\omega-y)}{(D+1)^{2}+(\omega-y)^{2}}.$
(40)

Finally, for $G$ in the limit of $x\rightarrow 0$, we have

$G(it)=\lim\limits_{x\rightarrow
0}\left(G_{1}-G_{2}\right)=+D\int\frac{d\omega\,g(\omega)}{D^{2}+(\omega-t)^{2}}-D\int_{-\infty}^{+\infty}\frac{d\omega\,g(\omega)}{(D+1)^{2}+(\omega-t)^{2}}-\\\
-\int\frac{d\omega\,g(\omega)}{(D+1)^{2}+(\omega-t)^{2}}+i\int\frac{d\omega\,g(\omega)(\omega-t)}{D^{2}+(\omega-t)^{2}}-i\int\frac{d\omega\,g(\omega)(\omega-t)}{(D+1)^{2}+(\omega-t)^{2}}.$
(41)

It defines the curve in complex plane $(x,y)$, where $t$ is the parameter.

## References

  * [1] Arkady Pikovsky, Michael Rosenblum, and Jürgen Kurths.  Synchronization: A Universal Concept in Nonlinear Sciences.  Cambridge University Press, 2001.  doi:10.1017/CBO9780511755743.
  * [2] Yoshiki Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In International symposium on mathematical problems in theoretical physics, pages 420–422. Springer, 1975.  doi:10.1007/BFb00133.
  * [3] Yoshiki Kuramoto.  Chemical oscillations, waves, and turbulence.  Springer, 1984.  doi:10.1007/978-3-642-69689-3.
  * [4] Hayato Chiba and Georgi S. Medvedev.  The mean field analysis of the Kuramoto model on graphs I. the mean field equation and transition point formulas.  Discrete and Continuous Dynamical Systems-Series A, 39(1):131–155, 2019.  doi:10.3934/dcds.2019006.
  * [5] Hayato Chiba and Georgi S. Medvedev.  The mean field analysis of the Kuramoto model on graphs II. asymptotic stability of the incoherent state, center manifold reduction, and bifurcations.  Discrete and Continuous Dynamical Systems, 39(7):3897–3921, 2019.
  * [6] Georgi S. Medvedev and Matthew S. Mizuhara.  Chimeras unfolded.  Journal of Statistical Physics, 186(3):46, 2022.  doi:10.1007/s10955-022-02881-y.
  * [7] Hayato Chiba and Georgi S. Medvedev.  Stability and bifurcation of mixing in the Kuramoto model with inertia.  SIAM Journal on Mathematical Analysis, 54(2):1797–1819, 2022.  doi:10.1137/21M1427000.
  * [8] Hayato Chiba, Georgi S. Medvedev, and Matthew S. Mizuhara.  Bifurcations and patterns in the Kuramoto model with inertia.  Journal of Nonlinear Science, 33(5):78, 2023.  doi:10.1007/s00332-023-09931-z.
  * [9] Oliver Penrose.  Electrostatic instabilities of a uniform non-Maxwellian plasma.  The Physics of Fluids, 3(2):258–265, 1960.  doi:10.1063/1.1706024.
  * [10] Hidetsugu Sakaguchi.  Cooperative phenomena in coupled oscillator systems under external fields.  Progress of theoretical physics, 79(1):39–46, 1988.  doi:10.1143/PTP.79.39.
  * [11] Carlo Lancellotti.  On the Vlasov limit for systems of nonlinearly coupled oscillators without noise.  Transport theory and statistical physics, 34(7):523–535, 2005.  doi:10.1080/00411450508951152.
  * [12] Steven H. Strogatz and Renato E. Mirollo.  Stability of incoherence in a population of coupled oscillators.  Journal of Statistical Physics, 63(3):613–635, 1991.  doi:10.1007/BF01029202.
  * [13] Hisa-Aki Tanaka, Allan J. Lichtenberg, and Shin’ichi Oishi.  First order phase transition resulting from finite inertia in coupled oscillator systems.  Physical Review Letters, 78(11):2104, 1997.  doi:10.1103/PhysRevLett.78.2104.
  * [14] M. Levi, Frank C. Hoppensteadt, and W. L. Miranker.  Dynamics of the Josephson junction.  Quarterly of Applied Mathematics, 36(2):167–198, 1978.  doi:10.1090/qam/484023.
  * [15] J. A. Acebrón, L. L. Bonilla, and R Spigler.  Synchronization in populations of globally coupled oscillators with inertial effects.  Physical Review E, 62(3):3437, 2000.  doi:10.1103/PhysRevE.62.3437.
  * [16] Simona Olmi, Adrian Navas, Stefano Boccaletti, and Alessandro Torcini.  Hysteretic transitions in the Kuramoto model with inertia.  Physical Review E, 90(4):042905, 2014.  doi:10.1103/PhysRevE.90.042905.
  * [17] Shamik Gupta, Alessandro Campa, and Stefano Ruffo.  Nonequilibrium first-order phase transition in coupled oscillator systems with inertia and noise.  Physical Review E, 89(2):022123, 2014.  doi:10.1103/PhysRevE.89.022123.
  * [18] Hayato Chiba.  A proof of the kuramoto conjecture for a bifurcation structure of the infinite-dimensional kuramoto model.  Ergodic Theory and Dynamical Systems, 35(3):762–834, 2015.  doi:10.1017/etds.2013.68.
  * [19] Helge Dietert.  Stability and bifurcation for the Kuramoto model.  Journal de Mathématiques Pures et Appliquées, 105(4):451–489, 2016.  doi:10.1016/j.matpur.2015.11.001.
  * [20] Thierry Dauxois, Stefano Ruffo, Ennio Arimondo, and Martin Wilkens.  Dynamics and Thermodynamics of Systems with Long Range Interactions.  Springer, 2002.  doi:10.1007/3-540-45835-2_1.
  * [21] Shinya Watanabe and Steven H Strogatz.  Integrability of a globally coupled oscillator array.  Physical Review Letters, 70(16):2391, 1993.  doi:10.1103/PhysRevLett.70.2391.
  * [22] Shinya Watanabe and Steven H Strogatz.  Constants of motion for superconducting Josephson arrays.  Physica D: Nonlinear Phenomena, 74(3-4):197–253, 1994.  doi:10.1016/0167-2789(94)90196-1.
  * [23] Seth A. Marvel, Renato E. Mirollo, and Steven H. Strogatz.  Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 19(4):043104, 2009.  doi:10.1063/1.3247089.
  * [24] Edward Ott and Thomas M. Antonsen.  Low dimensional behavior of large systems of globally coupled oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 18(3):037113, 2008.  doi:10.1063/1.2930766.
  * [25] Jan R. Engelbrecht and Renato Mirollo.  Is the Ott-Antonsen manifold attracting?  Physical Review Research, 2(2):023057, 2020.  doi:10.1103/PhysRevResearch.2.023057.
  * [26] Denis S. Goldobin and Anastasiya V. Dolmatova.  Ott-Antonsen ansatz truncation of a circular cumulant series.  Physical Review Research, 1(3):033139, 2019.  doi:10.1103/PhysRevResearch.1.033139.
  * [27] Evelina V Permyakova, Lyudmila S. Klimenko, Irina V. Tyulkina, and Denis S. Goldobin.  Four approaches for description of stochastic systems with small and finite inertia.  In Journal of Physics: Conference Series, page 012050. IOP Publishing, 2021.  doi:10.1088/1742-6596/1945/1/012050.
  * [28] Vladimir Vlasov, Yong Zou, and Tiago Pereira.  Explosive synchronization is discontinuous.  Physical Review E, 92(1):012904, 2015.  doi:10.1103/PhysRevE.92.012904.
  * [29] Artem Alexandrov.  Synchronization on star graph with noise.  Chaos, Solitons & Fractals, 167:113056, 2023.  doi:10.1016/j.chaos.2022.113056.
  * [30] Artem Alexandrov and Alexander Gorsky.  Information geometry and synchronization phase transition in the Kuramoto model.  Physical Review E, 107(4):044211, 2023.  doi:10.1103/PhysRevE.107.044211.
  * [31] Michael Breakspear, Stewart Heitmann, and Andreas Daffertshofer.  Generative models of cortical oscillations: neurobiological implications of the Kuramoto model.  Frontiers in human neuroscience, 4:190, 2010.  doi:10.3389/fnhum.2010.00190.
  * [32] Siyu Huo and Zonghua Liu.  Condensation of eigenmodes in functional brain network and its correlation to chimera state.  Communications Physics, 6(1):285, 2023.  doi:10.1038/s42005-023-01405-8.
  * [33] Hyoungkyu Kim, Joon-Young Moon, George A. Mashour, and UnCheol Lee.  Mechanisms of hysteresis in human brain networks during transitions of consciousness and unconsciousness: Theoretical principles and empirical evidence.  PLoS computational biology, 14(8):e1006424, 2018.  doi:10.1371/journal.pcbi.1006424.
  * [34] Miguel Romera, Philippe Talatchian, Sumito Tsunegi, Flavio Abreu Araujo, Vincent Cros, Paolo Bortolotti, Juan Trastoy, Kay Yakushiji, Akio Fukushima, Hitoshi Kubota, et al.  Vowel recognition with four coupled spin-torque nano-oscillators.  Nature, 563(7730):230–234, 2018.  doi:10.1038/s41586-018-0632-y.
  * [35] Ritwik Mondal, Levente Rózsa, Michael Farle, Peter M. Oppeneer, Ulrich Nowak, and Mikhail Cherkasskii.  Inertial effects in ultrafast spin dynamics.  Journal of Magnetism and Magnetic Materials, page 170830, 2023.  doi:10.1016/j.jmmm.2023.170830.