A Arenas et aL / Physics Reports 469(2008)93-153 the equation i= Kr sin(0i). On the other hand, the rest of the oscillators for which oil Kr holds, are drifting around the circle, sometimes accelerating and sometimes rotating at lower frequencies. Demanding some conditions for the stationary distribution of drifting oscillators with frequency ai and phases ei [27], a self-consistent equation for r can be derived as r=Kr/.(cose)g(@o)de where o= Kr sin (0). This equation admits a non-trivial solution, beyond which r >0. Eq. 8)is the Kuramoto mean field expression for the critical coupling at the onset of synchronization. Moreover, near the onset, the order parameter, r, obeys the usual square-root scaling law for mean field models, namely. (K-K) with B= 1/2. Numerical simulations of the model verified these results. The Kuramoto model(KM, from now on)approach to synchronization was a breakthrough for the understanding of synchronization in large populations of oscillators. Even in the simplest case of a mean field interaction, there are still unsolved problems that have resisted any analytical attempt. This is the case, e g. for finite populations of oscillators and some questions regarding global stability results 28 In what follows, we focus on another aspect of the models assumptions, namely that of the connection topology of real systems[14, 15], which usually do not show the all-to-all pattern of interconnections underneath the mean field approach. 3.1.2. Kuramoto model on complex networks To deal with the KM on complex topologies, it is necessary to reformulate eq ( 3)to include the connectivity =+∑qsi-)(=1,…N here oi is the coupling strength between pairs of connected oscillators and ai are the elements of the connectivity matrix. The original Kuramoto model is recovered by letting ai= 1, vi #j(all-to-all)and o k/N, vi, j The first problem when defining the KM in complex networks is how to state the interaction dynamics properly. In contrast with the mean field model, there are several ways to define how the connection topology enters in the governing uations of the dynamics. a good theory for Kuramoto oscillators in complex networks should be phenomenologically evant and provide formulas amenable to rigorous mathematical treatment. Therefore, such a theory should at least at the same time, should remain calculable in the thermodynamic limi k independently of the interaction dynamics, and For the original model, Eq (3), the coupling term on the right hand side of Eq (10 )is an intensive magnitude because the dependence on the size of the system cancels out. This independence on the number of oscillators n is achieved by choosing Ji= K/N. This prescription turns out to be essential for the analysis of the system in the thermodynamic limit N-0o in to become dependent on N. Therefore, in the thermodynamic limit, the coupling term tends to zero except for those nodes with a degree that scales with N. Note that the existence of such nodes is only possible in networks with power-law degree distributions[14,15, but this happens with a very small probability as k- ', with y > 2. In these cases, mean field solutions ndependent of N are recovered, with slight differences in the onset of synchronization of all-to-all and SF networks[31]. A second prescription consists in taking oi k/ki(where ki is the degree of node i)so that oi is a weighted interaction factor that also makes the right hand side of Eq. (10)intensive. This form has been used to solve the paradox o f heterogeneity [32 that states that the heterogeneity in the degree distribution, which often reduces the average distance between nodes, may suppress synchronization in networks of oscillators coupled symmetrically with uniform coupling strength. This result refers to the stability of the fully synchronized state, but not to the dependence of the order parameter on the coupling strength(where partially synchronized and unsynchronized states exist). Besides, the inclusion of weights in the interaction strongly affects the original KM dynamics in complex networks because it can impose a dynamic homogeneity that masks the real topological heterogeneity of the network. The prescription j= K/const, which may seem more appropriate, also causes some conceptual problems because the sum in the right hand side of Eq (10)could eventually diverge in the thermodynamic limit. The constant in the denominator could in principle be any quantity related to the topology such as the average connectivity of the graph. (k), or the maximum degree kmax. Its physical meaning is a re-scaling of the temporal scales involved in the dynamics. However, except for the case of oi =K/kmax, the other possible settings do not avoid the problems when N-o0. On the other hand, for a proper comparison of the results obtained for different complex topologies(e g SF or uniformly random), the global and local measures of coherence should be represented according to their respective time scales. Therefore, given two complex networks A and B with kmax = kA and kmax kB respectively, it follows that to make meaningful comparisons between
98 A. Arenas et al. / Physics Reports 469 (2008) 93–153 the equation ωi = Kr sin(θi). On the other hand, the rest of the oscillators for which |ωi | > Kr holds, are drifting around the circle, sometimes accelerating and sometimes rotating at lower frequencies. Demanding some conditions for the stationary distribution of drifting oscillators with frequency ωi and phases θi [27], a self-consistent equation for r can be derived as r = Kr Z π 2 − π 2 cos2 θ g(ω)dθ , where ω = Kr sin(θ ). This equation admits a non-trivial solution, Kc = 2 πg(0) (8) beyond which r > 0. Eq. (8) is the Kuramoto mean field expression for the critical coupling at the onset of synchronization. Moreover, near the onset, the order parameter, r, obeys the usual square-root scaling law for mean field models, namely, r ∼ (K − Kc ) β (9) with β = 1/2. Numerical simulations of the model verified these results. The Kuramoto model (KM, from now on) approach to synchronization was a breakthrough for the understanding of synchronization in large populations of oscillators. Even in the simplest case of a mean field interaction, there are still unsolved problems that have resisted any analytical attempt. This is the case, e.g. for finite populations of oscillators and some questions regarding global stability results [28]. In what follows, we focus on another aspect of the model’s assumptions, namely that of the connection topology of real systems [14,15], which usually do not show the all-to-all pattern of interconnections underneath the mean field approach. 3.1.2. Kuramoto model on complex networks To deal with the KM on complex topologies, it is necessary to reformulate Eq. (3) to include the connectivity θ˙ i = ωi + X j σijaij sin(θj − θi) (i = 1, . . . , N), (10) where σij is the coupling strength between pairs of connected oscillators and aij are the elements of the connectivity matrix. The original Kuramoto model is recovered by letting aij = 1, ∀i 6= j (all-to-all) and σij = K/N, ∀i, j. The first problem when defining the KM in complex networks is how to state the interaction dynamics properly. In contrast with the mean field model, there are several ways to define how the connection topology enters in the governing equations of the dynamics. A good theory for Kuramoto oscillators in complex networks should be phenomenologically relevant and provide formulas amenable to rigorous mathematical treatment. Therefore, such a theory should at least preserve the essential fact of treating the heterogeneity of the network independently of the interaction dynamics, and at the same time, should remain calculable in the thermodynamic limit. For the original model, Eq. (3), the coupling term on the right hand side of Eq. (10) is an intensive magnitude because the dependence on the size of the system cancels out. This independence on the number of oscillators N is achieved by choosing σij = K/N. This prescription turns out to be essential for the analysis of the system in the thermodynamic limit N → ∞ in the all-to-all case. However, choosing σij = K/N for the governing equations of the KM in a complex network makes them to become dependent on N. Therefore, in the thermodynamic limit, the coupling term tends to zero except for those nodes with a degree that scales with N. Note that the existence of such nodes is only possible in networks with power-law degree distributions [14,15], but this happens with a very small probability as k −γ , with γ > 2. In these cases, mean field solutions independent of N are recovered, with slight differences in the onset of synchronization of all-to-all and SF networks [31]. A second prescription consists in taking σij = K/ki (where ki is the degree of node i) so that σij is a weighted interaction factor that also makes the right hand side of Eq. (10) intensive. This form has been used to solve the paradox of heterogeneity [32] that states that the heterogeneity in the degree distribution, which often reduces the average distance between nodes, may suppress synchronization in networks of oscillators coupled symmetrically with uniform coupling strength. This result refers to the stability of the fully synchronized state, but not to the dependence of the order parameter on the coupling strength (where partially synchronized and unsynchronized states exist). Besides, the inclusion of weights in the interaction strongly affects the original KM dynamics in complex networks because it can impose a dynamic homogeneity that masks the real topological heterogeneity of the network. The prescription σij = K/const, which may seem more appropriate, also causes some conceptual problems because the sum in the right hand side of Eq. (10) could eventually diverge in the thermodynamic limit. The constant in the denominator could in principle be any quantity related to the topology, such as the average connectivity of the graph,hki, or the maximum degree kmax. Its physical meaning is a re-scaling of the temporal scales involved in the dynamics. However, except for the case of σij = K/kmax, the other possible settings do not avoid the problems when N → ∞. On the other hand, for a proper comparison of the results obtained for different complex topologies (e.g. SF or uniformly random), the global and local measures of coherence should be represented according to their respective time scales. Therefore, given two complex networks A and B with kmax = kA and kmax = kB respectively, it follows that to make meaningful comparisons between
A Arenas et aL/Physics Reports 469(2008)93-153 8 0.2 口N=104 ◇N=5x10 Fig. 2. Order parameter r(Eq(4))as a function of a for several BA networks of different sizes. Finite size scaling analysis shows that the onset of synchronization takes place at a critical value a =0.05(1). The inset is a zoom around ae From 35] observables, the equations of motion(Eq (10))should refer to the same time scales, i.e. d= ka/ka ke/kB= o. with this formulation in mind, Eq (10)reduces to B=m+∑q(-B)(=,……N however, that the same order parameter, Eq (4), is often used to describe the coherence of the synchronized state ZAaB ndently of the specific topology of the network. this allows us to study the dynamics of Eq (11on different topologies pare the results, and properly inspect the interplay between topology and dynamics in what concerns synchronization As we shall see, there are also several ways to define the order parameter that characterizes the global dynamics of the system, some of which were introduced to allow for analytical treatments at the onset of synchronization. We advance. 3.1.3. Onset of synchronization in complex networks Studies on synchronization in complex topologies where each node is considered to be a Kuramoto oscillator, were first reported for WS networks [33, 34 and Ba graphs 35, 36]. These works are mainly numerical explorations of the onset of synchronization, their main goal being the characterization of the critical coupling beyond which groups of nodes beating coherently first appear. In [34, the authors considered oscillators with intrinsic frequencies distributed according to a aussian distribution with unit variance arranged in a wS network with varying rewiring probability, p, and explored how the order parameter, Eq. (4), changes upon addition of long-range links. Moreover, they assumed a normalized coupling der variatio.K/(k), where (k)is the average degree of the graph. Numerical integration of the equations of motion(10) The results confirm that networks obtained from a regular ring by just rewiring a tiny fraction of links(p> 0)can be synchronized with a finite K. Moreover, in contrast with the arguments provided in [34 we notice that their results had been obtained for a fixed average degree and thus the Kuramoto's critical coupling cannot be recovered by simply taking p-1, which produces a random ER graph with a fixed minimum connectivity. This limit is recovered by letting(k)increase. Actually, numerical simulations of the same model in 33]showed that the Kuramoto limit is approached when the average connectivity grows. In[35] the same problem in BA networks is considered. The natural frequencies and the initial values of e; were randomly drawn from a uniform distribution in the interval (-1/2, 1/2)and(-, T), respectively. The global dynamics of the system, Eq (11). turns out to be qualitatively the same as for the original KM as shown in Fig. 2, where the dependence of the order parameter Eq (4)with o is shown for several system sizes. The existence of a critical point for the KM on SF networks came as a surprise. Admittedly, this is one of the few cases in which a dynamical process shows a critical behavior when the substrate is described by a power-law connectivity distribution with an exponent y <3[14, 15, 37). In principle it could be a finite size effect, but it turned out from numerical simulations that this was not the case. To determine the exact value of oc, one can make use of standard finite-size scaling analysis. At least two complementary strategies have been reported The first one allows bounding the critical point and is computationally more expensive. Consider a network of size N, for which no synchronization is attained below ac, wherer(t) decays to a small residual value of size O(1//N). Then, the critical point may be found by examining the N-dependence of r(o, N). In the sub-critical regime(o 0), the stationary value of r falls off as N-1 while for o >0, the order parameter reaches a stationary value as N oo(though still with O(1//N) fluctuations). Therefore, plots of r versus N allow us to
A. Arenas et al. / Physics Reports 469 (2008) 93–153 99 Fig. 2. Order parameter r (Eq. (4)) as a function of σ for several BA networks of different sizes. Finite size scaling analysis shows that the onset of synchronization takes place at a critical value σc = 0.05(1). The inset is a zoom around σc . From [35]. observables, the equations of motion (Eq. (10)) should refer to the same time scales, i.e. σij = KA/kA = KB/kB = σ. With this formulation in mind, Eq. (10) reduces to θ˙ i = ωi + σ X j aij sin(θj − θi) (i = 1, . . . , N), (11) independently of the specific topology of the network. This allows us to study the dynamics of Eq.(11) on different topologies, compare the results, and properly inspect the interplay between topology and dynamics in what concerns synchronization. As we shall see, there are also several ways to define the order parameter that characterizes the global dynamics of the system, some of which were introduced to allow for analytical treatments at the onset of synchronization. We advance, however, that the same order parameter, Eq. (4), is often used to describe the coherence of the synchronized state. 3.1.3. Onset of synchronization in complex networks Studies on synchronization in complex topologies where each node is considered to be a Kuramoto oscillator, were first reported for WS networks [33,34] and BA graphs [35,36]. These works are mainly numerical explorations of the onset of synchronization, their main goal being the characterization of the critical coupling beyond which groups of nodes beating coherently first appear. In [34], the authors considered oscillators with intrinsic frequencies distributed according to a Gaussian distribution with unit variance arranged in a WS network with varying rewiring probability, p, and explored how the order parameter, Eq. (4), changes upon addition of long-range links. Moreover, they assumed a normalized coupling strength σij = K/hki, where hki is the average degree of the graph. Numerical integration of the equations of motion (10) under variation of p shows that collective synchronization emerges even for very small values of the rewiring probability. The results confirm that networks obtained from a regular ring by just rewiring a tiny fraction of links (p & 0) can be synchronized with a finite K. Moreover, in contrast with the arguments provided in [34], we notice that their results had been obtained for a fixed average degree and thus the Kuramoto’s critical coupling cannot be recovered by simply taking p → 1, which produces a random ER graph with a fixed minimum connectivity. This limit is recovered by letting hkiincrease. Actually, numerical simulations of the same model in [33] showed that the Kuramoto limit is approached when the average connectivity grows. In [35] the same problem in BA networks is considered. The natural frequencies and the initial values of θi were randomly drawn from a uniform distribution in the interval(−1/2, 1/2) and (−π , π ), respectively. The global dynamics of the system, Eq. (11), turns out to be qualitatively the same as for the original KM as shown in Fig. 2, where the dependence of the order parameter Eq. (4) with σ is shown for several system sizes. The existence of a critical point for the KM on SF networks came as a surprise. Admittedly, this is one of the few cases in which a dynamical process shows a critical behavior when the substrate is described by a power-law connectivity distribution with an exponent γ ≤ 3 [14,15,37]. In principle it could be a finite size effect, but it turned out from numerical simulations that this was not the case. To determine the exact value of σc , one can make use of standard finite-size scaling analysis. At least two complementary strategies have been reported. The first one allows bounding the critical point and is computationally more expensive. Consider a network of size N, for which no synchronization is attained below σc , where r(t) decays to a small residual value of size O(1/ √ N). Then, the critical point may be found by examining the N-dependence of r(σ , N). In the sub-critical regime (σ < σc ), the stationary value of r falls off as N −1/2 , while for σ > σc , the order parameter reaches a stationary value as N → ∞ (though still with O(1/ √ N) fluctuations). Therefore, plots of r versus N allow us to
A Arenas et aL / Physics Reports 469(2008)93-153 locate the critical point o. Alternatively, a more accurate approach can be adopted Assume the scaling form for the order parameter[38: r=N (N(o-ac) where f(x) is a universal scaling function bounded as x ->+oo and a and v are two critical exponents to be determined. ofo for g s oc, the value of the functionf is independent of N, the estimation of oc can be done by plotting@ras a function a by-product, the method also allows us to calculate the two scaling exponents a and v, the latter can be obtained from the In[(dr/do )lo ]=(v-a)In n+const, once a is computed Following these scaling procedures, it was estimated a value for the critical coupling strength o 0.05(1)[35, 39, 40] Moreover, r M(o-o) when approaching the critical point from above with B= 0.46(2)indicating that the square-root behavior typical of the mean field version of the model(B= 1/2)seems to hold as well for BA networks. Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is worth briefly summarizing other numerical results that have explored how the critical coupling depends on other topological features of the underlying SF graph. Recent results have shed light on the influence of the topology of the local interactions on the route to, and the onset of, synchronization. In particular, the authors in[41-43 explored the Kuramoto dynamics on networks in which the degree distribution is kept fixed, while the clustering coefficient(C)and the average path length (e)of the graph change. The results suggest that the onset of synchronization is mainly determined by C, namely, networks with a high clustering coefficient promote synchronization at lower values of the coupling strength On the other hand, when the coupling is increased beyond the critical point, the effect of e dominates over C and the phase diagram is smoothed out(a sort of stretching), delaying the appearance of the fully synchronized state as the average shortest path length In a series of recent papers [31, 44-48 the onset of synchronization in large networks of coupled oscillators has been analyzed from a theoretical point of view under certain(sometimes strong) assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are available up to now. Moreover, the analytical approaches ponent y s 3, the critical coupling va anishes as n-o. in contrast to numerical simulations on BA model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity strongly hampers analytical solutions of the model Following [31]. consider the system in Eq (11), with a symmetric adjacency matrix i ji. Defining a local order parameter r as re=∑ae) (14) re(.)t stands for a time average, a new global order parameter to measure the macroscopic coherence is readily (15) Now, rewriting Eq (11)as a function of ri, yields, 6=0一 or sin(一中)一ah(t) (16) In Eq (16);(t)= Im(e- L, ai ((ei)t-e)) depends on time and contains time fluctuations. Assuming the terms in the previous sum to be statistically independent hi (t)is expected to be proportional to ki above the transition, where T ki. Therefore, except very close to the critical point, and assuming that the number of connections of each node is large enough(ki > 1 as to be able to neglect the time fluctuations entering h, i. e, hi <ri). the equation describing the dynamics of node i can be reduced to (17) I The reader can find the extension of the forthcoming formalism to directed networks in(44). This obviously restricts the range of real networks to which the approximation can be applied
100 A. Arenas et al. / Physics Reports 469 (2008) 93–153 locate the critical point σc . Alternatively, a more accurate approach can be adopted. Assume the scaling form for the order parameter [38]: r = N −α f(N ν (σ − σc )), (12) where f(x) is a universal scaling function bounded as x → ±∞ and α and ν are two critical exponents to be determined. Since at σ = σc , the value of the function f is independent of N, the estimation of σc can be done by plotting N α r as a function of σ for various sizes N and then finding the value of α that gives a well-defined crossing point, the critical coupling σc . As a by-product, the method also allows us to calculate the two scaling exponents α and ν, the latter can be obtained from the equality ln[(dr/dσ )|σc ] = (ν − α)ln N + const, (13) once α is computed. Following these scaling procedures, it was estimated a value for the critical coupling strength σc = 0.05(1) [35,39,40]. Moreover, r ∼ (σ − σc ) β when approaching the critical point from above with β = 0.46(2) indicating that the square-root behavior typical of the mean field version of the model (β = 1/2) seems to hold as well for BA networks. Before turning our attention to some theoretical attempts to tackle the onset of synchronization, it is worth briefly summarizing other numerical results that have explored how the critical coupling depends on other topological features of the underlying SF graph. Recent results have shed light on the influence of the topology of the local interactions on the route to, and the onset of, synchronization. In particular, the authors in [41–43] explored the Kuramoto dynamics on networks in which the degree distribution is kept fixed, while the clustering coefficient (C) and the average path length (`) of the graph change. The results suggest that the onset of synchronization is mainly determined by C, namely, networks with a high clustering coefficient promote synchronization at lower values of the coupling strength. On the other hand, when the coupling is increased beyond the critical point, the effect of ` dominates over C and the phase diagram is smoothed out (a sort of stretching), delaying the appearance of the fully synchronized state as the average shortest path length increases. In a series of recent papers [31,44–48], the onset of synchronization in large networks of coupled oscillators has been analyzed from a theoretical point of view under certain (sometimes strong) assumptions. Despite these efforts no exact analytical results for the KM on general complex networks are available up to now. Moreover, the analytical approaches predict that for uncorrelated SF networks with an exponent γ ≤ 3, the critical coupling vanishes as N → ∞, in contrast to numerical simulations on BA model networks. It appears that the strong heterogeneity of real networks and the finite average connectivity strongly hampers analytical solutions of the model. Following [31], consider the system in Eq. (11), with a symmetric1 adjacency matrix aij = aji. Defining a local order parameter ri as rie iφi = XN j=1 aijhe iθjit, (14) where h· · ·it stands for a time average, a new global order parameter to measure the macroscopic coherence is readily introduced as r = PN i=1 ri PN i=1 ki . (15) Now, rewriting Eq. (11) as a function of ri , yields, θ˙ i = ωi − σri sin(θi − φi) − σ hi(t). (16) In Eq. (16), hi(t) = Im{e −iθi PN j=1 aij(he iθjit − e iθj)} depends on time and contains time fluctuations. Assuming the terms in the previous sum to be statistically independent, hi(t) is expected to be proportional to √ ki above the transition, where ri ∼ ki . Therefore, except very close to the critical point, and assuming that the number of connections of each node is large enough2 (ki 1 as to be able to neglect the time fluctuations entering hi , i.e., hi ri), the equation describing the dynamics of node i can be reduced to θ˙ i = ωi − σri sin(θi − φi). (17) 1 The reader can find the extension of the forthcoming formalism to directed networks in [44]. 2 This obviously restricts the range of real networks to which the approximation can be applied
A. Arenas et aL Physics Reports 469(2008)93-153 Next, we look for stationary solutions of Eq (17), i.e. sin (0i -i)=o/or In particular, oscillators whose intrinsic frequency satisfies oil s or become locked. Then, as in the Kuramoto mean field model, there are two contributions(though in this case to the local order parameter), one from locked and the other from drifting oscillators such that (9-中) To move one step further, some assumptions are needed. Consider a graph such that the average degree of nearest neighbors is high(i. e if the neighbors of node i are well-connected ) Then it is reasonable to assume that these nodes are not affected by the intrinsic frequency of i. This is equivalent to assuming solutions (ri, i) that are, in a statistical sense, independent of the natural frequency oi. With this assumption, the second summand in Eq (18)can be neglected. Taking into account that the distribution g(o) is symmetric and centered at 32=0. after some algebra one is left with [31 n=∑0(-1-( 19) The critical coupling o is given by the solution of Eq (19)that yields the smallest o. It can be argued that it is obtained when cos(j -i)=1 in Eq (19), thus which is the main equation of the time average approximation(recall that time fluctuations have been neglected). Note however, that to obtain the critical coupling, one has to know the adjacency matrix as well as the particular values of o; for all i and then solve eq. (20)numerically for the ril. Finally the global order parameter defined in Eq (15)can be computed Even if the underlying graph satisfies the other aforementioned topological constraints, it seems unrealistic to require knowledge of the oils. A further approach, referred to as the frequency distribution approximation can be adopted. According to the assumption that ki >>1 for all i, or equivalently that the number of connections per node is large(a dense graph), one can also consider that the natural frequencies of the neighbors of node i follows the distribution g(o) Then, Eq (20)can be rewritten avoiding the dependence on the particular realization of oil to yield, n-2aC-(m)-业 (21) with x=o/(ori). This equation allows us to readily determine the order parameter r as a function of the network topology ai), the frequency distribution (g(o))and the control parameter(o ) On the other hand, Eq- (20)still does not provide approximating(xar])Ng(O)which is valid for small, but nonzero, values of r. Namely, when x uce explicit expressions for the order parameter and the critical coupling strength. To this end, one introduces a first-order r where K= 2/(g(O)is Kuramoto's critical coupling. Moreover, as the smallest value of o corresponds to oc, it follows that the critical coupling is related to both Ke and the largest eigenvalue amax of the adjacency matrix, yielding Eq (22)states that in complex networks, synchronization is first attained at a value of the coupling strength that inversely depends ong(O)and on the largest eigenvalue Amax of the adjacency matrix. Note that this equation also recovers Kuramoto's result when ai = 1, vi j. since Amax = N-1. It is worth stressing that although this method allows us to calculate oc analytically, it fails to explain why for uncorrelated random SF networks with y 3 and in the thermodynamic limit N-oo, the critical value remains finite. This disagreement comes from the fact that in these SF networks, Amax is proportional to the cutoff of the degree distribution, kmax which in turn scales with the system size. Putting the two dependencies together, one obtainsλmax~kax~Nx→∞asN→∞, thus predicting o=0 in the thermodynamic limit, in contrast to finite size scaling analysis for the critical coupling via numerical solution of the equations of motion. Note, however, that the difference may be due to the use of distinct order parameters. Moreover, even in the case of SF networks with y>3, max still diverges when we take the thermodynamic limit, so that o -0 as well. As we shall see soon, this is not the case when other approaches are adopted, at least for y>3
A. Arenas et al. / Physics Reports 469 (2008) 93–153 101 Next, we look for stationary solutions of Eq. (17), i.e. sin(θi−φi) = ωi/σri . In particular, oscillators whose intrinsic frequency satisfies |ωi | ≤ σri become locked. Then, as in the Kuramoto mean field model, there are two contributions (though in this case to the local order parameter), one from locked and the other from drifting oscillators such that ri = XN j=1 aijhe i(θj−φi ) it = X |ωj |≤σrj aije i(θj−φi ) + X |ωj |>σrj aijhe i(θj−φi ) it . (18) To move one step further, some assumptions are needed. Consider a graph such that the average degree of nearest neighbors is high (i.e. if the neighbors of node i are well-connected). Then it is reasonable to assume that these nodes are not affected by the intrinsic frequency of i. This is equivalent to assuming solutions (ri, φi) that are, in a statistical sense, independent of the natural frequency ωi . With this assumption, the second summand in Eq. (18) can be neglected. Taking into account that the distribution g(ω) is symmetric and centered at Ω = 0, after some algebra one is left with [31] ri = X |ωj |≤σrj aij cos(φj − φi) s 1 − ωj σrj 2 . (19) The critical coupling σc is given by the solution of Eq. (19) that yields the smallest σ. It can be argued that it is obtained when cos(φj − φi) = 1 in Eq. (19), thus ri = X |ωj |≤σrj aijs 1 − ωj σrj 2 , (20) which is the main equation of the time average approximation (recall that time fluctuations have been neglected). Note, however, that to obtain the critical coupling, one has to know the adjacency matrix as well as the particular values of ωi for all i and then solve Eq. (20) numerically for the {ri}. Finally, the global order parameter defined in Eq. (15) can be computed from ri . Even if the underlying graph satisfies the other aforementioned topological constraints, it seems unrealistic to require knowledge of the {ωi}’s. A further approach, referred to as the frequency distribution approximation can be adopted. According to the assumption that ki 1 for all i, or equivalently, that the number of connections per node is large (a dense graph), one can also consider that the natural frequencies of the neighbors of node i follows the distribution g(ω). Then, Eq. (20) can be rewritten avoiding the dependence on the particular realization of {ωi} to yield, ri = X j aij Z σrj −σrj g(ω)s 1 − ω σrj 2 dω = σ X j aijrj Z 1 −1 g(xσrj) p 1 − x 2dx, (21) with x = ω/(σrj). This equation allows us to readily determine the order parameter r as a function of the network topology (aij), the frequency distribution (g(ω)) and the control parameter (σ). On the other hand, Eq. (20) still does not provide explicit expressions for the order parameter and the critical coupling strength. To this end, one introduces a first-order approximation g(xσrj) ≈ g(0) which is valid for small, but nonzero, values of r. Namely, when rj → 0 + r 0 i = σ Kc X j aijr 0 j , where Kc = 2/(πg(0)) is Kuramoto’s critical coupling. Moreover, as the smallest value of σ corresponds to σc , it follows that the critical coupling is related to both Kc and the largest eigenvalue λmax of the adjacency matrix, yielding σc = Kc λmax . (22) Eq. (22) states that in complex networks, synchronization is first attained at a value of the coupling strength that inversely depends on g(0) and on the largest eigenvalue λmax of the adjacency matrix. Note that this equation also recovers Kuramoto’s result when aij = 1, ∀i 6= j, since λmax = N − 1. It is worth stressing that although this method allows us to calculate σc analytically, it fails to explain why for uncorrelated random SF networks with γ ≤ 3 and in the thermodynamic limit N → ∞, the critical value remains finite. This disagreement comes from the fact that in these SF networks, λmax is proportional to the cutoff of the degree distribution, kmax which in turn scales with the system size. Putting the two dependencies together, one obtains λmax ∼ k 1 2 max ∼ N 1 2(γ −1) → ∞as N → ∞, thus predicting σc = 0 in the thermodynamic limit, in contrast to finite size scaling analysis for the critical coupling via numerical solution of the equations of motion. Note, however, that the difference may be due to the use of distinct order parameters. Moreover, even in the case of SF networks with γ > 3, λmax still diverges when we take the thermodynamic limit, so that σc → 0 as well. As we shall see soon, this is not the case when other approaches are adopted, at least for γ > 3
A Arenas et aL / Physics Reports 469(2008)93-153 It is possible to go beyond the latter approximation and to determine the behavior of r near the critical point. In 31] a perturbative approach to higher orders of Eq. (21)is developed, which is valid for relatively homogeneous degree distributions(y>5). They showed that for(o/oc)-1-0+ k where n1=-Tg(O)Kc/16 and n (24) here u is the normalized eigenvector of the adjacency matrix corresponding to Amax and(u)=2i u4 The analytical insights discussed so far can also be reformulated in terms of a mean field approximation [31, 46-48] for complex networks. This approach(valid for large enough( k))considers that every oscillator is influenced by the local field created in its neighborhood, so that T is proportional to the degree of the nodes ki, i.e i ki. Assuming this is the case and introducing the order parameter r through k∠a(e吗) after summing over i and substituting ri rki in Eq (21)we obtain [31 k g(xr)√1-x2d The above relation, Eq (26), was independently derived in [46 who first studied analytically the problem of synchronization in complex networks, though using a different approach. Taking the continuum limit, Eq (26)becomes kP(k)dk=o/kP(k)dk/g(xo rk)v1-x2dx which for r-0+ verifies /k=o/koo一 og(O)T KP(k)dk, which leads to the condition for the onset of synchronization(r>0)as og(o)T/ P(dk> kP(k)dk, that is (0)(k2)-(k2 The mean field result, Eq (29), gives as a surprising result that the critical coupling o in complex networks is nothing else but the one corresponding to the all-to-all topology ke re-scaled by the ratio between the first two moments of the degree coupling strongly depends on the topology of the underlying graph. In particular, ac ->0 when the second momenta e distribution, regardless of the many differences between the patterns of interconnections. Precisely, it states that the critic distribution( k)diverges, which is the case for SF networks with y s 3. Note, that in contrast with the result in Eq (22). for y >3, the coupling strength does not vanish in the thermodynamic limit On the other hand, if the mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing y, the onset of synchronization occurs at smaller values of oc. Interestingly enough, the dependence gathered in Eq (29)has the same functional form for the critical points of other dynamical processes such as percolation and epidemic spreading processes [14, 15, 37]. while this result is still under numerical scrutiny, it would imply that the critical properties of many dynamical processes on complex networks are essentially determined by the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this last claim will be of extreme importance in physics, probably changing many preconceived ideas about the nature dynamical phenomena. The approach holds if the fourth of the degree distribution, ()=/i P(k)kdk remains finite when N00
102 A. Arenas et al. / Physics Reports 469 (2008) 93–153 It is possible to go beyond the latter approximation and to determine the behavior of r near the critical point. In [31] a perturbative approach to higher orders of Eq. (21) is developed, which is valid for relatively homogeneous degree distributions (γ > 5).3 They showed that for (σ /σc ) − 1 ∼ 0 + r 2 = η η1K 2 c σ σc − 1 σ σc 3 , (23) where η1 = −πg 00(0)Kc/16 and η = hui 2λ 2 max Nhki 2 hu 4 i , (24) where u is the normalized eigenvector of the adjacency matrix corresponding to λmax and hu 4 i = PN j u 4 j /N. The analytical insights discussed so far can also be reformulated in terms of a mean field approximation [31,46–48] for complex networks. This approach (valid for large enough hki) considers that every oscillator is influenced by the local field created in its neighborhood, so that ri is proportional to the degree of the nodes ki , i.e., ri ∼ ki . Assuming this is the case and introducing the order parameter r through r = ri ki = 1 ki XN j=1 aijhe iθjit , (25) after summing over i and substituting ri = rki in Eq. (21) we obtain [31] XN j kj = σ XN j k 2 j Z 1 −1 g(xσrkj) p 1 − x 2dx. (26) The above relation, Eq.(26), was independently derived in [46], who first studied analytically the problem of synchronization in complex networks, though using a different approach. Taking the continuum limit, Eq. (26) becomes Z kP(k)dk = σ Z k 2 P(k)dk Z 1 −1 g(xσrk) p 1 − x 2dx, (27) which for r → 0 + verifies Z kP(k)dk = σ Z k 2 P(k)dk Z 1 −1 g(0) p 1 − x 2dx = σg(0)π 2 Z k 2 P(k)dk, (28) which leads to the condition for the onset of synchronization (r > 0) as σg(0)π 2 Z k 2 P(k)dk > Z kP(k)dk, that is, σc = 2 πg(0) hki hk 2 i = Kc hki hk 2 i . (29) The mean field result, Eq. (29), gives as a surprising result that the critical coupling σc in complex networks is nothing else but the one corresponding to the all-to-all topology Kc re-scaled by the ratio between the first two moments of the degree distribution, regardless of the many differences between the patterns of interconnections. Precisely, it states that the critical coupling strongly depends on the topology of the underlying graph. In particular, σc → 0 when the second moment of the distribution hk 2 i diverges, which is the case for SF networks with γ ≤ 3. Note, that in contrast with the result in Eq. (22), for γ > 3, the coupling strength does not vanish in the thermodynamic limit. On the other hand, if the mean degree is kept fixed and the heterogeneity of the graph is increased by decreasing γ , the onset of synchronization occurs at smaller values of σc . Interestingly enough, the dependence gathered in Eq. (29) has the same functional form for the critical points of other dynamical processes such as percolation and epidemic spreading processes [14,15,37]. While this result is still under numerical scrutiny, it would imply that the critical properties of many dynamical processes on complex networks are essentially determined by the topology of the graph, no matter whether the dynamics is nonlinear or not. The corroboration of this last claim will be of extreme importance in physics, probably changing many preconceived ideas about the nature of dynamical phenomena. 3 The approach holds if the fourth moment of the degree distribution, hk 4 i = R ∞ 1 P(k)k 4dk remains finite when N → ∞