††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Collective dynamics of heterogeneously and nonlinearly coupled phase
oscillators

Can Xu  Institute of Systems Science and College of Information Science and
Engineering, Huaqiao University, Xiamen 361021, China  Xiaohuan Tang  School
of Physics and Electronic Engineering, Jiangsu Normal University, Xuzhou
221116, China  Huaping Lü  School of Physics and Electronic Engineering,
Jiangsu Normal University, Xuzhou 221116, China  Karin Alfaro-Bittner
Unmanned Systems Research Institute, Northwestern Polytechnical University,
Xi’an 710072, China  Departamento de Física, Universidad Técnica Federico
Santa María, Av. España 1680, Casilla 110V, Valparaíso, Chile  Stefano
Boccaletti  Unmanned Systems Research Institute, Northwestern Polytechnical
University, Xi’an 710072, China  Universidad Rey Juan Carlos, Calle Tulipán
s/n, 28933 Móstoles, Madrid, Spain  Moscow Institute of Physics and
Technology, 9 Institutskiy per., Dolgoprudny, 141701 Moscow, Russia  CNR -
Institute of Complex Systems, Via Madonna del Piano 10, I-50019 Sesto
Fiorentino, Italy  Matjaž Perc  Faculty of Natural Sciences and Mathematics,
University of Maribor, Koroška cesta 160, 2000 Maribor, Slovenia  Department
of Medical Research, China Medical University Hospital, China Medical
University, Taichung, Taiwan  Alma Mater Europaea ECM, Slovenska ulica 17,
2000 Maribor, Slovenia  Complexity Science Hub Vienna, Josefstädterstraße 39,
1080 Vienna, Austria  Shuguang Guan  Department of Physics, East China Normal
University, Shanghai 200241, China

###### Abstract

Coupled oscillators have been used to study synchronization in a wide range of
social, biological, and physical systems, including pedestrian-induced bridge
resonances, coordinated lighting up of firefly swarms, and enhanced output
peak intensity in synchronizing laser arrays. Here we advance this subject by
studying a variant of the Kuramoto model, where the coupling between the phase
oscillators is heterogeneous and nonlinear. In particular, the quenched
disorder in the coupling strength and the intrinsic frequencies are
correlated, and the coupling itself depends on the amplitude of the mean-field
of the system. We show that the interplay of these factors leads to a
fascinatingly rich collective dynamics, including explosive synchronization
transitions, hybrid transitions with hysteresis absence, abrupt irreversible
desynchronization transitions, and tiered phase transitions with or without a
vanishing onset. We develop an analytical treatment that enables us to
determine the observed equilibrium states of the system, as well as to explore
their asymptotic stability at various levels. Our research thus provides
theoretical foundations for a number of self-organized phenomena that may be
responsible for the emergence of collective rhythms in complex systems.

##  I Introduction

The Kuramoto model is a well-known and widely applicable paradigm for
theoretically describing and modeling collective behavior in large ensembles
of interacting units kuramoto1975in . Since its conception in 1975, the model
has been used to explore synchronization dynamics in numerous social,
biological, and physical systems, with examples ranging from arrays of
Josephson junctions marvel2009invar , flashing fireflies er1991an , cardiac
pacemaker cells osaka2017mod , to economic markets may2008comp . The model has
been particularly popular in physics and mathematics, because it often affords
analytical insights towards better understanding different types of collective
behavior and synchronization transitions that occur in systems consisting of a
large number of coupled nonlinear self-sustained oscillators strogatz2000from
; acebron2005the ; pikovsky2015dynamics ; Boccaletti2018 .

In its original version, the Kuramoto model elucidates synchronization at the
onset of a phase transition that emerges due to the interplay of the intrinsic
frequencies of individual oscillators and the global coupling among them.
Specifically, the heterogeneity of the natural frequencies of oscillators,
which are distributed randomly across the population according to a prescribed
probability density, tends to desynchronize the system, while the attracting
global coupling between the oscillators opposes this disorder and tends
towards synchrony. These two competing forces result in the rich dynamics of
the Kuramoto model that has been explored in great detail during the past
decades xu2020universal .

In addition to the heterogeneity in natural frequencies, recent work in
physics and network science has highlighted the importance of heterogeneous
coupling between dynamical units. Namely, the coupling strength between
individual oscillators need not be uniform, but rather oscillator-dependent
hong2011kuramoto ; latsenko2013stationary . Such heterogeneities can be
encoded by means of quenched random interactions, which often represent a
realistic setup for a wide range of systems and applications that go beyond
the traditional identical coupling daido1992 ; daido2000 ; stiller2000 ;
daido2005 . For example, such heterogeneous interactions have been used to
shed light on the collective dynamics of excitatory-inhibitory neurons, or to
describe opinion formation amongst conformists and contrarians in sociology
cb2003 ; ms2005 ; 2008 . In terms of theoretical results obtained with
heterogeneous coupling between phase oscillators, previous works have reported
the emergence of traveling waves, standing waves, $\pi$-states, glass states,
cluster states, as well as chimera states hong2011con ; hong2012mean ;
kloumann2104phase ; peron2021 ; jafari2021chimera .

Almost in parallel to coherent phases observed in network dynamics that emerge
from heterogeneous coupling, explosive synchronization transitions have also
received a lot of attention pazo2005 ; sb2016 ; souza2019 ; bick2021 . The
latter represent first-order-like phase transitions that switch abruptly
between incoherent and coherent states, and they are typically due to
correlations in network topology and oscillator dynamics gomze2011explosive ;
xu2015explosive ; xu2019syn . Since abrupt transitions between two stable
states underlie many real-life systems, the phenomenon has received ample
attention, for example to better understand bistable perception in the brain
dynamics wang2013brain  and the robustness and optimal functioning of large-
scale power grids skardal2015control . More recently, it has also been shown
that explosive synchronization transitions can be due to higher-order
interactions, where a link connects more than two oscillators, or due to
nonlinear coupling, where, for example, the coupling strength between
oscillators is proportional to the coherence of the system tt2011mul ;
komarov2015finite ; pikovsky2019low ; skardal2019abrupt ; millan2020exp .
Research has also shown that these two effects can actually be considered
equivalent xu2020bifurcation2020 ; xu2021spectrum ; xu2021stability .

But despite the wealth of findings and insights related to the Kuramoto model
and its variants, an unexplored question is what novel effects could different
coupling patterns have on the macroscopic dynamics of the system. Particularly
so if both heterogeneity and nonlinearity are present in the coupling among
phase oscillators. In this case, indeed little is known about the possible
collective behaviors and the phase transitions between them that may be
observed as a result. In this paper, we therefore consider a new variant of
the Kuramoto model, which incorporates heterogeneity and nonlinearity in the
coupling. More specifically, we consider correlations between the quenched
disorder in the coupling strength and the intrinsic frequencies, and we also
consider the coupling to be dependent on the amplitude of the mean-filed of
the population. As we will show, the interplay between these effects
significantly shapes the overall collective dynamics of the system, giving
rise to a number of fascinating phenomena of relevance for the emergence of
synchronization. Importantly, the studied model is analytically tractable, so
that our research provides theoretical foundations, and in fact a natural
mechanism, for the spontaneous emergence of various rhythmic states and the
bifurcations associated with them in complex systems, and it does so with
relatively small changes to system parameters.

The remainder of this paper is organized as follows. In Sec. II, we present
the new variant of the Kuramoto model, and we establish the self-consistency
approach to determine the long-term macroscopic dynamics. In Sec. III, we
carry out a linear stability analysis for the fully locked, two-cluster,
states in a special one-dimensional manifold, and we set up a necessary
stability criterion for their occurrence. In Sec. IV, we extend this stability
result to all perturbed directions in the phase space, and we explore the
characteristic dynamics of the system due to perturbations. In Sec. V, we
investigate the asymptotic stability of the steady states in the infinite
limit, and we obtain the associated eigenspectrum by using the Ott-Antonsen
reduction. In Sec. VI, we investigate the transitions by means of which the
system can move between different coherent states. Various types of phase
transition towards synchronization are identified, including the explosive
synchronization transition, the abrupt irreversible desynchronization
transition, and the tiered synchronization transition with or without a
vanishing onset. In Sec. VII, we summarize our main conclusions and briefly
discuss their possible wider implications.

##  II Mathematical model and its stationary solutions

As noted, we consider an extension of the Kuramoto model, whose governing
equations are

$\dot{\theta}_{i}=\omega_{i}+\frac{K_{i}f(R)}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}).$
(1)

Here $\theta_{i}$ is the phase of $i-$th oscillator with $i=1,...,N$.
$\\{\omega_{i}\\}$ are the natural frequencies selected from a probability
density function $g(\omega)$, which is assumed to be symmetric throughout the
paper, i.e., the mean frequency is set to zero and $g(\omega)=g(-\omega)$.

Eq. (1) differs from the classical Kuramoto model in that the uniform coupling
strength has been replaced by the heterogeneous coupling $K_{i}$ with a
feedback factor $f(R)$, which is a generic function of the amplitude of the
Kuramoto order parameter defined by

$Z(t)=R(t)e^{i\Theta(t)}=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}}.$ (2)

$R(t)$ measures the level of coherence of the system and $\Theta(t)$ denotes
the average phase of the population. $K_{i}$ and $f(R)$ account for the
inhomogeneity and nonlinearity of the coupling, respectively. Without loss of
generality, we set $f(R)>0$ (which can be made by rescaling the time) and
$K_{i}=K_{N-i}$. As we will show below, the long-term dynamics of this model
is significantly richer than the dynamics of the traditional Kuramoto mode
with linear homogeneous coupling.

We emphasize that the inherent nonlinearity, quenched random interactions, and
the large number of degrees of freedom make it difficult to understand the
quantitative dynamics of Eq. (1) with arbitrary choices of $K_{i}$ and $f(R)$.
We therefore consider only the case where the coupling is symmetric around
index $i$ of each oscillator. Aside from this constraint, however, we are free
to choose $K_{i}$ and $f(R)$ arbitrarily.

We first focus on a particular case, in which the inhomogeneities of natural
frequencies and the couplings are chosen deterministically rather than
randomly. Namely, we set $K_{i}=K\left|\omega_{i}\right|$, and
$f(R)=R^{\beta-1}$ with $K,\beta>0$. In this setting, the randomness is
intrinsic to the oscillators themselves rather than to the coupling between
them wang2011syn ; zhang2013exp ; bi2016coe ; hong2016pha ; xu2018origin ;
xu2016dyn ; xu2019universal ; xu2019bifur ; yuan2104dyn . Moreover, $\beta\geq
1$ ($\beta<1$) sets up a positive (negative) feedback between the coupling and
the coherence of the system rose2007self ; gf2007generalized ; zhang2015exp ;
hung2018 ; zou2020dyn ; gao2020reduction ; dai2020dis .

Before proceeding the analysis, we make a few comments of the model. As
mentioned above, the heterogeneous interactions have been used to shed light
on the intrinsic dynamics of glass states, or to describe opinion formation
amongst conformists and contrarians in sociology. On the other hand, the
natural frequencies, as a typical heterogeneity of the population, are the
intrinsic properties of oscillators themselves. Therefore, it is reasonable to
establish correlations between the coupling strength and natural frequencies
latsenko2013stationary ; iatsenko2014 ; hong2016pha . In all the quoted
studies, the coupling strength depends on the natural frequencies implicitly,
and the aim is to uncover the traveling waves and glass states in coupled
oscillator populations.

In our model, instead, the coupling strength depends on the natural
frequencies explicitly, i.e., we set $K_{i}=K|\omega_{i}|$ corresponding to a
frequency-weighted coupling (FWC). The FWC was proposed by Ref. wang2011syn ,
and then it was generalized in Ref. zhang2013exp ; bi2016coe . In all those
studies, it was shown that a FWC is equivalent to the frequency-degree
correlations in complex networks reported in Ref. gomze2011explosive . Both
models turn out to give rise to Explosive synchronization phenomena note3 .

Regarding the nonlinear coupling (NC), where $f(R)$ is a nonlinear function
with respect to $R$. Relevant examples include Ref. rose2007self ; bai2010 ,
where NC can be achieved by means of phase reduction in a system of Stuart-
Landau oscillators. For convenience, one may set $f(R)=R^{\beta-1}$ (a power-
law function), which serves as a typical NC. Similar discussions can be found
in Ref. gf2007generalized ; zou2020dyn . At the same time, recent works in
physics and neuroscience highlight the potential importance of higher-order
(or group) interactions, e.g., three-way or more, connections that may be
organized via higher-order simplexes or a simplicial complex. The effect of
such interactions on the dynamics thus represents an important topic of
research in the complex system community. More importantly, recent research
has also shown that the higher-order interactions are equivalent to NC
xu2021stability ; note4 .

The definition of the order parameter Eq. (2) allows us to rewrite Eq. (1) as

$\dot{\theta}_{i}=\omega_{i}+K|\omega_{i}|R^{\beta}\sin(\Theta-\theta_{i}).$
(3)

By doing so, Eq. (3) retains the mean-field character, and
$K|\omega_{i}|R^{\beta}$ displays an effective force acting on each phase
oscillator $\theta_{i}$.

Next, we briefly outline the self-consistency approach to analyze the
stationary dynamics of Eq. (1). The key idea for the self-consistent method is
to assume that, in the long-term evolution, the macroscopic quantity $R(t)$ is
a constant, and $\Theta(t)$ rotates uniformly, i.e., $\Theta(t)=\Omega
t+\Theta_{0}$. Indeed, we may set $\Omega=\Theta_{0}=0$ by entering into an
appropriate rotating frame and shifting initial conditions 2013pre ; note .
Then, Eq. (3) can be reformulated as

$\dot{\theta}_{i}=\omega_{i}-q|\omega_{i}|\sin\theta_{i},$ (4)

where the auxiliary parameter $q=KR^{\beta}$ is defined to simplify notation.

We can see from Eq. (4) that each individual oscillator exhibits two types of
long-term behavior depending on the magnitude of $q$. Specifically, if $q>1$,
all oscillators become phase-locked with

$\sin\theta_{i}=\frac{\omega_{i}}{q|\omega_{i}|},$ (5)

and

$\cos\theta_{i}=\sqrt{1-q^{-2}}.$ (6)

When $q<1$, all oscillators are drifting, i.e., they rotate nonuniformly on
the unit circle and form a stationary distribution given by

$\rho_{d}(\theta,\omega)=\frac{\rm{sgn}(\omega)\sqrt{\omega^{2}(1-q^{2})}}{2\pi(\omega-q|\omega|\sin\theta)}.$
(7)

The sign function $\rm sgn(\omega)$ appearing in Eq. (7) is needed to ensure
the positivity of the distribution.

In contrast to the traditional Kuramoto model and its previously studied
variants, we find that in our case the coexistence of the phase-locked and the
drifting populations can never occur. In other words, the system remains
either in the fully phase-locked or the totally drifting state for a fixed
value of $q$. The splitting of the system into the fully phase-locked or the
completely drifting state allows us to obtain an expression for the order
parameter $R$ self-consistently

$Z=R=\left\langle e^{i\theta}\right\rangle_{q<1},$ (8)

and

$Z=R=\left\langle e^{i\theta}\right\rangle_{q>1},$ (9)

where $\left\langle\cdot\right\rangle$ denotes the average over the
population.

Remarkably, the symmetry of the system implies that
$\left\langle\sin\theta\right\rangle_{q>1}=\left\langle
e^{i\theta}\right\rangle_{q<1}=0$. Therefore, for $q<1$, the only solution to
Eq. (8) is $R=q=0$, which corresponds to the incoherent state with the phases
scattered uniformly around the unit circle. Conveniently, for $q>1$, the
expression for $R$ becomes simply

$R=\cos\theta_{i},$ (10)

or

$\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}=F(q)=q^{-\frac{1}{\beta}}\sqrt{1-q^{-2}}.$
(11)

The derivations presented above reveal three important results. In the first
place, the presence of correlations between frequency and coupling results in
a universal phase transition towards synchronization, whose equilibrium is
independent of the frequency distribution and system size. Secondly, there are
no solutions to Eq. (11) for a sufficiently small value of $K$, since $F(q)$
remains finite for $q\in(1,+\infty)$. And thirdly, there exists a critical
coupling $K_{c}$, beyond which the two coherent state solutions exist located
at $q<q_{c}$ and $q>q_{c}$, respectively. The critical parameter $q_{c}$
satisfies the condition $F^{\prime}(q_{c})=0$, i.e.,

$q_{c}=\sqrt{\beta+1},$ (12)

which in turn determines the critical coupling $K_{c}$ and the critical order
parameter $R_{c}$, yielding

$K_{c}=\frac{(\beta+1)^{\frac{\beta+1}{2}}}{\beta^{\frac{\beta}{2}}},$ (13)

and

$R_{c}=\sqrt{\frac{\beta}{\beta+1}}.$ (14)

In what follows, we will prove that the solutions for $R>R_{c}$ and $q>q_{c}$
are attractive, and moreover, correspond exactly to the solutions along the
$F^{\prime}(q)<0$ branch.

##  III Stability in the one-dimensional manifold

We have shown by means of Eqs. (5) and (6) that, depending on the sign of the
natural frequencies and depending on the intensity of $q$, the oscillator
population splits into two populations in the long run at opposite locations.
These observations lead us to consider the following question: What is the
stability property of the steady states described by Eq. (11) in this
manifold? As we will show, this consideration is quite convenient, because the
low-dimensional description of the dynamics in Eq. (1) allows us to deduce the
associated critical points and gain some intuition about the stability of
phase-locked states in the presently studied modified Kuramoto model.

In the fully locked scenario, the system acts like two giant oscillators,
labeled as $\theta_{\pm}$, which correspond to natural frequencies
$\pm\omega$, respectively. In such a manifold, Eq. (1) reduces to

$\dot{\theta}_{+}=\omega+\frac{K|\omega|}{2}R^{\beta-1}\sin(\theta_{-}-\theta_{+})$
(15)

and

$\dot{\theta}_{-}=-\omega+\frac{K|\omega|}{2}R^{\beta-1}\sin(\theta_{+}-\theta_{-}).$
(16)

Without loss of generality, we assume that $\omega>0$, thus having
$\sin\theta_{\pm}=\pm q^{-1}$ and $\cos\theta_{\pm}=\sqrt{1-q^{-2}}$, and the
order parameter becoming

$R=\frac{1}{2}\Big{(}e^{i\theta_{+}}+e^{i\theta_{-}}\Big{)}.$ (17)

Since $\theta_{+}+\theta_{-}=0$, Eqs. (15) and (16) become one-dimensional if
we introduce the phase difference $\phi=\theta_{+}-\theta_{-}=2\theta_{+}$,
which then yields

$\dot{\phi}=\omega\big{[}2-K\Phi(\phi)\big{]},$ (18)

with

$\Phi(\phi)=\cos^{\beta-1}\frac{\phi}{2}\sin\phi.$ (19)

In this one-dimensional manifold, the dynamics is relatively simple to
analyze. The steady state solution of Eq. (18) requires that

$\frac{2}{K}=\Phi(\phi),\;\phi\in(0,\pi),$ (20)

which implies that the coupling strength $K\geq K_{c}$, such that

$\frac{2}{K}\leq\max\big{[}\Phi(\phi)\big{]}.$ (21)

We can then show that

$K_{c}=\frac{2}{\max\big{[}\Phi(\phi)\big{]}}=\frac{2}{\Phi(\phi_{0})}=\frac{(\beta+1)^{\frac{\beta+1}{2}}}{\beta^{\frac{\beta}{2}}},$
(22)

and that the critical order parameter $R_{c}$ is given by

$R_{c}=\cos\frac{\phi_{0}}{2}=\sqrt{\frac{\beta}{\beta+1}},$ (23)

which are exactly the expressions obtained previously in Eqs. (13) and (14).

In order to consider the stability of the solutions beyond $K_{c}$, we impose
a small perturbation $\delta\phi$ away from $\phi$ in Eq. (20) and neglect
higher order terms of $\delta\phi$. In doing so, we obtain

$\delta\dot{\phi}=-\omega K\Phi^{\prime}(\phi)\delta\phi,$ (24)

from which it follows that the characteristic value of the perturbation in
this direction is proportional to $-\Phi^{\prime}(\phi)$. This means further
that the stable (unstable) condition for $\phi$ is equivalent to
$\Phi^{\prime}(\phi)>0$ [$\Phi^{\prime}(\phi)<0$]. To corroborate this
conclusion, we note that

$\Phi^{\prime}(\phi)=\frac{1}{2}\cos^{\beta-1}\frac{\phi}{2}\Big{(}\cos^{2}\frac{\phi}{2}-\beta\sin^{2}\frac{\phi}{2}\Big{)},$
(25)

which can be simplified significantly by taking into account Eqs. (17) and
(20), thus obtaining

$\Phi^{\prime}(\phi)=\frac{1}{2}R^{\beta-1}\big{[}R^{2}(\beta+1)-\beta\big{]}.$
(26)

From this we find that $R>R_{c}$ leads to $\Phi^{\prime}(\phi)>0$, thereby
determining a stable branch of the solution, while $R<R_{c}$ leads to
$\Phi^{\prime}(\phi)<0$, corresponding to an unstable solution.

Based on the above treatment, we conclude that the steady state solution along
the branch $q>q_{c}$ (or $F^{\prime}(q)<0$) is stable at least in this special
perturbed direction. However, the stationary solution along the branch
$q<q_{c}$ (or $F^{\prime}(q)>0$) is unstable in all perturbed directions,
because it has a positive characteristic value of the perturbation even in the
synchronized manifold.

In the subsequent two sections, we will show that these stability criteria can
be extended to any perturbation direction, and this regardless of whether $N$
is finite or infinite.

##  IV Stability in the $N$-dimensional manifold

It is worth pointing out that what is unusual about the currently studied
system is that the dynamics of Eq. (4) naturally entrains the oscillators into
two groups with equal numbers and with opposite positions. To better
understand the nature of these fixed points, we restrict our attention to
their asymptotic stability in the $N$-dimensional manifold. We will show that
the solutions obtained with $R>R_{c}$ or $q>q_{c}$ are indeed stable in all
the perturbed directions, while the solutions obtained with $R<R_{c}$ or
$q<q_{c}$ are not.

For the locked states $(\theta_{1},\theta_{2}...\theta_{N})$, let
$\theta_{i}(t)=\theta_{i}+x_{i}(t)$ and $\mathbf{x}=(x_{1},x_{2}...x_{N})$ be
the small perturbations. The linearization of the dynamics described by Eq.
(1) is given by the following equation mirollo2005the

$\dot{\mathbf{x}}=\mathbf{J}\mathbf{x},$ (27)

where $\mathbf{J}$ is the Jacobian matrix, whose entries are

$\begin{split}J_{ij}=&\frac{\partial\dot{\theta}_{i}}{\partial\theta_{j}}=\frac{\partial}{\partial\theta_{j}}\bigg{[}\omega_{i}+\frac{K|\omega_{i}|}{N}R^{\beta-1}\sum_{k=1}^{N}\sin(\theta_{k}-\theta_{i})\bigg{]}\\\
=&\frac{\partial
R^{\beta-1}}{\partial\theta_{j}}\frac{K|\omega_{i}|}{N}\sum_{k=1}^{N}\sin(\theta_{k}-\theta_{i})\\\
&+\frac{K|\omega_{i}|}{N}R^{\beta-1}\sum_{k=1}^{N}\frac{\partial\sin(\theta_{k}-\theta_{i})}{\partial\theta_{j}}\\\
=&|\omega_{i}|\bigg{[}\frac{\beta
q}{NR}\sin\theta_{i}\sin\theta_{j}+\frac{q}{NR}\cos\theta_{i}\cos\theta_{j}-q\cos\theta_{i}\delta_{ij}\bigg{]}.\end{split}$
(28)

We note that the rotational invariance of the dynamics described by Eq. (1)
implies that $\mathbf{J}$ always has a trivial eigenvalue $\lambda=0$
corresponding to an eigenvector $(1,1...1)$. To verify that this indeed holds,
we point out that the rows of $\mathbf{J}$ all sum to $0$, i.e.,

$\sum_{j=1}^{N}J_{ij}=|\omega_{i}|q\cos\theta_{i}\bigg{(}\frac{\sum_{j=1}^{N}\cos\theta_{j}}{NR}-\sum_{j=1}^{N}\delta_{ij}\bigg{)}=0,$
(29)

where we have taken into account $\sum_{j=1}^{N}\sin\theta_{j}=0$ and
$\sum_{j=1}^{N}\cos\theta_{j}=NR$. Aside from $\lambda=0$, the stable
condition for the locked states requires that the remaining $N-1$ eigenvalues
are all negative.

Before a more in-depth exploration of the eigenvalues of $\mathbf{J}$, we here
provide an intuitive argument. We first express

$J_{ij}=|\omega_{i}|\bigg{(}\frac{\beta
q}{NR}s_{i}s_{j}+\frac{q}{NR}c_{i}c_{j}-qc_{i}\delta_{ij}\bigg{)},$ (30)

where $\mathbf{c}=(c_{1},c_{2}...c_{N})$ with $c_{i}=\cos\theta_{i}$ and
$\mathbf{s}=(s_{1},s_{2}...s_{N})$ with $s_{i}=\sin\theta_{i}$. Hence, the
stability condition is equivalent to $\mathbf{Jxx}<0$ for all
$\mathbf{x}\in\mathbb{R}^{N}$ with $\|\mathbf{x}\|\neq 0$, and $\mathbf{x}$ is
orthogonal to $(1,1...1)$. That is, $\sum_{i=1}^{N}x_{i}=0$, and we thus have

$\mathbf{Jxx}=\frac{\beta
q}{NR}\sum_{i=1}^{N}|\omega_{i}|s_{i}x_{i}\sum_{j=1}^{N}s_{j}x_{j}-q\sum_{i=1}^{N}|\omega_{i}|c_{i}x_{i}^{2}.$
(31)

Following Ref. mirollo2005the , it becomes apparent that $\mathbf{J}$ has at
most one positive eigenvalue. If there are two or more positive eigenvalues
for $\mathbf{J}$, the span of the corresponding eigenvectors defines at least
a two-dimensional subspace, on which $\mathbf{Jxx}>0$. On the other hand, one
can find a nonzero vector $\mathbf{x}$ on a space that is orthogonal to
$\mathbf{s}$, but from the identity above, we have $\mathbf{Jxx}<0$, which
thus leads to a contradiction.

To explore the eigenvalues of $\mathbf{J}$ more conveniently and in further
detail, we introduce a change of coordinates in order to partially diagonalize
the system. In particular, if $\mathbf{y}=(y_{1},y_{2}...y_{N})$ with
$y_{i}=\sqrt{|\omega_{i}|c_{i}}x_{i}$, we obtain

$\mathbf{Jxx}=\frac{\beta q}{NR}\mathbf{(u\cdot y)(v\cdot
y)}-q\|\mathbf{y}\|^{2}.$ (32)

The vectors $\mathbf{u}=(u_{1},u_{2}...u_{N})$ with
$u_{i}=\sqrt{\frac{|\omega_{i}|}{c_{i}}}s_{i}$, and
$\mathbf{v}=(v_{1},v_{2}...v_{N})$ with
$v_{i}=\frac{s_{i}}{\sqrt{|\omega_{i}|c_{i}}}$, are introduced to ease
notation, and their product is

$\mathbf{u\cdot
v}=\sum_{i=1}^{N}\frac{s_{i}^{2}}{c_{i}}=N\frac{q^{-2}}{\sqrt{1-q^{-2}}}.$
(33)

In this new coordinate system, we can express

$\mathbf{y}=a\mathbf{\frac{u}{\|u\|}}+b\mathbf{\frac{v}{\|v\|}}+\mathbf{y}^{\perp},$
(34)

with $a$ and $b$ being the projections to the unit vector
$\mathbf{\frac{u}{\|u\|}}$ and $\mathbf{\frac{v}{\|v\|}}$, respectively.
$\mathbf{y}^{\perp}$ is the remaining vector, such that
$\mathbf{y^{\perp}\cdot u=y^{\perp}\cdot v}=0$. We thus have

$\begin{split}\mathbf{Jxx}=&\big{(}a^{2}+b^{2}\big{)}\bigg{(}\frac{\beta
q^{-1}}{R^{2}}-q\bigg{)}-q\|\mathbf{y}^{\perp}\|^{2}\\\
&+ab\bigg{[}\frac{\beta q}{NR}\|\mathbf{u}\|\|\mathbf{v}\|+\bigg{(}\frac{\beta
q^{-3}}{R^{3}}-\frac{2q^{-1}}{R}\bigg{)}\frac{N}{\|\mathbf{u}\|\|\mathbf{v}\|}\bigg{]}.\end{split}$
(35)

In particular, for $b=0$, the term around $ab$ vanishes, and the identity
above becomes a quadratic form. Since the coefficient $a$ could be arbitrarily
small, the stability condition $\mathbf{Jxx}<0$ requires that

$\beta<R^{2}q^{2},$ (36)

which again recovers the necessary stability condition obtained in Eqs. (22)
and (23).

To proceed with extracting more details about the eigenvalues of $\mathbf{J}$,
we now compute its characteristic polynomial. We therefore express
$\mathbf{J}$ in the matrix form

$\mathbf{J=W(C-D)},$ (37)

where $\mathbf{W}$ and $\mathbf{D}$ are the diagonal matrices with the entries
$|\omega_{i}|$ and $qc_{i}$ along the diagonal, respectively. The matrix
$\mathbf{C}$ is given by

$C_{ij}=\frac{\beta q}{NR}s_{i}s_{j}+\frac{q}{NR}c_{i}c_{j}.$ (38)

The characteristic polynomial is

$\begin{split}P(\lambda)&=\rm det(\mathbf{J}-\lambda\mathbf{I})\\\ &=\rm
det\big{[}\mathbf{W(C-D)}-\lambda\mathbf{I}\big{]}\\\ &=\rm
det\big{[}\mathbf{W}(\mathbf{C-D-W}^{-1}\lambda)\big{]}\\\ &=\rm
det\Big{\\{}\mathbf{W}(\mathbf{D+W^{-1}}\lambda)\big{[}(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}-\mathbf{I}\big{]}\Big{\\}}\\\
&=\prod_{i=1}^{N}|\omega_{i}|\Big{(}qc_{i}+\frac{\lambda}{|\omega_{i}|}\Big{)}\rm
det\big{[}(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}-\mathbf{I}\big{]}.\end{split}$
(39)

The key task from here on is to determine the second term above. We recall
that the rank of $\mathbf{C}$ is only two, and $\rm
det\big{[}(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\big{]}\neq 0$. The same is
true for the matrix $(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}$.
Inspired by this property, we choose an orthogonal basis
$\big{(}\mathbf{\frac{c}{\|c\|},\frac{s}{\|s\|}}\big{)}$ and extend it into
$\mathbb{R}^{N}$. The remaining $N-2$ vectors span the kernel of
$(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}$. Therefore, the matrix
$(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}$ if restricted to the
orthogonal basis $(\mathbf{\frac{c}{\|c\|},\frac{s}{\|s\|}})$ reduces to a
$2\times 2$ matrix given by

$(\mathbf{D}+\lambda\mathbf{W}^{-1})^{-1}\mathbf{C}=\left(\begin{array}[]{cccc}Q_{11}(\lambda)&Q_{12}(\lambda)\\\\[8.0pt]
Q_{21}(\lambda)&Q_{22}(\lambda)\end{array}\right),$ (40)

where its entries are defined as

$\begin{split}Q_{11}(\lambda)&=\mathbf{\frac{c}{\|c\|}\cdot(D+\lambda
W^{-1})^{-1}C\cdot\frac{c}{\|c\|}}\\\
&=\frac{q}{NR}\sum_{i=1}^{N}\frac{|\omega_{i}|c_{i}^{2}}{\lambda+|\omega_{i}|qc_{i}},\end{split}$
(41) $\begin{split}Q_{12}(\lambda)&=\mathbf{\frac{c}{\|c\|}\cdot(D+\lambda
W^{-1})^{-1}C\cdot\frac{s}{\|s\|}}\\\ &=\frac{\beta
q^{-1}}{R\mathbf{\|c\|\|s\|}}\sum_{i=1}^{N}\frac{|\omega_{i}|s_{i}c_{i}}{\lambda+|\omega_{i}|qc_{i}},\end{split}$
(42) $\begin{split}Q_{21}(\lambda)&=\mathbf{\frac{s}{\|s\|}\cdot(D+\lambda
W^{-1})^{-1}C\cdot\frac{c}{\|c\|}}\\\
&=\frac{qR}{\mathbf{\|c\|\|s\|}}\sum_{i=1}^{N}\frac{|\omega_{i}|s_{i}c_{i}}{\lambda+|\omega_{i}|qc_{i}},\end{split}$
(43) $\begin{split}Q_{22}(\lambda)&=\mathbf{\frac{s}{\|s\|}\cdot(D+\lambda
W^{-1})^{-1}C\cdot\frac{s}{\|s\|}}\\\ &=\frac{\beta
q}{NR}\sum_{i=1}^{N}\frac{|\omega_{i}|s_{i}^{2}}{\lambda+|\omega_{i}|qc_{i}}.\end{split}$
(44)

We note further that the terms $Q_{12}(\lambda)$ and $Q_{21}(\lambda)$ are in
fact all zero, because the natural frequencies are grouped into plus-minus
pairs, i.e., $\omega_{i}=-\omega_{N-i}$. This symmetrical case significantly
simplifies the characteristic polynomial into

$P(\lambda)=\prod_{i=1}^{N}|\omega_{i}|\Big{(}qc_{i}+\frac{\lambda}{|\omega_{i}|}\Big{)}\big{[}1-Q_{11}(\lambda)\big{]}\big{[}1-Q_{22}(\lambda)\big{]}.$
(45)

The eigenvalues $\lambda$ corresponding to $P(\lambda)=0$ are therefore
determined by

$\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}=H_{c}(\lambda)=\frac{q^{\frac{\beta-1}{\beta}}}{N}\sum_{i=1}^{N}\frac{|\omega_{i}|(1-q^{-2})}{\lambda+|\omega_{i}|\sqrt{q^{2}-1}},$
(46)

and

$\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}=H_{s}(\lambda)=\frac{\beta
q^{\frac{\beta-1}{\beta}}}{N}\sum_{i=1}^{N}\frac{|\omega_{i}|q^{-2}}{\lambda+|\omega_{i}|\sqrt{q^{2}-1}}.$
(47)

For convenience, we order the natural frequencies so that
$|\omega_{1}|\leq|\omega_{2}|\leq...\leq|\omega_{N}|$.

The obtained eigenvalues merit three important comments. In the first place,
we note that $\lambda=0$ is a trivial solution of Eq. (46) because
$K^{-\frac{1}{\beta}}=H_{c}(0)$ is precisely the self-consistent equation
described by Eq. (11). As already mentioned, this property stems from the
rotational symmetry of the dynamical equation. Secondly, if
$\lambda<-|\omega_{N}|\sqrt{q^{2}-1}$, both functions $H_{c}(\lambda)$ and
$H_{s}(\lambda)$ are negative. This implies that the characteristic Eqs.
(46)-(47) have no roots in the region
$\lambda\in(-\infty,-|\omega_{N}|\sqrt{q^{2}-1})$. And lastly,
$\lim_{\lambda\xrightarrow{0^{\mp}}-|\omega_{i}|\sqrt{q^{2}-1}}H_{c(s)}(\lambda)=\pm\infty$,
which reflects that, away from their poles, the functions $H_{c}(\lambda)$ and
$H_{s}(\lambda)$ are strictly decreasing. Accordingly, the characteristic
equations must each have exactly one root between any two consecutive poles
$(-|\omega_{i+1}|\sqrt{q^{2}-1},-|\omega_{i}|\sqrt{q^{2}-1})$.

Since functions $H_{c}(\lambda)$ and $H_{s}(\lambda)$ are decreasing for
$\lambda\in(-|\omega_{1}|\sqrt{q^{2}-1},+\infty)$, and
$\lim_{\lambda\to+\infty}H_{c(s)}(\lambda)=0$, the remaining one eigenvalue is
completely determined by the difference $\Delta$ defined by

$\Delta=H_{s}(0)-\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}.$ (48)

Specifically, $\Delta<0$ implies that there is a negative root to Eq. (47) in
the region $\lambda\in(-|\omega_{1}|\sqrt{q^{2}-1},0)$. On the other hand, for
$\Delta>0$, a positive root to Eq. (47) exists in the region
$\lambda\in(0,+\infty)$. By combining this with the self-consistency equation,
we obtain the elegant identity

$\Delta=\beta qF^{\prime}(q),$ (49)

which directly relates the stability criterion for the locked states to the
shape of the self-consistent equation $F(q)$. Namely, if $F^{\prime}(q)>0$ the
steady state solutions are unstable, whereas if $F^{\prime}(q)<0$ the steady
state solutions are attractive and stable, which completes the proof.

##  V Stability in the infinite limit

Next, to better understand the dynamics that emerges from the above
considerations, we study the system in the thermodynamical limit. In
particular, we investigate two aspects of the asymptotic stability of the
equilibrium states, for which we reduce the system to a low-dimensional
equation that governs the long-term macroscopic dynamics. Under the limit
$N\to\infty$, the dynamical system is equivalent to the continuity equation
for the probability density function $\rho(\theta,\omega,t)$

$\frac{\partial\rho}{\partial t}+\frac{\partial(\rho v)}{\partial\theta}=0.$
(50)

Here the velocity $v(\theta,\omega,t)$ is given by

$v(\theta,\omega,t)=\omega+\frac{K|\omega|}{2i}R^{\beta-1}\big{(}Ze^{-i\theta}-\bar{Z}e^{i\theta}\big{)},$
(51)

and the complex order parameter $Z(t)$ becomes

$Z(t)=\int_{-\infty}^{+\infty}\int_{0}^{2\pi}e^{i\theta}\rho(\theta,\omega,t)g(\omega)d\omega
d\theta.$ (52)

By applying the Ott-Antonsen ansatz otto2008low ; otto2009long , the
probability density function can be expressed as

$\rho(\theta,\omega,t)=\frac{g(\omega)}{2\pi}\bigg{[}1+\sum_{n=1}^{\infty}\bar{\alpha}^{n}e^{in\theta}+\sum_{n=1}^{\infty}\alpha^{n}e^{-in\theta}\bigg{]},$
(53)

where the variable $\alpha(\omega,t)$ is a complex valued function with
$|\alpha(\omega,t)|\leq 1$, and the bar denotes the complex conjugate. This
consideration underlines the Poisson kernel that is dynamically invariant, as
long as $\alpha(\omega,t)$ evolves according to

$\frac{d\alpha}{dt}=i\omega\alpha+\frac{K|\omega|R^{\beta-1}}{2}\big{(}Z-\bar{Z}\alpha^{2}\big{)}.$
(54)

Meanwhile, $Z(t)$ reduces to

$Z(t)=\int_{-\infty}^{+\infty}g(\omega)\alpha(\omega,t)d\omega=\hat{g}\alpha(\omega,t),$
(55)

where we have defined the operator $\hat{g}$ to denote the integral with
respect to $\omega$. In this setting, the quantity $\alpha(\omega,t)$ can be
interpreted as the local order parameter that is formed by the oscillators
with frequencies near $\omega$.

With respect to Eq. (54), we are interested in the steady states
$\dot{\alpha}=0$, which yields two-composed solutions

$\alpha_{0}(\omega)=\left\\{\begin{array}[]{rl}i\rm
sgn(\omega)\frac{1-\sqrt{1-q^{2}}}{q},&q<1,\\\\[8.0pt]
\sqrt{1-q^{-2}}+i\frac{\omega}{q|\omega|},&q>1.\end{array}\right.$ (56)

We note that the first branch of Eq. (56) corresponds to the drifting
population, whereas the second branch corresponds to the phase-locked case.
Other possible steady state solutions of $\alpha_{0}(\omega)$ are ruled out by
the constraints $|\alpha_{0}(\omega)|\leq 1$ and $R\geq 0$ omel2013bifurcation
.

To perform the linear stability analysis of $\alpha_{0}(\omega)$, we introduce
a small perturbation to it. Namely, let

$\alpha(\omega,t)=\alpha_{0}(\omega)+\epsilon\eta(\omega,t),$ (57)

with $\epsilon$ being the perturbation magnitude and $\eta(\omega,t)$ the
perturbation function. The order parameter $Z(t)$ under perturbation becomes

$Z(t)=R+\epsilon\hat{g}\eta(\omega,t).$ (58)

Taking into account
$R^{\beta-1}=Z^{\frac{\beta-1}{2}}\cdot\bar{Z}^{\frac{\beta-1}{2}}$,
substituting Eqs. (57) and (58) into Eq. (54), and neglecting the higher order
terms of $\epsilon$, we get the linear evolution equation for the perturbation

$\frac{d\eta(\omega,t)}{dt}=c(\omega)\eta(\omega,t)+a(\omega)\hat{g}\eta(\omega,t)+b(\omega)\hat{g}\bar{\eta}(\omega,t),$
(59)

where

$c(\omega)=\left\\{\begin{array}[]{rl}i\omega\sqrt{1-q^{2}},&q<1,\\\\[8.0pt]
-|\omega|\sqrt{q^{2}-1},&q>1,\end{array}\right.$ (60)

and

$a(\omega)=\frac{|\omega|q}{2R}\bigg{(}\frac{\beta+1}{2}-\alpha_{0}^{2}\frac{\beta-1}{2}\bigg{)},$
(61)
$b(\omega)=\frac{|\omega|q}{2R}\bigg{(}\frac{\beta-1}{2}-\alpha_{0}^{2}\frac{\beta+1}{2}\bigg{)}.$
(62)

Since $\eta(\omega,t)$ is a complex valued function, we need to switch to a
system of two equations containing $\bar{\eta}(\omega,t)$. Thus, if
$\mathbf{V}=(\eta,\bar{\eta})$, the $\mathbf{V}$ evolves according to

$\frac{d\mathbf{V}}{dt}=\mathbf{MV+P}\hat{g}\mathbf{V},$ (63)

where the matrix $\mathbf{M}$ is a multiplication operator

$\mathbf{M}=\left(\begin{array}[]{cccc}c(\omega)&0\\\\[8.0pt]
0&\bar{c}(\omega)\end{array}\right),$ (64)

and $\mathbf{P}$ is the coefficient matrix given by

$\mathbf{P}=\left(\begin{array}[]{cccc}a(\omega)&b(\omega)\\\\[8.0pt]
\bar{b}(\omega)&\bar{a}(\omega)\end{array}\right).$ (65)

Figure 1: Phase diagrams of the time-averaged order parameter $\left\langle
R\right\rangle$ versus the global coupling strength $K$, as obtained for
different values of the exponent $\beta$. The panels show results for
$\beta=1$ (left column), $\beta<1$ (middle column), and $\beta>1$ (right
column). For (a1)-(a3) we have used $\beta=1$ and
$g(\omega)=\frac{\Delta}{2\pi}\Big{[}\frac{1}{(\omega-\omega_{0})^{2}+\Delta^{2}}+\frac{1}{(\omega+\omega_{0})^{2}+\Delta^{2}}\Big{]}$,
$K_{f}=4/\sqrt{1+(\omega_{0}/\Delta)^{2}}$, $K_{c}=2$, and $R_{c}=\sqrt{2}/2$.
In (a1) $\omega_{0}=\Delta=0.1$ and $K_{f}=2\sqrt{2}>K_{c}$, which leads to a
first-order-like phase transition. In (a2) $\omega_{0}=\sqrt{3}/10$,
$\Delta=0.1$, and $K_{f}=2=K_{c}$, which, on the other hand, leads to a hybrid
phase transition. In (a3) $\omega_{0}=\sqrt{15}/10$, $\Delta=0.1$, and
$K_{f}=1<K_{c}$, for which the system undergoes a tiered phase transition,
where the oscillatory state arises in the region $K\in(1,2)$. For (b1)-(b3) we
have used $\beta<1$, $g(\omega)=1/2$ with $\omega\in(-1,1)$, $K_{f}=0$, and
$K_{c}=(\beta+1)^{\frac{\beta+1}{2}}\beta^{-{\frac{\beta}{2}}}$. Under these
conditions the system undergoes a tiered phase transition with a vanishing
onset from the oscillatory state to the two-cluster state. Parameter values
were as follows: (b1) $\beta=0.2$, $K_{c}=1.31$, $R_{c}=0.41$; (b2)
$\beta=0.5$, $K_{c}=1.61$, $R_{c}=0.58$; (b3) $\beta=0.8$, $K_{c}=1.86$,
$R_{c}=0.67$. For (c1)-(c3) we have used
$\beta>1$,$g(\omega)=1/\sqrt{2\pi}e^{-\omega^{2}/2}$, $K_{f}=+\infty$, and
$K_{c}=(\beta+1)^{\frac{\beta+1}{2}}\beta^{-{\frac{\beta}{2}}}$, leading to an
irreversible abrupt desynchronization transition. Parameter values were as
follows: (c1) $\beta=1.5$, $K_{c}=2.32$, $R_{c}=0.77$; (c2) $\beta=2$,
$K_{c}=2.60$, $R_{c}=0.82$; (c3) $\beta=3$, $K_{c}=3.08$, $R_{c}=0.87$. For
the simulations, we have used $N=100000$, $t=950$ with time step $\Delta
t=0.01$, and $T=50$. The red circles and the blue crosses represent the
forward and backward continuations, respectively. The black lines depict
theoretical predictions obtained with the mean-field approach, whereby solid
and dash lines depict stable and unstable branches, respectively.

Finally, to conclude the linear stability analysis, let
$\lambda\in\mathbb{C}$, and

$\frac{d\mathbf{V}}{dt}=\mathbf{MV+P}\hat{g}\mathbf{V}=\lambda\mathbf{V}.$
(66)

The steady state solution $\mathbf{V}$ is then solved as

$\mathbf{V=\big{(}\lambda I-M\big{)}}^{-1}\mathbf{P}\hat{g}\mathbf{V}.$ (67)

By applying the integral operator $\hat{g}$ to both sides of Eq. (67), we have
the linear equations for the vector $\hat{g}\mathbf{V}$

$\big{(}\mathbf{T-I}\big{)}\hat{g}\mathbf{V}=0,$ (68)

where the matrix $\mathbf{T}$ is defined by

$\begin{split}\mathbf{T}&=\hat{g}\big{(}\lambda\mathbf{I}-\mathbf{M}\big{)}^{-1}\mathbf{P}\\\
&=\left(\begin{array}[]{cccc}\hat{g}\frac{a(\omega)}{\lambda-c(\omega)}&\hat{g}\frac{b(\omega)}{\lambda-c(\omega)}\\\\[8.0pt]
\hat{g}\frac{\bar{b}(\omega)}{\lambda-\bar{c}(\omega)}&\hat{g}\frac{\bar{a}(\omega)}{\lambda-\bar{c}(\omega)}\end{array}\right)\\\\[8.0pt]
&=\left(\begin{array}[]{cccc}T_{11}(\lambda)&T_{12}(\lambda)\\\\[8.0pt]
T_{21}(\lambda)&T_{22}(\lambda)\end{array}\right).\end{split}$ (69)

Therefore, the nontrivial solution of $\hat{g}\mathbf{V}$ in Eq. (68) requires
that

$\rm det\big{(}\mathbf{T}-\mathbf{I}\big{)}=0.$ (70)

For real $\lambda$, we have $T_{11}(\lambda)=\bar{T}_{22}(\lambda)$,
$T_{12}(\lambda)=\bar{T}_{21}(\lambda)$, and it is possible to show that, for
a symmetric $g(\omega)$, $T_{11}$ and $T_{12}$ are real. Hence,
$T_{11}(\lambda)=T_{22}(\lambda)$ and $T_{12}(\lambda)=T_{21}(\lambda)$. The
characteristic equation (70) reduces to

$\big{(}1-T_{11}-T_{12}\big{)}\big{(}1-T_{11}+T_{12}\big{)}=0$ (71)

or equivalently to mirollo2007the

$1=\hat{g}\frac{a(\omega)+b(\omega)}{\lambda-c(\omega)},$ (72)
$1=\hat{g}\frac{a(\omega)-b(\omega)}{\lambda-c(\omega)}.$ (73)

These equations require the consideration of two different scenarios. Firstly,
for $q<1$, the equations reduce to

$\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}=\hat{g}\frac{\lambda\beta|\omega|q^{-\frac{\beta+1}{\beta}}\big{(}1-\sqrt{1-q^{2}}\big{)}}{\lambda^{2}+\omega^{2}(1-q^{2})},$
(74)

and

$\bigg{(}\frac{1}{K}\bigg{)}^{\frac{1}{\beta}}=\hat{g}\frac{\lambda|\omega|q^{-\frac{\beta+1}{\beta}}\big{(}q^{2}-1+\sqrt{1-q^{2}}\big{)}}{\lambda^{2}+\omega^{2}(1-q^{2})}.$
(75)

According to Eq. (8), the only solution to the self-consistency equation is
$R=q=0$. For $\beta<1$, $\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(1-\sqrt{1-q^{2}})=\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(q^{2}-1+\sqrt{1-q^{2}})=\lim\limits_{q\to
0}\frac{1}{2}q^{\frac{\beta-1}{\beta}}=+\infty$, which implies that the
incoherent state loses its stability once $K>0$. This in fact corresponds to
the so-called vanishing onset. Moreover, for $\beta>1$, $\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(1-\sqrt{1-q^{2}})=\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(q^{2}-1+\sqrt{1-q^{2}})=\lim\limits_{q\to
0}\frac{1}{2}q^{\frac{\beta-1}{\beta}}=0$. In order to balance Eqs. (74) and
(75), $\lambda=\pm i\omega$ must correspond to a continuous spectrum, implying
that the incoherent state remains neutrally, or asymptotically, stable for all
$K>0$. Particularly, for $\beta=1$ $\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(1-\sqrt{1-q^{2}})=\lim\limits_{q\to
0}q^{-\frac{\beta+1}{\beta}}(q^{2}-1+\sqrt{1-q^{2}})=\lim\limits_{q\to
0}\frac{1}{2}q^{\frac{\beta-1}{\beta}}=\frac{1}{2}$, so that thus Eqs. (74)
and (75) degenerate to the same characteristic equation strogatz1991stability

$\frac{2}{K}=\hat{g}\frac{\lambda|\omega|}{\lambda^{2}+\omega^{2}},$ (76)

which is consistent with the results obtained in Ref. xu2018origin . The
incoherent state loses its stability at a critical point $K_{f}$ which is
determined by the condition $\rm{Re}(\lambda)\to 0^{+}$.

Secondly, and lastly, for $q>1$, Eqs. (74) and (75) correspond precisely to
the limit formula of Eqs. (46) and (47), respectively, where
$\hat{g}\to\frac{1}{N}\sum_{i=1}^{N}$, which thus completes the proof.

Figure 2: Microscopic properties of the different oscillatory states. The top
columns show the instantaneous phases in dependence on the natural
frequencies, while the bottom columns show the effective frequencies
$\left\langle\omega_{i}\right\rangle=1/T\int_{t}^{t+T}\dot{\theta}_{i}(t)dt$
in dependence on the natural frequencies. The inserts show the time series of
the order parameter. The following parameters have been used: in (a1) and (a2)
$\beta=0.5$, $g(\omega)=1/2$ with $\omega\in(-1,1)$, $K=1.5$, in (b1) and (b2)
$\beta=0.8$, $g(\omega)=1/2$ with $\omega\in(-1,1)$, $K=1.7$, in (c1) and (c2)
$\beta=1.0$,
$g(\omega)=\frac{\Delta}{2\pi}\Big{[}\frac{1}{(\omega-\omega_{0})^{2}+\Delta^{2}}+\frac{1}{(\omega+\omega_{0})^{2}+\Delta^{2}}\Big{]}$,
(i.e., $g(\omega)$ is a bimodal Lorentzian distribution with $\pm\omega_{0}$
being the locations of the two peaks of the distribution, and $\Delta$ is the
half width of the distribution ) with $\omega_{0}=0.3,\Delta=0.1$, $K=1.6$,
and in (d1) and (d2) $\beta=1.0$, the bimodal Lorentzian distribution
$g(\omega)=\frac{\Delta}{2\pi}\Big{[}\frac{1}{(\omega-\omega_{0})^{2}+\Delta^{2}}+\frac{1}{(\omega+\omega_{0})^{2}+\Delta^{2}}\Big{]}$,
with $\omega_{0}=0.3,\Delta=0.1$, $K=1.8$. For the simulations, we have used
the same values for $N$, $t$, $\Delta t$, and $T$ as in Fig. 1.

##  VI Numerical results

To show how the heterogeneous and nonlinear coupling affects the phase
coherence, we study the system numerically. Figs. 1 and 2 illustrate the
macroscopical and microscopical synchronization properties, respectively. When
the exponent $\beta$ is fixed, and the natural frequencies are drawn
independently from a symmetric distribution $g(\omega)$, we can indeed observe
rich collective dynamics from the still relatively simple extension of the
traditional Kuramoto model. In Fig. 1, we show the time-averaged amplitude of
the order parameter, $\left\langle
R\right\rangle=\frac{1}{T}\int_{t}^{t+T}R(t)dt$ with $T$ being the time-
averaged window, as a function of the global coupling strength $K$. In the
simulations, both directions for synchronization transitions are considered,
i.e., $K$ is first increased adiabatically from $0$ and then decreased. The
presented simulation results reveal that the inhomogeneity and nonlinearity in
the coupling combines together to affect the overall collective dynamics of
the system, leading to a number of interesting dynamical phenomena.

In particular, for $\beta=1$ (left column), the system displays three types of
phase transitions to synchronization depending on the difference between
$K_{f}$ and $K_{c}$. Specifically, if $K_{f}>K_{c}$ [Fig. 1(a1)], the system
exhibits a first-order (explosive) phase transition to synchronization. As the
system transitions from the completely incoherent state ($\left\langle
R\right\rangle\sim 0$) to the fully locked state ($\left\langle
R\right\rangle\sim 1$) at the critical coupling $K_{f}$ under forward
continuation, when $K$ is gradually decreased, another abrupt transition from
synchronization to incoherence occurs at $K_{c}$. In the hysteretic region,
$K\in(K_{c},K_{f})$, the system admits a bistability, where the stable
incoherent state and synchronized state coexist. If $K_{f}=K_{c}$ [Fig.
1(a2)], the system undergoes a hybrid phase transition, which is characterized
by a vanishing hysteresis loop. Interestingly, if $K_{f}<K_{c}$ [Fig. 1(a3)],
a tiered phase transition to synchronization occurs, where the system converts
from the incoherent state to the fully synchronized state, which is mediated
by an oscillatory state that emerges when $K\in(K_{f},K_{c})$.

When the exponent $\beta\neq 1$, the phase diagrams change considerably. For
$\beta<1$ (middle column), we have shown that $K_{f}=0<K_{c}$, so that the
system displays a tiered phase transition to synchronization with a vanishing
onset, where the incoherent state is replaced by an oscillatory state the
amplitude of the mean-field is time-varying that arises in the region
$K\in(0,K_{c})$. However, for $\beta>1$ (right column), $K_{f}=+\infty>K_{c}$,
an irreversible abrupt desynchronization transition is observed, featuring a
hysteresis area, whose width is infinite. The system experiences a
discontinuous transition from the fully locked state to the incoherent state
at $K_{c}$ in the backward continuation, where $\left\langle R\right\rangle$
jumps from $R_{c}$ to $0$ directly. In contrast, there is no corresponding
counterpart of this abrupt transition when going from the incoherent to the
coherent state.

In Fig. 2, we highlight the oscillatory states that emerge for the
intermediate region of coupling strengths note1 . We show the instantaneous
phases, effective frequencies, and the time series of $R(t)$. These results
illustrate that, in the oscillatory state, the oscillators form coexisting
synchronized and drifting groups, or spontaneously partition into different
clusters according to their natural frequencies. Moreover, the effective
frequencies remain entrained in each cluster, but they do not synchronize with
one another. As a result, the instantaneous phases show a correlation at a
fixed time, thereby causing the order parameter to behave periodically. For
further details, we refer to the results shown in Fig. 2 and the corresponding
description of the parameter values in the figure caption.

##  VII Discussion

In summary, we have studied an extension of the Kuramoto model, where the
randomness of the coupling strength is related to the natural frequencies of
the coupled oscillators, and the mean-field is modified to depend on the
global coherence of the system. In doing so, we have thus jointly considered
heterogeneity and nonlinearity in the coupling of the Kuramoto model, which
has led us to uncover a number of fascinating dynamical states and transitions
associated with their emergence. These include coherent and incoherent states,
two-cluster states, and oscillatory states, which can emerge via various
transitions among them, including explosive synchronization transitions,
hybrid transitions with hysteresis absence, abrupt irreversible
desynchronization transitions, and tiered phase transitions with or without a
vanishing onset.

Despite the wealth of different dynamical states and the transitions between
them, our model still proved to be analytically tractable, thus allowing us to
provide firm theoretical foundations for the observed collective behavior. In
particular, we have provided a rigorous stability analysis of the equilibrium
states by studying their eigenspecturm properties in the finite and infinite
size limit. Moreover, we have shown that the critical points that correspond
to bifurcations of steady states and the stable conditions for their
occurrence in the phase space can be obtained by means of the matrix theory
and the Ott-Antonsen reduction.

Although the Kuramoto model and its many variants have been studied
extensively in the past decades, we have shown that the joint consideration of
heterogeneity and nonlinearity in the coupling among oscillators affords new
insights into the emergent collective dynamics that could be of relevance for
complex systems where such conditions apply. Previous research has already
highlighted the importance of heterogeneous interactions for the dynamics of
excitatory-inhibitory neurons, and for the emergence of consensus amongst
conformists and contrarians in different social settings cb2003 ; ms2005 ;
2008 . Heterogeneity and nonlinearity are also often related to plasticity,
which leads to changes in the connection strength among individuals to achieve
efficient global states. This may have implications in neuroscience as well as
in the social sciences, where such interaction scenarios often unfold. More
generally, our research reveals how an alternative coupling scheme can lead to
the emergence of diverse rhythmic states in complex systems that exhibit
synchronization transitions with relatively small changes to system
parameters. Since such complexity is rarely analytically tractable, we hope
that our model will prove to be inspirational for future research along these
lines and find applicability also in realistic complex systems.

###### Acknowledgements.

This work is supported by the National Natural Science Foundation of China
(Grants No. 11905068, 11875132, and 11835003), the Scientific Research Funds
of Huaqiao University (Grant No. ZQN-810), the Natural Science Foundation of
Shanghai (Grants No. 18ZR1411800), and the Slovenian Research Agency (Grant
Nos. P1-0403 and J1-2457).

## References

  * (1) Y. Kuramoto, in _International Symposium on Mathematical Problems in Theoretical Physics_ , edited by H. Araki, Lecture Notes in Physics No. 30 (Springer, New York, 1975), p. 420.
  * (2) S. A. Marvela and S. H. Strogatz, Invariant submanifold for series arrays of Josephson junctions, Chaos 19, 013132 (2009).
  * (3) B. Ermentrout, An adaptive model for synchrony in the firefly Pteroptyx malaccae, J. Math. Biol. 29, 571 (1991).
  * (4) M. Osaka, Modified Kuramoto phase model for simulating acrdiac pacemaker cell synchronization, Appl. Math. 8, 1227 (2017).
  * (5) R. M. May, S. A. Levin, and G. Sugihara, Complex systems: Ecology for bankers, Nature 451, 893 (2008).
  * (6) S. H. Strogatz, From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators, Physica. D 143, 1 (2000).
  * (7) J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77, 137 (2005).
  * (8) A. Pikovsky and M. Rosenblum, Dynamics of globally coupled oscillators: Progress and perspectives, Chaos 25, 097616 (2015).
  * (9) Boccaletti S, Pisarchik A N, Del Genio C I, Amann A. Synchronization: from coupled systems to complex networks. Cambridge: Cambridge University Press, 2018.
  * (10) C. Xu, X. Wang, and P. S. Skardal, Universal scaling and phase transitions of coupled phase oscillator populations, Phys. Rev. E 102, 042310 (2020).
  * (11) H. Hong and S. H. Strogatz, Kuramoto model of coupled oscillators with positive and negative coupling parameters: an example of conformist and contrarian oscillators, Phys. Rev. Lett. 106, 054102 (2011).
  * (12) D. Iatsenko, S. Petkoski, P. V. E. McClintock, and A. Stefanovska, Stationary and traveling wave states of the Kuramoto model with an arbitrary distribution of frequencies and coupling strengths, Phys. Rev. Lett. 110, 064101 (2013).
  * (13) H. Daido, Quasientrainment and slow relaxation in a population of oscillators with random and frustrated interactions, Phys. Rev. Lett. 68, 1073 (1992).
  * (14) H. Daido, Algebraic relaxation of an order parameter in randomly coupled limit-cycle oscillators, Phys. Rev. E 61, 2145 (2000).
  * (15) J. C. Stiller and G. Radons, Self-averaging of an order parameter in randomly coupled limit-cycle oscillators, Phys. Rev. E 61, 2148 (2000).
  * (16) D. H. Zanette, Synchronization and frustration in oscillator networks with attractive and repulsive interactions, Europhys. Lett. 72, 190 (2005).
  * (17) C. Börgers and N. Kopell, Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity, Neural Comput. 15, 509 (2003).
  * (18) M. S. de la Lama, J. M. López, and H. S. Wio, Spontaneous emergence of contrarian-like behaviour in an opinion spreading model, Europhys. Lett. 72, 851 (2005).
  * (19) G. H. Paissan and D. H. Zanette, Synchronization of phase oscillators with heterogeneous coupling: A solvable case, Physica D 237, 818 (2008).
  * (20) H. Hong and S. H. Strogatz, Conformists and contrarians in a Kuramoto model with identical natural frequencies, Phys. Rev. E 84, 046202 (2011).
  * (21) H. Hong and S. H. Strogatz, Mean-field behavior in coupled oscillators with attractive and repulsive interactions, Phys. Rev. E 85, 056210 (2012).
  * (22) I. M. Kloumann, I. M. Lizarraga, and S. H. Strogatz, Phase diagram for the Kuramoto model with van Hemmen interactions, Phys. Rev. E 89, 012904 (2014).
  * (23) T. Peron, Discordant synchronization patterns on directed networks of identical phase oscillators with attractive and repulsive couplings, Phys. Rev. E 103, 042210 (2021).
  * (24) F. Parastesh, S. Jafari, H. Azarnoush, Z. Shahriari, Z. Wang, S. Boccaletti, and M. Perc, Chimeras, Phys. Rep. 898, 1-114 (2021).
  * (25) D. Pazó, Thermodynamic limit of the first-order phase transition in the Kuramoto model, Phys. Rev. E 72, 046211 (2005).
  * (26) S. Boccaletti, J. A. Almendral, S. Guan, I. Leyva, Z. Liu, I. Sendia-Nadal, Z. Wang, and Y. Zou, Explosive transitions in complex networks structure and dynamics: Percolation and synchronization, Phys. Rep, 660, 1 (2016).
  * (27) R. D’ Souza, J. Gómez-Gardeñes, J. Nagler, and A. Arenas, Explosive phenomena in complex networks, Adv. Phys. 68, 123 (2019).
  * (28) C. Kuehn and C. Bick, A universal route to explosive phenomena, Sci. Adv. 7, eabe3824 (2021).
  * (29) J. Gómez-Gardeñes, S. Gómez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Phys. Rev. Lett. 106, 128701 (2011).
  * (30) C. Xu, J. Gao, Y. Sun, X. Huang, and Z. Zheng, Explosive or continuous: Incoherent state determines the route to synchronization, Sci. Rep. 5, 12039 (2015).
  * (31) C. Xu, J. Gao, S. Boccaletti, Z. Zheng, and S. Guan, Synchronization in starlike networks of phase oscillators, Phys. Rev. E 100, 012212 (2019).
  * (32) M. Wang, D. Arteaga, and B. J. He, Brain mechanisms for simple perception and bistable perception, Proc. Natl. Acad. Sci. 110, 3350 (2013).
  * (33) P. S. Skardal and A. Arenas, Control of coupled oscillator networks with application to microgrid technologies, Sci. Adv. 1, 1500339 (2015).
  * (34) T. Tanaka and T. Aoyagi, Multistable attractors in a network of phase oscillators with three-body interactions, Phys. Rev. Lett. 106, 224101 (2011).
  * (35) M. Komarov and A. Pikovsky, Finite-size-induced transitions to synchrony in oscillator ensembles with nonlinear global coupling, Phys. Rev. E 92, 020901(R) (2015).
  * (36) A. Pikovsky, Low-dimensional dynamics for higher-order harmonic, globally coupled phase-oscillator ensembles, Phys. Rev. E 100, 062210 (2019).
  * (37) P. S. Skardal and A. Arenas, Abrupt desynchronization and extensive multistability in globally coupled oscillator simplexes, Phys. Rev. Lett. 122, 248301 (2019).
  * (38) A. P. Millán, J. J. Torres, and G. Bianconi, Explosive higher-order Kuramoto dynamics on simplicial complexes, Phys. Rev. Lett. 124, 218301 (2020).
  * (39) C. Xu, X. Wang, and P. S. Skardal, Bifurcation analysis and structural stability of simplicial oscillator populations, Phys. Rev. Research 2, 023281 (2020).
  * (40) C. Xu and P. S. Skardal, Spectrum of extensive multiclusters in the Kuramoto model with higher-order interactions, Phys. Rev. Research 3, 013013 (2021).
  * (41) C. Xu, X. Wang, Z. Zheng, and Z. Cai, Stability and bifurcation of collective dynamics in phase oscillator populations with general coupling, Phys. Rev. E 103, 032307 (2021).
  * (42) H. Wang and X. Li, Synchronization and chimera states of frequency-weighted Kuramoto-oscillator networks, Phys. Rev. E 83, 066214 (2011).
  * (43) X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive synchronization in a general complex network, Phys. Rev. E 88, 010802(R) (2013).
  * (44) D. Yuan, M. Zhang, and J. Zhong, Dynamics of the Kuramoto model in the presence of correlation between distributions of frequencies and coupling strengths, Phys. Rev. E 89, 012910 (2014).
  * (45) H. Bi, X. Hu, S. Boccaletti, X. Wang, Y. Zou, Z. Liu, and S. Guan, Coexistence of quantized, time dependent, clusters in globally coupled oscillators, Phys. Rev. Lett. 117, 204101 (2016).
  * (46) H. Hong, K. P. O’Keeffe, and S. H. Strogatz, Phase coherence induced by correlated disorder, Phys. Rev. E 93, 022219 (2016).
  * (47) C. Xu, J. Gao, H. Xiang, W. Jia, S. Guan, and Z. Zheng, Dynamics of phase oscillators with generalized frequency-weighted coupling, Phys. Rev. E 94, 062204 (2016).
  * (48) C. Xu, S. Boccaletti, S. Guan, and Z. Zheng, Origin of Bellerophon states in globally coupled phase oscillators, Phys. Rev. E 98, 050202(R) (2018).
  * (49) C. Xu, S. Boccaletti, Z. Zheng, and S. Guan, Universal phase transitions to synchronization in Kuramoto-like models with heterogeneous coupling, New J. Phys. 21, 113018 (2019).
  * (50) C. Xu and Z. Zheng, Bifurcation of the collective oscillatory state in phase oscillators with heterogeneity coupling, Nonlin. Dyn. 98, 2365 (2019).
  * (51) M. Rosenblum and A. Pikovsky, Self-organized quasiperiodicity in oscillator ensembles with global nonlinear coupling, Phys. Rev. Lett. 98, 064101 (2007).
  * (52) G. Filatrella, N. F. Pedersen, and K. Wiesenfeld, Generalized coupling in the Kuramoto model. Phys. Rev. E 75, 017201 (2007).
  * (53) X. Zhang, S. Boccaletti, S. Guan, and Z. Liu, Explosive synchronization in adaptive and multilayer networks, Phys. Rev. Lett. 114, 038701 (2015).
  * (54) H. Ly, H. Kuo, and K. Wu, Kuramoto model of weakly conformed oscillators, Phys. Rev. E 98, 052230 (2018).
  * (55) W. Zou and J. Wang, Dynamics of the generalized Kuramoto model with nonlinear coupling: Bifurcation and stability, Phys. Rev. E 102, 012219 (2020).
  * (56) J. Gao, and K. Efstathiou, Reduction of oscillator dynamics on complex networks to dynamics on complete graphs through virtual frequencies, Phys. Rev. E 101, 022302 (2020).
  * (57) X. Dai, X. Li, H. Guo, D. Jia, M. Perc, P. Manshour, Z. Wang, and S. Boccaletti, Discontinuous transitions and rhythmic states in the D-dimensional Kuramoto model induced by a positive feedback with the global order parameter, Phys. Rev. Lett. 125, 194101 (2020).
  * (58) D. Iatsenko, P. V. E. McClintock, and A. Stefanovska, Glassy states and super-relaxation in populations of coupled phase oscillators, Nat. Commun. 5, 4118 (2014).
  * (59) We emphasize that the form of the coupling that we are considering strongly deviates from the classical Kuramoto model, the huge body of literature quoted above points actually to the fact that a such form of coupling (or similar ones which include a dependence on the oscillators’ natural frequencies) are in fact of physical interest, and may be considered as a reasonable hypothesis for the description of collective dynamics observed in real-world systems.
  * (60) Y. Baibolatov, M. Rosenblum, Z. Zhanabaev, and Arkady Pikovsky, Complex dynamics of an oscillator ensemble with uniformly distributed natural frequencies and global nonlinear coupling, Phys. Rev. E 82, 016212 (2010).
  * (61) We remark that conformity is a well-established social psychology phenomenon demonstrated through experiments, which showed that people could conform to a majority to gain the social approval of others. Individual behavioral responses that vary with the others’ behavior are also commonly seen in animals. For example, social insects such as Temnothorax ants and honeybees select their new nest site based on a popularity competition of recruitment between potential nests. These examples suggest that individual behavioral response to conform or synchronize with others depends on the degree of consensus or synchronization of the group at that moment. From the point of view of coupled phase oscillator models, one may use NC to better model conformity behavior.
  * (62) S. Petkoski, D. Iatsenko, L. Basnarkov, and A. Stefanovska, Mean-field and mean-ensemble frequencies of a system of coupled oscillators, Phys. Rev. E 87, 032908 (2013).
  * (63) In principle, a nonzero $\Omega$ corresponding to the traveling waves should be considered for a complete self-consistent analysis. However, in this paper we do not discuss this scenario. On the one hand, the symmetry of the system ensures that $\Omega=0$ holds for the self-consistent equations. On the other hand, numerical simulations do not find the evidence of the traveling waves.
  * (64) R. E. Mirollo and S. H. Strogatz, The spectrum of the locked state for the Kuramoto model of coupled oscillators, Physica D 205, 249 (2005).
  * (65) E. Ott and T. M. Antonsen, Low dimensional behavior of large systems of globally coupled oscillators, Chaos 18, 037113 (2008).
  * (66) E. Ott and T. M. Antonsen, Long time evolution of phase oscillator systems, Chaos 19, 023117 (2009).
  * (67) O. E. Omel’chenko and M. Wolfrum, Bifurcations in the Sakaguchi-Kuramoto model, Physica D (Amsterdam, Neth.) 263, 74 (2013).
  * (68) R. Mirollo and S. H. Strogatz, The spectrum of the partially locked state for the Kuramoto model, J. Nonlinear Sci. 17, 309 (2007).
  * (69) S. H. Strogatz and R. E. Mirollo, Stability of incoherence in a population of coupled oscillators, J. Stat. Phys. 63, 613 (1991).
  * (70) Throughout this paper, the oscillatory states represent the collective states, in which the amplitude of the order parameter is time-dependent.