![FIG. 3. Numerical calculation of critical coupling Kc as a function of number of oscillators N of SF networks with γ = 2.25 under different normalizations N . Dashed line marks the result Kc = 2/[πg(0)]. Natural frequencies distributed according to g(ω) = 1/π (ω2 + 1). Each point is an average over 100 network realizations. Error bars are smaller than symbols.](/figures/fig-3-numerical-calculation-of-critical-coupling-kc-as-a-1qbothu2.webp)
![FIG. 1. Comparison between the estimations of Kc by susceptibilities χ [Eq. (14)] and χr [Eq. (15)]. Networks generated according to the UCM with degree distribution P(k) ∼ k−γ , with γ = 2.25 and kmin = 5. Natural frequencies are assigned according to Eq. (13). Each point is an average over 100 network realizations. Error bars are smaller than symbols.](/figures/fig-1-comparison-between-the-estimations-of-kc-by-3rpd7135.webp)
PHYSICAL REVIEW E 100, 042302 (2019)
Onset of synchronization of Kuramoto oscillators in scale-free networks
Thomas Peron,
1,4,*
Bruno Messias F. de Resende,
2
Angélica S. Mata,
3
Francisco A. Rodrigues,
1
and Yamir Moreno
4,5,6
1
Institute of Mathematics and Computer Science, University of São Paulo, São Carlos, São Paulo 13566-590, Brazil
2
São Carlos Institute of Physics, University of São Paulo, São Carlos, São Paulo 13566-590, Brazil
3
Departamento de Física, Universidade Federal de Lavras, 37200-000 Lavras, Minas Gerais, Brazil
4
Institute for Biocomputation and Physics of Complex Systems (BIFI), University of Zaragoza, E-Zaragoza 50018, Spain
5
Department of Theoretical Physics, University of Zaragoza, E-Zaragoza 50009, Spain
6
ISI Foundation, I-10126 Torino, Italy
(Received 10 May 2019; published 8 October 2019)
Despite the great attention devoted to the study of phase oscillators on complex networks in the last two
decades, it remains unclear whether scale-free networks exhibit a nonzero critical coupling strength for the
onset of synchronization in the thermodynamic limit. Here, we systematically compare predictions from the
heterogeneous degree mean-field (HMF) and the quenched mean-field (QMF) approaches to extensive numerical
simulations on large networks. We provide compelling evidence that the critical coupling vanishes as the
number of oscillators increases for scale-free networks characterized by a power-law degree distribution with
an exponent 2 <γ 3, in line with what has been observed for other dynamical processes in such networks.
For γ>3, we show that the critical coupling remains finite, in agreement with HMF calculations and highlight
phenomenological differences between critical properties of phase oscillators and epidemic models on scale-free
networks. Finally, we also discuss at length a key choice when studying synchronization phenomena in complex
networks, namely, how to normalize the coupling between oscillators.
DOI: 10.1103/PhysRevE.100.042302
I. INTRODUCTION
Synchronization processes are pervasively observed in a
wide range of physical, chemical, technological, and biolog-
ical systems [1]. These phenomena can, to a great extent, be
described by models of coupled phase oscillators. Arguably,
one of the most studied models in this context is the one pro-
posed by Kuramoto [2], which in the last decade was exten-
sively investigated when the oscillators are placed on complex
networks (see [3,4] and references therein). A key question
addressed in these studies is how the heterogeneous connec-
tivity pattern impacts on the onset of synchronization—or, in
other words, how the critical coupling strength required for
the emergence of collective motion is affected by the network
topology.
The relationship between structure and synchronous dy-
namics has been studied in many scenarios: From homoge-
neous and unclustered networks to heterogeneous and modu-
lar ones, in addition to variations of phase oscillator models
including correlations between intrinsic dynamics and local
topology [3,4]. Yet, despite the notorious advances achieved
over the past years, fundamental questions regarding the
collective dynamics of large ensembles of oscillators still
remain elusive. One of these problems is whether the critical
coupling strength for the onset of synchronization remains
finite in the thermodynamic limit for scale-free (SF) networks
characterized by a power-law degree distribution with an
exponent 2 <γ 3. Another important question concerns
*
thomaskaue@gmail.com
the very definition of the coupling strength in the dynamical
equations. This paper will address both challenges.
The above questions were already pointed out in the first
work that dealt with the dynamics of Kuramoto oscillators
on heterogeneous scale-free structures [5]. There, the authors
remarked on the supposed finite magnitude of the critical
coupling and highlighted the apparent contrast of the Ku-
ramoto dynamics with epidemic spreading and percolation—
processes which were already known to exhibit vanishing
critical points in the thermodynamic limit for SF topologies.
Subsequent theoretical approaches [6,7] estimated via mean-
field approximations that, in the absence of degree-degree
correlations, the critical coupling should converge to zero
as the number of oscillators tends to infinity—similarly to
what happens for other dynamical processes on networks [8].
However, later investigations reported significant deviations
between predictions of mean-field theories and numerical
simulations [9], casting further doubts on the validity of
the classical result on the nonexistence of a synchronization
threshold [3,4,10].
One clear difficulty for a precise estimation of the onset of
synchronous motion is, naturally, the sizes of the simulated
networks. Indeed, the first hypotheses on the existence or
absence of a critical point in the Kuramoto dynamics were
supported by numerical experiments considering populations
with sizes of the order of up to 10
4
oscillators [5–7,9]—a
value that potentially limits the accuracy of finite-size analysis
and calculations, especially in what concerns the detection
of the onset of synchronization for highly heterogeneous
structures. It is noteworthy to mention, though, that recent
contributions (see, e.g. [4,11–13]) have investigated finite-size
2470-0045/2019/100(4)/042302(8) 042302-1 ©2019 American Physical Society
THOMAS PERON et al. PHYSICAL REVIEW E 100, 042302 (2019)
effects of the dynamics and reported excellent agreement
between simulation and mean-field theories. However, most
of those analyses have focused on Erd
˝
os-Rènyi (ER) random
graphs and SF networks with degree exponent γ>3, situ-
ations in which the critical coupling is expected to be finite
(according to the heterogeneous degree mean-field approx-
imations) [4]. Of particular interest is a recent contribution
[14], where the authors investigated finite-size effects of ER
graphs, reaching networks with very large sizes (up to N =
2
27
nodes).
Another possible source of disagreement between mean-
field theories and numerical simulations in the estimation of
the critical coupling strength is the consideration of different
definitions of order parameters [3,8]. Very recently, Yook and
Kim [15] performed a thorough comparison between the clas-
sical Kuramoto order parameter and the order parameter ac-
counting for heterogeneous degree distributions. The authors
verified that, indeed, the definition of the order parameter
crucially affects the assessment of the asynchronous state
in highly heterogeneous SF networks. However, although
simulations with networks of size up to 10
7
oscillators were
carried out, it is not clear from [15] how the transition point
behaves as the network size increases. Therefore, the ques-
tion regarding whether or not there is a well-defined critical
coupling for the onset of synchronization in SF networks with
2 <γ 3 has remained without a concluding answer.
In order to address this problem, and given the difficulties
in performing very large numerical simulations, here we adopt
an alternative approach: We perform a systematic compar-
ison between simulations and the results derived using the
heterogeneous degree mean-field (HMF) and the quenched
mean-field (QMF) formulations in networks of sizes up to
N = 3 ×10
6
nodes. We show that the critical coupling pre-
dicted by both the HMF and the QMF agrees with the values
measured in numerical experiments for networks with power-
law exponent γ 3, hence providing stronger evidence that
the critical coupling of such systems vanishes in the thermo-
dynamic limit. For SF networks whose degree distribution has
a finite second statistical moment, we find that the onset of
synchronization remains constant in the thermodynamic limit.
Furthermore, we highlight differences between the critical
behavior of synchronization dynamics and that found in the
disease spreading process. In particular, we verify that HMF
correctly predicts a finite critical threshold in the thermody-
namic limit for γ>3, in contrast to results obtained in the
context of epidemic dynamics [16,17].
Additionally, we also revisit another issue debated in
early works on phase oscillator models on networks: How
to define and properly normalize the coupling function in
the dynamical equations. In particular, we verify that some
choices previously considered as appropriate for SF networks
actually induce undesired dependencies on the system’s size,
including the increase of the onset of synchronization as
networks become larger, and an infinite coupling strength that
locks low degree nodes in the thermodynamic limit of highly
heterogeneous networks. The rest of the paper is organized as
follows: In Sec. II, we provide a brief review of the mean-
field approximations to treat coupled oscillators in heteroge-
neous networks. In Sec. III, we compare the estimations by
mean-field theories with numerical simulations. Section IV is
devoted to the discussion on the coupling normalization. We
give our conclusions in Sec. V.
II. MEAN-FIELD THEORIES FOR PHASE OSCILLATORS
IN HETEROGENEOUS NETWORKS
In this section, we provide a brief review of the main
analytical approximations used to deal with ensembles of
phase oscillators in heterogeneous networks. The Kuramoto
model consists of the following system of equations [3,4]:
˙
θ
i
(t ) = ω
i
+ K
N
j=1
A
ij
sin(θ
j
− θ
i
), (1)
where θ
i
and ω
i
are the phase and natural frequency of the ith
oscillator, respectively; K is the coupling strength, and A is the
adjacency matrix, with A
ij
= 1 if nodes i and j are connected,
and 0 otherwise.
In order to assess the overall synchrony of an ensemble of
oscillators, Kuramoto [2] introduced the order parameter,
Re
i
=
1
N
N
j=1
e
iθ
j
(t )
, (2)
where R and are the magnitude and phase of the centroid
associated with the N points e
iθ
j
(t )
in the complex plane,
respectively. If phases are uniformly distributed over [0, 2π],
it follows that R ≈ 0, whereas R ≈ 1 if oscillators rotate
grouped into a synchronous cluster.
In contrast to the case of globally connected populations,
the original analytical treatment via a self-consistent analysis
by Kuramoto [2] cannot be directly extended to the network
case. The reason for this relies on the fact that Eq. (1)is
not exactly decoupled by a global order parameter. Instead,
an exact decoupling is only achieved by defining local order
parameters as
r
i
e
iψ
i
(t )
=
N
j=1
A
ij
e
iθ
j
(t )
, (3)
which leads to
˙
θ
i
(t ) = ω
i
+ Kr
i
sin(ψ
i
− θ
i
). (4)
In this paper, we consider the oscillators frequencies ω
i
to be distributed according to a smooth and unimodal distri-
bution g(ω) centered at ω = 0. By inserting the fixed point
solution (
˙
θ
i
(t ) = 0) of the equation above into Eq. (3), and
performing a self-consistent analysis of the resulting equation,
one arrives at the critical coupling given by [4,9]
K
QMF
c
=
2
πg(0)
1
max
, (5)
where
max
is the largest eigenvalue of A. The latter result
was first derived in [9], with what the authors called pertur-
bation theory of the Kuramoto model on complex networks.
Henceforth, we refer to Eq. (5) as the QMF critical coupling
strength, owing to the similarity with epidemic thresholds
derived with techniques that preserve the quenched structure
of the network [18]. To gain further insights on the predictions
of Eq. (5) to the dynamics on SF networks, we recall the
042302-2
ONSET OF SYNCHRONIZATION OF KURAMOTO … PHYSICAL REVIEW E 100, 042302 (2019)
result [19],
max
∼
⎧
⎨
⎩
k
2
k
if
k
2
k
>
√
k
max
ln(N ),
x
√
k
max
if
√
k
max
>
k
2
k
ln
2
(N ),
(6)
where k
max
is the maximum degree of the network. In uncor-
related SF networks, k
max
scales as k
max
∼ N
1/2
if 2 <γ
3, and k
max
∼ N
1/(γ −1)
,forγ>3. By noticing further that
k
2
/k∼k
3−γ
max
√
k
max
, we then estimate [18]
K
QMF
c
2
πg(0)
×
k
k
2
if 2 <γ <5/2,
1
√
k
max
if γ>5/2.
(7)
Therefore, according to the QMF approach, the critical cou-
pling K
c
should vanish in the thermodynamic limit as k
max
diverges, even if k
2
remains finite (i.e., the case when
γ>3).
Another way of modeling synchronization processes on
networks is by virtue of the annealed network approximation
[4]. It consists of replacing the elements of the adjacency
matrix A
ij
by its ensemble average
˜
A
ij
, which corresponds to
the probability that two nodes, i and j, are connected in the
configuration model; that is,
˜
A
ij
=
k
i
k
j
Nk
, (8)
where k
i
is the degree of node i. Substituting Eq. (8)into
Eq. (1) yields
˙
θ
i
(t ) = ω
i
+
Kk
i
Nk
j
k
j
sin(θ
j
− θ
i
). (9)
The previous equation motivates the definition of the follow-
ing order parameter:
re
iψ (t )
=
1
Nk
N
j=1
k
j
e
iθ
j
(t )
. (10)
Equation (8) is equivalent to the so-called heterogeneous
degree mean-field approximation (HMF) [4] and leads to the
definition of the order parameter in Eq. (10). Essentially, in the
HMF approximation, one assumes that the network topology
(initially fully represented by the adjacency matrix A)is
abstracted in the degree distribution P(k); that is, nodes are
coarse grained according to their degrees and the oscillators
become statistically equivalent, differing only by the parame-
ters k
i
and ω
i
.
By decoupling Eq. (9) with Eq. (10), and performing a self-
consistent analysis of the equations, one can show that the
onset of synchronization within the annealed approximation
occurs when [4,6]
K > K
HMF
c
=
2
πg(0)
k
k
2
, (11)
where k
n
=
k
k
n
P(k)isthenth moment of the degree
distribution P(k).
As previously mentioned, it has been recently shown [15]
that the traditional order parameter [Eq. (2)] and the one
introduced by the HMF approximation yield different results
when assessing the synchronization of networks. In particular,
the latter overestimates the level of coherence among the
oscillators in the asynchronous regime for SF networks. This
effect is particularly evident in networks having hubs whose
degree scales with O(N ); however, discrepancies between R
and r are also likely to emerge for networks with power-law
exponent γ>3[15]. Therefore, in this paper, we evaluate the
onset synchronization numerically using the standard order
parameter R in Eq. (2).
Our goal is to systematically investigate the behavior
of the onset of synchronization as the size of SF networks
increases, comparing the theoretical predictions provided by
the current mean-field approaches. Seeking to keep the source
of fluctuations across network realizations to a minimum, we
assign natural frequencies deterministically according to [12]
i
N
−
1
2N
=
ω
i
−∞
g(ω)dω, (12)
which for the Lorentzian distribution g(ω) =
π
1
ω
2
+
2
yields
ω
i
= tan
iπ
N
−
(N + 1)π
2N
, i = 1,...,N. (13)
In this way, we generate a set of quasiuniformly spaced
frequencies, removing, thus, the disorder introduced by
different realizations of frequencies [12].
III. CRITICAL COUPLING OF UNCORRELATED
SCALE-FREE NETWORKS
All networks analyzed in this section were generated
following the uncorrelated configuration model (UCM) [20]
considering a power-law degree distribution P(k) ∼ k
−γ
with
the cutoff k
max
∼ N
1/2
for γ 3, and k
max
∼ N
1/(γ −1)
for
γ>3. Furthermore, in order to avoid sample-sample fluctu-
ations on k
max
, for each value of N,wefixedk
max
=k
max
across the network realizations. Simulations were performed
on graphical processing units (GPUs) and by using the Heun’s
method with time steps adapted according to the value of N .
Our implementation uses
TENSORFLOW and is available in the
Python package
STDOG [21].
Typically, the critical coupling strength of finite networks
can be estimated numerically via detecting the divergent peak
of the susceptibility,
χ = N
R
2
t
−R
2
t
, (14)
where ···
t
denotes a temporal average. However, we em-
ploy the modified susceptibility defined as [16]
χ
r
= N
R
2
t
−R
2
t
R
t
. (15)
As with the definition in Eq. (14), the modified susceptibil-
ity also exhibits a peak at K = K
c
. Nonetheless, analogous
forms of χ
r
have been shown to be better suited to detect
transition points in epidemic spreading and contact processes
in networks with diverging k
2
[16,22,23]. Thus, motivated
by those results, we extend this measure for the detection of
onset of the synchronous state. Our choice is confirmed by
the numerical results presented in Fig. 1.Forγ = 2.25, the
critical points estimated via Eq. (15) are in better agreement
with HMF and QMF theories than that estimated via Eq. (14),
042302-3
THOMAS PERON et al. PHYSICAL REVIEW E 100, 042302 (2019)
FIG. 1. Comparison between the estimations of K
c
by suscepti-
bilities χ [Eq. (14)] and χ
r
[Eq. (15)]. Networks generated according
to the UCM with degree distribution P(k) ∼ k
−γ
,withγ = 2.25 and
k
min
= 5. Natural frequencies are assigned according to Eq. (13).
Each point is an average over 100 network realizations. Error bars
are smaller than symbols.
especially for low values of N. Similar results are found for
different values of γ . Thus, we henceforth detect the critical
points via χ
r
.
Let us now analyze how the mean-field theories perform
in comparison with simulations for the different regimes of
γ .First,forγ<5/2, as discussed in the previous section,
both HMF and QMF predict a vanishing K
c
, which should
scale with k/k
2
. Indeed, as it is seen in Fig. 2(a),for
γ = 2.25, both theories predict quite accurately the onset of
synchronization.
Discrepancies between the approximations appear when
γ>5/2. To be precise, in this regime, HMF yields K
c
∼
k/k
2
, while QMF gives K
c
∼ k
−1/2
max
. As depicted in
Fig. 2(b), the mean-field theories provide a satisfactory
approximation of the synchronization thresholds for networks
with γ = 2.7. Note that, although QMF contains in its for-
mulation the whole information about the network topology,
it performs slightly worse than HMF (see inset). Similar
dependencies with the system size are found for epidemic
thresholds in SF networks with 5/2 <γ <3[16,17].
For γ = 3.5 [Fig. 2(c)], we observe that the numerical
calculation of K
c
converges to a constant value as N increases,
in agreement with the HMF prediction, whereas QMF theory
clearly fails in capturing the onset of synchronization. That
is, while simulations show that large SF networks in this case
exhibit a finite synchronization threshold, QMF reveals a van-
ishing K
c
. Furthermore, it is interesting to point out discrep-
ancies between synchronization and the epidemic spreading
described by the susceptible-infected-susceptible (SIS) model
[24,25] regarding the dependence on the system size for
γ>3. In contrast to the finite onset of synchronization seen in
Fig. 2(c), epidemic thresholds of the SIS model are known to
decay as N increases for γ>3[16,17]. In fact, Chatterjee and
Durrett [26] proved rigorously that, for uncorrelated random
networks with a power-law degree distribution P(k) ∼ k
−γ
with any γ , the SIS model presents an unstable absorbing
phase in the thermodynamic limit, resulting in a null epi-
demic threshold. Afterwards, Boguñá et al. [27] physically
interpreted this proof with a semianalytical approach taking
into account a long-range reinfection mechanism between
hubs and found a vanishing epidemic threshold including for
γ>3.
Actually, the behavior of the SIS model is distinct and more
intricate than other dynamical processes that also present a
phase transition from active to inactive states. This epidemic
(a)
(c)
(b)
FIG. 2. Critical coupling K
c
against network size N for UCM networks with power-law exponent (a) γ = 2.25, (b) γ = 2.7, (c) γ = 3.5.
All networks have k
min
= 5. Insets in (a) and (b) depict the difference between numerical estimation of K
c
and mean-field theories. Each point
is an average over 100 network realizations. Error bars are smaller than symbols.
042302-4
ONSET OF SYNCHRONIZATION OF KURAMOTO … PHYSICAL REVIEW E 100, 042302 (2019)
model is governed by mutual activation of hubs. Outliers,
a small amount of vertices with connectivity much larger
than the other nodes of the network, can sustain localized
epidemics for long times. This phenomenon causes a double-
peaked shape in the susceptibility curve [16,17] and the emer-
gence of Griffiths effects [28]—in this case, QMF captures the
peak associated with the activation of the largest hub in the
network [22]. Surprisingly, simulations with networks with
N = 10
7
(not shown here) did not reveal signs of multiple
peaks in susceptibility curves of Kuramoto oscillators. How-
ever, with the aim of understanding the nature of the threshold
in epidemic models on uncorrelated random networks, recent
works [23,29,30] showed that different epidemic models such
as, for instance, susceptible-infected-recovered-susceptibility
[31], contact process [32], the generalized SIS model with
weighted infection rates [33], and other alterations of the SIS
model [30], have a finite threshold in the thermodynamic
limit. This behavior is related to standard phase transitions
given by collective activation processes involving essentially
the whole network, as observed in the synchronization phe-
nomenon of the Kuramoto oscillators.
At last, the results in Fig. 2 point to a different scenario
as the one in [9] regarding the accuracy of mean-field theo-
ries. More precisely, in Fig. 2 we see that HMF exhibits an
excellent agreement with numerical simulations for γ = 2.25
and 2.7, and a qualitative agreement with the scaling with
N for γ = 3.5. Conversely, Ref. [9] found that HMF agrees
best with the numerical results obtained for γ>3, while
significantly deviating from simulations for 2 <γ <5/2; i.e.,
the opposite situation observed in Fig. 2. These discrepancies
are possibly due to structural correlations induced by the large
artificial cutoffs (k
max
∼ N) and the relative small size of the
SF networks (N ∼ 10
3
) considered in [9].
IV. COUPLING NORMALIZATION
The size dependence of the onset of synchronization on
the system’s size brings back to attention a topic intensively
debated in early studies of network synchronization [3,4],
namely, the choice for the normalization of the coupling
function. When dealing with phase oscillators on networks,
it is a common practice to let the oscillators interact through
unnormalized couplings, as in Eq. (1). The reason for this
resides in the fact that the definition of the coupling is not
as straightforward as for the model on fully connected graphs.
In the latter scenario, the number of neighbors of a given node
scales linearly with N; it thus suffices to set the K/N to assure
that the coupling is an intensive quantity.
The connectivity in real and synthetic networks, on the
other hand, scales differently with the number of oscillators,
making the definition of the coupling function to be not unique
and, therefore, motivating the formulation of the equations of
motion as Eq. (1). Nevertheless, the lack of an appropriate nor-
malization has several major consequences to the collective
dynamics of Kuramoto oscillators: (i) the vanishing character
of K
c
in the thermodynamic of limit for SF networks, as seen
in the previous section; (ii) the difficulty in comparing the
dynamics of networks with different connectivity patterns [3],
and (iii) the second term in the right-hand side of Eq. (1)
diverges in the thermodynamic limit for networks in which
FIG. 3. Numerical calculation of critical coupling K
c
as a func-
tion of number of oscillators N of SF networks with γ = 2.25
under different normalizations N . Dashed line marks the result
K
c
= 2/[πg(0)]. Natural frequencies distributed according to g(ω) =
1/π (ω
2
+ 1). Each point is an average over 100 network realizations.
Error bars are smaller than symbols.
the maximum degree is not bounded when N →∞.Inthis
section, we compare the impact of different prescriptions for
the coupling function in large heterogeneous networks in the
light of the latter points.
Let us now consider the phase equations defined as
˙
θ
i
= ω
i
+
K
N
i
N
j=1
A
ij
sin(θ
j
− θ
i
), (16)
where N
i
is the normalization constant of node i. Reasonable
choices for N
i
would then be quantities that are related to
the network topology. One of these prescriptions discussed in
previous works is N
i
= k
max
∀i [3,4]. It makes the summation
to be an intensive quantity, since it prevents this term from
diverging in highly heterogeneous networks. However, by
repeating the analysis of the previous section for SF net-
works with N = k
max
, we observe that the new normalization
yields a critical coupling that depends on the system’s size
(Fig. 3). This result is easily understood by noticing that,
under the HMF approximation and considering P(k) ∼ k
−γ
,
K
c
is rescaled to
K
c
=
2
πg(0)
k
k
2
k
max
∼ N
γ
2
−1
, (17)
explaining why the onset of synchronization increases in
this case. Curiously, this conclusion is not evident from
early works [3,4]. Phenomenologically, the previous result
can be understood by noticing that this normalization also
makes the coupling
K
N
i
→ 0 when N →∞for all nodes with
bounded connectivity in the thermodynamic limit. Thus, as
these degree-bounded nodes are effectively decoupled of the
hubs, one should expect K
c
→∞when N →∞.
On the other hand, if a K
c
that is independent of the
system’s size is sought, then a natural choice would be
to rescale the coupling according to N =k
2
/k. Indeed,
observing the corresponding result in Fig. 3, it looks like as
042302-5