Universal scaling and phase transitions of coupled phase oscillator populations
aa r X i v : . [ n li n . AO ] J u l Universal scaling and phase transitions of coupled phase oscillator populations
Can Xu, ∗ Xuebin Wang, and Per Sebastian Skardal † Institute of Systems Science and College of Information Science and Engineering, Huaqiao University, Xiamen 361021, China Department of Mathematics, Trinity College, Hartford, Connecticut 06106, USA
The Kuramoto model, which serves as a paradigm for investigating synchronization phenomenonof oscillatory system, is known to exhibit second-order, i.e., continuous, phase transitions in themacroscopic order parameter. Here, we generalize a number of classical results by presenting ageneral framework for capturing, analytically, the critical scaling of the order parameter at the onsetof synchronization. Using a self-consistent approach and constructing a characteristic function, weidentify various phase transitions toward synchrony and establish scaling relations describing theasymptotic dependence of the order parameter on coupling strength near the critical point. We findthat the geometric properties of the characteristic function, which depends on the natural frequencydistribution, determines the scaling properties of order parameter above the criticality.
I. INTRODUCTION
Synchronization of large ensembles of interacting unitsis a universal phenomenon that arises in a variety natu-ral processes. Typical examples include the flashing offireflies [1], colonies of yeast cells [2], pacemaker cellsin the heart [3], and neural activity [6]. Additional ex-amples where synchronization plays a prominent role inengineered systems include power grid dynamics [4] andJosephson junction arrays [5]. Exploring the routes thesesystems take towards consensus and self-organization in-variably begins at the onset of synchronization, and thebehavior of the system thereafter remains an active areaof research [7].A remarkable prototype for studying synchronizationissue is the Kuramoto model which displays such syn-chronization transitions [8]. In particular, when the cou-pling strength between oscillators is increased a transi-tion from incoherence to partial synchrony, as measuredby the classical Kuramoto order parameter, takes placein a manner analogous to phase transitions observed inother physical systems [9]. In this analogy, the phasetransition between incoherence and synchronization canbe characterized by a supercritical bifurcation of the or-der parameter. Because of its analytical tractability, theKuramoto model together with its variant versions wasextensively investigated and a great deal of progress havebeen made including the discovering of low-dimensionaldescription for the order parameter [10, 11], and theidentification of various coherent states on the way tosynchronization[12–14]. However, the scaling propertiesof the order parameter at the onset of synchronization,i.e., at the phase transition, remains not completely un-derstood.In this paper, we provide a general framework for re-vealing the inner-relation between the phase transitionand the critical behavior, i.e., scaling properties, of the ∗ [email protected] † [email protected] order parameter in the Kuramoto model. Using a self-consistent approach, various phase transitions towardssynchrony and the associated scaling behaviors of theorder parameter near the onset are established in a uni-versal form. Depending on the properties of the naturalfrequency distributions, we show that the bifurcation ofcollective dynamics involves three manners including su-percritical, marginal, and subcritical that correspond tocontinuous (i.e., second-order), hybrid, and tiered phasetransitions, respectively. In particular, we demonstratethat the structural property of characteristic function inthe vicinity of its critical point determine the scaling ex-ponents as well as asymptotic coefficients of the order pa-rameter with leading and sub-leading terms. Our studyserves as a promising strategy in unveiling the rich scal-ing behaviors and phase transitions in coupled oscillatornetworks.The remainder of this paper is organized as follows.In Sec. II we review the self-consistent equation for thestationary solution of the coupled phase oscillator systemand introduce the characteristic function. In Sec. III weuse the characteristic function to derive a general frame-work for uncovering the scaling behavior of the order pa-rameter at the onset of synchronization. In Sec. IV wegive several examples that illustrate the utility of thisnew framework. Finally, in Sec. V we conclude with adiscussion of our results. II. SELF-CONSISTENT APPROACH
We begin by considering the classical Kuramoto modelin which a system of oscillators evolve according to˙ θ i = ω i + KN N X j =1 sin( θ j − θ i ) , i = 1 , . . . , N, (1)where θ i is the phase of oscillator i , N is the size of sys-tem, ω i is the natural frequency of oscillator i , which isassumed to be drawn from the probability density func-tion g ( ω ) with mean ω , and K > ω of g ( ω ) to zero by entering a suitablerotating reference frame, and we will further assume that g ( ω ) is an even function throughout this paper. As weshall see, the symmetry assumption is essential for thebifurcation and scaling analysis below. We will also con-sider the continuum limit N → ∞ where the populationof oscillators may be described accurately by a densityfunction.Next, a complex order parameter (mean-field) Z ( t ) isneeded to quantified phase transition in Kuramoto modeldefined as Z ( t ) = R ( t ) exp { i Ψ( t ) } = 1 N N X j =1 exp { iθ j ( t ) } . (2)The amplitude R ( t ) ∈ [0 ,
1] measures the coherence of thephases, and Ψ( t ) denotes the average phase of the ensem-ble. For the self-consistent analysis, we are interested inthe steady state, which implies that R is a constant andΨ rotates uniformly with Ψ( t ) = Ω t + Ψ . We howeverremark that a suitable global phase shift of initial condi-tions [under which the dynamics of Eq. (1) are invariant]and the assumed symmetry of g ( ω ) makes Ψ and Ω tobe zero, respectively. As a result, the governing equationfor each oscillator is simplified to the mean-field form˙ θ i = ω i − q sin θ i , (3)with q = KR . Depending on the natural frequencies, os-cillators in Eq. (3) can be divided into those locked by themean-field with sin θ i = ω i /q and cos θ i = p − ( ω i /q ) ( | ω i | ≤ q ) and the unlocked (drifting) ones rotating non-uniformly with period T i = 2 π ( ω i − q ) − / ( | ω i | > q ).Then, the order parameter is given by Z = R = h e iθ i lock + h e iθ i drift , (4)where h·i denotes the average over the populations. Notethat, the symmetric assumptions further reduce to h sin θ i lock = h e iθ i drift = 0 , (5)whether N is infinite or not. Thus, the self-consistentequation R = h cos θ i lock in thermodynamic limit N → ∞ may be written1 K = F ( q ) = 1 q Z | ω | Proceeding with our analysis, we consider a small per-turbation 0 < δq ≪ q c , or equivalently considering K = K c + δK , R = R c + δR with 0 < δK ≪ < δR ≪ 1. We then consider theexpanding of F ( q ) in powers of δq up to the sub-leadingorder, F ( q ) = F ( q c ) + A ( δq ) µ + B ( δq ) v . (8)Based on the property of F ( q ) in the neighborhood of q c , µ and v are either integers or fractions with 0 < µ < v and A and B are the corresponding expansion coeffi-cients. Substituting Eq. (8) into the self-consistent equa-tion, i.e., Eq. (6), we have δK = − AK c δq µ + BK c δq v AK c δq µ + BK c δq v . (9)Furthermore, expanding Eq. (9) up to the sub-leadingorder of δq leads to δK = X ( δq ) α + Y ( δq ) β , (10)with 0 < α < β . The relation between ( A, B, µ, v ) and( X, Y, α, β ) is determined by the specifics of the system,in particular the distribution of natural frequencies g ( ω ).To gain better insight into δq in Eq. (10), we make thefollowing ansatz: δq = X − /α ( δK ) /α + H ( δK ) ǫ . (11)The first term of δq corresponds to a leading solutionwhereas the second term is assumed to be the perturbedsolution. Inserting Eq. (11) into Eq. (10), we get δK = δK + αHX /α ( δK ) ǫ +1 − /α + Y X − β/α ( δK ) β/α + βY HX /α − β/α ( δK ) ǫ + β/α − /α . (12)Since ǫ > /α , then H and ǫ can be determined self-consistently by balancing the sub-leading term whichgives ǫ = β − α + 1 α , (13)and H = − Y X − ( β +1) /α α . (14)The corresponding deviation of the order parameter R close to phase transition point becomes δR = δqK c − q c K c δK. (15)The asymptotic behavior of the order parameter near theonset of synchronization can then be described in thefollowing universal form, δR = P ( δK ) η + Q ( δK ) ξ . (16)For simplicity, we define a parameter array χ =( η, ξ, P, Q ) for the critical exponents and asymptotic co-efficients of the order parameter near K c , in which η and ξ represent the leading and sub-leading scaling ex-ponents, P and Q denote the corresponding asymptoticcoefficients, respectively. For the case of q c = 0, the sec-ond term of right hand side (r.h.s.) of Eq. (15) vanishes.Substituting Eq. (11) into Eq. (15) and combining it withEq. (16), we obtain the parameter array χ = ( α − , ǫ, X − /α K − c , HK − c ) . (17)Whereas q c > 0, Eq. (15) shows that there exists a fixedsub-leading term for the order parameter with the scalingexponent being 1. Hence, χ should be discussed in threedistinct cases according to ǫ in Eq. (11). For ǫ < χ is the same as Eq. (17). For ǫ = 1,the second term of r.h.s. of Eqs. (11) and (15) should beconsidered at the same time, χ = ( α − , , X − /α K − c , HK − c − q c K − c ) . (18)Likewise if ǫ > 1, the second term of r.h.s. of Eqs. (11)becomes a negligible higher-order term compared withEq. (15), then the parameter array is determined as χ = ( α − , , X − /α K − c , − q c K − c ) . (19) IV. EXAMPLES FOR PTS AND SCALINGPROPERTIESA. Continuous phase transitions As a demonstration of this approach, we first consider g ( ω ) to be the following form, g n ( ω ) = nπ sin (cid:16) π n (cid:17) γ n − ω n + γ n , (20)with γ > n = 1 , , . . . [see Fig. 1(a)]. As studiedin [16], g n ( ω ) represents a family of rational functionsin which the low-dimensional dynamics for the long-termevolution of the order parameter were obtained by meansof Ott-Antonsen ansatz. The unimodality of g n ( ω ) en-sures that the incoherent state loses its stability via a -3 -2 -1 0 1 2 30.00.10.20.30.40.5 -1.5-1.0-0.5 0.0 0.5 1.0 1.50.00.20.40.60.81.0-2 -1 0 1 20.00.10.20.30.4 -3 -2 -1 0 1 2 30.00.20.40.60.8 ( ) n g ( )b ( ) m g ( )c ( )a ( ) k g ( )g ( )d FIG. 1: Sketch diagram of frequency distributions, (a-c) arethe natural frequency distributions with γ = 1 and (d) are thevirtual frequency distributions. (a) g n ( ω ) of Eq.(20), n = 1(black solid line), n = 2 (red dot line), n = 3 (green dash line),and n = 10 (blue dash-dot line); (b) g m ( ω ) of Eq.(26), m = 1(black solid line), m = 2 (red dot line), m = 3 (green dashline), and m = 10 (blue dash-dot line); (c) g k ( ω ) of Eq.(30)with g (0) = C , k = 1 (black solid line), k = 2 (red dot line), k = 3 (green dash line), and k = 10 (blue dash-dot line);(d) virtual frequency distributions of heterogenous couplingsystem, ˆ g ( ω ) = | ω | exp {− ω / } / g ( ω ) = | ω | π ( ω +1) (red dot line), ˆ g ( ω ) = | ω | [1 − ( | ω | − Θ( | ω | − ω ∈ [ − , 2] (green dash line), and ˆ g ( ω ) = 2 | ω | (1 − | ω | ), ω ∈ [ − , 1] (blue dash-dot line). real eigenvalue crossing the origin, which implies that K a = 2 / [ πg n (0)]. Additionally, the characteristic func-tion reduces to F ( q ) = Z − g n ( qx ) p − x dx, (21)which is strictly decreasing in q ∈ [0 , ∞ ) for each n [seeFig. 2(a)]. Hence, a continuous phase transition takesplace at q c = 0 with K c = F (0) − = K a [see Fig. 3(a)]corresponding to a trans-critical bifurcation in the pa-rameter space.For the scaling analysis, we note that g n ( ω ) near itsmaximum ω = 0 can be expanded as g n ( ω ) = 1 πγ n sin (cid:16) π n (cid:17) ∞ X m =0 ( − m ( ω/γ ) mn , (22)indicating that the derivative F ( k ) (0) = 0 holds only for k = 2 mn . Power series expansion of F ( q ) around zeroleads to( µ, v, A, B ) = (2 n, n, F (2 n ) (0)(2 n )! , F (4 n ) (0)(4 n )! ) , (23)then the corresponding values of parameters in Eq. (10)become ( α, β, X, Y ) = (2 n, n, − AK c , A K c ) . (24) ( )F s s qq ( )F q( )F q ( )F q q FIG. 2: Characteristic functions are generated by the naturalfrequency distributions (a-c) and virtual distributions (d). (a) F ( q ) marked with different colors corresponds to each g n ( ω )of fig. 1(a) respectively; (b) F ( q ) marked with different col-ors corresponds to each g m ( ω ) of fig. 1(b) respectively; (c) F ( q ) marked with different colors corresponds to each g k ( ω )of fig. 1(c) respectively; (d) F ( s ) marked with different colorscorresponds to each ˆ g ( ω ) of fig. 1(d) respectively. The order parameter R near K c obeys the scaling lawwith the critical exponents ( η, ξ ) = ( n , n +12 n ) and theasymptotic coefficients ( P, Q ) are evaluated directly fromEq. (17).To verify its validity, we take n = 1 as an example. Inthis case, g ( ω ) is a Lorentizan distribution with K a = K c = 2 γ , and the parameter array is given by χ = ( 12 , , q − πg (2)1 (0) K c , − q − πg (2)1 (0) K c ) . (25)Using the self-consistent equation, the exact solutionof the order parameter with K > K c is R = p − γ/K . Taylor expansion of R near K c yields δR = ( γ ) / ( δK ) / − ( γ ) / ( δK ) / , which coincideswith χ perfectly [see Figs. 4(a) and 5(a)]. The resultsobtained above can be applied to arbitrary g ( ω ) whichhas integer-order of nonzero derivative at ω = 0 (like g n ( ω ) ∼ exp {− ω n /γ n } ). Therefore, we conclude thatfor the Kuramoto model with continuous phase transi-tion, the asymptotic behavior of the order parameter nearphase transition point is only determined by the deriv-able property (concave-convex) of the natural frequencydistribution at its maximum.We now extend the analysis above to a family of uni-modal polynomial distributions defined in the finite sup-port, g m ( ω ) = m + 12 m γ m − | ω | m γ m +1 , ω ∈ [ − γ, γ ] , (26)with γ > m > F ( q ) should be discussed in two distinct cases. If q ≤ γ , we have F ( q ) = F (0) + Eq m , (27)where F (0) = π ( m +1)4 mγ and E = − ( m +1) mγ m +1 R x m √ − x dx .If q > γ , we get F ( q ) = m + 1 mqγ m +1 Z γ ( γ m − ω m ) p − ω /q dω. (28)Clearly, F ( q ) is continuous and monotonically decreasingin q ∈ [0 , ∞ ) for each m [see Fig. 2(b)]. Hence, the form ofphase transition remains the same as above, and the sys-tem undergoes a supercritical bifurcation with exchangeof the stability between the incoherence and partial syn-chronization at K a = K c = F (0) − . The explicit formulaof F ( q ) makes scaling analysis easier in comparison withthe series form Eq. (8), solving the self-consistent equa-tion yields R = (cid:18) − K − K c EK c K m +1 (cid:19) /m . (29)Expanding Eq. (29) near K c , the critical parameters forthe asymptotic behavior of R are ( η, ξ ) = (1 /m, /m ), P = ( − E ) − m K − m +2 m c , and Q = − m +1 m ( − E ) − m K − m +2 m c [see Figs. 4(b) and 5(b)]. B. A hybrid phase transition The discussion above indicates that the phase tran-sition corresponding to the onset of synchronization issecond-order as long as there is a single maximum inthe frequency distribution, and that the index of g ( ω ) ( n or m ) controls the critical exponent of the scaling law.Increasing the index of g ( ω ), the continuous phase tran-sition becomes steeper and steeper. A question here ishow the order parameter in the vicinity of the criticalcoupling behaves as if the index tends to infinity. Be-cause lim n →∞ g n ( ω ) = lim m →∞ g m ( ω ) = γ Θ( γ − | ω | ) (Θ is aheaviside function), we have that g ∞ ( ω ) is uniform. Forthe purpose of scaling analysis, we further extend thelimit case to a family of distributions having a plateau atthe maximum, which are defined as g k ( ω ) = g (0) − C ( | ω | − γ ) k Θ( | ω | − γ ) (30)with γ, k > 0. According to the constants g (0) and C ,the distribution is defined in a finite or infinite supportwith the restriction that g k ( ω ) is non-negative and nor-malizable [see Fig. 1(c)]. As shown in refs. [17, 18], g k ( ω )represents a family of unimodal functions with a plateausection in the middle. In particular, k = 0 and C = g (0)degenerate to uniform distribution.The characteristic function F ( q ) corresponding to g k ( ω ) should be discussed in three distinct cases, if0 ≤ q ≤ γ , we have F ( q ) = πg (0)2 ; (31) ( )b R ( )a R FIG. 3: Phase diagram of the order parameter with differentnatural frequency distributions. (a) Continuous phase tran-sition corresponding to g n ( ω ) in Eq. (20) with γ = 1, n = 1(red circle), n = 2 (green triangle), n = 3 (blue diamond), and n = 10 (pink star). (b) Hybrid phase transition correspond-ing to g k ( ω ) in Eq. (30) and the parameters are the same asfig. 1(c) with k = 0 (red circle), k = 1 (green triangle), k = 3(blue diamond), and k = 10 (pink star). The symbols repre-sent numerical simulations with N = 10 and the solid linescorrespond to the solutions of the self-consistent equation. and if γ < q ≤ ω b , F ( q ) = 2 q "Z γ g (0) s − ω q dω + Z qγ G k ( ω, q ) dω ;(32)and if q > ω b , F ( q ) = 2 q "Z γ g (0) s − ω q dω + Z ω b γ G k ( ω, q ) dω . (33)Here, we denote with G k ( ω, q ) = [ g (0) − C ( ω − γ ) k ] p − ω /q and ω b is the boundary frequency suchthat g k ( ω ) = 0 for ω ≥ ω b . Remarkably, F ( q ) is a con-stant for q < γ and monotonously decreases in q ∈ [ γ, ∞ )[see Fig. 2(c)]. Consequently, the plateau structure of F ( q ) informs the form of the phase transition. By in-creasing the coupling strength, the system undergoes anabrupt transition at q c = γ with K a = K c = πg (0) ,where the order parameter R jumps from 0 to a value R c = πg (0)2 γ [see Fig. 3(b)]. In contrast to a conven-tional first order phase transition [19], the reversibly dis-continuous transition is termed the hybrid phase tran-sition [20], where the hysteresis region is replaced by avertical line at K c characterizing an infinite number ofhidden metastable states formed by the oscillators with | ω i | < γ [21].The property of F ( q ) in the neighbourhood of the flatregion is crucial for the scaling behavior of the orderparameter. Changing variables with δq = q − γ and y = ω − γ leads to F ( δq ) = πg (0)2 + Z δq y k f ( δq, y ) dy, (34)where we have defined the function f ( δq, y ) = − Cδq + γ s − (cid:18) y + γδq + γ (cid:19) . (35) -8 -6 -4 -2-4-3-2-10 -5 -4 -3 -2 -1-4-3-2-10-6 -4 -2-4-3-2-10 -6 -4 -2-3-2-1 n=1 log ( ) l og ( R ) ( )c ( )b( )a m=1 l og ( R ) log ( ) k=0 l og ( R ) log ( ) ( )d in-coupling l og ( R ) log ( ) FIG. 4: Scaling behavior of the order parameter with leadingorder. (a) g n ( ω ) of Eq. (20) with n = 1 and γ = 1, (b) g m ( ω )of Eq. (26) with m = 1 and γ = 1, (c) g k ( ω ) of Eq. (30) with k = 0 and γ = 0 . 5, (d) in-coupling case with ˆ g ( ω ) of Eq. (45).The circles represent numerical simulations with N = 10 andthe solid lines correspond to analytical prediction. -3 -2 -1-5-4-3-2-1 -3 -2 -1-4-20-4 -3 -2 -1 0-4-3-2-10 -5 -4 -3 -2 -1-6-4-2 in-couplingk=0 m=1n=1 l og ( P K - R ) log ( ) ( )a ( )b l og ( P K - R ) log ( )c l og ( P K - R ) log ( ) ( )d l og ( P K - R ) log ( ) FIG. 5: Scaling behavior of the order parameter correspond-ing to fig. 4 with sub-leading order and the correspondingparameters are the same as fig. 4. After some calculations, the dominant term of deviationof F ( δq ) with each order reads F ( n ) ( δq ) ∼ Z δq y k ∂ n f∂ ( δq ) n dy ∼ Z δq y k [1 − ( y + γδq + γ ) ] − ( n − ) dy, (36)with n = 1 , , . . . . In the limit δq → F (1) (0) = F (2) (0) = . . . = F ( k +1) (0) = 0, whereas F ( k +2) ( δq ) ∼ ( δq ) − / . To avoid the divergence of the expansion coef-ficients, Taylor expansion of F ( δq ) at δq = 0 should takethe following form, F ( δq ) = 1 K c + A ( δq ) k + + B ( δq ) k + , (37)and the coefficients are given by ( A, B ) = ( d k +3 Fdx k +3 | x = √ δq , d k +5 Fdx k +5 | x = √ δq ). Using Eq.(8), we have( α, β, X, Y ) = ( k + 32 , k + 52 , − AK c , − BK c ) . (38)Hence, the critical scaling exponents for the asymptoticorder parameter near K c are( η, ξ ) = ( 22 k + 3 , min { , k + 3 } ) . (39)Depending on the value of k , the asymptotic coefficients( P, Q ) are determined from the standard procedure [seeeqs.(18, 19)] .A typical example for the illustration is to consider k = 0 corresponding to the uniform distribution of g ( ω ), in which g (0) = C = 1 / (2 γ ). The criticalpoint for the hybrid phase transition is ( K c , R c ) =(4 γ/π, π/ µ, v ) = ( α, β ) = (3 / , / 2) and( A, B ) = ( − √ γ − / , √ γ − / ). The parameter array χ = (2 / , , ( π γ ) / , − π γ ) [see Figs. 4(c) and 5(c)],which recovers the result obtained in [22]. C. A tiered phase transition As we have shown above, the natural frequency dis-tribution of Kuramoto model controls the structure ofcharacteristic function, thereby determining the type ofphase transition as well as the scaling behavior of orderparameter near its onset. To better understanding thisproperty, we generalize Kuramoto model by consideringthe nonuniform coupling with˙ θ i = ω i + 1 N N X j =1 K ij sin( θ j − θ i ) , (40)where K ij is the element of the coupling matrix. In con-trast to the uniform coupling in Eq. (1), K ij accountsfor a class of heterogeneity coupling of the oscillator sys-tem, such as the network topology [23], non-local cou-pling [24], and excitation and inhibition coupling [25].As demonstrated in Ref. [26], the oscillator dynamicsfor certain type of heterogeneity can be reduced to Ku-ramoto model with uniform coupling through a remark-able transformation in the parameter space. The inter-play between phase transition and the critical behaviorof the order parameter can be understood in terms ofa virtual frequency distribution. To this end, we adopt the coupling scheme introduced in refs. [27–29] establish-ing the frequency-coupling correlation between oscillatorsnamely K ij = K | ω j |h| ω |i (out-coupling) and K ij = K | ω i | (in-coupling).Introducing the weighted order parameter W e i Φ = 1 N N X j =1 | ω j |h| ω |i e iθ j , (41)the self-consistent equation of the out-coupling holdsfor replacing ( q, g ( ω )) with ( s, ˆ g ( ω )) that are defined as s = KW and ˆ g ( ω ) = | ω | g ( ω ) h| ω |i . Here ˆ g ( ω ) is called thevirtual frequency distribution, since it plays a equiva-lent role of the natural frequency distribution in the Ku-ramoto model. This distribution combines the hetero-geneity of natural frequencies and coupling, thus pro-viding an intuitive understanding of phase transition forsynchronized dynamics. As depicted in Fig. 1(d), ˆ g ( ω )is typically bimodal for a symmetrical unimodal g ( ω )with ˆ g (0) = ˆ g ( ±∞ ) = 0. Consequently, F ( s ) exhibitsa parabolic-like shape with F (0) = F ( ∞ ) = 0 [seeFig. 2(d)], and its maximum at s c with F ′ ( s c ) = 0.In this setup, the bifurcation mechanism for phasetransition in oscillator system becomes clear. On the onehand, ˆ g (0) = 0 implies that the incoherent state losesits stability via Hopf bifurcation at the critical coupling K a = π ˆ g (Ω c ) [30, 31], where Ω c is the imaginary partof the eigenvalue of the linear stability analysis about R = W = 0. On the other hand, the characteristicfunction F ( s ) indicates that a saddle node bifurcationtakes place at K c = F ( s c ) − , at which a number of os-cillators lock their phases forming a macroscopic ordercharacterized by a non-zero W c = s c F ( s c ). As the re-sult, K a > K c corresponds to the explosive synchroniza-tion observed previously [27], whereas K a < K c refers tothe tiered phase transition where the system goes fromthe incoherence to the partial synchrony mediated by anoscillatory state (standing wave). The oscillatory statedisplaces a time-dependent coherent behavior emergingin the proximity of the critical point, where the tran-sition from the incoherence to synchrony converts fromexplosive to a continuous phase transition. Such an in-termediate state is born in the Hopf bifurcation way anddisappears in a saddle-node infinite periodic or homo-clinic bifurcation manner [32, 33].Within this framework, the power series of F ( s )at s c gives ( µ, v, A, B ) = (2 , , F (2) ( s c ) / , F (3) ( s c ) / α, β, X, Y ) =(2 , , − F (2) ( s c ) / [2 F ( s c )] , − F (3) ( s c ) / [6 F ( s c )]). Thus,the asymptotic behavior of the weighted order param-eter W near the tiered phase transition point follows χ = , , √ F ( s c ) p − F (2) ( s c ) , F ( s c ) F (3) ( s c )3[ F (2) ( s c )] − s c F ( s c ) ! . (42)Likewise, the order parameter R near K c obeys R = R c + P R ( δK ) + Q R ( δK ) , (43)with R c = Z − g ( s c x ) s c p − x dx. (44)The coefficients are P R = √ F ( s c ) I/ p − F (2) ( s c )and Q R = F (3) ( s c ) F ( s c ) I/ [ √ F (2) ( s c )] with I = R − [ g ( s c x ) + g (1) ( s c x ) xs c ] √ − x dx .As an analytical illustration, considering the in-coupling case, the virtual frequency distribution degen-erates to a universal formˆ g ( ω ) = 12 [ δ ( ω − 1) + δ ( ω + 1)] (45)(bimodal distribution with vanishing width), which leadsto a simple characteristic function F ( q ) = q p − q − .The critical point corresponding to the saddle node bifur-cation is q c = √ K c , R c ) = (2 , / √ R near K c is χ = ( 12 , , √ , − √ 216 ) . (46)The analytical expression for the order parameter with K ≥ R = √ r q − K [see Figs. 4(d) and5(d)], which coincides with Eq. (46) perfectly by impos-ing a power series expansion at K = K c . V. CONCLUSION In summary, we have developed an analytical descrip-tion for the critical behavior of the order parameter nearthe onset synchronization transition in the Kuramotomodel. Based on the structures of characteristic func-tions in the self-consistent equation, various phase tran-sitions were identified corresponding to different bifurca-tions for the emergence of collective dynamics in phasespace. More importantly, this analysis provides a univer-sal framework for uncovering the scaling properties of or-der parameter near the onset of synchronization. In par-ticular, the leading and sub-leading critical exponents aswell as asymptotic coefficients are all determined withinthis general framework. Our study provides a new wayfor exploring the synchronization transition in coupledoscillator systems, which could deepen the understand-ing of the mechanism of a phase transition in coupleddynamical networks. ACKNOWLEDGMENTS This work is supported by the National Natural Sci-ence Foundation of China (Grants No. 11905068) and theScientific Research Funds of Huaqiao University (GrantNo. ZQN-810). [1] B. Ermentrout, An adaptive model for synchrony in thefirefly Pteroptyx malaccae, J. Math. Biol. , 571 (1991).[2] P. Richard, B.M. Bakker, B. Teusink, K.V. Dam, andH.V. Westerhoff, Acetaldehyde Mediates the Synchro-nization of Sustained Glycolytic Oscillations in Popula-tions of Yeast Cells, Eur. J. Biochem. , 238 (1996).[3] D. Taylor, E. Ott, and J. G. Restrepo, Spontaneous syn-chronization of coupled oscillator systems with frequencyadaptation, Phys. Rev. E 81, 046214 (2010).[4] G. Filatrella, A.H. Nielsen, and N.F. Pedersen, Analysisof a power grid using a Kuramoto-like model, Eur. Phys.J. B , 485 (2008).[5] S. Benz and C. Burroughs, Coherent emission from two-dimensional Josephson junction arrays, Appl. Phys. Lett. , 2162 (1991).[6] G. Buzs´aki and A. Draguhn, Neuronal Oscillations inCortical Networks, Science , 1926 (2004).[7] I.Z. Kiss, Y. Zhai, and J.L. Hudson, Emerging Coherencein a Population of Chemical Oscillators, Science ,1676 (2002).[8] Y. Kuramoto, in International Symposium on Mathemat-ical Problems in Theoretical Physics , edited by H. Araki,Lecture Notes in Physics No. 30 (Springer, New York,1975), p. 420.[9] J.A. Acebr´on, L.L. Bonilla, C.J. P´erez Vicente, F. Ritort,and R. Spigler, The Kuramoto model: A simple paradigmfor synchronization phenomena. Rev. Mod. Phys. , 137 (2005).[10] E. Ott and T. M. Antonsen, Low dimensional behavior oflarge systems of globally coupled oscillators, Chaos ,037113 (2008).[11] E. Ott and T. M. Antonsen, Long time evolution of phaseoscillator systems, Chaos , 023117 (2009).[12] D. M. Abrams and S. H. Strogatz, Chimera States forCoupled Oscillators, Phys. Rev. Lett. , 174102 (2004).[13] D. M. Abrams, R. Mirollo, S. H. Strogatz, and D. A.Wiley, Solvable Model for Chimera States of CoupledOscillators, Phys. Rev. Lett. , 084103 (2008).[14] C. Xu, X. Wang, and P. S. Skardal, Bifurcation analysisand structural stability of simplicial oscillator popula-tions, Phy. Rev. Research , 023281 (2020).[15] S.H. Strogatz and R.E. Mirollo, Stability of incoherencein a population of coupled oscillators, J. Stat. Phys. ,613 (1991).[16] P. S. Skardal, Low-dimensional dynamics of the Ku-ramoto model with rational frequency distributions,Phys. Rev. E, , 022207 (2018).[17] L. Basnarkov, and V. Urumov, Phase transitions in theKuramoto model, Phys. Rev. E, , 057201 (2007).[18] L. Basnarkov, and V. Urumov, Kuramoto model withasymmetric distribution of natural frequencies, Phys.Rev. E, , 011113 (2008).[19] J. G˜omez-Garde˜nes, S. G´omez, A. Arenas, and Y.Moreno, Explosive synchronization transitions in scale- free networks, Phys. Rev. Lett. , 128701 (2011).[20] J. Park, and B. Kahng, Metastable state en route totraveling-wave synchronization state, Phys. Rev. E, ,020203 (2018).[21] R.E. Mirollo, and S.H. Strogatz, The spectrum of thelocked state for the Kuramoto model of coupled oscilla-tors, Physica D , 249-266 (2005).[22] D. Paz´o, Thermodynamic limit of the first-order phasetransition in the Kuramoto model, Phys. Rev. E, ,046211 (2005).[23] S. Yoon, M. Sorbaro Sindaci, A. V. Goltsev, and J. F.F. Mendes, Critical behavior of the relaxation rate, thesusceptibility, and a pair correlation function in the Ku-ramoto model on scale-free networks, Phys. Rev. E, ,032814 (2015).[24] Y. Kuramoto, and B. Dorjsuren , Coexistence of coher-ence and incoherence in nonlocally coupled phase oscil-lators. arXiv preprint cond-mat/0210694 (2002).[25] H. Hong and S. H. Strogatz, Kuramoto Model of CoupledOscillators with Positive and Negative Coupling Param-eters: An Example of Conformist and Contrarian Oscil-lators, Phys. Rev. Lett. , 054102 (2011).[26] J. Gao, and K. Efstathiou, Reduction of oscillator dy-namics on complex networks to dynamics on completegraphs through virtual frequencies, Phys. Rev. E, ,022302 (2020). [27] X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive syn-chronization in a general complex network, Phys. Rev.E, 88, 010802 (2013).[28] H. Bi, X. Hu, S. Boccaletti, X. Wang, Y. Zou, Z. Liu,and S. Guan, Coexistence of quantized, time dependent,clustersing lobally coupled oscillators, Phys. Rev. Lett. , 204101 (2016).[29] C. Xu, S. Boccaletti, S. Guan, and Z. Zheng, Origin ofBellerophon states in globally coupled phase oscillators,Phys. Rev. E, , 050202(R) (2018).[30] C. Xu, J. Gao, H. Xiang, W. Jia, S. Guan, and Z. Zheng,Dynamics of phase oscillators with generalized frequency-weighted coupling, Phys. Rev. E, , 062204 (2016).[31] C. Xu, S. Boccaletti, Z. Zheng, and S. Guan, Univer-sal phase transitions to synchronization in Kuramoto-likemodels with heterogeneous coupling, New J. Phys. ,113018 (2019).[32] E. A. Martens, E. Barreto, S. H. Strogatz, E. Ott, P.So, and T. M. Antonsen, Exact results for the Kuramotomodel with a bimodal frequency distribution, Phys. Rev.E, , 026204 (2009).[33] D, Paz´o, and E. Montbri´o, Existence of hysteresis in theKuramoto model with bimodal frequency distributions,Phys. Rev. E,80