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 heterogenous couplings.

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 couplings 31,32 . 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 in 33 ) to the model is conducted, which reveals the local bifurcation mechanism of the incoherent state near the critical point 34 . 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 model 28,30 , in which the dynamics of phase oscillators are governed by the following equations where K denotes the coupling strength, and ω i is the natural frequency of the ith oscillator. Without loss of generality, the natural frequencies {ω i } 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 c . 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 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  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 ρ θ ω θ ( , , ) t d 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).
is a normalization constant. It is easy to obtain that for drifting oscillators 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, Equation (16) is called as the phase balance equation 35,36 . 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 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 c for the onset of synchronization reads where Ω c is the critical mean-field frequency. Thus, following the analysis of Eq. (19), we can conclude that Ω c = 0 means K c → ∞, which is not supported by numerical simulation. By Taylor expansion of Eq. (16), we find that Ω c satisfies the following balance equation c 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 Ω c , and the critical coupling strength K c 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 simulations 30 .
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 i 0 2 and the velocity is given by  From Eq. (21), it is easy to verify that the order parameter η(t) is the integral of Z 1 (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 ≡ Z 0 n for = , ,  n 1 2 , 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 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 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., 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 axis 37 . 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 c ) the discrete spectrum is empty, but as K increases, discrete eigenvalues emerge with real part x > 0 for K > K c . Imposing the critical condition x → 0 + for Eq. (30), once again we obtain the critical coupling strength as where y j are determined by the Eq. (31) with the limit x → 0 + . Evidently, Ω c 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 j means that we choose the jth root y j which makes the product ( ) y g y is maximal, so that K c 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 physics 38 Here C 0 is a constant related to the initial value and it is convenient to set C 0 = 1. Equation (35) represents a closed form for the dependence of η p (t) on the coupling strength K. Unfortunately, it is difficult to get the expression of η p (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, η | ( )| t p . 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)], ( ) = ( − )/ R t t t 2 1 cos 2 for the triangle distribution [ Fig. 1(d)], and R(t) = e −t 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, ( ) ∝ δ R t e t 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. 2 is linearly unstable 30 . 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 b = 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 c , the incoherent state becomes unstable with the emergence of a pair of complex conjugated eigenvalues λ c = ± iΩ c ; 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 0 30 . Thus, if K c > 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 c . 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 theory 34 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 ρ θ ω ( , , ) t obeys the Fokker-Planck equation where D is the strength of noise. Similarly, imposing the small perturbation to the incoherent state, i.e., ρ µ θ ω = + ( , , ) π t 1 2 , we obtain the following perturbed equation . 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 c . As a result, a traveling wave solution and a standing wave solution emerge above K c . 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 distribution [40][41][42] . 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 c the incoherent state loses its stability. Meanwhile, non-stationary R(t) emerges with a hysteresis near K c [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.

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. , ω ∈ (− , ) 1 1 . 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.