Universal phase transitions to synchronization in Kuramoto-like models with heterogeneous coupling

We reveal a class of universal phase transitions to synchronization in Kuramoto-like models with both in- and out-coupling heterogeneity. By analogy with metastable states, an oscillatory state occurs as a high-order coherent phase accompanying explosive synchronization in the system. The critical points of synchronization transition and the stationary solutions are obtained analytically, by the use of mean-field theory. In particular, the stable conditions for the emergence of phase-locked states are determined analytically, consistently with the analysis based on the Ott–Antonsen manifold. We demonstrate that the in- or out-coupling heterogeneity have influence on both the dynamical properties (eigen’spectrum) and the synchronizability of the system.


Introduction
Synchronization of interacting units is a macroscopic self-organized behavior, and is ubiquitous in nature and human society. Examples range from physics, chemistry, engineering to social science [1], such as the flashing of fireflies, power grids, Josephon junction arrays, and neurons in human brain, etc [2,3]. Synchronous evolutions are furthermore at the basis of cooperative functioning of many biological and man'made systems, and therefore revealing the mechanisms behind the setting and maintenance of synchronization is a very important issue [4].
The first prototypic model for investigating synchronization in coupled-phase oscillator was proposed by Winfree [5]. Later, Kuramoto refined the model, and made it mathematically tractable [6]. The classical Kuramoto model describes an ensemble of phase oscillators interacting via a sinusoidal function of the phase differences, and shows that a large population of oscillators can in fact synchronize despite the diversity in the individuals' natural frequencies. Such a pioneering work has inspired extensive studies on synchronization in the following four decades [7,8]. Typically, synchronization in Kuramoto-like models turns out to be a continuous process, i.e. it follows a second-order phase transition. Recently, however, first-order-like synchronization transitions (the so called explosive synchronization) has also been revealed [9].
Besides the diversity in natural frequencies, heterogeneity in the coupling strength is also an important characteristic that may dominate in many real systems [10,11], inducing the emergence of various coherent states in the route to synchronization, such as cardio-respiratory synchronization [12,13] and network physiology studying networks between different organ systems [14][15][16]. Furthermore, non-local coupling is known as relevant for the formation of chimera states [17], and Kuramoto models with both positive and negative coupling can exhibit π states, travelling-wave states, standing-wave states, and glass states, among others [18,19]. On the other hand, network structure leads to nontrivial phase transition, such as explosive synchronization [9] and cluster synchronization [20]. Recently, [21] shows that a metastable state appears in the route to the travelling-wave state near the hybrid phase transition, which has potential applications in the recover of human consciousness.
In this paper, we provide a comprehensive analysis for a system of coupled-phase oscillators by considering both in-and out-coupling heterogeneity. We characterize the dynamical behaviors in the system, and obtain the critical points of synchronization transition analytically. Remarkably, we show that the in-and out-coupling schemes have no influence on the critical point where the incoherent state becomes unstable, whereas synchronization ability is changed dramatically. A rigorous stability analysis for the phase-locked state is provided from different levels. We give a general condition for the stability of the phase-locked state based on the properties of eigen spectrum and the associated eigenvectors which are clarified. We reveal a universal route toward synchronization, i.e. from the incoherent state to the oscillatory state, then from the oscillatory state to the phase-locked state. Such route occurs even for purely attractive coupling, which is a signature that is different from the metastable state appearing near the hybrid phase transition.
The paper is organized as follow: in section 2 we introduce the dynamical model and the mean-field theory. In section 3 we focus on the stability analysis of the phase-locked state and its eigen-spectrum properties in detail. Section 4 provides a discussion about the tiered phase transition to the oscillatory state. Finally, conclusions are drawn in the last section.

Dynamical model and mean-field theory
Let us consider an ensemble of N intercating pase oscillators. The model we consider here has the form where θ i is the instantaneous phase of the ith oscillator, dot denotes temporal derivative, and ω i is the natural frequency of the ith oscillator, which is taken from a specific distribution g(ω). N is the number of oscillators, and K ij is the coupling value between the pair of oscillators i and j. In particular, K ij =κ/N (κ>0) stands for the global and uniform coupling that is typically used in the Kuramoto model. In order to account for heterogeneity in the coupling, we set , which leads to a frequency-weighted Kuramoto model with κ>0 [22,23]. In the present study, we consider a uniform distribution where Δ is the half-width of the distribution [24].
For the sake of convenience, we pay attention to the case of heterogeneity in the out-coupling, where k w = | | K N ij j , i.e. the weighted factor is inside the summation in equation (1). To characterize the collective behavior in the model is defined as the order parameter of the system, and 0 denotes the incoherent state, where all oscillators are desynchronized and distribute uniformly on the unit circle. By applying the mean-field theory or linear stability analysis, one can obtain the critical coupling κ c,1 for the onset of synchronization as where Ω c is the critical mean-field frequency (the imaginary part of the eigenvalue), which satisfies the balanced principal-value integral equation [25] ò w w w For the uniform distribution in equation ( . This suggests that the critical points for the emergence of synchronization in both out-coupling and in-coupling are the same [22]. Next, we seek for solutions of the coherent state. The dynamical equation (1) can be rewritten in the meanfield form as In the long time limit (  ¥ t ) the stationary solution of equation (7) indicates that D is independent of time and Θ=Ωt. Considering the symmetry of the system, we assume a rotating frequency Ω=0. Hence, the steady state for equation (7) is for the phase-locked oscillators. The drifting oscillators (with w k > | | D i ) cannot be entrained by the mean-field and they rotate non-uniformly on the unit circle. It has been proved that the drifting oscillators have no contribution to the order parameter Letting N ℓ be the number of phase-locked oscillators, there are ℓ 2 N possibilities of configuration q { } cos i except for those D<0 in equation (10). However, if we perform an independent perturbation δθ i for each locked phase θ i , while keeping the others invariant, the corresponding eigenvalue equation is dq ldq The necessary condition for stable coherent state D>0 requires all q > cos 0 i . As  ¥ N , the order parameter D should be re-defined in its integral form Throughout the paper, we set Δ=1 for κ D<1. The solutions are D=3κ −2 (κ>3) and R=3π/(4κ) which correspond to partial synchronization. For κD>1, all oscillators become phase-locked. Setting κD=q, . Figure 1(a) plots the schematic function f (q) (  q 1), and one can easily see that there is no intersection between 1/κ and f (q) when κ is small enough. As κ increases, the line 1/κ is tangent to the curve f (q) at (q c , f (q c )) and there are two solutions in the interval k Î -- The situation for the in-coupling is relatively simple due to the symmetry. The system splits into two groups once κR>1, where sign of the natural frequency. The coherent solutions are , and they are independent on N and g(ω).

Stability of the phase-locked states
The above analysis shows that there are multi-branch solutions of the coherent state. Therefore, determining their stability becomes an important theoretical task. As we know, the incoherent state is neutrally stable due to the absence of eigenvalues for the linear operator in the region κ<κ c,1 , while its stability is further determined by the resonant pole calculated through analytical continuation [25]. The natural question is how the eigenspectrum for the phase-locked state appears. In the following, we provide a detailed stability analysis of the coherent state for finite size N. We emphasize that the results can be straightforwardly extended to the case  ¥ N . Choosing a pair of symmetric oscillators q p n , with opposite natural frequencies w p n , , the governing equation where D 1 represents the local order parameter formed by q p n , , , and D 0 =D−D 1 stands for the order parameter formed by the other N−2 oscillators. Considering a special perturbation dq p n , that makes D 0 invariant, the linearized equation for dq p n , is governed by dq The first eigenvalue of the Jacobian matrix is l k q = -D cos p 1 , which is always negative. The corresponding eigenvector satisfies δθ p =δθ n representing identical perturbations, and it is stable. The second eigenvalue is , and the associated eigenvector satisfies δθ p =−δθ n characterizing the perturbations in the reverse direction. As we will see, the two eigenvectors are the basis of the eigendirections for general perturbations. When κ D<1, λ 2 >0 in the limit w k  | | D p , and this implies that oscillators q p n , near the phase-locked boundary ( w k  | | D p ) would first lose their stabilities under reverse perturbation (see figure 1(b)). Hence, one can conclude that the configuration where synchronous and drifting oscillators coexist is unstable.
A rigorous analysis needs to account for all perturbations δθ i (i=1, 2, K, N) [26], then the linearized equation for phase-locked state is governed by where δθ is a vector (δθ 1 , δθ 2 , K, δθ N ) and J is the Jacobian matrix with entries = , which implies that (1, 1, K,1) is a trivial eigenvector corresponding to the eigenvalue λ 1 =0. This property is associated with the rotation invariance of the Kuramoto model equation (1). Since J is a real and symmetric matrix, the stability condition for the phase-locked state is to have all other N−1 eigenvalues  l 0 i (i=2, 3, K, N). All detailed information about the eigen-spectrum comes from the characteristic equation of J. To compute the characteristic polynomial, we express where C and W are the diagonal matrix with entries q cos i and w | | i along the diagonal, respectively, and M is a real-symmetric matrix with elements Apart from the poles −κDc i , the functions Q c (λ) and Q s (λ) are strictly decreasing, and so the polynomial p (λ)=0 must have exactly two roots (Q c (λ)=1 and Q s (λ)=1) between any two consecutive poles. Obviously, Q c (0)=1 corresponds to the rotation invariance and the only remaining eigenvalue Q s (λ)=1 determines the stability of the phase-locked state, so the stable condition (λ<0) is Q s (0)<1 which yields In the continuous limit, equation (23) is equivalent to the following inequality which requires ¢ < ( ) f q 0. Therefore, the eigen-spectrum of J is made up of three pieces (see figures 2 and 3). The first part, widely distributed in the interval of poles, merges into the continuous spectrum ( The second part is a trivial eigenvalue λ=0 due to the rational symmetry of the Kuramoto model, and the third part is a separated eigenvalue located in between those two, which turns out to be the nontrivial part of the discrete spectrum as  ¥ N . For qä [1, q c To better understand the property of the eigen-spectrum, we concentrate on the continuous limit  ¥ N . In this way, the state of the system is described by a real density function ρ(θ, ω, t) which satisfies the Figure 2. The eigen-spectrum of the Jacobian matrix J (see equation (17)). Panels (a) and (c) refer to q=1.01 (q<q c ), for which the largest eigenvalue λ is positive; panels (b) and (d) refer to q=2.0 (q>q c ), for which the largest eigenvalue λ vanishes. The frequencies are evenly spaced: w = -+ - normalization condition. Meanwhile, ρ(θ, ω, t) is 2π-periodic with respect to θ and can be expanded into Fourier series as 2 , e , 25 n n n i where the Fourier coefficients α 0 =1 and * a a = -n n . Considering the Ott-Antonsen ansatz, ρ(θ, ω, t) takes the form of Poisson kernel α n =α n [27][28][29][30]. Equation (2) is equivalent to the continuous equation Apart from the incoherent state α(ω)=0, the phase-locked state a w = q w (8) and (9)) is another fixed point of equation (26). Linearizing equation (26) around this steady state, one gets da wda k a da wa d lda Combining the expression of D with equation (27), the self-consistent equation for λ is which is consistent with Q s (λ)=1 in the finite dimensional case. Note that l k w = --( ) D 2 2 just corresponds to continuous spectrum of the linear operator, and the discrete eigenvalue is only determined by equation (29). On the one hand, when κD is large enough, λ must be negative to balance equation (29) (stable solution). On the other hand, when κD is small (k w  { } D max ), λ should be positive to avoid singularity of equation (29) (unstable solution). As λ is a continuous function of κD, the critical point κ c and D c for the emergence of phase-locked state can be obtained by imposing l  + 0 . Figure 3. The eigen-spectrum of the Jacobian matrix ¢ J (see equation (30)). Panels (a) and (c) refer to q=1.1 (q<q c ), the largest eigenvalue being λ>0. Panels (b) and (d) refer instead to q=3.0 (q>q c ), the largest eigenvalue being λ=0. The frequencies are evenly spaced: w = -+ - The stability analysis for in-coupling case follows a similar procedure, where the Jacobian matrix ¢ J has the following elements The eigenvalues of ¢ J are the roots of The structure of the eigen-spectrum has the same form as equation (22). There It is easily checked that the stable condition holds for θ(R + )<π/4, and it is violated since θ(R − ) > π/4. We emphasize that in contrast to the out-coupling case, the discrete eigenvalue for R + is empty as  ¥ N . Since the existence interval k qw -( | | ) R cos , 0 1 for λ(R + ) gets shorter and shorter, w , the interval vanishes and λ is absent.
Let us move now to describe the eigenvectors associated to the eigenvalues. The natural frequencies are reordered as and w w = -i i 2 2 1 . For an even N, R N can be split into two subspace V even and V odd , such that x ä V even , . Thus, dim(V even )=dim(V odd )=N/2, andV V even odd , consequently, the vector c ä V even and s ä V odd . If x ä V even , then for the out-coupling case If one defines y=Wx and y ä V even , then y satisfies My c y c.  Figure 4 illustrates the phase diagram of the system by direct numerical simulation with N=100 000. The blue solid lines are the theoretical predictions indicating stability (q>q c ). The green dashed lines represent a branch of unstable solution (q<q c ). Both the in-and out-coupling schemes exhibit a similar route to synchronization, namely, a tiered phase transition appears from the incoherent state to the oscillatory state, then to the phaselocked state. For the incoherent state, all oscillators are disordered and distribute uniformly on the unit circle, and two coupling schemes show the same critical point k p = 4 2 c,1 . Interestingly, these two coupling schemes reveal different properties of the coherent state in spite of  ¥ N . The critical point κ c,2 for the emergence of phase-locked state is postponed for the out-coupling case. In this sense, we can say that the whole synchronization ability in the out-coupling case is weaker than the in-coupling case (see figure 4(d)). The differences of synchronization ability between these two cases may provide significant strategies for synchronization optimization in real system. Figure 5 plots several typical features of the oscillatory state for both cases. In contrast to the metastable state reported in [21], where the system (the order parameter) experiences a long-lasting fluctuation near the critical point of the travelling-wave state (a signature of a hybrid phase transition), here we find that the oscillatory state behaves as a stable periodic-T vibration of the order parameters Z 1,2 without any fluctuation (figures 5(e) and (f)). In addition, the oscillatory state takes place in a particular regime (κ c,1 <κ<κ c,2 ) [31 -35], where the incoherent state loses its stability and phase-locked state is empty. As a type of time-dependent clustering state, a number of phases appear locked according to their natural frequencies (see figures 5(a) and (b)). Consequently, the effective frequencies averaged in long time limit correlate in a way that they converge to a common value w p á ñ =  T 2 i [36] (figures 5(c) and (d)). It should be pointed out that the instantaneous velocities of all oscillators within each group are distinct, which leads to the oscillation of the order parameters [37]. We argue that the onset of such state is related to the Hopf bifurcation in the phase space characterized by nonzero Ω c , and it disappears at κ c,2 accompanying saddle-node bifurcation of the phase-locked state.

Conclusion
In summary, we extended the Kuramoto model, and accounted for both in-and out-coupling heterogeneity. We reveal a type of universal phase transition from the incoherent state to the phase-locked state via an oscillatory  The left panels refer to the out-coupling case (κ=2.5) and the right panels refer to the in-coupling case (κ=1.9). Panels (a) and (b) report the instantaneous phases θ i versus ω i , while panels (c) and (d) report the effective frequencies w á ñ i versus ω i . Panels (e) and (f) are time series of the order parameters Z 1,2 , where Z 1,2 are purely real. Pink represents Z 1 and green represents Z 2 . state without any repulsive coupling. The critical points for the occurrence of the coherent states and stable solutions of the phase-locked states are obtained analytically, and the predictions have been perfectly supported by numerical simulations. A simple stable condition for the formation of the phase-locked state is provided in terms of matrix analysis theory, and in full consistency with the Ott-Antonsen analysis in the limit  ¥ N . This work not only enhances our understanding of dynamical phase transitions in coupled oscillators with general heterogenous coupling, but also sheds light on synchronization control and optimization.