Classification of bifurcation diagrams in coupled phase-oscillator models with asymmetric natural frequency distributions

TL;DR

This paper develops a theoretical framework to classify five types of bifurcation diagrams in coupled phase-oscillator models with asymmetric natural frequency distributions, expanding understanding of synchronization phenomena.

Contribution

It introduces a general classification method for bifurcation diagrams, including new types caused by asymmetry, using linear analysis and amplitude equations.

Findings

Five bifurcation diagram types identified

Theoretical classification matches numerical results

New bifurcation behaviors explained by asymmetry

Abstract

Synchronization among rhythmic elements is modeled by coupled phase-oscillators each of which has the so-called natural frequency. A symmetric natural frequency distribution induces a continuous or discontinuous synchronization transition from the nonsynchronized state, for instance. It has been numerically reported that asymmetry in the natural frequency distribution brings new types of bifurcation diagram having, in the order parameter, oscillation or a discontinuous jump which emerges from a partially synchronized state. We propose a theoretical classification method of five types of bifurcation diagrams including the new ones, paying attention to the generality of the theory. The oscillation and the jump from partially synchronized states are discussed respectively by the linear analysis around the nonsynchronized state and by extending the amplitude equation up to the third leading term. The theoretical classification is examined by comparing with numerically obtained one.

Tables

Table 1. Table 1: Techniques we use to determine the domains. λ 1 subscript 𝜆 1 \lambda_{1} corresponds to the linear theory, and c 3 subscript 𝑐 3 c_{3} and c 5 subscript 𝑐 5 c_{5} correspond to the nonlinear theory.

	${\lambda }_1$	$c_3$	$c_5$	Phenomenology
	✓	✓	✓	✓
	✓
	✓			✓
	✓	✓	✓	✓
	✓	✓	✓

Equations

$$\frac{\mathrm{d}{\theta }_j}{\mathrm{d}t}={\omega }_j-\frac{K}{N}\sum _{k=1}^N\sin \left({\theta }_j-{\theta }_k\right),(j=1,\cdots ,N).$$

$$g(\omega )=\frac{C}{[(\omega -\Omega {)}^2+{\gamma }_1^2][(\omega +\Omega {)}^2+{\gamma }_2^2]}$$

$$C=\frac{{\gamma }_1{\gamma }_2[({\gamma }_1+{\gamma }_2{)}^2+4{\Omega }^2]}{\pi ({\gamma }_1+{\gamma }_2)}$$

$${\int }_{-\infty }^{\infty }g(\omega )\mathrm{d}\omega =1.$$

$$g(-\omega )\propto \frac{1/{\tau }^4}{[(\omega +\Omega /\tau {)}^2+({\gamma }_1/\tau {)}^2][(\omega -\Omega /\tau {)}^2+({\gamma }_2/\tau {)}^2]}.$$

$$z=r{e}^{i\psi }=\frac{1}{N}\sum _{j=1}^N{e}^{i{\theta }_j},(r,\psi \in \mathbb{R}).$$

$$\begin{array}{rl}& \frac{\partial F}{\partial t}+\frac{\partial }{\partial \theta }(v[F]F)=0,\\ & v[F]=\omega -K{\int }_{-\infty }^{\infty }\mathrm{d}{\omega }^{\prime }{\int }_{-\pi }^{\pi }\mathrm{d}{\theta }^{\prime }\sin (\theta -{\theta }^{\prime })F({\theta }^{\prime },{\omega }^{\prime },t),\end{array}$$

$${\int }_{-\pi }^{\pi }F(\theta ,\omega ,t)\mathrm{d}\theta =g(\omega ).$$

$$z=r{e}^{i\psi }={\int }_{-\infty }^{\infty }\mathrm{d}\omega {\int }_{-\pi }^{\pi }\mathrm{d}\theta F(\theta ,\omega ,t)e^{i\theta }.$$

$$\begin{array}{rl}& \frac{\mathrm{d}{z}_1}{\mathrm{d}t}=(i\Omega -{\gamma }_1)z_1-\frac{K}{2}(z^{\ast }{z}_1^2-z),\\ & \frac{\mathrm{d}{z}_2}{\mathrm{d}t}=-(i\Omega +{\gamma }_2)z_2-\frac{K}{2}(z^{\ast }{z}_2^2-z).\end{array}$$

$$z=k_1{z}_1+k_2{z}_2$$

$$\begin{array}{rl}& k_1=\frac{{\gamma }_2[2\Omega -i({\gamma }_1+{\gamma }_2)]}{({\gamma }_1+{\gamma }_2)[2\Omega +i({\gamma }_1-{\gamma }_2)]},\\ & k_2=\frac{{\gamma }_1[2\Omega +i({\gamma }_1+{\gamma }_2)]}{({\gamma }_1+{\gamma }_2)[2\Omega +i({\gamma }_1-{\gamma }_2)]}.\end{array}$$

$$f^0(\omega )=\frac{g(\omega )}{2\pi }.$$

$$\frac{\partial }{\partial \theta }(v[f^0]f^0)=0$$

$$F(\theta ,\omega ,t)=f^0(\omega )+f(\theta ,\omega ,t).$$

$$\frac{\partial f}{\partial t}=\mathcal{L}f+\mathcal{N}[f],$$

$$\mathcal{L}f=-\omega \frac{\partial f}{\partial \theta }-K{f}^0\frac{\partial }{\partial \theta }\mathrm{Im}\left(z(t)e^{-i\theta }\right)$$

$$\mathcal{N}[f]=-K\frac{\partial }{\partial \theta }\left[f\mathrm{Im}\left(z(t)e^{-i\theta }\right)\right]$$

$$z(t)={\int }_{-\infty }^{\infty }\mathrm{d}\omega {\int }_{-\pi }^{\pi }\mathrm{d}\theta e^{i\theta }f(\theta ,\omega ,t).$$

$$f(\theta ,\omega ,t)=\sum _{k\in \mathbb{Z}}{\tilde{f}}_k(\omega ,t)e^{ik\theta }.$$

$${\Lambda }_{\pm 1}(\lambda )=1-\frac{K}{2}{\int }_{-\infty }^{\infty }\frac{g(\omega )}{\lambda \pm i\omega }\mathrm{d}\omega .$$

$${\Lambda }_{-1}({\lambda }^{\ast })={\Lambda }_1^{\ast }(\lambda )$$

$$D_{-1}({\lambda }^{\ast })=D_1^{\ast }(\lambda )$$

$$f(\theta ,\omega ,t)=A(t)\Psi (\theta ,\omega )+A^{\ast }(t){\Psi }^{\ast }(\theta ,\omega )+H(\theta ,\omega ,A,A^{\ast }).$$

$$z=2\pi A^{\ast }+O(|A{|}^3),$$

$$\left({\mathcal{L}}^{†}{f}_1,f_2\right)=\left(f_1,\mathcal{L}{f}_2\right),\mathrm{for}\forall f_1,f_2,$$

$$\left(f_1,f_2\right)={\int }_{-\infty }^{\infty }\mathrm{d}\omega {\int }_{-\pi }^{\pi }\mathrm{d}\theta f_1^{\ast }(\theta ,\omega )f_2(\theta ,\omega ).$$

$$\begin{array}{rl}& \left(\tilde{\Psi },\Psi \right)=1,\left(\tilde{\Psi },{\Psi }^{\ast }\right)=0,\\ & \left({\tilde{\Psi }}^{\ast },\Psi \right)=0,\left({\tilde{\Psi }}^{\ast },{\Psi }^{\ast }\right)=1,\end{array}$$

$$\left(\tilde{\Psi },H\right)=\left({\tilde{\Psi }}^{\ast },H\right)=0$$

$$\frac{\mathrm{d}A}{\mathrm{d}t}={\lambda }_{1}A+\left(\tilde{\Psi },\mathcal{N}[f]\right)$$

Full text

I Introduction

Synchronization among rhythmic elements is observed in various fields of nature: metronomes pantaleone2002 , flashing fireflies smith1935 ; buck1968 , frog choruses aihara2014 , and Josephson junction arrays wiesenfeld1996 ; wiesenfeld1998 . The synchronization has been theoretically studied through coupled phase-oscillator models winfree1967 ; kuramoto1975 . The Kuramoto model kuramoto1975 , a paradigmatic model, describes the bifurcation from the nonsynchronized state to partially synchronized states and the type of bifurcation depends on the distribution of the so-called natural frequencies each of which is assigned for each phase-oscillator. The bifurcation is continuous if the natural frequency distribution is symmetric and unimodal kuramoto1975 ; strogatz2000 ; chiba2015 , but the continuity breaks if the distribution is flat at the peak basnarkov-urumov-07 ; bastian-18 or bimodal martens2009 . A bimodal distribution also yields temporally oscillating states martens2009 .

A significant part of previous studies assumes symmetry of natural frequency distributions, and bifurcation is found as the synchronization transition referring to the transition from the nonsynchronized state to partially synchronized states. Interestingly, symmetry breaking brings new types of bifurcations that emerge not from the nonsynchronized state but from partially synchronized states and are followed by a discontinuous jump or oscillation of the order parameter terada2017 . The new types of bifurcations have been observed by performing numerical simulations, and this paper aims to propose a theoretical explanation of the new types of bifurcations. For analyzing bifurcations, we have three theoretical methods based on the large population limit and the equation of continuity.

The first method is the self-consistent equation for the order parameter of system kuramoto1975 ; strogatz2000 ; pazo-05 ; basnarkov-urumov-07 ; park-kahng-19 . A stationary solution to the equation of continuity is formally written by using an unknown value of the order parameter and the value is determined self-consistently. The self-consistent strategy is powerful to obtain stationary values of the order parameter, but the oscillating states are not captured. Moreover, writing down stationary states is not straightforward when a coupling function between a pair of oscillators has several sinusoidal functions komarov-pikovsky-13 ; komarov-pikovsky-14 .

The second method is the Ott-Antonsen ansatz oa2008 ; oa2009 . This ansatz reduces the equation of continuity to a simpler partial differential equations. This strategy is applicable for the Kuramoto model and its variances lu-etal-16 ; akao-19 , and is extremely useful when we use a single Lorentzian as the natural frequency distribution, thereby resulting a real two-dimensional reduced system. The dimension of the reduced system becomes higher as the distribution becomes complicated martens2009 ; terada2017 ; bastian-18 , and searching for stationary and oscillating states becomes difficult accordingly.

The third method is the amplitude equation, which describes dynamics projected onto the unstable manifold of the nonsynchronized state crawford1994 . In contrast to the previous two methods, the amplitude equation is widely useful in the Kuramoto model crawford1994 , in a coupled phase-oscillator model with a generic coupling function crawford1995 , in the time-delayed Kuramoto model metivier-19 , in an extended Kuramoto model with inertia barre2016 , and in Hamiltonian systems crawford1994b ; crawford1995b ; barre-metivier-yamaguchi-16 . Other advantages of the amplitude equation are that one-dimensional dynamics is sufficient to judge the continuity of bifurcation and that any natural frequency distributions are equally tractable.

The first and second methods profit from a single sinusoidal coupling function, but such a simple coupling function is not always the case hansel-mato-meunier-93 ; crawford1995 ; kiss-zhai-hudson-05 ; kiss-zhai-hudson-06 . We, therefore, adopt the third method, the amplitude equation, to explain the new types of bifurcations. One basic block of the amplitude equation is the linear analysis around the nonsynchronized state, and the linear analysis gives, for instance, eigenvalues. We will demonstrate that the linear analysis is also useful to explain a mechanism for the emergence of an oscillating state.

The method proposed in this paper is in principle applicable to general systems beyond the Kuramoto model and its variances, but we investigate the Kuramoto model to illustrate the usefulness of the method. One reason is that the new types of bifurcations are reported in the Kuramoto model, and this choice makes it straightforward to compare the theory with the previous work terada2017 . Another reason is that, as we discussed above, the reduction by the Ott-Antonsen ansatz permits us to obtain accurate classification, which helps to examine the theory.

This paper is organized as follows. In Sec. II, we introduce the Kuramoto model and its large population limit written by the equation of continuity. The considering family of the natural frequency distribution, which is characterized by two parameters, is also exhibited. In Sec. III, we numerically divide the parameter plane into five domains corresponding to types of bifurcation diagrams and prepare the reference parameter plane to examine the theory. This numerical search is performed by using the Ott-Antonsen reduction. We stress that this reduction is used only for obtaining the reference parameter plane and is not used in our theory. In Sec. IV, linear and nonlinear analyses of the equation of continuity are shortly reviewed. With the aid of these analyses, we propose ideas to identify the domains on the parameter 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,

[TABLE]

The real constant $K>0$ is the coupling constant, $\theta_{j}\in(-\pi,\pi]=\mathbb{S}^{1}$ and $\omega_{j}\in\mathbb{R}$ are respectively the phase and the natural frequency of the $j$th oscillator. The natural frequencies $\{\omega_{j}\}$ obey a probability distribution function $g(\omega)$. By using the parameters $\gamma_{1},\gamma_{2}>0$ and $\Omega\geq 0$, we introduce a family of $g(\omega)$ as

[TABLE]

to systematically consider unimodal and bimodal, and symmetric and asymmetric distributions terada2017 ; terada2018 . The family of $g(\omega)$ consists of rational functions, which is useful to draw the reference parameter plane by using the Ott-Antonsen ansatz, but is not crucial in our methodology. The normalization constant

[TABLE]

is determined from the condition

[TABLE]

The three parameters $\gamma_{1},\gamma_{2},$ and $\Omega$ are assumed to be positive. By scaling the variables $t,\omega_{j},K,$ and $\gamma_{1}$, we may set $\gamma_{2}=1$ without loss of generality, and the family of $g(\omega)$ is characterized by a point on the parameter plane $(\gamma_{1},\Omega)$. The parameter plane is further restricted into the region of $\gamma_{1}\leq 1$ by making the one-to-one correspondence between the regions of $\gamma_{1}>1$ and $\gamma_{1}<1$ as follows. The equation of motion, (1), is invariant under the transformations of $\theta_{j}\to-\theta_{j}$ and $\omega_{j}\to-\omega_{j}$, and the latter induces the transformation of $g(\omega)\to g(-\omega)$. The equation of motion is also invariant under the change of time scale $t\to t/\tau$ with $\omega_{j}\to\tau\omega_{j}$ and $K\to\tau K$. These changes do not modify type of the bifurcation diagram for $\tau>0$ and the transformed natural frequency distribution is

[TABLE]

The point $(\gamma_{1},\Omega)$, therefore, corresponds to the point $(1/\gamma_{1},\Omega/\gamma_{1})$ by selecting $\tau=\gamma_{1}$ under $\gamma_{2}=1$. The line $\gamma_{1}=1$ gives a family of symmetric distributions.

The complex order parameter $z$ is defined by

[TABLE]

The absolute value $r$ of $z$ measures the extent of synchronization. If all the oscillators distribute uniformly on $\mathbb{S}^{1}$, then $r\simeq 0$. If $r\simeq 1$, the majority of oscillators gathers around a point on $\mathbb{S}^{1}$.

In the limit of large population $N\to\infty$, by the conservation of the number of oscillators, the equation of motion (1) can be written in the equation of continuity lancellotti2004 ,

[TABLE]

where $F(\theta,\omega,t)$ is the probability distribution function of $\theta$ and $\omega$ at the time $t$. In other words, $F(\theta,\omega,t)d\theta~{}\mathrm{d}\omega$ represents the fraction of oscillators having phases between $\theta$ and $\theta+d\theta$ and natural frequencies between $\omega$ and $\omega+~{}\mathrm{d}\omega$ at the time $t$. From the normalization condition $\int Fd\theta~{}\mathrm{d}\omega=1$, we have

[TABLE]

In this limit the order parameter is expressed by

[TABLE]

III Numerical Classification on Parameter plane

The aim of this section is to classify numerically the parameter plane $(\gamma_{1},\Omega)$ into five types of bifurcation diagrams terada2018 before developing the theory. Direct $N$-body simulations of the model (1) include finite-size fluctuations, which make it difficult to judge the types of bifurcation diagrams. For eliminating the finite-size fluctuations, we perform the Ott-Antonsen reduction oa2008 ; oa2009 , which reduces the equation of continuity to a real four-dimensional differential equation for the family of $g(\omega)$, (2).

III.1 Ott-Antonsen reduction

The reduced equation by using the Ott-Antonsen ansatz is expressed as terada2017

[TABLE]

The two variables $z_{1}$ and $z_{2}$ are complex and the complex order parameter $z$ is written as

[TABLE]

with the complex constants $k_{1}$ and $k_{2}$ defined by

[TABLE]

$z^{\ast}$ represents the complex conjugate of $z$. In the later computations we set $\gamma_{2}=1$ without loss of generality as we mentioned in the previous section II, but we kept $\gamma_{2}$ free to show the dependence explicitly. See Appendix A for details of the reduction.

III.2 Bifurcation diagrams

Numerical integration of the reduced system (10) is performed by using the fourth-order Runge-Kutta algorithm with the time step $\Delta t=0.01$. For a given set of $(\gamma_{1},\Omega)$, we start from $K=0$ and increase the value up to $K=10$ with the step $\Delta K=0.0001$. This increasing process is called the forward process. At each $K$, the time average and standard deviation of the order parameter are taken in the time interval $t\in[4500,5000]$ to avoid a transient regime. The final state at $t=5000$ is used as the initial state at the successive value of $K$. If the final state is the origin, then we shift the initial state from the origin to $z_{1}=z_{2}=0.01$ to escape from the trivial stationary state. After arriving $K=10$, we decrease the value of $K$ from $K=10$ to $K=0$ to check existence of a hysteresis, which reveals a discontinuous bifurcation. This decreasing process is called the backward process.

The parameter plane $(\gamma_{1},\Omega)$ is classified into five domains as shown in Fig. 1 terada2018 , which domains correspond to the five types of bifurcation diagrams reported in Fig. 2. The symmetry line $\gamma_{1}=1$, in particular, is classified into the three intervals included in the domains A, B, and C. The separating points $\Omega=1$ and $\sqrt{3}$ are obtained by applying the amplitude equation and the eigenvalue analysis respectively, which are explained later. The asymmetry region $\gamma_{1}<1$ includes the new domains D and E with two known domains A and B. The goal of this paper is to reproduce the parameter plane theoretically for understanding mechanisms yielding the new domains D and E.

IV Theoretical analyses of equation of continuity

We shortly review linear and nonlinear analyses of the equation of continuity (7) around the nonsynchronized state $f^{0}$, which is explicitly written as

[TABLE]

It is straightforward to check stationarity of $f^{0}$ as

[TABLE]

from the fact $v[f^{0}]=\omega$. We expand the equation of continuity by substituting

[TABLE]

The perturbation $f$ is governed by the equation

[TABLE]

where the linear part is

[TABLE]

and the nonlinear part is

[TABLE]

with the order parameter

[TABLE]

Note that the equation (16) is an exact transform from the equation of continuity (7). The linear part will be used to explain the oscillating state emerging from a partially synchronized state as well as a basic block of the amplitude equation.

IV.1 Linear analysis

The perturbation $f$ is a periodic function of $\theta$ and is expanded into the Fourier series as

[TABLE]

The linear analysis of (17) can be performed independently in each Fourier mode $k$ because the nonsynchronized state $f^{0}(\omega)$ does not depend on $\theta$. The Fourier modes $k\neq\pm 1$ give only rotations, and instability comes from the modes $k=\pm 1$. The eigenvalues for the modes $k=\pm 1$ are obtained as roots of the spectral functions

[TABLE]

See Appendix B for derivations. If the real part of a root is positive, the eigenvalue induces instability. We call such an eigenvalue as an unstable eigenvalue, which is the target of the amplitude equation introduced in the next subsection IV.2.

We define the synchronization transition point $K_{\rm c}$ at which the nonsynchronized state $f^{0}$ changes stability. This definition suggests that the eigenvalue having the largest real part must be on the imaginary axis at $K_{\rm c}$. A pure imaginary eigenvalue, however, induces a singularity in the integrands of the spectral functions. To avoid the singularity, we perform the analytic continuation of $\Lambda_{\pm 1}(\lambda)$ and denote the continued functions as $D_{\pm 1}(\lambda)$. See Appendix C for the continuation.

We give three remarks. First, the continuation does not modify the spectral functions in the region ${\rm Re}(\lambda)>0$, and a root of $D_{\pm 1}(\lambda)$ whose real part is positive is also a root of $\Lambda_{\pm 1}(\lambda)$. Second, a root of $D_{\pm 1}(\lambda)$, however, may not be an eigenvalue if it is on the region ${\rm Re}(\lambda)\leq 0$, and the root is called a fake eigenvalue accordingly. Third, the relations

[TABLE]

and

[TABLE]

hold, where $\Lambda_{1}^{\ast}(\lambda)$ is, for instance, the complex conjugate of $\Lambda_{1}(\lambda)$. These relations imply that $\lambda^{\ast}$ is a (fake) eigenvalue if $\lambda$ is.

For the family of natural frequency distributions (2), the equation $D_{1}(\lambda)=0$ leads a quadratic equation of $\lambda$ and the two roots are denoted by $\lambda_{1}$ and $\lambda_{2}$ ordered as ${\rm Re}(\lambda_{1})\geq{\rm Re}(\lambda_{2})$. In this paper, the two (fake) eigenvalues $\lambda_{1}$ and $\lambda_{2}$ are called the first and the second (fake) eigenvalues respectively. See Appendix D for the quadratic equation of $\lambda$, the determination of the critical point $K_{\rm c}$, and the number of unstable eigenvalues.

IV.2 Amplitude equation

We assume that the linear operator $\mathcal{L}$ has only one pair of unstable eigenvalues, $\lambda_{1}$ and $\lambda_{1}^{\ast}$ coming from the relation (23), and that the corresponding eigenfunctions are $\Psi(\theta,\omega)$ and $\Psi^{\ast}(\theta,\omega)$ respectively. To derive the amplitude equation, we expand $f$ into

[TABLE]

The amplitude $A$ relates to the order parameter $z$ as

[TABLE]

and the asymptotic value of $z$ is approximately obtained by considering temporal evolution of the amplitude $A(t)$. The function $H$ represents the unstable manifold of the nonsynchronized state $f^{0}$. In other words, $H$ represents the height of the unstable manifold from the eigenspace ${\rm Span}\{\Psi,\Psi^{\ast}\}$. We assume that the unstable manifold is tangent to the eigenspace ${\rm Span}\{\Psi,\Psi^{\ast}\}$ at $A=A^{\ast}=0$ and $H=O(|A|^{2})$ accordingly.

For deriving the amplitude equation, we introduce the adjoint linear operator $\mathcal{L}^{\dagger}$ of $\mathcal{L}$ defined by

[TABLE]

where the inner product $(\cdot,\cdot)$ is defined by

[TABLE]

Let $\widetilde{\Psi}$ and $\widetilde{\Psi}^{\ast}$ be the eigenfunctions of $\mathcal{L}^{\dagger}$ corresponding to the eigenvalues $\lambda_{1}^{\ast}$ and $\lambda_{1}$ respectively. We can choose the eigenfunctions satisfying the relations

[TABLE]

without loss of generality. We also have

[TABLE]

because $H$ does not belong to the eigenspace ${\rm Span}\{\Psi,\Psi^{\ast}\}$. See Appendix E for the explicit expressions of the adjoint operator and the eigenfunctions.

Substituting the expansion (24) into the equation of continuity (7) and using the relations (28) and (29), we have the amplitude equation as

[TABLE]

and the equation for the unstable manifold $H$ as

[TABLE]

These equations can be solved perturbatively for sufficiently small $|A|$, and the right-hand-side of (30) is expanded as

[TABLE]

where the eigenvalue $\lambda_{1}$ and the coefficients $c_{3},c_{5},\cdots$ depend on the coupling constant $K$. Note that the right-hand-side of (32) has only odd order terms. See Appendix F for derivations and the explicit forms of coefficients.

The complex amplitude equation (32) can be reduced to a real equation written as

[TABLE]

where $\sigma=|A|^{2}\geq 0$ and

[TABLE]

We will search for stationary solutions with which the right-hand-side of the amplitude equation (33) is zero. The nonsynchronized state, for instance, corresponds to $\sigma=0$, which always satisfies the stationary condition.

The amplitude equation is useful to determine the continuity of the synchronization transition from the nonsynchronized state $f^{0}$ to a partially synchronized state. Around the critical point $K_{\rm c}$, the order parameter must be small if the transition is continuous and we can use the truncated equation

[TABLE]

This equation has a nontrivial solution for ${\rm Re}(c_{3})<0$ from the instability condition ${\rm Re}(\lambda_{1})>0$, while no real solution for ${\rm Re}(c_{3})>0$ suggests a discontinuous transition. Assuming continuity of $c_{3}(K)$ with respect to $K$ and taking the limit $K\to K_{\rm c}+0$, we say that the synchronization transition is continuous if ${\rm Re}(c_{3}(K_{\rm c}))<0$, and is discontinuous if ${\rm Re}(c_{3}(K_{\rm c}))>0$ barre2016 . The separating point $(\gamma_{1},\Omega)=(1,1)$ between the domains A and B on the line $\gamma_{1}=1$ is obtained by searching for the point which satisfies ${\rm Re}(c_{3}(K_{\rm c}))=0$.

V Criteria for determining domains

Looking back the five bifurcation diagrams exhibited in Fig. 2, we have two elements to characterize the diagrams, which are oscillation and a jump of the order parameter. We first discuss mechanisms of the oscillation and of the jump in Secs. V.1 and V.2 respectively with the aid of the linear and nonlinear analyses reviewed in the previous section IV. After that, we propose a procedure to divide the parameter plane into the five domains in Sec. V.3.

V.1 Oscillation of order parameter

The asymptotically periodic oscillation of the order parameter has been also observed in the Hamiltonian mean-field model, and the oscillation comes from existence of two small clusters running with different velocities morita-kaneko-06 . Inspired by this work, we propose the following idea; The first unstable eigenvalue induces a partially synchronized state and the oscillation is excited by appearance of a second unstable eigenvalue.

The order parameter is computed as

[TABLE]

the time evolution of the order parameter $z(t)$ is hence described by the eigenvalues arisen from the Fourier $-1$ mode. As defined in the end of Sec. IV.1, the two fake eigenvalues, the roots of $D_{-1}(\lambda)$, are denoted by $\lambda_{1}^{\ast}$ and $\lambda_{2}^{\ast}$, where ${\rm Re}(\lambda_{1}^{\ast})\geq{\rm Re}(\lambda_{2}^{\ast})$.

Let us increase the value of the coupling constant $K$ from the nonsynchronized region $K<K_{\rm c}$. Beyond the critical point $K_{\rm c}$, the first fake eigenvalue $\lambda_{1}^{\ast}$ becomes unstable and the resonant oscillators, whose natural frequencies are close to ${\rm Im}(\lambda_{1}^{\ast})$, form a cluster. For instance, if $g(\omega)$ is symmetric and unimodal, the resonant frequency is zero and oscillators around $\omega=0$ form a synchronized cluster. This small cluster corresponds to the continuous synchronization transition in the bifurcation diagram (see Fig. 2 (d)). Further increasing $K$, for suitable pairs of $(\gamma_{1},\Omega)$, the second fake eigenvalue $\lambda_{2}^{\ast}$ also becomes unstable and forms the second cluster. The two resonant frequencies differ in general, ${\rm Im}(\lambda_{1}^{\ast})\neq{\rm Im}(\lambda_{2}^{\ast})$, and the order parameter oscillates; The order parameter is large when the two clusters are close on $\mathbb{S}^{1}$, and is small when the clusters are in the antiphase positions each other. The symmetry of $g(\omega)$, in particular, results in inducing simultaneous destabilization of the two fake eigenvalues $\lambda_{1}^{\ast}$ and $\lambda_{2}^{\ast}$ at $K_{\rm c}$ and yielding the bifurcation diagram Fig. 2 (c).

Summarizing, existence of two unstable eigenvalues, which mean two pairs of the unstable eigenvalues by counting the roots of $D_{1}(\lambda)$ in addition to the roots of $D_{-1}(\lambda)$, suggests the oscillation of the order parameter. To support this mechanism, we investigate $K$ dependence of the fake eigenvalues $\lambda_{1}^{\ast}$ and $\lambda_{2}^{\ast}$ at points chosen from the domains C, D, and E of the parameter plane $(\gamma_{1},\Omega)$.

A symmetric case is examined in Fig. 3 for $(\gamma_{1},\Omega)=(1.0,3.0)$ which belongs to the domain C. As we expected, two eigenvalues become unstable at the same $K_{\rm c}$ with different imaginary parts. In Fig. 1 the separating point $\Omega=\sqrt{3}$ between the domains B and C on the symmetry line $\gamma_{1}=1$ is obtained by checking existence of the two unstable eigenvalues.

For the point $(\gamma_{1},\Omega)=(0.8,3.0)$ belonging to the domain D, the simultaneity of destabilization breaks as shown in Fig. 4 owing to asymmetry of $g(\omega)$. A correspondence between the two unstable eigenvalues and the two clusters is exhibited in Fig. 5, where $N$-body simulations of the Kuramoto model (1) were performed by using the fourth-order Runge-Kutta algorithm with the time step $\Delta t=0.01$. See Appendix G for the procedure to determine the synchronized intervals of $\omega$ reported as red bars in Fig. 5. The second unstable eigenvalue and the second cluster emerge at almost same values of $K$, and this observation is consistent with the discussion above. The scenario holds even when the first cluster is further developed as shown in Fig. 6 for the point $(\gamma_{1},\Omega)=(0.6,4.0)$, although the linear analysis is performed around the nonsynchronized state. It is worth noting that the second fake eigenvalue approaches to zero but does not become unstable at the point $(\gamma_{1},\Omega)=(0.9,1.8)$ close to the domain D but belonging to the domain E as reported in Fig. 7.

We note that the discussion above does not always hold because the second unstable eigenvalue does not always yield the second cluster. At the point $(\gamma_{1},\Omega)=(0.8,2.6)$ belonging to the domain E, at $K=5.032$, the second unstable eigenvalue emerges but its imaginary part is not sufficiently far from the grown first cluster, and the second virtual cluster is absorbed by the first cluster without emerging as shown in Fig. 8. One more condition must be therefore added to characterize the oscillation;

[TABLE]

for some $s$ barre-yamaguchi-09 . See Appendix H for the further investigation on this condition.

We give another remark on multi-cluster states; The maximum number of clusters is not necessarily two for the bimodal natural frequency distribution. Multi-cluster states called the Bellerophon states are reported in bi2016 ; li2019 for symmetric natural frequency distributions. The symmetry is, however, not essential to the multi-cluster states as an example is shown in Fig. 9 for an asymmetric distribution with $(\gamma_{1},\Omega)=(0.82,4.5)$.

V.2 Jump of order parameter

A jump of the order parameter emerges from $r=0$ (domain B) or from $r>0$ (domain E) . We expect that the jump from $r=0$ is identified by ${\rm Re}(c_{3}(K_{\rm c}))>0$ as discussed in Sec. IV.2. This criterion is successfully used in a generalized model barre2016 , but for symmetric natural frequency distributions. We will verify this criterion for asymmetric ones with $\gamma<1$.

To characterize the jump from $r>0$, a typical bifurcation diagram of the type E is schematically shown in Fig. 10. We focus on the saddle-node bifurcation point $K_{\rm Q}$, at which $G(\sigma)$ in the amplitude equation (33) has degenerated nontrivial roots. The degeneracy can be approximately captured by $G_{2}(\sigma)$, the truncation of $G(\sigma)$ up to the second order of $\sigma$, as the vanishing discriminant.

Let us discuss validity of the expansion (32) around the saddle-node bifurcation point $K_{\rm Q}$. The point $K_{\rm Q}$ is approximated by $K_{\rm Q^{\prime}}$ obtained from $G_{2}(\sigma)$. The value of $\sigma$ at $K_{\rm Q}$ is also approximated as

[TABLE]

and gives the equality

[TABLE]

This equality suggests that the expansion (32) up to the third term is not excellent but does not break around the bifurcation point $K_{\rm Q}$. We will discuss validity of (32) later from the view point of smallness of the order parameter.

We remark on the upper stable branch appearing at another saddle-node bifurcation point $K_{\rm P}$. Coexistence of the three branches between $K_{\rm P}$ and $K_{\rm Q}$ is a distinguishing feature of the type E, but we focus on $K_{\rm Q}$ because of the two advantages: (i) Capturing the coexistence in the amplitude equation (33) requires the third order term of $G(\sigma)$, the coefficient $c_{7}$, which is hard to compute. (ii) The upper branch has a larger value of $\sigma$ than $\sigma_{\ast}^{\prime}$, and validity of the expansion (32) becomes worse accordingly.

V.3 Theoretical division of the parameter plane

In the basis of the discussions done in Secs. V.1 and V.2, we divide the parameter plane $(\gamma_{1},\Omega)$ into the five domains by the following flow. For preparation, we compute the critical point $K_{\rm c}$ for a given pair of parameters $(\gamma_{1},\Omega)$ as shown in Appendix D.

We first focus on the critical point $K=K_{\rm c}$. The discontinuous jump appears at $K_{\rm c}$ only in the type B among the considered five bifurcation diagrams (see Fig. 2). The point $(\gamma_{1},\Omega)$ must belong to the domain B if ${\rm Re}(c_{3}(K_{\rm c}))>0$ holds, which suggests a jump from $r=0$. When this condition does not hold, we check if there are two unstable eigenvalues at $K_{\rm c}$. If this condition is satisfied, oscillation of the order parameter starts from $K_{\rm c}$ and the point $(\gamma_{1},\Omega)$ is considered to be in the domain C.

When both checks at $K_{\rm c}$ are negative, we increase $K$ from $K_{\rm c}$ and examine the following two propositions:

[TABLE]

and

[TABLE]

where the discriminant $\Delta(K)$ of $G_{2}(\sigma)$ is defined by

[TABLE]

Note that the degenerated root of $G_{2}(\sigma)$, $-c_{3}/(2c_{5})$, is positive under the assumptions $c_{3}(K_{\rm J})<0$ and $c_{5}(K_{\rm J})>0$. If both the propositions (40) and (41) are false, we decide that the point $(\gamma_{1},\Omega)$ belongs to the domain A. If the oscillation proposition (40) is true but the jump proposition (41) is false, the point must be in the domain D. If the oscillation proposition (40) is false but the jump proposition (41) is true, the point must be in the domain E. If both the propositions are true, there is a competition between $K_{\rm O}$ and $K_{\rm J}$; $K_{\rm O}<K_{\rm J}$ suggests the domain D and $K_{\rm O}>K_{\rm J}$ suggests the domain E. In Table 1, we summarize techniques which are used to determine the domains

The flow proposed above provides three theoretical lines, I, II, and III which divide the parameter plane $(\gamma_{1},\Omega)$ into the five domains, as reported in Fig. 11. the upper side of the line I is the theoretical domain D, the inside of the line II is the theoretical domain B, and the lines I, II, and III enclose the theoretical domain E. The theoretical lines reproduce qualitatively the numerically obtained domains.

For the line I, our result is fairly consistent with numerical ones, but the domain D is overestimated. Our method says that the emergence of two clusters is associated with the emergence of two unstable eigenvalues. The second eigenvalue, however, does not always give rise to the second cluster as we discussed in Sec. V.1. We have to reduce the domain D by introducing an additional condition that $|{\rm Im}(\lambda_{1}^{\ast})-{\rm Im}(\lambda_{2}^{\ast})|$ is sufficiently large.

The line II is perfect because the order parameter is small around the critical point of the synchronization transition and the perturbatively obtained amplitude equation is valid for such a small amplitude. This criterion had been used for symmetric $g(\omega)$, but we have confirmed that it is also powerful for asymmetric ones. The line II also reveals existence of discontinuous synchronization transition for unimodal but asymmetric distributions. The existence has been pointed out numerically in a previous study terada2018 , but our theoretical analysis ensures it.

The line III is not perfect. We used a perturbative method to derive the amplitude equation (32), assuming that the order parameter is sufficiently small at the jumping point $K_{\rm Q}$. This assumption becomes worse as the point $(\gamma_{1},\Omega)$ approaches to the boundary of the domain E where the width of the hysteresis vanishes as reported in Fig. 12. It is thus expectable that the boundary of the domain E is not perfectly described by the theoretical line III. Nevertheless, the theoretical line III captures the domain E qualitatively and is useful to have a strong candidate of the domain E.

VI Conclusion and Discussions

We have proposed a theoretical method to classify bifurcation diagrams in the Kuramoto model by using the amplitude equation with the aid of the linear analysis around the nonsynchronized reference state. The amplitude equation, obtained perturbatively, is usually stopped up to the second leading term because it is sufficient to judge the continuity of the synchronization transition, bifurcation from the nonsynchronized state. The premise, however, is broken by introducing asymmetry in the natural frequency distribution because the asymmetry induces bifurcations from partially synchronized states. We have extended the amplitude equation up to the third leading term and successfully captured a discontinuous bifurcation after a continuous one.

For bifurcation diagrams having oscillating states, we have focused on unstable eigenvalues of the nonsynchronized state. Roughly speaking, one unstable eigenvalue corresponds to one cluster formation, and hence existence of two unstable eigenvalues suggest appearance of two clusters rotating with different speeds. This idea is qualitatively in good agreement with numerical results, but overestimates the domain having the oscillation. We have to reduce the domain by adding one more condition which requires a large discrepancy of rotating speeds of the two clusters in order to avoid the absorption of the second virtual cluster by the first grown cluster. The additional condition has to be investigated. We remark that a similar condition was discussed in a Hamiltonian system barre-yamaguchi-09 .

We mainly used the amplitude equation to predict the bifurcation diagram. This method is a perturbation technique and is not powerful if the considering point is far from the critical point. This restriction prevents us from getting precise theoretical lines for $\Omega$ sufficiently large, for instance, because the oscillation/jump emerging from a partially synchronized state are not close to the critical point. On the contrary, for symmetric and bimodal natural frequency distributions, our method successfully identified the oscillating states, the domain C, solely by a linear analysis with phenomenology, while they have been previously analyzed by a nonlinear analysis, the amplitude equation crawford1994 . Our analysis hence captures an essence of the oscillating states.

It is worth noting that a cluster is created by the resonance between natural frequencies an unstable eigenvalue, but the number of clusters, $n_{c}$, is not necessarily the same with the number of unstable eigenvalues, $n_{u}$. In our computations, $n_{c}=1$ when $n_{u}=1$ but $n_{c}\geq 2$ when $n_{u}=2$. Such a multi-cluster state with $n_{c}>2$, called the Bellerophon state, was reported previously bi2016 ; li2019 but the mechanism is an open question. Noting that two or more clusters induce oscillation of the order parameter and following the scenario that a resonance makes a cluster, we expect that the multi-cluster state with $n_{c}>2$ is realized hierarchically by resonances between natural frequencies and the frequencies of the order parameter oscillation. This scenario has to be investigated.

Finally, we discuss applications of our theory to general coupled oscillator models. Introducing the phase-lag parameter in the coupling function sakaguchi-kuramoto-86 , we have a phase diagram having a continuous synchronization transition followed by a discontinuous jump as the type E omelchenko-wolfrum-12 . Analyzing this system is a straightforward application of our theory. Studying time delay yeung-strogatz-99 ; montbrio-pazo-schmidt-06 is also an interesting application. If the coupling function becomes general, the coefficients of the amplitude equation may have divergences as the coupling constant goes to the critical value as shown in a collisionless plasma crawford1995b . For instance, the coefficient ${\rm Re}(c_{3})$ may be proportional to a power of $1/{\rm Re}(\lambda_{1})$. Nevertheless, the amplitude equation captures the asymptotic scaling of the order parameter against the strength of instability in plasmas, and we may expect that the amplitude equation still has power to predict the bifurcation diagram. It is another future problem to check the precision of the proposed criterion under such divergent coefficients.

Acknowledgements.

Y.Y.Y. acknowledges the support of JSPS KAKENHI Grant No. 16K05472.

Appendix A Ott-Antonsen ansatz and reduction

We derive the reduced system (10) from the equation of continuity (7) by using the Ott-Antonsen ansatz. The Ott-Antonsen ansatz assumes the form of the probability distribution $F(\theta,\omega,t)$ as

[TABLE]

with $\left|a(\omega,t)\right|<1$ to ensure the convergence of the series. By substituting (43) to (7), we have

[TABLE]

where the order parameter is written as

[TABLE]

The integration of the right-hand side can be computed by using the residue theorem because $|g(\omega)a^{\ast}(\omega,t)|<g(\omega)$ holds and $g(\omega)$ of (2) decays faster than $1/|\omega|$ when $|\omega|\to\infty$. The residue theorem gives

[TABLE]

where the constants $k_{1}$ and $k_{2}$ are defined by (12). Introducing the complex variables $z_{1}$ and $z_{2}$ as

[TABLE]

the order parameter is expressed as

[TABLE]

Setting $\omega=\Omega+i\gamma_{1}$ and $\omega=-\Omega+i\gamma_{2}$ in (44), we have dynamics of $z_{1}$ and $z_{2}$ as the equation (10).

Appendix B Spectral functions and eigenfunctions

are obtained as roots of the spectral functions. We derive the spectral functions and the eigenfunction for an eigenvalue. The linear operator $\mathcal{L}$ is expanded into the Fourier series as

[TABLE]

where $\widetilde{f}_{k}$ are the Fourier components of $f$ defined by (20). The linear operator $\mathcal{L}_{k}$ is defined by

[TABLE]

The symbol $\delta_{k,l}$ represents the Kronecker delta.

Once we have an eigenfunction $\psi(\omega)$ of the linear operator $\mathcal{L}_{k}$ satisfying

[TABLE]

we have an eigenfunction $\Psi(\theta,\omega)$ of the linear operator $\mathcal{L}$ as

[TABLE]

which induces

[TABLE]

We therefore discuss eigenvalues and eigenfunctions of the Fourier expanded linear operator $\mathcal{L}_{k}$. The mode $k\neq\pm 1$ has $\mathcal{L}_{k}=-ik\omega$ and has only rotations. From now on, we focus on the modes $k=\pm 1$.

The equation (51) is explicitly rewritten as

[TABLE]

Let $\lambda+ik\omega\neq 0$. Multiplying $(\lambda+ik\omega)^{-1}$ and integrating over $\omega$, we have

[TABLE]

where the spectral functions $\Lambda_{\pm 1}(\omega)$ are defined in (21). The vanishing integral, $\int_{-\infty}^{\infty}\psi(\omega)~{}\mathrm{d}\omega=0$, implies $\psi(\omega)\equiv 0$ for any $\omega\in\mathbb{R}$ from (54) and the assumption $\lambda+ik\omega\neq 0$, but this is not compatible with the assumption that $\psi$ is an eigenfunction. Consequently, the eigenvalue $\lambda$ must satisfy the equation $\Lambda_{k}(\lambda)=0$. The nonvanishing integral may be assumed to be $\int_{-\infty}^{\infty}\psi(\omega)~{}\mathrm{d}\omega=1$ without loss of generality, and the eigenfunction $\psi$ is expressed as

[TABLE]

Appendix C Analytic continuation

We introduce the analytic continuation of the spectral functions $\Lambda_{\pm 1}(\lambda)$, (21). We start from a point $\lambda$ whose real part is positive, ${\rm Re}(\lambda)>0$, and decrease the real part. The condition ${\rm Re}(\lambda)>0$ of the starting point comes from the Laplace transform strogatz1992

[TABLE]

which analyzes temporal evolution and is defined in the positive real half-plane ${\rm Re}(s)>0$ to ensure convergence of the integral. The variable $s$ corresponds to $\lambda$.

When $\lambda$ passes the imaginary axis, the integrands of the spectral functions $\Lambda_{\pm 1}(\lambda)$, (21), meet the singularities at $\omega=\pm i\lambda$. To avoid this singularity at $\omega=i\lambda$ ($\omega=-i\lambda$), we modify continuously the integral contour from the real axis by adding a small lower (upper) half-circle around the singular point. This modification brings a residue part to the integral. Continuing this modification for ${\rm Re}(\lambda)<0$, we have analytically continued functions of $\Lambda_{\pm 1}(\lambda)$ as

[TABLE]

where

[TABLE]

and PV represents the Cauchy principal value. The roots of $D_{\pm 1}(\lambda)$ and of $\Lambda_{\pm 1}(\lambda)$ are identical from their definitions if ${\rm Re}(\lambda)>0$, and hence such roots of $D_{\pm 1}(\lambda)$ are eigenvalues. They do not coincide, however, for ${\rm Re}(\lambda)\leq 0$ in general. A root of $D_{\pm 1}(\lambda)$ with ${\rm Re}(\lambda)\leq 0$ is not an eigenvalue but is called a resonance pole, a Landau pole, or a fake eigenvalue ogawa2013 .

The difference between $\Lambda_{1}(\lambda)$ and $D_{1}(\lambda)$ is demonstrated by considering a Lorentzian natural frequency distribution,

[TABLE]

Straightforward computations give

[TABLE]

and

[TABLE]

$\lambda=K/2-\gamma$ is a root of $D_{1}(\lambda)$ but is not of $\Lambda_{1}(\lambda)$ if $K\leq K_{\rm c}$ and ${\rm Re}(\lambda)\leq 0$ accordingly, where $K_{\rm c}=2\gamma$ is the synchronization transition point. Similarly, we can reproduce the critical point $K_{\rm c}=2/[\pi g(0)]$ for a symmetric unimodal $g(\omega)$ by computing the roots of $D_{\pm 1}(\lambda)$. The critical point $K_{\rm c}$ for the considering family (2) is given in Appendix D.

Appendix D Fake eigenvalues and synchronization transition point $K_{\rm c}$

The fake eigenvalues are roots of the continued spectrum functions, (58). They move continuously on the complex plane, and the synchronization transition point $K_{\rm c}$ is computed as the point where one of the fake eigenvalue passes the imaginary axis from left to right with the smallest coupling constant $K>0$. The fake eigenvalues are analyzed in Appendix D.1 for the natural frequency distribution (2). We will show that the first fake eigenvalue, which has the largest real part, passes the imaginary axis at most once time. After that, the equation to determine $K_{\rm c}$ will be derived in Appendix D.2 with an analysis of the movement of the second fake eigenvalue.

D.1 Fake eigenvalues

Substituting (2) to (62), we obtain the explicit form of $D_{1}(\lambda)$ as

[TABLE]

where

[TABLE]

The equation $D_{1}(\lambda)=0$ induces the quadratic equation

[TABLE]

where

[TABLE]

To write down the two solutions to (70), we introduce the complex variable

[TABLE]

The argument $\theta$ is in the interval $[\pi,2\pi]$ from $a_{\rm I}\leq 0$. The two fake eigenvalues $\lambda_{1}$ and $\lambda_{2}$, which satisfy ${\rm Re}(\lambda_{1})\geq{\rm Re}(\lambda_{2})$, are written as

[TABLE]

as $\cos(\theta/2)\leq 0$. Note that the signs of the imaginary parts are ${\rm Im}(\lambda_{1})\leq 0$ and ${\rm Im}(\lambda_{2})\geq 0$ from $\sin(\theta/2)\geq 0$ and that they are consistent with Figs. 3-8. We remark that the frequency of the order parameter corresponds to the imaginary part of not $\lambda_{1}$ but $\lambda_{1}^{\ast}$.

The synchronization transition point $K_{\rm c}$ is determined by the equation ${\rm Re}(\lambda_{1}(K_{\rm c}))=0$. Let us show that the solution is at most one. Using the relation

[TABLE]

and the definitions of $\rho$ and $\theta$,

[TABLE]

we have

[TABLE]

The functions $b(K),$ $a_{\rm I}^{2}(K)$, and $b^{2}(K)+4a_{\rm R}(K)$ are increasing functions of $K$ for $K>0$ from

[TABLE]

These increasing functions imply that the real part ${\rm Re}(\lambda_{1}(K))$ is also an increasing function of $K$ for $K>0$ and takes zero at most once time.

D.2 Synchronization transition point $K_{\rm c}$

We show that there exists the unique solution to the equation ${\rm Re}(\lambda_{1}(K))=0$, which determines the synchronization transition point $K_{\rm c}$. At $K_{\rm c}$, the fake eigenvalue $\lambda$ must be pure imaginary of the form $\lambda=i\lambda_{\rm I}~{}(\lambda_{\rm I}\in\mathbb{R})$. Substituting this form with $K=K_{\rm c}$ into the quadratic equation (70), we have

[TABLE]

The real part reads

[TABLE]

and the imaginary part reads

[TABLE]

Eliminating $\lambda_{\rm I}$ from (79) and (80), we have the cubic equation to determine $K_{\rm c}$ as

[TABLE]

The number of real solutions to (81) is one or three, and all the real solutions are positive because the left-hand side of (81) is always negative for $K_{\rm c}<0$. The synchronization transition point is determined as the smallest real solution to (81).

To investigate the number of unstable eigenvalues, we introduce the discriminant $\Delta_{3}$ for the cubic equation (81). In general, the discriminant is defined as

[TABLE]

for the cubic equation

[TABLE]

The number of real solutions is one for $\Delta_{3}<0$ and is three for $\Delta_{3}>0$.

In the case $\Delta_{3}<0$ the first fake eigenvalue $\lambda_{1}$ passes the imaginary axis at the unique real solution $K_{\rm c}$ and no other passing occurs. In the case $\Delta_{3}>0$ we have three different real solutions of $K_{\rm c}^{(1)}$, $K_{\rm c}^{(2)}$, and $K_{\rm c}^{(2)}$, where $K_{\rm c}^{(1)}<K_{\rm c}^{(2)}<K_{\rm c}^{(3)}$. By the definition, the first fake eigenvalue $\lambda_{1}$ passes the imaginary axis at $K_{\rm c}=K_{\rm c}^{(1)}$ and it does not pass the imaginary axis any more as shown in the previous subsection D.1. As a result, the other two solutions $K_{\rm c}^{(2)}$ and $K_{\rm c}^{(3)}$ are realized by the second fake eigenvalue $\lambda_{2}$; It passes the imaginary axis from left to right at $K_{\rm c}^{(2)}$, and from right to left at $K_{\rm c}^{(3)}$, as observed in Fig. 4.

Appendix E Adjoint linear operator

The idea of deriving the amplitude equation presented in Sec. IV.2 is to project the full dynamics onto the reduced space spanned by some eigenfunctions of the linear operator $\mathcal{L}$, (17). The projection operator is defined by eigenfunctions of the adjoint operator of $\mathcal{L}$. In this Appendix, a necessary eigenfunction of the adjoint operator $\mathcal{L}^{\dagger}$ is presented.

The adjoint operator $\mathcal{L}^{\dagger}$ acts on $f$ as

[TABLE]

where

[TABLE]

As done for the linear operator $\mathcal{L}$, we expand $\mathcal{L}^{\dagger}$ into the Fourier series as

[TABLE]

The linear operator $\mathcal{L}^{\dagger}_{k}$ is defined by

[TABLE]

We focus on the modes $k=\pm 1$.

Let $\mu$ be an eigenvalue of $\mathcal{L}_{k}^{\dagger}$ and $\widetilde{\psi}(\omega)$ be the corresponding eigenfunction which satisfies

[TABLE]

This equation brings

[TABLE]

Repeating the same discussion done in Appendix B, the eigenvalue $\mu$ must satisfy

[TABLE]

This equation implies that $\mu^{\ast}$ is an eigenvalue of $\mathcal{L}^{\dagger}_{k}$ if $\mu$ is an eigenvalue of $\mathcal{L}_{k}$.

Let $\lambda$ be an eigenvalue of $\mathcal{L}_{1}$ and $\psi$ be the corresponding eigenfunction (56). From the above discussion, $\mathcal{L}_{1}^{\dagger}$ has an eigenvalue $\lambda^{\ast}$ and the corresponding eigenfunction is computed as

[TABLE]

As the case of $\mathcal{L}$, $\lambda^{\ast}$ is also an eigenvalue of $\mathcal{L}^{\dagger}$ and the corresponding eigenfunction is

[TABLE]

which satisfies the normalization condition

[TABLE]

Appendix F Derivation of the amplitude equation

The coefficients of the amplitude equation (32) are obtained perturbatively with an expression of the unstable manifold. We give explicit forms of the coefficients.

F.1 Derivations of equations for $A$ and $H$

We assume that $\lambda$ is the unique unstable eigenvalue of $\mathcal{L}_{1}$ and $\psi(\omega)$ is the corresponding eigenfunction. The relation $\Lambda_{-1}(\lambda^{\ast})=\Lambda_{1}^{\ast}(\lambda)$ implies that $\lambda^{\ast}$ is an eigenvalue of $\mathcal{L}_{-1}$ and $\psi^{\ast}(\omega)$ is the corresponding eigenfunction. The linear operator $\mathcal{L}$, therefore, has two unstable eigenvalues $\lambda$ and $\lambda^{\ast}$, and the corresponding eigenfunctions respectively

[TABLE]

Using these eigenfunctions, we expand the perturbation $f$ into the form of (24),

[TABLE]

We assume that the unstable manifold $H$ is tangent to the unstable eigenspace ${\rm Span}(\Psi,\Psi^{\ast})$ at $A=A^{\ast}=0$. Substituting this expansion into the equation of continuity, we have

[TABLE]

In (96), taking the inner product with $\widetilde{\Psi}$, we obtain the equation for $A$ as

[TABLE]

Extracting this equality and its complex conjugate from (96), we have the equation for $H$ as

[TABLE]

The left-hand side of (98) is read as

[TABLE]

The right-hand side of (97) starts from the linear term $\lambda A$, while one of (98) starts from the quadratic terms of $A$ and $A^{\ast}$ by the tangency assumption of the unstable manifold, $H=O(|A|^{2})$. We, therefore, solve the two equations (97) and (98) perturbatively by assuming that $|A|$ is small.

F.2 Fourier series expansion of $H$

Before going to the Taylor series expansion with respect to $|A|$, we expand $H$ into the Fourier series as

[TABLE]

The rotational symmetry of the system yields crawford1994

[TABLE]

where $\sigma=|A|^{2}$.

Now we expand $h_{n}$ into the Taylor series. For sufficiently small $\sigma$ we expand

[TABLE]

Algebraic computations gives the equation of $A$ as

[TABLE]

where

[TABLE]

and

[TABLE]

For simplicity of notation, we introduce a symbol

[TABLE]

We compute $\langle\widetilde{\psi},h_{2,0}\rangle,~{}\langle\langle h_{1,0}\rangle\rangle,~{}\langle\widetilde{\psi},h_{2,1}\rangle,~{}\langle\langle h_{1,1}\rangle\rangle$ and $\langle\widetilde{\psi},h_{2,2}\rangle$ by expanding (98) into the Fourier series. Four Fourier modes of $n=1,2,3$ and $4$ are sufficient to compute $c_{3}$, $c_{5}$, and $c_{7}$ as shown in Fig. 13.

F.3 Functions $h_{n,m}$

The leading order of the Fourier mode $n$, $h_{n,0}(\omega)$, is expressed as

[TABLE]

The other $5$ functions appearing in Fig. 13 solve the following equations:

[TABLE]

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

F.4 Inner products and integrals in the coefficients

The expression (107) and the eigenfunction $\widetilde{\psi}$, (91), give the inner product $\langle\widetilde{\psi},h_{2,0}\rangle$ as

[TABLE]

To express the other quantities, we introduce two integral operators of

[TABLE]

and

[TABLE]

Remarking that $\langle\langle h_{1,0}\rangle\rangle$ and $\langle\langle h_{1,1}\rangle\rangle$ are solved self-consistently and using the relation

[TABLE]

we have the necessary quantities as

[TABLE]

[TABLE]

[TABLE]

and

[TABLE]

where all the integral operators, $\mathcal{I}_{m}$ and $\mathcal{J}_{m}$, are evaluated at $\lambda$.

Appendix G Procedure to determine the synchronized intervals of $\omega$

This Appendix summarizes the procedure to determine the intervals of $\omega$, which we call the synchronized intervals, in which oscillators synchronize. The synchronized intervals are identified through dividing the $\omega$-axis into bins of the width $w$. The bin $\alpha$ defined by $[\alpha w,(\alpha+1)w)$ induces the set $P_{\alpha}$ consisting of the indices of the oscillators included in this bin as

[TABLE]

We set $w=0.1$ here.

The oscillators in a bin gather around a point on the $(\theta,\omega)$ plane, when the bin is in the synchronized intervals. In the Kuramoto model, the gathering point with $\omega=(\alpha+1/2)w$ is unique for the bin $\alpha$, and the standard deviation of the phases $\{\theta_{i}\}_{i\in P_{\alpha}}$, denoted by $s_{\alpha}$, must be sufficiently small. In contrast, the standard deviation $s_{\alpha}$ must be large when the bin is not in the synchronized intervals. We note that the standard deviation $s_{\alpha}$ is defined by

[TABLE]

where $|P_{\alpha}|$ is the number of elements of the set $P_{\alpha}$ and $\mu_{\alpha}$ is the mean in the bin $\alpha$ as

[TABLE]

The standard deviation and the mean are computed after taking mod $2\pi$ for each $\theta_{i}$.

We have to be careful for picking up all the bins belonging to the synchronized intervals because the standard deviation $s_{\alpha}$ becomes large if the gathering point is close to [math] or $2\pi$. To overcome this difficulty, we prepare the shifted phases

[TABLE]

and compute the standard deviation $\widetilde{s}_{\alpha}$ for $\{\widetilde{\theta}_{i}\}_{i\in P_{\alpha}}$. By using $s_{\alpha}$ and $\widetilde{s}_{\alpha}$, we conclude that the bin $\alpha$ is in the synchronized intervals if

[TABLE]

holds with $s_{\ast}=0.1$.

Appendix H Criterion of the imaginary difference (37)

The domain D is characterized by oscillation of the order parameter. As discussed in Sec. V.1, the oscillation must have one more condition to satisfy, the imaginary difference (37), to ensure emergence of the second cluster without being absorbed by the first cluster. We will show a trial to introduce this criterion to shave the theoretically obtained domain D (see Fig. 11).

For the Kuramoto model with a unimodal and symmetric natural frequency distribution, only one unstable eigenvalue for $K>K_{\mathrm{c}}$ exists, resulting in only one cluster and the continuous transition. The half width of the cluster along the $\omega$ axis is calculated as $Kr$ strogatz2000 . Inspired by this estimation, we introduce the difference criterion as

[TABLE]

where $K_{2}$ is the emerging point of the second unstable eigenvalue. We estimate $r(K_{2})$, the value of order parameter at $K=K_{2}$, from the solution of the amplitude equation up to the third order, that is

[TABLE]

In Fig. 14, a domain D restricted by the criterion (126) is compared with the original theoretical domain D. The criterion shaves a restricted area of the domain D, which is close to the original theoretical line I. The criterion still needs to be investigated.