On the onset of synchronization of Kuramoto oscillators in scale-free networks

TL;DR

This paper investigates the onset of synchronization in Kuramoto oscillators on scale-free networks, revealing that the critical coupling strength vanishes for certain degree exponents as network size grows, contrasting with finite values for others.

Contribution

It compares mean-field approaches with numerical simulations to clarify how the critical coupling depends on the degree distribution exponent in large scale-free networks.

Findings

Critical coupling vanishes for 2<γ≤3 as network size increases.

Critical coupling remains finite for γ>3.

Differences between synchronization and epidemic thresholds are highlighted.

Abstract

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 < \gamma \leq 3$, in line with what has been observed for other dynamical processes in such networks. For $\gamma > 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.

Equations

$${\dot{\theta }}_i(t)={\omega }_i+K\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$$R{e}^{\mathrm{i}\psi (t)}=\frac{1}{N}\sum _{j=1}^N{e}^{\mathrm{i}{\theta }_j(t)},$$

$$r_i{e}^{\mathrm{i}{\psi }_i(t)}=\sum _{j=1}^N{A}_{ij}{e}^{\mathrm{i}{\theta }_j(t)},$$

$${\dot{\theta }}_i(t)={\omega }_i+K{r}_i\sin ({\psi }_i-{\theta }_i).$$

$$K_c^{\mathrm{QMF}}=\frac{2}{\pi g(0)}\frac{1}{{\Lambda }_{\max }},$$

$${\Lambda }_{\max }\sim \{\begin{array}{ll}\frac{⟨k^2⟩}{⟨k⟩}& \text{mbox}if\frac{⟨k^2⟩}{⟨k⟩}>\sqrt{k_{\max }}\ln (N),\\ \sqrt{k_{\max }}& \text{mbox}if\sqrt{k_{\max }}>\frac{⟨k^2⟩}{⟨k⟩}{\ln }^2(N),\end{array}$$

$$K_c^{\mathrm{QMF}}\simeq \frac{2}{\pi g(0)}\times \{\begin{array}{ll}\frac{⟨k⟩}{⟨k^2⟩}& \text{mbox}if2<\gamma <5/2,\\ \frac{1}{\sqrt{k_{\max }}}& \text{mbox}if\gamma >5/2.\end{array}$$

$${\tilde{A}}_{ij}=\frac{k_i{k}_j}{N⟨k⟩},$$

$${\dot{\theta }}_i(t)={\omega }_i+\frac{K{k}_i}{N⟨k⟩}\sum _j{k}_j\sin ({\theta }_j-{\theta }_i).$$

$$r{e}^{\mathrm{i}\psi (t)}=\frac{1}{N⟨k⟩}\sum _{j=1}^N{k}_j{e}^{\mathrm{i}{\theta }_j(t)}.$$

$$K>K_c^{\mathrm{HMF}}=\frac{2}{\pi g(0)}\frac{⟨k⟩}{⟨k^2⟩},$$

$$\frac{i}{N}-\frac{1}{2N}={\int }_{-\infty }^{{\omega }_i}g(\omega )d\omega ,$$

$${\omega }_i=\Delta \tan \left[\frac{i\pi }{N}-\frac{(N+1)\pi }{2N}\right],i=1,\dots ,N.$$

$$\chi =N(⟨R^2{⟩}_t-⟨R{⟩}_t^2),$$

$${\chi }_r=N\frac{({⟨R^2⟩}_t-{⟨R⟩}_t^2)}{⟨R{⟩}_t},$$

$${\dot{\theta }}_i={\omega }_i+\frac{K}{{\mathcal{N}}_i}\sum _{j=1}^N{A}_{ij}\sin ({\theta }_j-{\theta }_i),$$

$$K_c=\frac{2}{\pi g(0)}\frac{⟨k⟩}{⟨k^2⟩}{k}_{\max }\sim N^{\frac{\gamma }{2}-1},$$

Full text

I Introduction

Synchronization processes are pervasively observed in a wide range of physical, chemical, technological and biological systems Pikovsky et al. (2003). 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 proposed by Kuramoto Acebrón et al. (2005), which in the last decade was extensively investigated when the oscillators are placed on complex networks (see Arenas et al. (2008); Rodrigues et al. (2016) and references therein). A key question addressed in these studies is how the heterogeneous connectivity 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 dynamics has been studied in many scenarios: from homogeneous and unclustered networks to heterogeneous and modular ones, in addition to variations of phase oscillator models including correlations between intrinsic dynamics and local topology Arenas et al. (2008); Rodrigues et al. (2016). 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<\gamma\leq 3$. Another important question concerns 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 Moreno and Pacheco (2004). There, the authors remarked on the supposed finite magnitude of the critical coupling and highlighted the apparent contrast of the Kuramoto 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 Ichinomiya (2004); Lee (2005) 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 Boccaletti et al. (2006). However, later investigations reported significant deviations between predictions of mean-field theories and numerical simulations Restrepo et al. (2005), casting further doubts on the validity of the classical result on the nonexistence of a synchronization threshold Rodrigues et al. (2016); Arenas et al. (2008); Dorogovtsev (2010).

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 Moreno and Pacheco (2004); Ichinomiya (2004); Lee (2005); Restrepo et al. (2005) – 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. Rodrigues et al. (2016); Hong et al. (2007, 2013); Um et al. (2014)) have investigated finite-size effects of the dynamics and reported excellent agreement between simulation and mean-field theories. However, most of those analyses have focused on Erdős-Rènyi (ER) random graphs and SF networks with degree exponent $\gamma>3$, situations in which the critical coupling is expected to be finite (according to the heterogeneous degree mean-field approximations) Rodrigues et al. (2016). Of particular interest is a recent contribution Juhász et al. (2019), 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 Arenas et al. (2008); Boccaletti et al. (2006). Very recently, Yook and Kim Yook and Kim (2018) performed a thorough comparison between the classical Kuramoto order parameter and the order parameter accounting 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 Yook and Kim (2018) how the transition point behaves as the network size increases. Therefore, the question regarding whether or not there is a well defined critical coupling for the onset of synchronization in SF networks with $2<\gamma\leq 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 comparison 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\times 10^{6}$ nodes. We show that the critical coupling predicted by both the HMF and the QMF agrees with the values measured in numerical experiments for networks with power-law exponent $\gamma\leq 3$, hence providing stronger evidence that the critical coupling of such systems vanishes in the thermodynamic 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 thermodynamic limit for $\gamma>3$, in contrast to results obtained in the context of epidemic dynamics Ferreira et al. (2012); Mata and Ferreira (2013).

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 dependences 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 heterogeneous 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 Arenas et al. (2008); Rodrigues et al. (2016)

[TABLE]

where $\theta_{i}$ and $\omega_{i}$ are the phase and natural frequency of the $i$-th oscillator, respectively; $K$ is the coupling strength, and $\mathbf{A}$ is the adjacency matrix, with $A_{ij}=1$ if nodes $i$ and $j$ are connected, and [math] otherwise.

In order to assess the overall synchrony of an ensemble of oscillators, Kuramoto Acebrón et al. (2005) introduced the order parameter

[TABLE]

where $R$ and $\psi$ are the magnitude and phase of the centroid associated to the $N$ points $e^{\mathrm{i}\theta_{j}(t)}$ in the complex plane, respectively. If phases are uniformly distributed over $[0,2\pi]$, it follows that $R\approx 0$, whereas $R\approx 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 Acebrón et al. (2005) cannot be directly extended to the network case. The reason for this relies in the fact that Eqs. 1 are not exactly decoupled by a global order parameter. Instead, an exact decoupling is only achieved by defining local order parameters as

[TABLE]

which leads to

[TABLE]

In this paper, we consider the oscillators frequencies $\omega_{i}$ to be distributed according to a smooth and unimodal distribution $g(\omega)$ centered at $\omega=0$. By inserting the fixed point solution $(\dot{\theta}_{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 Restrepo et al. (2005); Rodrigues et al. (2016)

[TABLE]

where $\Lambda_{\max}$ is the largest eigenvalue of $\mathbf{A}$. The latter result was first derived in Restrepo et al. (2005), with what the authors called perturbation 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 Castellano and Pastor-Satorras (2010). To gain further insights on the predictions of Eq. 5 to the dynamics on SF networks, we recall the result Chung et al. (2004)

[TABLE]

where $k_{\max}$ is the maximum degree or the network. In uncorrelated SF networks, $k_{\max}$ scales as $k_{\max}\sim N^{1/2}$ if $2<\gamma\leq 3$, and $k_{\max}\sim N^{1/(\gamma-1)}$, for $\gamma>3$. By noticing further that $\langle k^{2}\rangle/\langle k\rangle\sim k_{\max}^{3-\gamma}\ll\sqrt{k_{\max}}$, we then estimate Castellano and Pastor-Satorras (2010)

[TABLE]

Therefore, according to the QMF approach, the critical coupling $K_{c}$ should vanish in the thermodynamic limit as $k_{\max}$ diverges, even if $\langle k^{2}\rangle$ remains finite (i.e., the case when $\gamma>3$).

Another way of modeling synchronization processes on networks is by virtue of the annealed network approximation Rodrigues et al. (2016). It consists of replacing the elements of the adjacency matrix $A_{ij}$ by its ensemble average $\tilde{A}_{ij}$, which corresponds to the probability that two nodes, $i$ and $j$, are connected in the configuration model; that is,

[TABLE]

where $k_{i}$ is the degree of node $i$. Substituting Eq. 8 into Eq. 1 yields

[TABLE]

The previous equation motivates the definition of the following order parameter

[TABLE]

Equation 8 is equivalent to the so-called heterogeneous degree mean-field approximation (HMF) Rodrigues et al. (2016) 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 $\mathbf{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 parameters $k_{i}$ and $\omega_{i}$.

By decoupling Eqs. 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 Rodrigues et al. (2016); Ichinomiya (2004)

[TABLE]

where $\langle k^{n}\rangle=\sum_{k}k^{n}P(k)$ is the $n-$th moment of the degree distribution $P(k)$.

As previously mentioned, it has been recently shown Yook and Kim (2018) 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 $\mathcal{O}(N)$; however, discrepancies between $R$ and $r$ are also likely to emerge for networks with power-law exponent $\gamma>3$ Yook and Kim (2018). 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 Hong et al. (2013)

[TABLE]

which for the Lorentzian distribution $g(\omega)=\frac{\Delta}{\pi}\frac{1}{\omega^{2}+\Delta^{2}}$ yields

[TABLE]

In this way, we generate a set of quasi-uniformly spaced frequencies, removing, thus, the disorder introduced by different realizations of frequencies Hong et al. (2013).

III Critical coupling of uncorrelated scale-free networks

All networks analyzed in this section were generated following the uncorrelated configuration model (UCM) Catanzaro et al. (2005) considering a power-law degree distribution $P(k)\sim k^{-\gamma}$ with the cutoff $k_{\max}\sim N^{1/2}$ for $\gamma\leq 3$, and $k_{\max}\sim N^{1/(\gamma-1)}$, for $\gamma>3$. Furthermore, in order to avoid sample-sample fluctuations on $k_{\max}$, for each value of $N$, we fixed $k_{\max}=\langle k_{\max}\rangle$ across the network realizations.

Typically, the critical coupling strength of finite networks can be estimated numerically via detecting the divergent peak of the susceptibility

[TABLE]

where $\langle\cdots\rangle_{t}$ denotes a temporal average. However, we employ the modified susceptibility defined as Ferreira et al. (2012)

[TABLE]

As with the definition in Eq. 14, the modified susceptibility also exhibits a peak at $K=K_{c}$. Nonetheless, analogous forms of $\chi_{r}$ have been shown to be better suited to detect transition points in epidemic spreading and contact processes in networks with diverging $\langle k^{2}\rangle$ Ferreira et al. (2012); Mata and Ferreira (2015); Mata et al. (2014). 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 $\gamma=2.25$, the critical points estimated via Eq. 15 are in better agreement with HMF and QMF theories than that estimated via Eq. 14, especially for low values of $N$. Similar results are found for different values of $\gamma$. Thus, we henceforth detect the critical points via $\chi_{r}$.

Let us now analyze how the mean-field theories perform in comparison with simulations for the different regimes of $\gamma$. First, for $\gamma<5/2$, as discussed in the previous section, both HMF and QMF predict a vanishing $K_{c}$, which should scale with $\langle k\rangle/\langle k^{2}\rangle$. Indeed, as it is seen in Fig. 2(a), for $\gamma=2.25$, both theories predict quite accurately the onset of synchronization.

Discrepancies between the approximations appear when $\gamma>5/2$. To be precise, in this regime, HMF yields $K_{c}\sim\langle k\rangle/\langle k^{2}\rangle$, while QMF gives $K_{c}\sim k_{\max}^{-1/2}$. As depicted in Fig. 2(b), the mean-field theories provide a satisfactory approximation of the synchronization thresholds for networks with $\gamma=2.7$. Note that, although QMF contains in its formulation the whole information about the network topology, it performs slightly worse than HMF (see inset). Similar dependences with the system size are found for epidemic thresholds in SF networks with $5/2<\gamma<3$ Mata and Ferreira (2013); Ferreira et al. (2012).

For $\gamma=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 vanishing $K_{c}$. Furthermore, it is interesting to point out discrepancies between synchronization and the epidemic spreading described by the susceptible-infected-susceptible (SIS) model Harris (1974); Diekmann and Heesterbeek (2000) regarding the dependence on the system size for $\gamma>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 $\gamma>3$ Mata and Ferreira (2013); Ferreira et al. (2012). In fact, Chatterjee and Durrett Chatterjee and Durrett (2009) proved rigorously that, for uncorrelated random networks with a power-law degree distribution $P(k)\sim k^{-\gamma}$ with any $\gamma$, the SIS model presents an unstable absorbing phase in the thermodynamic limit, resulting in a null epidemic threshold. Afterwards, Boguñá et al. Boguñá et al. (2013) physically interpreted this proof with a semi-analytical approach taking into account a long-range reinfection mechanism between hubs and found a vanishing epidemic threshold including for $\gamma>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 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 Mata and Ferreira (2013); Ferreira et al. (2012) and the emergence of Griffiths effects Cota et al. (2016)–in this case, QMF captures the peak associated to the activation of the largest hub in the network Mata and Ferreira (2015). Surprisingly, simulations with networks with $N=10^{7}$ (not shown here) did not reveal signs of multiple peaks in susceptibility curves of Kuramoto oscillators. However, with the aim of understanding the nature of the threshold in epidemic models on uncorrelated random networks, recent works Mata et al. (2014); Ferreira et al. (2016); Cota et al. (2018) showed that different epidemic models such as, for instance, susceptible-infected-recovered-susceptibility Anderson and May (1992), contact process Castellano and Pastor-Satorras (2006), the generalized SIS model with weighted infection rates Karsai et al. (2006) and other alterations of the SIS model Cota et al. (2018), 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 phenomenon of the Kuramoto oscillators.

At last, the results in Fig. 2 point to a different scenario as the one in Restrepo et al. (2005) regarding the accuracy of mean-field theories. More precisely, in Fig. 2 we see that HMF exhibits an excellent agreement with numerical simulations for $\gamma=2.25$ and $2.7$, and a qualitative agreement with the scaling with $N$ for $\gamma=3.5$. Conversely, Ref. Restrepo et al. (2005) found that HMF agrees best with the numerical results obtained for $\gamma>3$, while significantly deviating from simulations for $2<\gamma<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}\sim N$) and the relative small size of the SF networks ($N\sim 10^{3}$) considered in Restrepo et al. (2005).

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 Arenas et al. (2008); Rodrigues et al. (2016), 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 normalization 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 Arenas et al. (2008), and (iii) the second term in the r.h.s. of Eq. 1 diverges in the thermodynamic limit for networks in which the maximum degree is not bounded when $N\rightarrow\infty$. In this 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

[TABLE]

where $\mathcal{N}_{i}$ is the normalization constant of node $i$. Reasonable choices for $\mathcal{N}_{i}$ would then be quantities that are related to the network topology. One of these prescriptions discussed in previous works is $\mathcal{N}_{i}=k_{\max}$ $\forall i$ Arenas et al. (2008); Rodrigues et al. (2016). 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 networks with $\mathcal{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)\sim k^{-\gamma}$, $K_{c}$ is rescaled to

[TABLE]

explaning why the onset of synchronization increases in this case. Curiously, this conclusion is not evident from early works Arenas et al. (2008); Rodrigues et al. (2016). Phenomenologically, the previous result can be understood by noticing that this normalization also makes the coupling $\frac{K}{\mathcal{N}_{i}}\rightarrow 0$ when $N\rightarrow\infty$ 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}\rightarrow\infty$ when $N\rightarrow\infty$.

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 $\mathcal{N}=\langle k^{2}\rangle/\langle k\rangle$. Indeed, observing the corresponding result in Fig. 3, it looks like as if the problem of finding the appropriate normalization has been solved: as $N\rightarrow\infty$, $K_{c}$ converges to $2/\pi g(0)$, which is the same value encountered for the fully connected graph. Nevertheless, while this choice leads to a finite onset of synchronization in the thermodynamic limit – and moreover sets the same $K_{c}$ for all heterogeneous networks – , it imposes a vanishing coupling strength to low connected nodes. In other words, for infinitely large networks, such nodes will require an infinite $K$ in order to lock in synchrony with mean-field. This effect is evident in Fig. 4, where we see that even though SF networks with $\mathcal{N}=\langle k^{2}\rangle/\langle k\rangle$ have similar $K_{c}$, the level of synchronization for $K>K_{c}$ decreases as $N$ gets larger. The solution for the problems of having a vanishing critical coupling and a diverging normalization for poorly connected nodes seems to be the choice $\mathcal{N}_{i}=k_{i}$. However, this comes with the price of washing out from the dynamics effects that are intrinsic to the network topology, since the normalization acts as an average over the contribution of the nearest neighbors Arenas et al. (2008); Rodrigues et al. (2016). For instance, as seen in Fig. 4, the synchronization curves of large ER and SF networks become qualitatively equivalent under $\mathcal{N}_{i}=k_{i}$. What would then be an appropriate normalization for the coupling function for the Kuramoto dynamics on networks? It turns out that, if differences between the network structures must be highlighted, the most natural choice is the classical normalization $\mathcal{N}=1$, at the expense of having a vanishing $K_{c}$ for large networks with diverging $\langle k^{2}\rangle$.

V Conclusion

In this paper, we have analyzed the onset of synchronization of Kuramoto phase oscillators in scale-free networks. First, we revisited a long-standing problem about the dynamics of Kuramoto oscillators coupled in heterogeneous topologies, namely, whether there is a nonzero critical coupling for the onset of synchronization in SF networks. The debate around this question arose already in the early days of network science and, although there has been a substantial amount of work on network dynamics, this question has been seldom addressed in the last years. For SF networks with $\gamma<3$, our extensive simulations showed that QMF and HMF solutions turned out to be equivalent in estimating the critical coupling strength. Specifically, both theories predicted a vanishing critical value for the onset of synchronization. For $\gamma>3$, on the other hand, the HMF correctly predicted the finite threshold in the thermodynamic limit, whereas the QMF erroneously estimated a decaying critical coupling.

We pointed out that this is a noticeable difference between the critical properties of synchronization of phase oscillators and the SIS dynamics. In particular, concerning the latter dynamics, experimental evidence reveals Ferreira et al. (2012); Mata and Ferreira (2013) that critical thresholds of SF networks with $\gamma>3$ decrease as $N\rightarrow\infty$, although with a different scaling as yielded by QMF. Nevertheless, in this regime of $\gamma$, the latter approximation estimates correctly secondary peaks in susceptibility curves associated to localization effects due to the epidemic activation of the largest hub in the network – a phenomenon for which we have not observed a counterpart in the synchronization dynamics of large SF networks. Synchronization thresholds, on the other hand, present the same behavior as observed in most dynamical processes that exhibit a phase transition from active to inactive states, such as contact process and SIRS model Ferreira et al. (2016); Cota et al. (2018). In fact, the phase transition observed in the Kuramoto model is a standard phase transition associated to a collective phenomenon (i.e., the activation of the entire network), whereas the phase transition in the SIS model is related to a mutual reinfection of hubs. Therefore, although synchronization and SIS epidemic thresholds behave similarly in SF networks with $\gamma\leq 3$, fundamental differences between the critical properties of these dynamics emerge for $\gamma>3$. Future investigations should test if other types of localization effects and multiple transitions de Arruda et al. (2017); Colomer-de Simón and Boguñá (2014) can be detectable in ensembles of phase oscillators.

In addition, we have also discussed the influence of different normalization choices in the long term dynamics of large networks. We pointed out that choices previously considered to be suitable for the dynamics of highly heterogeneous networks actually have major drawbacks. Regarding the normalization by the maximum degree, while it prevents the hub’s interaction to diverge, it yields a critical coupling that grows with the network’s size – a fact that remained unnoticed in previous works. Normalizing the coupling function by the ratio $\langle k^{2}\rangle/\langle k\rangle$ circumvents the inconvenient of a size-dependent threshold. However, as for the case $\mathcal{N}_{i}=k_{\max}$, it ends up establishing a divergent normalization for low connected nodes, which requires them to have an infinite coupling strength to lock with the mean-field in the thermodynamic limit. The alternative then to these drawbacks is the choice $\mathcal{N}_{i}=k_{i}$, which, as discussed here and in previous texts Arenas et al. (2008); Rodrigues et al. (2016), removes from the dynamics the contribution from network topology, making networks with significantly different structures to exhibit similar synchronous dynamics. While this prescription could be appropriate in cases in which the focus of the analysis is not on the role played by the network topology in the dynamics (e.g. Jörg et al. (2014)), it seems counterintuitive that large heterogeneous networks should synchronize similarly as homogeneous ones. This scenario, therefore, points back to the conclusion that the most natural choice for the interaction between oscillators is the classical unnormalized couplings.

Acknowledgments

TP acknowledges FAPESP (Grants No. 2016/23827-6 and 2018/15589-3). BM thanks CAPES for financial support. FAR acknowledges the Leverhulme Trust, CNPq (Grant No. 305940/2010-4) and FAPESP (Grants No. 2016/25682-5 and grants 2013/07375-0) for the financial support given to this research. AM acknowledges FAPEMIG and CNPq (Grant No. 423185/2018-7). YM acknowledges partial support from the Government of Aragón, Spain through grant E36-17R (FENOL), by MINECO and FEDER funds (FIS2017-87519-P) and by Intesa Sanpaolo Innovation Center. This research was carried out using the computational resources of the Center for Mathematical Sciences Applied to Industry (CeMEAI) funded by FAPESP (grant 2013/07375-0). The funders had no role in study design, data collection and analysis, or preparation of the manuscript.