Synchronization of phase oscillators with frequency-weighted coupling.

Recently, the first-order synchronization transition has been studied in systems of coupled phase oscillators. In this paper, we propose a framework to investigate the synchronization in the frequency-weighted Kuramoto model with all-to-all couplings. A rigorous mean-field analysis is implemented to predict the possible steady states. Furthermore, a detailed linear stability analysis proves that the incoherent state is only neutrally stable below the synchronization threshold. Nevertheless, interestingly, the amplitude of the order parameter decays exponentially (at least for short time) in this regime, resembling the Landau damping effect in plasma physics. Moreover, the explicit expression for the critical coupling strength is determined by both the mean-field method and linear operator theory. The mechanism of bifurcation for the incoherent state near the critical point is further revealed by the amplitude expansion theory, which shows that the oscillating standing wave state could also occur in this model for certain frequency distributions. Our theoretical analysis and numerical results are consistent with each other, which can help us understand the synchronization transition in general networks with heterogeneous couplings.

Synchronization in dynamical systems of coupled oscillators is one important issue in the frontier of nonlinear dynamics and complex systems. This study provides insights for understanding the collective behaviors in many fields, such as the power grids, the flashing of fireflies, the rhythm of pacemaker cells of the heart, and even some social phenomena1234. Theoretically, the classical Kuramoto model with its generalizations turn out to be paradigms for synchronization problem, which have inspired a wealth of works because of both their simplicity for mathematical treatment and their relevance to practice56. A latest review of Kuramoto model in complex network is presented in7.

Recently, the first-order synchronization transition in networked Kuramoto-like oscillators has attracted much attention. For instance, it has been shown that the positive correlation of frequency-degree in the scale-free network, or a particular realization of frequency distribution of oscillators in an all-to-all network, or certain special couplings among oscillators, etc, would cause a discontinuous phase transition to synchronization891011121314151617181920212223242526. In particular, our recent work27 analytically investigated the mechanism of the first-order phase transition on star network. We revealed that the structural relationship between the incoherent state and the synchronous state leads to different routes to the transition of synchronization. Furthermore, it has been shown that the generalized Kuramoto model with frequency-weighted coupling can generate first-order synchronization transition in general networks2829. In ref. 30, the critical coupling strength for both forward and backward transitions, as well as the stability of the two-cluster coherent state, have been further determined analytically for typical frequency distributions.

In this paper, we present a complete framework to investigate the synchronization in the frequency-weighted Kuramoto model with all-to-all couplings. It includes three separate analyses from different angles, which together presents a global picture for our understanding of the synchronization in the model. First, a rigorous mean-field analysis is implemented where the possible steady states of the model are predicted, such as the incoherent state, the two-cluster synchronous state, and the traveling wave state. It is shown that in this model the mean-field frequency is not necessarily equal to 0. Instead, the non-vanishing mean-field frequency plays a crucial role in determining the critical coupling strength. Second, a detailed linear stability analysis of the incoherent state is performed. Also, the exact expression for the critical coupling strength is obtained, which is consistent with the results of the mean-field analysis, and keeps the same form for general heterogenous couplings3132. Furthermore, it has been proved that the linearized operator has no discrete spectrum when the coupling strength is below a threshold. This implies that in this model the incoherent state is only neutrally stable below the synchronization threshold. Interestingly, numerical simulations demonstrate that in this neutrally stable regime predicted by the linear theory, the perturbed order parameter decays to zero and its decaying envelope follows exponential form for short time. Finally, a nonlinear center-manifold reduction (see the recent development of this theory in33) to the model is conducted, which reveals the local bifurcation mechanism of the incoherent state near the critical point34. As expected, the non-stationary standing wave state could also exist in this model with certain frequency distributions. Extensive numerical simulations have been carried out to verify our theoretical analyses. In the following, we report our main results, both theoretically and numerically.

Results

The mean-field theory

We start by considering the frequency-weighted Kuramoto model2830, in which the dynamics of phase oscillators are governed by the following equations

where K denotes the coupling strength, and ω is the natural frequency of the ith oscillator. Without loss of generality, the natural frequencies {ω} satisfy certain density function g(ω) that is assumed to be symmetric and centered at 0 throughout the paper. The most important characteristic of this generalized Kuramoto model is to introduce a frequency weight to the coupling, which leads to heterogeneous interactions in networks. Equation (1) exhibits a transition to synchronization as the coupling strength K increases above a critical threshold K. Typically, the collective behavior in Eq. (1) can be characterized by the order parameter

Here, η is the average complex amplitude of all oscillators on the unit circle. R is the magnitude of complex amplitude characterizing the level of synchronization, and Θ is the phase of the mean-field corresponding to the peak of the distribution of phases. When K is small enough, R(t) ≈ 0 characterizing the incoherent state in which the phases of oscillators are almost randomly distributed. As K increases, usually a cluster of phase-locked oscillators appear, characterized by an order parameter 0 < R(t) < 1. Then the system is in the synchronous (coherent) state where the phase-locked oscillators coexist with the phase-drifting ones.

One central issue in the study of synchronization is to identify all the possible asymptotic coherent states of the system as the coupling strength K increases. To this end, the self-consistence method turns out to be effective. In the following, we conduct theoretical analysis to Eq. (1) based on this method.

Substituting Eq. (2) into Eq. (1), we obtain the dynamical equation of the mean-field form

We assume that the mean-field phase Θ rotates uniformly with frequency Ω, i.e., Θ(t) = Ωt + Θ(0). Without loss of generality, Θ(0) = 0 after an appropriate time shift. In the rotating frame with frequency Ω, we introduce the phase difference

and Eq. (3) can be transformed into

in the rotating frame. It should be pointed out that Ω = 0 when the coupling form is uniform and g(ω) is even and unimodal. However, for more general cases, any asymmetry of the system, such as asymmetric frequency distribution, or asymmetric coupling function (phase lag or time delay), or asymmetric coupling strength (heterogeneous coupling or time varying coupling), etc, will cause Ω ≠ 0. Therefore, there is no guarantee that Ω in Eq. (5) is necessarily equal to 0. Actually, there are several macroscopic characteristic frequencies for all oscillators, for example, the average frequency of oscillators , the mean-field frequency Ω, and the mean-ensemble frequency (or effective frequency) 3536.

Since we are interested in the steady coherent states of the system, Eq. (5) should be discussed in two situations corresponding to the phase-locked oscillators and the drifting ones, respectively. On the one hand, when , Eq. (5) has solution of fixed point, i.e., , which leads to

corresponding to the phase-locked oscillators entrained by the mean-field. On the other hand, for those drifting oscillators, . Taking into account both the phase-locked and the drifting oscillators, the order parameter in Eq. (2) can be rewritten as

where H(x) is the Heaviside function. In the thermodynamical limit N → ∞, the summation over the frequency should be replaced by the integration. As a result, the contribution of the phase-looked oscillators to the order parameter R reads

In contrast to the phase-locked oscillators, the drifting oscillators could not be entrained by the mean-field. In the thermodynamic limit N → ∞, Eq. (1) is equivalent to the following continuity equation as a consequence of the conservation of the number of oscillators, i.e.,

Here gives the fraction of oscillators of natural frequency ω which lie between θ and (θ + dθ) at time t with the appropriate normalization condition

and 2π period in θ. Then the stationary distribution of the drifting oscillators in the rotating frame could be obtained explicitly as (∂ρ/∂t = 0).

where

is a normalization constant. It is easy to obtain that for drifting oscillators

and

Equation (13) shows that the drifting oscillators have no contributions to the real part of R. However, their contributions to the imaginary part of R should not be neglected. Substituting Eqs (8, 9, 10, 11, 12, 13, 14) into Eq. (7), the closed form of self-consistence equations take the following form. For the real part of R,

and for the imaginary part of it

Equation (16) is called as the phase balance equation3536. Equations (15 and 16) together provide a closed equation for the dependence of the magnitude R and the frequency Ω of the mean field on K.

We notice that Ω = 0 is always a trivial solution of Eq. (16), but it may not be the only solution. There may be more than one value for Ω that satisfies the phase balance equation. Considering g(ω) = g(−ω), a pair of Ω with opposite sign might emerge. Define α = KR ≥ 0 and x = (ω − Ω)/ω, Ω ≠ 0, Eq. (15) can be expressed as

For the case of α > 1, to avoid divergency of Eq. (17), the only choice is Ω = 0, and Eq. (15) is reduced as

which exactly corresponds to the two-cluster synchronous state in ref. 30. For the case of α < 1, the solution of Eqs (16) and (17) can be solved numerically, which corresponds to the traveling wave state. In such a state, the mean-field amplitude R keeps stationary, whereas the mean-field frequency Ω differs from the mean of the natural frequencies. In particular, in the limit case α → 0+, the critical coupling K for the onset of synchronization reads

where Ω is the critical mean-field frequency. Thus, following the analysis of Eq. (19), we can conclude that Ω = 0 means K → ∞, which is not supported by numerical simulation. By Taylor expansion of Eq. (16), we find that Ω satisfies the following balance equation

where the symbol P means the principal-value integration within the real line. As an example, Table 1 shows the balance eq. (20), the critical mean-field frequency Ω, and the critical coupling strength K with respect to different frequency distributions where the θ is the Heaviside function, a and g equals to 1 respectively. All these analytical results were supported by the previous numerical simulations30.

Table 1

Summary of the frequency distributions, the balance equations, the critical mean-field frequencies Ω , and the critical coupling strength K .

frequency distribution	balance equation	Ωc	Kc
		0,
		0, ±0.40	2.65
		0, ±0.47	2.34
		0, ±γ	4

From top to bottom: the uniform distribution, the triangle distribution, the parabolic distribution, and the Lorentzian distribution.

The linear stability analysis

The analysis of the mean-field theory above reveals four macroscopic steady states, including the incoherent state (R = 0), the traveling wave state (Ω ≠ 0, 0 < α < 1), and the two-cluster synchronous states (Ω = 0, α > 1), respectively. However, a thorough stability analysis to every possible solution has not been performed due to the limitation of the mean-field method. In the following, we conduct a detailed linear stability analysis to the incoherent state because its instability usually signals the onset of synchronization. In particular, we will show that the critical coupling strength for synchronization can be alternatively obtained via the linear operator theory.

The continuum limit of the order parameter, i.e., Eq. (2), is rewritten as

and the velocity is given by

where “*” denotes the complex conjugate of η(t). Let

be the nth Fourier coefficient of , then Z0(t, ω) = 1 and Z satisfies the following differential equations

From Eq. (21), it is easy to verify that the order parameter η(t) is the integral of Z1(t, ω) with the frequency density function g(ω), and the higher Fourier harmonics have no contribution to the order parameter. Since the incoherent state corresponds to the trivial solution for , to study its stability we can consider the evolution of a perturbation away from the incoherent state. In this spirit, Eq. (24) can be linearized around the origin as

and

Here is the operator defined as

Where q(ω) is a function in the weighted Lebesgue space, and is defined as

From Eq. (26), it is obvious that the higher Fourier harmonics are neutrally stable to perturbation. Hence, the key is to study the spectrum of Eq. (25). Following ref. 37, Eq. (25) has continuous spectrum on the whole imaginary axis. For the discrete spectrum, we assume that the perturbation of the first Fourier coefficient has the form . Then the self-consistent eigenvalue equation Eq. (25) for the operator takes the form

where λ is the complex eigenvalue of except for those points iω. Notice that Eq. (29) relates implicitly the coupling strength K with the eigenvalue λ. Since the real part of the eigenvalue λ determines the stability of the incoherent state, we rewrite Eq. (29) into two equations by letting λ = x + iy, i.e.,

and

From Eq. (30), we see that x, i.e., the real part of λ, can never be negative, otherwise K < 0, which makes no physical sense. Hence, the incoherent state in model (1) cannot be linearly stable. In fact, it is neutrally stable due to the existence of continuous spectrum on the imaginary axis37. Furthermore, if the coupling strength K > 0 but sufficiently small, we have proved that the eigenvalue λ does not exist (the details are included in the supplementary material).

The analysis above reveals that the linearized operator has continuous spectrum iω lying on the whole imaginary axis (with real part equals to 0) for all K, and it may also have discrete spectrum (eigenvalues) depending on K. When K is small (K < K) the discrete spectrum is empty, but as K increases, discrete eigenvalues emerge with real part x > 0 for K > K. Imposing the critical condition x → 0+ for Eq. (30), once again we obtain the critical coupling strength as

where y are determined by the Eq. (31) with the limit x → 0+. Evidently, Ω is the imaginary part of the eigenvalues of operator at the boundary of stability. Generally Eq. (31) may have more than one root with x → 0+. sup means that we choose the jth root y which makes the product is maximal, so that K corresponds to the foremost critical point for the onset of synchronization.

According to the above linear stability analysis, the incoherent state of model (1) is only neutrally stable below the synchronization threshold. However, interestingly, we find that in this regime the perturbed order parameter η(t) actually decays to zero in the long time limit (t → ∞). This phenomenon was first found in the classical Kuramoto model, and was revealed to be analogous to the famous Landau damping in plasma physics38. To investigate the Landau damping effect in our model, we rewrite Eq. (25) as

which could be solved explicitly, i.e.,

Substituting Eq. (34) into the expression of η(t), we obtain the perturbed order parameter as

Here C0 is a constant related to the initial value and it is convenient to set C0 = 1. Equation (35) represents a closed form for the dependence of η(t) on the coupling strength K. Unfortunately, it is difficult to get the expression of η(t) analytically for the present model. However, we still can obtain useful information via direct numerical simulations. In Fig. 1, the numerical solutions of Eq. (35) are illustrated for different values of K and typical frequency distributions g(ω). Generally, we observe the decay of R(t), i.e, . Depending on g(ω) and K, the scenarios turn out to be different. We notice that when the coupling constant is absent, Eq. (35) can be solved analytically. Specifically, for example, when K = 0, R(t) = sin t/t for the uniform distribution [Fig. 1(a)], for the triangle distribution [Fig. 1(d)], and R(t) = e− for the Lorentzian distribution [Fig. 1(g)]. In these cases, the decay phenomena strongly depend on the form of g(ω). However, with the increasing of K the situation changes. It is observed that the order parameter decays in a way with significant oscillation. Nevertheless, its envelope follows the form of exponential decay, namely, for a short time, where δ is the decay exponent. While the general dynamical mechanism of this decay is still an open issue, ref. 39 pointed out that this exponential decay of order parameter in the neutrally stable regime is closely related to the resonance pole on the left-half complex plane, and the decaying rate δ is the real part of it.

Figure 1

Different scenarios of the decay of R(t) with different frequency distributions and coupling strength below the critical threshold (K < K).

(a–c) Uniform distribution , . K = 0, 1.6, 1.78, respectively. (d–f) Triangle distribution , . K = 0, 2.2, 2.6, respectively. (g–i) Lorentzian distribution , . K = 0, 3.7, 3.9, respectively. The red solid lines denote the fitted curves of the envelopes which all satisfy the exponential form e.

For the two-cluster synchronous states Eq. (18), previous analysis has shown that is linearly stable, and is linearly unstable30. For the traveling wave state in the range (0 < α < 1), its stability can only be studied through numerical simulations. We have conducted extensive simulations by choosing different initial conditions for phase oscillators. We even specially choose a proper initial condition to make the system artificially locate onto the traveling wave state. It is found that in all these cases the system evolves to the two-cluster synchronous state as long as K > K = 2 (the subscript b denotes the backward transition point). Thus the numerical results give evidence that the traveling wave state predicted by the mean-field theory turns out to be unstable in the current model.

The bifurcation analysis

The above stability analysis leads to the conclusion that near the finite critical coupling K, the incoherent state becomes unstable with the emergence of a pair of complex conjugated eigenvalues λ = ±iΩ; meanwhile the traveling wave solution is unstable. Moreover, due to the absolute coupling, Eq. (1) always has the two-cluster synchronous solution [Eq. (18)] when K > 2. This is independent of the specific form of g(ω) as long as g(ω) is symmetric and centered at 030. Thus, if K > 2, the first-order synchronization transition would take place. However, the mechanism underlying the instability of the incoherent state is still unclear, for example, the bifurcation type and the local stability of the traveling wave solution near K. These information is crucial for us to get a global picture of the synchronization transition in the dynamical system. Generally, the dynamic behavior near the critical point can be investigated through the local bifurcation theory. For this purpose, we refer to the framework of nonlinear analysis developed in ref. 34 to reveal the local bifurcation type for the incoherent state of model (1).

The main idea of the theory34 is that when the perturbed equation of the incoherent state satisfies O(2) symmetry, the center manifold reduction could be applied to obtain the amplitude equations for both steady state and limit cycle. Moreover, in order to avoid dealing with the continuous spectrum, the Gaussian white noise is added and eventually the noise magnitude is extended to zero for all calculations. Following this treatment, now the evolution of the density function obeys the Fokker-Planck equation

where D is the strength of noise. Similarly, imposing the small perturbation to the incoherent state, i.e., , we obtain the following perturbed equation

Here, is lineared equation for μ, i.e.,

and is the nonlinear term, i.e.,

where . Then we can derive the normal form of amplitude equation for both the steady state and the Hopf bifurcation in the frequency-weighted model based on the center manifold assumption. Since the complete calculation is tedious, we put the details into the Supplementary Material for interested readers. In the following we only report the main results.

It is found that for the frequency-weighted model, the system undergoes Hopf bifurcation near the critical coupling K. As a result, a traveling wave solution and a standing wave solution emerge above K. The standard standing wave solution consists of two counter-rotating clusters of phase-locked oscillators. Thus the order parameter η(t) plots a limit cycle on the complex plane. Previously, such state has been found in the classical Kuramoto model with symmetric bimodal frequency distribution404142. It should be pointed out that the mean-field theory fails to predict such state due to the fact that neither the distribution function nor R(t) are stationary in any rotating frame for such a state.

As an example to illustrate our results, we focus on the case of uniform distribution g(ω) = 1/2. The exact expression for critical coupling strength is . The nonlinear analysis shows that the bifurcation for the traveling wave solution is supercritical and unstable (which is consistent with the mean-field theory). In addition, the standing wave solution is subcritical. This implies that a hysteresis would occur by taking the high order terms of the amplitude equation into account. Numerical evidence suggests that above K the incoherent state loses its stability. Meanwhile, non-stationary R(t) emerges with a hysteresis near K [branch 3 in Fig. 2]. As the coupling strength increases, and it eventually vanishes at K = 2 via a discontinuous transition (with very small hysteresis loop) to the two-cluster synchronous state. We have also conducted calculations for other typical frequency distributions, such as the triangle, the Lorentzian, and the parabolic. The results show that the bifurcations for the standing wave solution are all subcritical and the traveling wave solution are all unstable locally. Moreover, the direction of bifurcation for the traveling wave state supports the numerical solution of the mean-field equation. It should be pointed out that the stable branch of subcritical bifurcation for both states are not observed numerically. One possible reason for this is that their basins of attraction might be so small in such a high-dimensional phase space that most of the initial conditions eventually lead to the stable two-cluster synchronous state as long as K > 2.

Figure 2

Characterizing various coherent states in the phase diagram.

R [long time average of R(t)] vs. K for the uniform frequency distribution , . Branches 1 − −5 are the incoherent state, the (unstable) traveling wave state predicted by the mean-field theory, the standing wave state, the unstable and the stable two-cluster synchronous states, respectively. The blue and red lines denote the forward and the backward transitions, respectively. In both directions, K is changed adiabatically in simulations. There is a hysteresis region of the standing wave solution within K = 1.725 − 1.8. In the simulations oscillators number N = 50,000, and a fourth-order Runge-Kutta integration method with time step 0.01 is used.

Discussion

To summarize, we investigated the synchronization transition in the frequency-weighted Kuramoto model with all-to-all couplings. Theoretically, mean-field analysis, linear stability analysis, and bifurcation analysis have been carried out to obtain insights. Together with the numerical simulations, our study presented the following main results. First, we predicted the possible steady states in this model, including the incoherent state, the two-cluster synchronous state, the traveling wave state, and the standing wave state. Second, the critical coupling strength for synchronization transition has been obtained analytically. Third, we proved that in this model the incoherent state is only neutrally stable below the synchronization threshold. However, in this regime, the perturbed order parameter decays exponentially to zero for short time. Finally, the amplitude equations near the bifurcation point have been derived based on the center-manifold reduction, which predicted that the non-stationary standing wave state could also exist in this model. This work provided a complete framework to deal with the frequency-weighted Kuramoto model, and the obtained results will enhance our understandings of the first-order synchronization transition in networks.

Additional Information

How to cite this article: Xu, C. et al. Synchronization of phase oscillators with frequency-weighted coupling. Sci. Rep. 6, 21926; doi: 10.1038/srep21926 (2016).