Classification of bifurcation diagrams in coupled phase-oscillator models with asymmetric natural frequency distributions
aa r X i v : . [ n li n . C D ] N ov Classification of bifurcation diagrams in coupled phase-oscillator models withasymmetric natural frequency distributions
Ryosuke Yoneda and Yoshiyuki Y. Yamaguchi
Graduate School of Informatics, Kyoto University, 606-8501 Kyoto, Japan
Synchronization among rhythmic elements is modeled by coupled phase-oscillators each of whichhas the so-called natural frequency. A symmetric natural frequency distribution induces a continuousor discontinuous synchronization transition from the nonsynchronized state, for instance. It hasbeen numerically reported that asymmetry in the natural frequency distribution brings new typesof bifurcation diagram having, in the order parameter, oscillation or a discontinuous jump whichemerges from a partially synchronized state. We propose a theoretical classification method of fivetypes of bifurcation diagrams including the new ones, paying attention to the generality of thetheory. The oscillation and the jump from partially synchronized states are discussed respectivelyby the linear analysis around the nonsynchronized state and by extending the amplitude equation upto the third leading term. The theoretical classification is examined by comparing with numericallyobtained one.
I. INTRODUCTION
Synchronization among rhythmic elements is observedin various fields of nature: metronomes [1], flashing fire-flies [2, 3], frog choruses [4], and Josephson junction ar-rays [5, 6]. The synchronization has been theoreticallystudied through coupled phase-oscillator models [7, 8].The Kuramoto model [8], a paradigmatic model, de-scribes the bifurcation from the nonsynchronized stateto partially synchronized states and the type of bifurca-tion depends on the distribution of the so-called natu-ral frequencies each of which is assigned for each phase-oscillator. The bifurcation is continuous if the naturalfrequency distribution is symmetric and unimodal [8–10],but the continuity breaks if the distribution is flat at thepeak [11, 12] or bimodal [13]. A bimodal distributionalso yields temporally oscillating states [13].A significant part of previous studies assumes symme-try of natural frequency distributions, and bifurcation isfound as the synchronization transition referring to thetransition from the nonsynchronized state to partiallysynchronized states. Interestingly, symmetry breakingbrings new types of bifurcations that emerge not fromthe nonsynchronized state but from partially synchro-nized states and are followed by a discontinuous jump oroscillation of the order parameter [14]. The new types ofbifurcations have been observed by performing numericalsimulations, and this paper aims to propose a theoreticalexplanation of the new types of bifurcations. For an-alyzing bifurcations, we have three theoretical methodsbased on the large population limit and the equation ofcontinuity.The first method is the self-consistent equation for theorder parameter of system [8, 9, 11, 15, 16]. A stationarysolution to the equation of continuity is formally writ-ten by using an unknown value of the order parameterand the value is determined self-consistently. The self-consistent strategy is powerful to obtain stationary valuesof the order parameter, but the oscillating states are notcaptured. Moreover, writing down stationary states is not straightforward when a coupling function between apair of oscillators has several sinusoidal functions [17, 18].The second method is the Ott-Antonsen ansatz [19,20]. This ansatz reduces the equation of continuity toa simpler partial differential equations. This strategyis applicable for the Kuramoto model and its variances[21, 22], and is extremely useful when we use a singleLorentzian as the natural frequency distribution, therebyresulting a real two-dimensional reduced system. The di-mension of the reduced system becomes higher as thedistribution becomes complicated [12–14], and searchingfor stationary and oscillating states becomes difficult ac-cordingly.The third method is the amplitude equation, which de-scribes dynamics projected onto the unstable manifold ofthe nonsynchronized state [23]. In contrast to the previ-ous two methods, the amplitude equation is widely usefulin the Kuramoto model [23], in a coupled phase-oscillatormodel with a generic coupling function [24], in the time-delayed Kuramoto model [25], in an extended Kuramotomodel with inertia [26], and in Hamiltonian systems [27–29]. Other advantages of the amplitude equation are thatone-dimensional dynamics is sufficient to judge the con-tinuity of bifurcation and that any natural frequency dis-tributions are equally tractable.The first and second methods profit from a single si-nusoidal coupling function, but such a simple couplingfunction is not always the case [24, 30–32]. We, there-fore, adopt the third method, the amplitude equation, toexplain the new types of bifurcations. One basic blockof the amplitude equation is the linear analysis aroundthe nonsynchronized state, and the linear analysis gives,for instance, eigenvalues. We will demonstrate that thelinear analysis is also useful to explain a mechanism forthe emergence of an oscillating state.The method proposed in this paper is in principle ap-plicable to general systems beyond the Kuramoto modeland its variances, but we investigate the Kuramoto modelto illustrate the usefulness of the method. One reason isthat the new types of bifurcations are reported in the Ku-ramoto model, and this choice makes it straightforwardto compare the theory with the previous work [14]. An-other reason is that, as we discussed above, the reductionby the Ott-Antonsen ansatz permits us to obtain accu-rate classification, which helps to examine the theory.This paper is organized as follows. In Sec. II, we in-troduce the Kuramoto model and its large populationlimit written by the equation of continuity. The consid-ering family of the natural frequency distribution, whichis characterized by two parameters, is also exhibited. InSec. III, we numerically divide the parameter plane intofive domains corresponding to types of bifurcation dia-grams and prepare the reference parameter plane to ex-amine the theory. This numerical search is performed byusing the Ott-Antonsen reduction. We stress that thisreduction is used only for obtaining the reference param-eter plane and is not used in our theory. In Sec. IV,linear and nonlinear analyses of the equation of continu-ity are shortly reviewed. With the aid of these analyses,we propose ideas to identify the domains on the parame-ter plane and report the theoretical consequence in Sec.V. Finally, we summarize this paper in Sec. VI.
II. MODEL
The Kuramoto model is expressed by the N -dimensional ordinary differential equation,d θ j d t = ω j − KN N X k =1 sin ( θ j − θ k ) , ( j = 1 , · · · , N ) . (1)The real constant K > θ j ∈ ( − π, π ] = S and ω j ∈ R are respectively the phase andthe natural frequency of the j th oscillator. The naturalfrequencies { ω j } obey a probability distribution function g ( ω ). By using the parameters γ , γ > ≥
0, weintroduce a family of g ( ω ) as g ( ω ) = C [( ω − Ω) + γ ][( ω + Ω) + γ ] (2)to systematically consider unimodal and bimodal, andsymmetric and asymmetric distributions [14, 33]. Thefamily of g ( ω ) consists of rational functions, which is use-ful to draw the reference parameter plane by using theOtt-Antonsen ansatz, but is not crucial in our method-ology. The normalization constant C = γ γ [( γ + γ ) + 4Ω ] π ( γ + γ ) (3)is determined from the condition Z ∞−∞ g ( ω )d ω = 1 . (4)The three parameters γ , γ , and Ω are assumed to bepositive. By scaling the variables t, ω j , K, and γ , wemay set γ = 1 without loss of generality, and the fam-ily of g ( ω ) is characterized by a point on the parameter plane ( γ , Ω). The parameter plane is further restrictedinto the region of γ ≤ γ > γ < θ j → − θ j and ω j → − ω j , and thelatter induces the transformation of g ( ω ) → g ( − ω ). Theequation of motion is also invariant under the change oftime scale t → t/τ with ω j → τ ω j and K → τ K . Thesechanges do not modify type of the bifurcation diagramfor τ > g ( − ω ) ∝ /τ [( ω + Ω /τ ) + ( γ /τ ) ][( ω − Ω /τ ) + ( γ /τ ) ] . (5)The point ( γ , Ω), therefore, corresponds to the point(1 /γ , Ω /γ ) by selecting τ = γ under γ = 1. The line γ = 1 gives a family of symmetric distributions.The complex order parameter z is defined by z = re iψ = 1 N N X j =1 e iθ j , ( r, ψ ∈ R ) . (6)The absolute value r of z measures the extent of synchro-nization. If all the oscillators distribute uniformly on S ,then r ≃
0. If r ≃
1, the majority of oscillators gathersaround a point on S .In the limit of large population N → ∞ , by the con-servation of the number of oscillators, the equation ofmotion (1) can be written in the equation of continuity[34], ∂F∂t + ∂∂θ ( v [ F ] F ) = 0 ,v [ F ] = ω − K Z ∞−∞ d ω ′ Z π − π d θ ′ sin( θ − θ ′ ) F ( θ ′ , ω ′ , t ) , (7)where F ( θ, ω, t ) is the probability distribution functionof θ and ω at the time t . In other words, F ( θ, ω, t ) dθ d ω represents the fraction of oscillators having phases be-tween θ and θ + dθ and natural frequencies between ω and ω + d ω at the time t . From the normalization con-dition R F dθ d ω = 1, we have Z π − π F ( θ, ω, t )d θ = g ( ω ) . (8)In this limit the order parameter is expressed by z = re iψ = Z ∞−∞ d ω Z π − π d θ F ( θ, ω, t ) e iθ . (9) III. NUMERICAL CLASSIFICATION ONPARAMETER PLANE
The aim of this section is to classify numerically theparameter plane ( γ , Ω) into five types of bifurcation di-agrams [33] before developing the theory. Direct N -bodysimulations of the model (1) include finite-size fluctua-tions, which make it difficult to judge the types of bifur-cation diagrams. For eliminating the finite-size fluctua-tions, we perform the Ott-Antonsen reduction [19, 20],which reduces the equation of continuity to a real four-dimensional differential equation for the family of g ( ω ),(2). A. Ott-Antonsen reduction
The reduced equation by using the Ott-Antonsenansatz is expressed as [14]d z d t = ( i Ω − γ ) z − K z ∗ z − z ) , d z d t = − ( i Ω + γ ) z − K z ∗ z − z ) . (10)The two variables z and z are complex and the complexorder parameter z is written as z = k z + k z (11)with the complex constants k and k defined by k = γ [2Ω − i ( γ + γ )]( γ + γ )[2Ω + i ( γ − γ )] ,k = γ [2Ω + i ( γ + γ )]( γ + γ )[2Ω + i ( γ − γ )] . (12) z ∗ represents the complex conjugate of z . In the latercomputations we set γ = 1 without loss of generality aswe mentioned in the previous section II, but we kept γ free to show the dependence explicitly. See Appendix Afor details of the reduction. B. Bifurcation diagrams
Numerical integration of the reduced system (10) isperformed by using the fourth-order Runge-Kutta algo-rithm with the time step ∆ t = 0 .
01. For a given set of( γ , Ω), we start from K = 0 and increase the value up to K = 10 with the step ∆ K = 0 . K , the timeaverage and standard deviation of the order parameterare taken in the time interval t ∈ [4500 , t = 5000 is used asthe initial state at the successive value of K . If the finalstate is the origin, then we shift the initial state fromthe origin to z = z = 0 .
01 to escape from the trivialstationary state. After arriving K = 10, we decrease thevalue of K from K = 10 to K = 0 to check existence ofa hysteresis, which reveals a discontinuous bifurcation.This decreasing process is called the backward process.The parameter plane ( γ , Ω) is classified into five do-mains as shown in Fig. 1 [33], which domains corre-spond to the five types of bifurcation diagrams reported . . . . . . γ . . . . . . . . . . Ω A E BD C a bcde FIG. 1. Classification of the parameter plane ( γ , Ω) into thefive domains: A (filled orange triangles), B (open blue cir-cles), C (open magenta squares), D (filled magenta squares),and E (filled blue circles). This classification is obtained byperforming numerical simulations of the reduced system (10).The gray line represents the borderline between the unimodaland bimodal natural frequency distributions g ( ω ). The bifur-cation diagrams at the five points marked by the crosses arereported in Fig. 2. in Fig. 2. The symmetry line γ = 1, in particular,is classified into the three intervals included in the do-mains A, B, and C. The separating points Ω = 1 and √ γ < IV. THEORETICAL ANALYSES OF EQUATIONOF CONTINUITY
We shortly review linear and nonlinear analyses of theequation of continuity (7) around the nonsynchronizedstate f , which is explicitly written as f ( ω ) = g ( ω )2 π . (13)It is straightforward to check stationarity of f as ∂∂θ ( v [ f ] f ) = 0 (14)from the fact v [ f ] = ω . We expand the equation ofcontinuity by substituting F ( θ, ω, t ) = f ( ω ) + f ( θ, ω, t ) . (15)The perturbation f is governed by the equation ∂f∂t = L f + N [ f ] , (16) . . . r (a1) − . . . g ( ω ) (a2) . . . (b1) − . . . (b2) . . (c1) − . . (c2) . . (d1) − . . . (d2) K . . . (e1) − ω . . . (e2) FIG. 2. Five bifurcation diagrams at the five points markedon the parameter plane in Fig. 1 and the corresponding nat-ural frequency distributions g ( ω ). The values of parame-ters ( γ , Ω) are (a) (1 . , . . , . . , . . , . . , . g ( ω ) against ω . where the linear part is L f = − ω ∂f∂θ − Kf ∂∂θ Im (cid:0) z ( t ) e − iθ (cid:1) (17)and the nonlinear part is N [ f ] = − K ∂∂θ (cid:2) f Im (cid:0) z ( t ) e − iθ (cid:1)(cid:3) (18) with the order parameter z ( t ) = Z ∞−∞ d ω Z π − π d θe iθ f ( θ, ω, t ) . (19)Note that the equation (16) is an exact transform fromthe equation of continuity (7). The linear part will beused to explain the oscillating state emerging from a par-tially synchronized state as well as a basic block of theamplitude equation. A. Linear analysis
The perturbation f is a periodic function of θ and isexpanded into the Fourier series as f ( θ, ω, t ) = X k ∈ Z e f k ( ω, t ) e ikθ . (20)The linear analysis of (17) can be performed indepen-dently in each Fourier mode k because the nonsynchro-nized state f ( ω ) does not depend on θ . The Fouriermodes k = ± k = ±
1. The eigenvalues for the modes k = ± ± ( λ ) = 1 − K Z ∞−∞ g ( ω ) λ ± iω d ω. (21)See Appendix B for derivations. If the real part of aroot is positive, the eigenvalue induces instability. Wecall such an eigenvalue as an unstable eigenvalue, whichis the target of the amplitude equation introduced in thenext subsection IV B.We define the synchronization transition point K c atwhich the nonsynchronized state f changes stability.This definition suggests that the eigenvalue having thelargest real part must be on the imaginary axis at K c . Apure imaginary eigenvalue, however, induces a singular-ity in the integrands of the spectral functions. To avoidthe singularity, we perform the analytic continuation ofΛ ± ( λ ) and denote the continued functions as D ± ( λ ).See Appendix C for the continuation.We give three remarks. First, the continuation doesnot modify the spectral functions in the region Re( λ ) > D ± ( λ ) whose real part is positive is also aroot of Λ ± ( λ ). Second, a root of D ± ( λ ), however, maynot be an eigenvalue if it is on the region Re( λ ) ≤
0, andthe root is called a fake eigenvalue accordingly. Third,the relations Λ − ( λ ∗ ) = Λ ∗ ( λ ) (22)and D − ( λ ∗ ) = D ∗ ( λ ) (23)hold, where Λ ∗ ( λ ) is, for instance, the complex conju-gate of Λ ( λ ). These relations imply that λ ∗ is a (fake)eigenvalue if λ is.For the family of natural frequency distributions (2),the equation D ( λ ) = 0 leads a quadratic equation of λ and the two roots are denoted by λ and λ ordered asRe( λ ) ≥ Re( λ ). In this paper, the two (fake) eigen-values λ and λ are called the first and the second(fake) eigenvalues respectively. See Appendix D for thequadratic equation of λ , the determination of the criticalpoint K c , and the number of unstable eigenvalues. B. Amplitude equation
We assume that the linear operator L has only onepair of unstable eigenvalues, λ and λ ∗ coming from therelation (23), and that the corresponding eigenfunctionsare Ψ( θ, ω ) and Ψ ∗ ( θ, ω ) respectively. To derive the am-plitude equation, we expand f into f ( θ, ω, t ) = A ( t )Ψ( θ, ω ) + A ∗ ( t )Ψ ∗ ( θ, ω ) + H ( θ, ω, A, A ∗ ) . (24)The amplitude A relates to the order parameter z as z = 2 πA ∗ + O ( | A | ) , (25)and the asymptotic value of z is approximately obtainedby considering temporal evolution of the amplitude A ( t ).The function H represents the unstable manifold of thenonsynchronized state f . In other words, H representsthe height of the unstable manifold from the eigenspaceSpan { Ψ , Ψ ∗ } . We assume that the unstable manifold istangent to the eigenspace Span { Ψ , Ψ ∗ } at A = A ∗ = 0and H = O ( | A | ) accordingly.For deriving the amplitude equation, we introduce theadjoint linear operator L † of L defined by (cid:0) L † f , f (cid:1) = ( f , L f ) , for ∀ f , f , (26)where the inner product ( · , · ) is defined by( f , f ) = Z ∞−∞ d ω Z π − π d θf ∗ ( θ, ω ) f ( θ, ω ) . (27)Let e Ψ and e Ψ ∗ be the eigenfunctions of L † correspondingto the eigenvalues λ ∗ and λ respectively. We can choosethe eigenfunctions satisfying the relations (cid:16) e Ψ , Ψ (cid:17) = 1 , (cid:16) e Ψ , Ψ ∗ (cid:17) = 0 , (cid:16) e Ψ ∗ , Ψ (cid:17) = 0 , (cid:16) e Ψ ∗ , Ψ ∗ (cid:17) = 1 , (28)without loss of generality. We also have (cid:16) e Ψ , H (cid:17) = (cid:16) e Ψ ∗ , H (cid:17) = 0 (29)because H does not belong to the eigenspaceSpan { Ψ , Ψ ∗ } . See Appendix E for the explicit expres-sions of the adjoint operator and the eigenfunctions. Substituting the expansion (24) into the equation ofcontinuity (7) and using the relations (28) and (29), wehave the amplitude equation asd A d t = λ A + (cid:16) e Ψ , N [ f ] (cid:17) (30)and the equation for the unstable manifold H asd H d t = L H + N [ f ] − h(cid:16) e Ψ , N [ f ] (cid:17) Ψ + (cid:16) e Ψ ∗ , N [ f ] (cid:17) Ψ ∗ i . (31)These equations can be solved perturbatively for suffi-ciently small | A | , and the right-hand-side of (30) is ex-panded asd A d t = λ A + c A | A | + c A | A | + c A | A | + · · · , (32)where the eigenvalue λ and the coefficients c , c , · · · depend on the coupling constant K . Note that the right-hand-side of (32) has only odd order terms. See Ap-pendix F for derivations and the explicit forms of coeffi-cients.The complex amplitude equation (32) can be reducedto a real equation written asd σ d t = 2 σG ( σ ) , (33)where σ = | A | ≥ G ( σ ) = Re( λ ) + Re( c ) σ + Re( c ) σ + · · · . (34)We will search for stationary solutions with which theright-hand-side of the amplitude equation (33) is zero.The nonsynchronized state, for instance, corresponds to σ = 0, which always satisfies the stationary condition.The amplitude equation is useful to determine the con-tinuity of the synchronization transition from the non-synchronized state f to a partially synchronized state.Around the critical point K c , the order parameter mustbe small if the transition is continuous and we can usethe truncated equationRe( λ ) + Re( c ) σ = 0 . (35)This equation has a nontrivial solution for Re( c ) < λ ) >
0, while no realsolution for Re( c ) > c ( K ) with respect to K and taking the limit K → K c + 0, we say that the syn-chronization transition is continuous if Re( c ( K c )) < c ( K c )) > γ , Ω) = (1 ,
1) between the domains A and Bon the line γ = 1 is obtained by searching for the pointwhich satisfies Re( c ( K c )) = 0. V. CRITERIA FOR DETERMINING DOMAINS
Looking back the five bifurcation diagrams exhibitedin Fig. 2, we have two elements to characterize the dia-grams, which are oscillation and a jump of the order pa-rameter. We first discuss mechanisms of the oscillationand of the jump in Secs. V A and V B respectively withthe aid of the linear and nonlinear analyses reviewed inthe previous section IV. After that, we propose a proce-dure to divide the parameter plane into the five domainsin Sec. V C.
A. Oscillation of order parameter
The asymptotically periodic oscillation of the order pa-rameter has been also observed in the Hamiltonian mean-field model, and the oscillation comes from existence oftwo small clusters running with different velocities [35].Inspired by this work, we propose the following idea; Thefirst unstable eigenvalue induces a partially synchronizedstate and the oscillation is excited by appearance of asecond unstable eigenvalue.The order parameter is computed as z ( t ) = 2 π Z ∞−∞ e f − ( ω, t )d ω, (36)the time evolution of the order parameter z ( t ) is hencedescribed by the eigenvalues arisen from the Fourier − D − ( λ ), are denoted by λ ∗ and λ ∗ , where Re( λ ∗ ) ≥ Re( λ ∗ ).Let us increase the value of the coupling constant K from the nonsynchronized region K < K c . Beyond thecritical point K c , the first fake eigenvalue λ ∗ becomesunstable and the resonant oscillators, whose natural fre-quencies are close to Im( λ ∗ ), form a cluster. For in-stance, if g ( ω ) is symmetric and unimodal, the resonantfrequency is zero and oscillators around ω = 0 form a syn-chronized cluster. This small cluster corresponds to thecontinuous synchronization transition in the bifurcationdiagram (see Fig. 2 (d)). Further increasing K , for suit-able pairs of ( γ , Ω), the second fake eigenvalue λ ∗ alsobecomes unstable and forms the second cluster. The tworesonant frequencies differ in general, Im( λ ∗ ) = Im( λ ∗ ),and the order parameter oscillates; The order parame-ter is large when the two clusters are close on S , andis small when the clusters are in the antiphase positionseach other. The symmetry of g ( ω ), in particular, resultsin inducing simultaneous destabilization of the two fakeeigenvalues λ ∗ and λ ∗ at K c and yielding the bifurcationdiagram Fig. 2 (c).Summarizing, existence of two unstable eigenvalues,which mean two pairs of the unstable eigenvalues bycounting the roots of D ( λ ) in addition to the roots of D − ( λ ), suggests the oscillation of the order parameter.To support this mechanism, we investigate K dependence of the fake eigenvalues λ ∗ and λ ∗ at points chosen fromthe domains C, D, and E of the parameter plane ( γ , Ω).A symmetric case is examined in Fig. 3 for ( γ , Ω) =(1 . , .
0) which belongs to the domain C. As we expected,two eigenvalues become unstable at the same K c withdifferent imaginary parts. In Fig. 1 the separating pointΩ = √ γ = 1 is obtained by checking existence of the twounstable eigenvalues. K − − − K c Re( λ ∗ ( K ))Im( λ ∗ ( K ))Re( λ ∗ ( K ))Im( λ ∗ ( K )) FIG. 3. K dependence of the fake eigenvalues λ ∗ and λ ∗ for ( γ , Ω) = (1 . , . λ ∗ ( K ))(red solid line) and Re( λ ∗ ( K )) (red dashed line) collapse for K ≤ √ − ≃ . λ ∗ and λ ∗ become unstable simultaneously at K c = 4. For the point ( γ , Ω) = (0 . , .
0) belonging to the do-main D, the simultaneity of destabilization breaks asshown in Fig. 4 owing to asymmetry of g ( ω ). A cor-respondence between the two unstable eigenvalues andthe two clusters is exhibited in Fig. 5, where N -bodysimulations of the Kuramoto model (1) were performedby using the fourth-order Runge-Kutta algorithm withthe time step ∆ t = 0 .
01. See Appendix G for the pro-cedure to determine the synchronized intervals of ω re-ported as red bars in Fig. 5. The second unstable eigen-value and the second cluster emerge at almost same val-ues of K , and this observation is consistent with thediscussion above. The scenario holds even when thefirst cluster is further developed as shown in Fig. 6for the point ( γ , Ω) = (0 . , . γ , Ω) = (0 . , .
8) close to the domain D but be-longing to the domain E as reported in Fig. 7.We note that the discussion above does not always holdbecause the second unstable eigenvalue does not alwaysyield the second cluster. At the point ( γ , Ω) = (0 . , . K = 5 . K − − − − K c Re( λ ∗ ( K ))Im( λ ∗ ( K ))Re( λ ∗ ( K ))Im( λ ∗ ( K )) FIG. 4. K dependence of the fake eigenvalues λ ∗ and λ ∗ for( γ , Ω) = (0 . , .
0) belonging to the domain D. The first fakeeigenvalue λ ∗ becomes unstable at the critical point K c =2 . λ ∗ becomes unstableat a larger K = 4 . . . . . . . . . K − − − − ω Im( λ ∗ (2 . λ ∗ (4 . . . g ( ω ) FIG. 5. The natural frequency distribution and the emerg-ing clusters for ( γ , Ω) = (0 . , . N -body system (1)with N = 10 , are found in the ranges of red bars. The sec-ond fake eigenvalue λ ∗ becomes unstable at K = 4 . K = 6 .
04, which is plotted with gray dash-dotted line, andthis leads to the jump in the order parameter, as seen inFig. 2(d1). stable eigenvalue emerges but its imaginary part is notsufficiently far from the grown first cluster, and the sec-ond virtual cluster is absorbed by the first cluster withoutemerging as shown in Fig. 8. One more condition must . . . . . . . K − − − − ω Im( λ ∗ (1 . λ ∗ (5 . . . g ( ω ) FIG. 6. Same with Fig. 5 but for ( γ , Ω) = (0 . , .
0) belong-ing to the domain D. The second fake eigenvalue λ ∗ becomesunstable at K = 5 . K = 7 . K − − K c Re( λ ∗ ( K ))Im( λ ∗ ( K ))Re( λ ∗ ( K ))Im( λ ∗ ( K )) FIG. 7. Same with Fig. 4 but for ( γ , Ω) = (0 . , .
8) belong-ing to the domain E. The real part of the fake eigenvalue λ approaches to zero, but does not become positive. be therefore added to characterize the oscillation; | Im( λ ∗ ) − Im( λ ∗ ) | > s (37)for some s [36]. See Appendix H for the further investi-gation on this condition.We give another remark on multi-cluster states; Themaximum number of clusters is not necessarily twofor the bimodal natural frequency distribution. Multi-cluster states called the Bellerophon states are reportedin [37, 38] for symmetric natural frequency distributions.The symmetry is, however, not essential to the multi- . . . . . . . . K − − − ω Im( λ ∗ (2 . λ ∗ (5 . . . g ( ω ) FIG. 8. Same with Fig. 5 but for ( γ , Ω) = (0 . , .
6) be-longing to the domain E. The second fake eigenvalue becomesunstable at K = 5 . cluster states as an example is shown in Fig. 9 for anasymmetric distribution with ( γ , Ω) = (0 . , . K − . − . − . − . . . . . . ω Im( λ ∗ (2 . λ ∗ (4 . . . g ( ω ) FIG. 9. Same with Fig. 5 but for ( γ , Ω) = (0 . , .
5) belong-ing to the domain E. The second fake eigenvalue λ ∗ becomesunstable at K = 4 . g ( ω ). This multi-cluster state is referred to as theBellerophon state [37, 38]. B. Jump of order parameter
A jump of the order parameter emerges from r = 0(domain B) or from r > r = 0 is identified by Re( c ( K c )) > γ < r >
0, a typical bifur-cation diagram of the type E is schematically shown inFig. 10. We focus on the saddle-node bifurcation point K Q , at which G ( σ ) in the amplitude equation (33) hasdegenerated nontrivial roots. The degeneracy can be ap-proximately captured by G ( σ ), the truncation of G ( σ )up to the second order of σ , as the vanishing discrimi-nant.Let us discuss validity of the expansion (32) aroundthe saddle-node bifurcation point K Q . The point K Q isapproximated by K Q ′ obtained from G ( σ ). The valueof σ at K Q is also approximated as σ ′∗ = − Re( c ( K Q ′ ))Re( c ( K Q ′ )) (38)and gives the equality | Re( c ( K Q ′ ))( σ ′∗ ) || Re( c ( K Q ′ )) σ ′∗ | = 12 . (39)This equality suggests that the expansion (32) up to thethird term is not excellent but does not break aroundthe bifurcation point K Q . We will discuss validity of(32) later from the view point of smallness of the orderparameter.We remark on the upper stable branch appearing at an-other saddle-node bifurcation point K P . Coexistence ofthe three branches between K P and K Q is a distinguish-ing feature of the type E, but we focus on K Q because ofthe two advantages: (i) Capturing the coexistence in theamplitude equation (33) requires the third order term of G ( σ ), the coefficient c , which is hard to compute. (ii)The upper branch has a larger value of σ than σ ′∗ , andvalidity of the expansion (32) becomes worse accordingly. C. Theoretical division of the parameter plane
In the basis of the discussions done in Secs. V A andV B, we divide the parameter plane ( γ , Ω) into the fivedomains by the following flow. For preparation, we com-pute the critical point K c for a given pair of parameters( γ , Ω) as shown in Appendix D.We first focus on the critical point K = K c . Thediscontinuous jump appears at K c only in the typeB among the considered five bifurcation diagrams (seeFig. 2). The point ( γ , Ω) must belong to the domainB if Re( c ( K c )) > r = 0. When this condition does not hold, we check ifthere are two unstable eigenvalues at K c . If this condi-tion is satisfied, oscillation of the order parameter startsfrom K c and the point ( γ , Ω) is considered to be in thedomain C. σ G ( σ ) K K Q K c K P r FIG. 10. Schematic picture for the bifurcation diagram oftype E. The red solid lines denote the stable branches, andthe blue dashed lines denote the unstable branches. Insetsare schematic graphs of G ( σ ). When both checks at K c are negative, we increase K from K c and examine the following two propositions:(Oscillation) ∃ K O ( > K c ) s . t . Re( λ ( K O )) > , (40)and (Jump) ∃ K J ( > K c ) s . t . ∆( K J ) = 0 , (41)where the discriminant ∆( K ) of G ( σ ) is defined by∆( K ) = Re( c ( K )) − λ ( K ))Re( c ( K )) . (42)Note that the degenerated root of G ( σ ), − c / (2 c ), ispositive under the assumptions c ( K J ) < c ( K J ) >
0. If both the propositions (40) and (41) are false, wedecide that the point ( γ , Ω) belongs to the domain A.If the oscillation proposition (40) is true but the jumpproposition (41) is false, the point must be in the domainD. If the oscillation proposition (40) is false but the jumpproposition (41) is true, the point must be in the domainE. If both the propositions are true, there is a competitionbetween K O and K J ; K O < K J suggests the domain Dand K O > K J suggests the domain E. In Table I, wesummarize techniques which are used to determine thedomainsThe flow proposed above provides three theoreticallines, I, II, and III which divide the parameter plane( γ , Ω) into the five domains, as reported in Fig. 11. theupper side of the line I is the theoretical domain D, theinside of the line II is the theoretical domain B, and thelines I, II, and III enclose the theoretical domain E. Thetheoretical lines reproduce qualitatively the numericallyobtained domains.For the line I, our result is fairly consistent with nu-merical ones, but the domain D is overestimated. Ourmethod says that the emergence of two clusters is asso-ciated with the emergence of two unstable eigenvalues.
TABLE I. Techniques we use to determine the domains. λ corresponds to the linear theory, and c and c correspond tothe nonlinear theory. λ c c Phenomenology A X X X X B X C X X D X X X X E X X X
FIG. 11. The theoretical lines I(red broken line), II(greensolid line), and III(blue dash-dotted line) are overlapped onFig. 1. Line I denotes the borderline between the domains Dand A ∪ B ∪ E. Line II denotes the borderline between thedomains B and A ∪ C ∪ E. Line III denotes the borderlinebetween the domains E and A. The gray line represents theborderline between the unimodal and bimodal natural fre-quency distributions g ( ω ). The second eigenvalue, however, does not always giverise to the second cluster as we discussed in Sec. V A.We have to reduce the domain D by introducing an ad-ditional condition that | Im( λ ∗ ) − Im( λ ∗ ) | is sufficientlylarge.The line II is perfect because the order parameteris small around the critical point of the synchroniza-tion transition and the perturbatively obtained ampli-tude equation is valid for such a small amplitude. Thiscriterion had been used for symmetric g ( ω ), but we haveconfirmed that it is also powerful for asymmetric ones.0 .
86 0 .
88 0 .
90 0 .
92 0 .
94 0 .
96 0 .
98 1 . γ . . . . . . r Q . . . . . . K Q − K P , K Q − K c KK Q K c K P r KK Q K c K P r FIG. 12. The absolute value r of the order parameter at K Q (blue solid line) and the value K Q − K P (red dashed line)for Ω = 1 .
7. As γ decreases, r Q , the smaller value of r at K = K Q , increases monotonically, and the value K Q − K P decreases monotonically. The red dash-dotted line representsthe value K Q − K c , which is drawn only for γ giving K P < K c ;The upper stable branch emerges before the lower unstablebranch emerges. The line II also reveals existence of discontinuous syn-chronization transition for unimodal but asymmetric dis-tributions. The existence has been pointed out numeri-cally in a previous study [33], but our theoretical analysisensures it.The line III is not perfect. We used a perturbativemethod to derive the amplitude equation (32), assumingthat the order parameter is sufficiently small at the jump-ing point K Q . This assumption becomes worse as thepoint ( γ , Ω) approaches to the boundary of the domainE where the width of the hysteresis vanishes as reportedin Fig. 12. It is thus expectable that the boundary ofthe domain E is not perfectly described by the theoreti-cal line III. Nevertheless, the theoretical line III capturesthe domain E qualitatively and is useful to have a strongcandidate of the domain E.
VI. CONCLUSION AND DISCUSSIONS
We have proposed a theoretical method to classifybifurcation diagrams in the Kuramoto model by usingthe amplitude equation with the aid of the linear anal-ysis around the nonsynchronized reference state. Theamplitude equation, obtained perturbatively, is usuallystopped up to the second leading term because it is suffi-cient to judge the continuity of the synchronization tran-sition, bifurcation from the nonsynchronized state. Thepremise, however, is broken by introducing asymmetryin the natural frequency distribution because the asym-metry induces bifurcations from partially synchronizedstates. We have extended the amplitude equation up tothe third leading term and successfully captured a dis-continuous bifurcation after a continuous one. For bifurcation diagrams having oscillating states, wehave focused on unstable eigenvalues of the nonsynchro-nized state. Roughly speaking, one unstable eigenvaluecorresponds to one cluster formation, and hence exis-tence of two unstable eigenvalues suggest appearance oftwo clusters rotating with different speeds. This idea isqualitatively in good agreement with numerical results,but overestimates the domain having the oscillation. Wehave to reduce the domain by adding one more conditionwhich requires a large discrepancy of rotating speeds ofthe two clusters in order to avoid the absorption of thesecond virtual cluster by the first grown cluster. Theadditional condition has to be investigated. We remarkthat a similar condition was discussed in a Hamiltoniansystem [36].We mainly used the amplitude equation to predict thebifurcation diagram. This method is a perturbation tech-nique and is not powerful if the considering point is farfrom the critical point. This restriction prevents us fromgetting precise theoretical lines for Ω sufficiently large,for instance, because the oscillation/jump emerging froma partially synchronized state are not close to the criti-cal point. On the contrary, for symmetric and bimodalnatural frequency distributions, our method successfullyidentified the oscillating states, the domain C, solely bya linear analysis with phenomenology, while they havebeen previously analyzed by a nonlinear analysis, theamplitude equation [23]. Our analysis hence capturesan essence of the oscillating states.It is worth noting that a cluster is created by the res-onance between natural frequencies an unstable eigen-value, but the number of clusters, n c , is not necessarilythe same with the number of unstable eigenvalues, n u . Inour computations, n c = 1 when n u = 1 but n c ≥ n u = 2. Such a multi-cluster state with n c >
2, calledthe Bellerophon state, was reported previously [37, 38]but the mechanism is an open question. Noting that twoor more clusters induce oscillation of the order parameterand following the scenario that a resonance makes a clus-ter, we expect that the multi-cluster state with n c > c ) may be proportional to apower of 1 / Re( λ ). Nevertheless, the amplitude equationcaptures the asymptotic scaling of the order parameteragainst the strength of instability in plasmas, and we1may expect that the amplitude equation still has powerto predict the bifurcation diagram. It is another futureproblem to check the precision of the proposed criterionunder such divergent coefficients. ACKNOWLEDGMENTS
Y.Y.Y. acknowledges the support of JSPS KAKENHIGrant No. 16K05472.
Appendix A: Ott-Antonsen ansatz and reduction
We derive the reduced system (10) from the equationof continuity (7) by using the Ott-Antonsen ansatz. TheOtt-Antonsen ansatz assumes the form of the probabilitydistribution F ( θ, ω, t ) as F = g ( ω )2 π ( ∞ X k =1 (cid:2) a ( ω, t ) k e ikθ + a ∗ ( ω, t ) k e − ikθ (cid:3)) (A1)with | a ( ω, t ) | < ∂a∂t + iωa + K (cid:0) a z − z ∗ (cid:1) = 0 , (A2)where the order parameter is written as z = Z ∞−∞ g ( ω ) a ∗ ( ω, t )d ω. (A3)The integration of the right-hand side can be computedby using the residue theorem because | g ( ω ) a ∗ ( ω, t ) | We introduce the analytic continuation of the spectralfunctions Λ ± ( λ ), (21). We start from a point λ whose2real part is positive, Re( λ ) > 0, and decrease the realpart. The condition Re( λ ) > b f ( s ) = Z ∞ f ( t ) e − st d t, (C1)which analyzes temporal evolution and is defined in thepositive real half-plane Re( s ) > s corresponds to λ .When λ passes the imaginary axis, the integrands ofthe spectral functions Λ ± ( λ ), (21), meet the singularitiesat ω = ± iλ . To avoid this singularity at ω = iλ ( ω = − iλ ), we modify continuously the integral contour fromthe real axis by adding a small lower (upper) half-circlearound the singular point. This modification brings aresidue part to the integral. Continuing this modificationfor Re( λ ) < 0, we have analytically continued functionsof Λ ± ( λ ) as D ± ( λ ) = 1 − K I ± ( λ ) , (C2)where I ± ( λ ) = Z ∞−∞ g ( ω ) λ ± iω d ω, Re( λ ) > Z ∞−∞ g ( ω ) λ ± iω d ω + πg ( ± iλ ) , Re( λ ) = 0 Z ∞−∞ g ( ω ) λ ± iω d ω + 2 πg ( ± iλ ) , Re( λ ) < D ± ( λ ) and of Λ ± ( λ ) are identical from their defi-nitions if Re( λ ) > 0, and hence such roots of D ± ( λ )are eigenvalues. They do not coincide, however, forRe( λ ) ≤ D ± ( λ ) with Re( λ ) ≤ ( λ ) and D ( λ ) is demon-strated by considering a Lorentzian natural frequencydistribution, g ( ω ) = γπ ω + γ . (C4)Straightforward computations giveΛ ( λ ) = − K λ + γ ) , Re( λ ) > − K λ − γ ) , Re( λ ) < D ( λ ) = 1 − K λ + γ ) , λ ∈ C . (C6) λ = K/ − γ is a root of D ( λ ) but is not of Λ ( λ ) if K ≤ K c and Re( λ ) ≤ K c = 2 γ is the synchronization transition point. Similarly, we can repro-duce the critical point K c = 2 / [ πg (0)] for a symmetricunimodal g ( ω ) by computing the roots of D ± ( λ ). Thecritical point K c for the considering family (2) is givenin Appendix D. Appendix D: Fake eigenvalues and synchronizationtransition point K c The fake eigenvalues are roots of the continued spec-trum functions, (C2). They move continuously on thecomplex plane, and the synchronization transition point K c is computed as the point where one of the fake eigen-value passes the imaginary axis from left to right with thesmallest coupling constant K > 0. The fake eigenvaluesare analyzed in Appendix D 1 for the natural frequencydistribution (2). We will show that the first fake eigen-value, which has the largest real part, passes the imagi-nary axis at most once time. After that, the equation todetermine K c will be derived in Appendix D 2 with ananalysis of the movement of the second fake eigenvalue. 1. Fake eigenvalues Substituting (2) to (C3), we obtain the explicit formof D ( λ ) as D ( λ ) = 1 − K γ + γ + λ + ( γ + ) + i Ω γ − ( λ + γ + i Ω)( λ + γ − i Ω) (D1)where γ + = γ + γ > , γ − = γ − γ ≤ . (D2)The equation D ( λ ) = 0 induces the quadratic equation λ − b ( K ) λ − a R ( K ) − ia I ( K ) = 0 (D3)where b ( K ) = K − γ + ,a R ( K ) = K γ + − ( γ γ + Ω ) ,a I ( K ) = Ω γ − (cid:18) K γ + (cid:19) ≤ . (D4)To write down the two solutions to (D3), we introducethe complex variable x = b + 4 a R + i a I = ρe iθ , ρ, θ ∈ R . (D5)The argument θ is in the interval [ π, π ] from a I ≤ 0. Thetwo fake eigenvalues λ and λ , which satisfy Re( λ ) ≥ Re( λ ), are written as2 λ = b − √ ρ cos θ − i √ ρ sin θ , λ = b + √ ρ cos θ i √ ρ sin θ , (D6)3as cos( θ/ ≤ 0. Note that the signs of the imaginaryparts are Im( λ ) ≤ λ ) ≥ θ/ ≥ λ but λ ∗ .The synchronization transition point K c is determinedby the equation Re( λ ( K c )) = 0. Let us show that thesolution is at most one. Using the relationcos θ − r θ ρ and θ , ρ = p ( b + 4 a R ) + (4 a I ) , cos θ = b + 4 a R ρ , (D8)we have2Re( λ ( K )) = b + 1 √ qp ( b + 4 a R ) + (4 a I ) + b + 4 a R . (D9)The functions b ( K ) , a ( K ), and b ( K ) + 4 a R ( K ) are in-creasing functions of K for K > b + 4 a R = (cid:18) K γ + (cid:19) − γ γ + Ω ) . (D10)These increasing functions imply that the real partRe( λ ( K )) is also an increasing function of K for K > 2. Synchronization transition point K c We show that there exists the unique solution to theequation Re( λ ( K )) = 0, which determines the synchro-nization transition point K c . At K c , the fake eigenvalue λ must be pure imaginary of the form λ = iλ I ( λ I ∈ R ).Substituting this form with K = K c into the quadraticequation (D3), we have − λ − ib ( K c ) λ I − a R ( K c ) − ia I ( K c ) = 0 . (D11)The real part reads λ + a R ( K c ) = 0 (D12)and the imaginary part reads b ( K c ) λ I + a I ( K c ) = 0 . (D13)Eliminating λ I from (D12) and (D13), we have the cubicequation to determine K c as( γ + ) K − (cid:8) γ + ) + γ γ (cid:2) ( γ + ) + 4Ω (cid:3)(cid:9) K + 4 γ + (cid:2) ( γ + ) + 2 γ γ ( γ + ) + 4Ω (cid:0) γ + γ (cid:1)(cid:3) K c − γ γ ( γ + ) (cid:2) ( γ + ) + 4Ω (cid:3) = 0 . (D14) The number of real solutions to (D14) is one or three, andall the real solutions are positive because the left-handside of (D14) is always negative for K c < 0. The synchro-nization transition point is determined as the smallestreal solution to (D14).To investigate the number of unstable eigenvalues, weintroduce the discriminant ∆ for the cubic equation(D14). In general, the discriminant is defined as∆ = b c − a d − ac − b d + 18 abcd (D15)for the cubic equation ax + bx + cx + d = 0 . (D16)The number of real solutions is one for ∆ < > < λ passesthe imaginary axis at the unique real solution K c and noother passing occurs. In the case ∆ > K (1)c , K (2)c , and K (2)c , where K (1)c < K (2)c < K (3)c . By the definition, the first fakeeigenvalue λ passes the imaginary axis at K c = K (1)c andit does not pass the imaginary axis any more as shown inthe previous subsection D 1. As a result, the other twosolutions K (2)c and K (3)c are realized by the second fakeeigenvalue λ ; It passes the imaginary axis from left toright at K (2)c , and from right to left at K (3)c , as observedin Fig. 4. Appendix E: Adjoint linear operator The idea of deriving the amplitude equation presentedin Sec. IV B is to project the full dynamics onto the re-duced space spanned by some eigenfunctions of the linearoperator L , (17). The projection operator is defined byeigenfunctions of the adjoint operator of L . In this Ap-pendix, a necessary eigenfunction of the adjoint operator L † is presented.The adjoint operator L † acts on f as L † f = ω ∂f∂θ + K (cid:0) r [ f f ] e − iθ + r − [ f f ] e iθ (cid:1) , (E1)where r n [ f ] = Z ∞−∞ d ω Z π − π d θ f ( θ, ω, t ) e inθ . (E2)As done for the linear operator L , we expand L † into theFourier series as L † f = X k ∈ Z L † k e f k ( ω, t ) e ikθ . (E3)The linear operator L † k is defined by L † k e f k = ikω e f k + K δ k, + δ k, − ] Z ∞−∞ e f k ( ω, t ) g ( ω )d ω. (E4)4We focus on the modes k = ± µ be an eigenvalue of L † k and e ψ ( ω ) be the corre-sponding eigenfunction which satisfies L † k e ψ = µ e ψ. (E5)This equation brings( µ − ikω ) e ψ ( ω ) = K Z ∞−∞ e ψ ( ω ) g ( ω )d ω. (E6)Repeating the same discussion done in Appendix B, theeigenvalue µ must satisfyΛ ∗ k ( µ ∗ ) = 0 . (E7)This equation implies that µ ∗ is an eigenvalue of L † k if µ is an eigenvalue of L k .Let λ be an eigenvalue of L and ψ be the correspond-ing eigenfunction (B8). From the above discussion, L † has an eigenvalue λ ∗ and the corresponding eigenfunc-tion is computed as e ψ ( ω ) = 1[Λ ′ ( λ )] ∗ λ ∗ − iω . (E8)As the case of L , λ ∗ is also an eigenvalue of L † and thecorresponding eigenfunction is e Ψ( θ, ω ) = e ψ ( ω )2 π e iθ , (E9)which satisfies the normalization condition (cid:16) e Ψ , Ψ (cid:17) = 1 . (E10) Appendix F: Derivation of the amplitude equation The coefficients of the amplitude equation (32) are ob-tained perturbatively with an expression of the unstablemanifold. We give explicit forms of the coefficients. 1. Derivations of equations for A and H We assume that λ is the unique unstable eigenvalueof L and ψ ( ω ) is the corresponding eigenfunction. Therelation Λ − ( λ ∗ ) = Λ ∗ ( λ ) implies that λ ∗ is an eigen-value of L − and ψ ∗ ( ω ) is the corresponding eigenfunc-tion. The linear operator L , therefore, has two unstableeigenvalues λ and λ ∗ , and the corresponding eigenfunc-tions respectivelyΨ( θ, ω ) = ψ ( ω ) e iθ , Ψ ∗ ( θ, ω ) = ψ ∗ ( ω ) e − iθ . (F1)Using these eigenfunctions, we expand the perturbation f into the form of (24), f = A ( t )Ψ + A ∗ ( t )Ψ ∗ + H ( θ, ω, A, A ∗ ) . (F2) We assume that the unstable manifold H is tangent tothe unstable eigenspace Span(Ψ , Ψ ∗ ) at A = A ∗ = 0.Substituting this expansion into the equation of continu-ity, we haved A d t Ψ + d A ∗ d t Ψ ∗ + d H d t = λA Ψ + λ ∗ A ∗ Ψ ∗ + L H + N [ f ] . (F3)In (F3), taking the inner product with e Ψ, we obtainthe equation for A asd A d t = λA + (cid:16) e Ψ , N [ f ] (cid:17) . (F4)Extracting this equality and its complex conjugate from(F3), we have the equation for H asd H d t = L H + N [ f ] − h(cid:16) e Ψ , N [ f ] (cid:17) Ψ + (cid:16) e Ψ ∗ , N [ f ] (cid:17) Ψ ∗ i . (F5)The left-hand side of (F5) is read asd H d t = ∂H∂A d A d t + ∂H∂A ∗ d A ∗ d t . (F6)The right-hand side of (F4) starts from the linear term λA , while one of (F5) starts from the quadratic terms of A and A ∗ by the tangency assumption of the unstablemanifold, H = O ( | A | ). We, therefore, solve the twoequations (F4) and (F5) perturbatively by assuming that | A | is small. 2. Fourier series expansion of H Before going to the Taylor series expansion with re-spect to | A | , we expand H into the Fourier series as H = X n ∈ Z H n ( ω, A, A ∗ ) e inθ . (F7)The rotational symmetry of the system yields [23] H n ( ω, A, A ∗ ) = n = 0 Aσh ( σ, ω ) n = 1 A n h n ( σ, ω ) n ≥ , (F8)where σ = | A | .Now we expand h n into the Taylor series. For suffi-ciently small σ we expand h n ( σ, ω ) = h n, ( ω ) + σh n, ( ω ) + · · · . (F9)Algebraic computations gives the equation of A asd A d t = λA + c Aσ + c Aσ + c Aσ + · · · (F10)5 c c c h , ✟✟✙ ❄ h , h , h , ✲ ✛❄ ❄ ❄✟✟✙❍❍❥ ✟✟✙❍❍❥ h , h , h , h , ✲ ✛ ✛ FIG. 13. Hierarchy of h n,m . Each line represents h n,m whichnewly appear to compute the coefficient written in the mostleft column. Arrows indicate dependency, although some ar-rows are omitted for a graphical reason. The functions withunderlines are directly used in (F11). where c = − πK h e ψ, h , i ,c = − πK (cid:20) h e ψ, h , i (cid:18)Z h , d ω (cid:19) ∗ + h e ψ, h , i (cid:21) ,c = − πK (cid:20)(cid:18)Z h , d ω (cid:19) ∗ h e ψ, h , i + h e ψ, h , i + (cid:18)Z h , d ω (cid:19) ∗ h e ψ, h , i (cid:21) , (F11)and h f , f i = Z ∞−∞ f ∗ f d ω. (F12)For simplicity of notation, we introduce a symbol hh f ii = Z ∞−∞ f d ω. (F13)We compute h e ψ, h , i , hh h , ii , h e ψ, h , i , hh h , ii and h e ψ, h , i by expanding (F5) into the Fourier series. FourFourier modes of n = 1 , , c , c , and c as shown in Fig. 13. 3. Functions h n,m The leading order of the Fourier mode n , h n, ( ω ), isexpressed as h n, ( ω ) = π n − K n g ( ω )( λ + iω ) n , (2 ≤ n ≤ . (F14)The other 5 functions appearing in Fig. 13 solve the fol-lowing equations:(2 λ + λ ∗ + iω ) h , = K g ( ω ) hh h , ii− πKh , − c ψ, (F15)(3 λ + λ ∗ + 2 iω ) h , = 2 πK ( hh h , ii ψ + h , − h , ) − c h , , (F16) (3 λ + 2 λ ∗ + iω ) h , = K g ( ω ) hh h , ii− πK ( hh h , ii ∗ h , + h , ) − (2 c + c ∗ ) h , − c ψ (F17)2(2 λ + λ ∗ + iω ) h , = − (3 c + c ∗ ) h , − c h , + 2 πK ( hh h , ii ψ + h , + hh h , ii h , − h , − hh h , ii ∗ h , ) , (F18)and(4 λ + λ ∗ + 3 iω ) h , = − c h , + 3 πK ( hh h , ii h , + h , − h , ) . (F19) 4. Inner products and integrals in the coefficients The expression (F14) and the eigenfunction e ψ , (E8),give the inner product h e ψ, h , i as h e ψ, h , i = − πK D ′′ ( λ ) D ( λ ) . (F20)To express the other quantities, we introduce two integraloperators of I m [ f ]( λ ) = Z ∞−∞ fmλ + ( m − λ ∗ + iω d ω (F21)and J m [ f ]( λ ) = Z ∞−∞ f ( λ + iω )[ mλ + λ ∗ + ( m − iω ] d ω. (F22)Remarking that hh h , ii and hh h , ii are solved self-consistently and using the relation1 − K I m [ g ]( λ ) = D ( mλ + ( m − λ ∗ ) , (F23)we have the necessary quantities as hh h , ii = − πK I [ h , ] + c I [ ψ ] D (2 λ + λ ∗ ) , (F24) h e ψ, h , i = − c D ′ ( λ ) J [ h , ]+ 2 πKD ′ ( λ ) ( hh h , iiJ [ ψ ] + J [ h , ] − J [ h , ]) , (F25) hh h , ii = − (2 c + c ∗ ) I [ h , ] + c I [ ψ ] D (3 λ + 2 λ ∗ ) − πKD (3 λ + 2 λ ∗ ) ( hh h , ii ∗ I [ h , ] + I [ h , ]) (F26)6and h e ψ, h , i = − c + c ∗ D ′ ( λ ) J [ h , ] − c D ′ ( λ ) J [ h , ]+ πKD ′ ( λ ) ( hh h , iiJ [ ψ ] + J [ h , ] + hh h , iiJ [ h , ] −J [ h , ] − hh h , ii ∗ J [ h , ]) , (F27)where all the integral operators, I m and J m , are evalu-ated at λ . Appendix G: Procedure to determine thesynchronized intervals of ω This Appendix summarizes the procedure to determinethe intervals of ω , which we call the synchronized inter-vals, in which oscillators synchronize. The synchronizedintervals are identified through dividing the ω -axis intobins of the width w . The bin α defined by [ αw, ( α + 1) w )induces the set P α consisting of the indices of the oscil-lators included in this bin as P α = { i | i ∈ { , · · · , N } , ω i ∈ [ αw, ( α + 1) w ) } . (G1)We set w = 0 . θ, ω ) plane, when the bin is in the synchronized inter-vals. In the Kuramoto model, the gathering point with ω = ( α + 1 / w is unique for the bin α , and the standarddeviation of the phases { θ i } i ∈ P α , denoted by s α , must besufficiently small. In contrast, the standard deviation s α must be large when the bin is not in the synchronizedintervals. We note that the standard deviation s α is de-fined by s α = s | P α | X i ∈ P α ( θ i − µ α ) , (G2)where | P α | is the number of elements of the set P α and µ α is the mean in the bin α as µ α = 1 | P α | X i ∈ P α θ i . (G3)The standard deviation and the mean are computed aftertaking mod 2 π for each θ i .We have to be careful for picking up all the bins be-longing to the synchronized intervals because the stan-dard deviation s α becomes large if the gathering point isclose to 0 or 2 π . To overcome this difficulty, we preparethe shifted phases e θ i = θ i + π (mod 2 π ) (G4)and compute the standard deviation e s α for { e θ i } i ∈ P α . Byusing s α and e s α , we conclude that the bin α is in the synchronized intervals ifmin { s α , e s α } ≤ s ∗ (G5)holds with s ∗ = 0 . . . . . . . γ . . . . . . . . . . Ω FIG. 14. Comparison among a restricted domain D by thecriterion (H1) (blue solid line), the original theoretical domainD (red broken line, which is the same with the line I reportedin Fig. 11), and numerically obtained domain D (magentasquares), together with domain C (open magenta squares). Appendix H: Criterion of the imaginary difference(37) The domain D is characterized by oscillation of the or-der parameter. As discussed in Sec. V A, the oscillationmust have one more condition to satisfy, the imaginarydifference (37), to ensure emergence of the second clusterwithout being absorbed by the first cluster. We will showa trial to introduce this criterion to shave the theoreti-cally obtained domain D (see Fig. 11).For the Kuramoto model with a unimodal and sym-metric natural frequency distribution, only one unstableeigenvalue for K > K c exists, resulting in only one clus-ter and the continuous transition. The half width of thecluster along the ω axis is calculated as Kr [9]. Inspiredby this estimation, we introduce the difference criterionas | Im( λ ∗ ) − Im( λ ∗ ) | > K r ( K ) (H1)where K is the emerging point of the second unstableeigenvalue. We estimate r ( K ), the value of order pa-rameter at K = K , from the solution of the amplitudeequation up to the third order, that is r ( K ) = 2 π s − Re( λ ( K ))Re( c ( K )) . (H2)In Fig. 14, a domain D restricted by the criterion (H1)is compared with the original theoretical domain D. Thecriterion shaves a restricted area of the domain D, whichis close to the original theoretical line I. The criterionstill needs to be investigated.7 [1] J. Pantaleone, Synchronization of metronomes, Am. J.Phys. , 992 (2002).[2] H. M. Smith, Synchronous flashing of fireflies, Science , 151 (1935).[3] J. Buck and E. Buck, Mechanism of rhythmic syn-chronous flashing of fireflies, Science , 1319 (1968).[4] I. Aihara, T. Mizumoto, T. Otsuka, H. Awano, K. Nagira,H. G. Okuno and K. Aihara, Spatio-Temporal Dynam-ics in Collective Frog Choruses Examined by Mathemat-ical Modeling and Field Observations, Sci. Rep. , 3891(2014).[5] K. Wiesenfeld, P. Colet and S. H. Strogatz, Synchroniza-tion Transitions in a Disordered Josephson Series Array, Phys. Rev. Lett. , 404 (1996).[6] K. Wiesenfeld, P. Colet and S. H. Strogatz, Frequencylocking in Josephson arrays: Connection with the Ku-ramoto model, Phys. Rev. E , 1563 (1998).[7] A. T. Winfree, Biological Rhythms and the Behavior ofPopulations of Coupled Oscillators, J. Theoret. Biol. ,15 (1967).[8] Y. Kuramoto, Self-entertainment of a population of cou-pled non-linear oscillators, International Symposium onMathematical Problems in Theoretical Physics, LectureNotes in Physics , 420 (1975).[9] S. H. Strogatz, From Kuramoto to Crawford: exploringthe onset of synchronization in populations of coupledoscillators, Physica , 1 (2000).[10] H. Chiba, A proof of the Kuramoto conjecture for a bi-furcation structure of the infinite-dimensional Kuramotomodel, Ergod. Theo. Dyn. Syst. , 762 (2015).[11] L. Basnarkov and V. Urumov, Phase transitions in theKuramoto model, Phys. Rev. E , 057201 (2007).[12] B. Pietras, N. Deschle, and A. Daffertshofer, First-orderphase transitions in the Kuramoto model with com-pact bimodal frequency distributions, Phys. Rev. E ,062219 (2018).[13] 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).[14] Y. Terada, K. Ito, T. Aoyagi, and Y. Y. Yamaguchi, Non-standard transition in the Kuramoto model: a role ofasymmetry in natural frequency distributions, J. Stat.Mech. (2017) 0134403.[15] D. Paz´o, Thermodynamic limit of the first-order phasetransition in the Kuramoto model, Phys. Rev. E ,046211 (2005).[16] J. Park and B. Kahng, Synchronization transitionsthrough metastable state on structured networks,arXiv:1901.02123.[17] M. Komarov and A. Pikovsky, Multiplicity of SingularSynchronous States in the Kuramoto Model of CoupledOscillators, Phys. Rev. Lett. , 204101 (2013).[18] M. Komarov and A. Pikovsky, The Kuramoto model ofcoupled oscillators with a bi-harmonic coupling function, Physica D , 18 (2014).[19] E. Ott and T. M. Antonsen, Low dimensional behavior oflarge systems of globally coupled oscillators, Chaos ,037113 (2008).[20] E. Ott and T. M. Antonsen, Long time evolution of phaseoscillator systems, Chaos , 023117 (2009). [21] Z. Lu, K. Klein-Carde˜na, S. Lee, T. M. Antonsen, M.Girvan, and E. Ott, Resynchronization of circadian os-cillators and the east-west asymmetry of jet-lag, Chaos , 094811 (2016).[22] A. Akao, S. Shirasaka, Y. Jimbo, B. Ermentrout,and K. Kotani, Theta-gamma cross-frequency cou-pling enables covariance between distant brain regions,arXiv:1903.12155.[23] J. D. Crawford, Amplitude Expansions for Instabilitiesin Populations of Globally-Coupled Oscillators, J. Stat.Phys. , 1047 (1994).[24] J. D. Crawford, Scaling and Singularities in the Entrain-ment of Globally Coupled Oscillators, Phys. Rev. Lett. , 4341 (1995).[25] D. M´etivier and S. Gupta, Bifurcations in the Time-Delayed Kuramoto Model of Coupled Oscillators: ExactResults, J. Stat. Phys. , (2019).[26] J. Barr´e and D. M´etivier, Bifurcations and Singulari-ties for Coupled Oscillators with Inertia and Frustration, Phys. Rev. Lett. , 214102 (2016).[27] J. D. Crawford, Universal Trapping Scaling on the Un-stable Manifold for a Collisionless Electrostatic Mode, Phys. Rev. Lett. , 656 (1994).[28] J. D. Crawford, Amplitude equations for electrostaticwaves: Universal singular behavior in the limit of weakinstability, Phys. Plasmas , 97 (1995).[29] J. Barr´e, D. M´etivier, and Y. Y. Yamaguchi, Trappingscaling for bifurcations in the Vlasov systems, Phys. Rev.E , 042207 (2016).[30] D. Hansel, G. Mato and C. Meunier, Phase Dynamicsfor Weakly Coupled Hodgkin-Huxley Neurons, Europhys.Lett. , 367 (1993).[31] I. Z. Kiss, Y. Zhai and J. L. Hudson, Predicting Mu-tual Entrainment of Oscillators with Experiment-BasedPhase Models, Phys. Rev. Lett. , 248301 (2005).[32] I. Z. Kiss, Y. Zhai and J. L. Hudson, Characteristicsof Cluster Formation in a Population of Globally Cou-pled Electrochemical Oscillators: An Experiment-BasedPhase Model Approach, Prog. Theor. Phys. Suppl. ,99 (2006).[33] Y. Terada, K. Ito, R. Yoneda, T. Aoyagi and Y. Y. Yam-aguchi, A role of asymmetry in linear response of globallycoupled oscillator systems, arXiv:1802.08383.[34] C. Lancellotti, On the Vlasov Limit for Systems of Non-linearly Coupled Oscillators without Noise, Transp. The-ory Stat. , 523 (2004).[35] H. Morita and K. Kaneko, Collective Oscillation in aHamiltonian System, Phys. Rev. Lett. , 050602 (2006).[36] J. Barr´e and Y. Y. Yamaguchi, Small traveling clusters inattractive and repulsive Hamiltonian mean-field models, Phys. Rev. E , 036208 (2009).[37] H. Bi, X. Hu, S. Boccaletti, X. Wang, Y. Zou, Z. Liu, andS. Guan, Coexistence of Quantized, Time Dependent,Clusters in Globally Coupled Oscillators, Phys. Rev. Lett. , 204101 (2016).[38] X. Li, T. Qiu, S. Boccaletti, I. Sendi-Nadal, Z. Liu andS. Guan, Synchronization clusters emerge as the result ofa global coupling among classical phase oscillators, NewJ. Phys. , 053002 (2019).[39] H. Sakaguchi and Y. Kuramoto, A Soluble Active Ro-tator Model Showing Phase Transitions via Mutual En- trainment, Prog. Theor. Phys. , 576 (1986).[40] O. E. Omel’chenko and M. Wlfrum, Nonuniversal Tran-sitions to Synchrony in the Sakaguchi-Kuramoto Model, Phys. Rev. Lett. , 164101 (2012).[41] M. K. S. Yeung and S. H. Strogatz, Time Delay in theKuramoto Model of Coupled Oscillators, Phys. Rev. Lett. , 648 (1999).[42] E. Montbri´o, D. Paz´o and J. Schmidt, Time delay in theKuramoto model with bimodal frequency distribution, Phys. Rev. E , 056201 (2006). [43] S. H. Strogatz, R. E. Mirollo and P. C. Matthews, Cou-pled Nonlinear Oscillators below the SynchronizationThreshold: Relaxation be Generalized Landau Damping, Phys. Rev. Lett. , 2730 (1992).[44] S. Ogawa, Spectral and formal stability criteria of spa-tially inhomogeneous stationary solutions to the Vlasovequation for the Hamiltonian mean-field model, Phys.Rev. E87