Collective dynamics of identical phase oscillators with high-order coupling Can Xu,1 Hairong Xiang,1 Jian Gao,1 and Zhigang Zheng2, ∗ 1 Department of Physics and the Beijing-Hong Kong-Singapore Joint Center for Nonlinear and Complex Systems (Beijing), Beijing Normal University, Beijing 100875, China 2 College of Information Science and Engineering, Huaqiao University, Xiamen 361021, China arXiv:1604.04816v1 [nlin.AO] 17 Apr 2016 (Dated: November 24, 2021) Abstract In this paper, we propose a framework to investigate the collective dynamics in ensembles of globally coupled phase oscillators when higher-order modes dominate the coupling. The spatiotemporal properties of the attractors in various regions of parameter space are analyzed. Furthermore, a detailed linear stability analysis proves that the stationary symmetric distribution is only neutrally stable in the marginal regime which stems from the generalized time-reversal symmetry. Moreover, the critical parameters of the transition among various regimes are determined analytically by both the Ott-Antonsen method and linear stability analysis, the transient dynamics are further revealed in terms of the characteristic curves method. Finally, for the more general initial condition the symmetric dynamics could be reduced to a rigorous three-dimensional manifold which shows that the neutrally stable chaos could also occur in this model for particular parameter. Our theoretical analysis and numerical results are consistent with each other, which can help us understand the dynamical properties in general system with higher-order harmonics couplings. ∗ Electronic address: zgzheng@hqu.edu.cn 1 Large system of coupled oscillators occur in a wide variety of situation throughout the nature which has attracted much attention from the scientific community during the last decades [1]. Examples are including the flashing of fireflies [2], electrochemical and spin-toque oscillators [3, 4], pedestrians on footbridges [5], applauding person in a large audience [6] and others. Understanding the cooperative dynamical properties of such system is therefore of considerable theoretical and experiment interest. Indeed, in the weak interaction limit, the dynamics of limit-cycle oscillators could be effectively described in terms of their phase variable θ, while the most famous case is the Kuramoto model [7] which stands out as the classical paradigm for studying the spontaneous emergence collective synchronization in such system [8]. The form of phase equation obeys θ̇i = ωi + N X j=1 Γ(θj − θi ) , (1) where θi denotes the phase of the ith oscillators, ωi is its natural frequency. Γ(θ) is 2π-periodic function representing the interaction between units. The simple choice of Γ(θ) = (K/N ) sin θ leads to the classical Kuramoto model. Along the past decades the Kuramoto model with its generalizations have inspired and simulated a wealth of extensive studies because of both their simplicity for mathematical treatment and their relevance to practice [9–11]. In particular, when Γ(θ) includes higher harmonics the system exhibits nontrivial dynamical features which were reported in the recent works [12–15]. In many realistic systems, the higher harmonic term (especial the second) always plays an essential role in the interaction and even dominates the coupling function, such as the Huygens pendulum system [16], the neuronal oscillators and genetic networks system [17], the globally coupled photochemical oscillators system [18], etc. In contrast to the pervious discussions of Kuramoto model containing higher order coupling. In the present work, we study phase equations of the following form θ̇i = ω − λ sin 2θi + N σ X sin 2θj , i = 1, · · · , N, N j=1 (2) where ω is the frequency of identical oscillators, λ and σ are the coupling strength respectively. There are various motivations for the study of this model, for examples, in the Josephson junction arrays [19, 20] the dynamics of a single element is determinated by a double well potential and therefore the strong effects caused by the second harmonics is important. Another example is the star-like model with a central element while the case of single harmonic term was considered in the paper [21–29]. Similar to mean-field coupling in the Kuramoto model, the interaction term in Eq.(2) is a driving force that does not depends on the phase of driven oscillators and is equal on every oscillator. In this paper, we present a complete framework to investigate the collective dynamics of globally coupled phase oscillator when higher-order modes dominate the coupling function. It includes several aspects, which together presents a global picture for the understanding of the dynamical properties in this system. First, we use the Ott-Antonsen method [30] to obtain the lowdimensional description of symmetric dynamics, where various states are illustrated schematically in the phase-diagram (Fig.1), they are double-center, single-center, center-synchrony coexistence, and synchrony regime respectively. Furthermore, a detailed linear stability analysis which is based 2 on the self-consistent method is implemented and both the boundary curves and eigenvalues of steady states are obtained analytically which are consistent with the Ott-Antonsen method. Additionally, it has been proved that the linearized operator for the stationary symmetric distribution has infinite many pure imaginary eigenvalues which implies the stationary symmetric distribution is neutrally stable to perturbation in all the directions. Second, the two-cluster synchrony state is determined which is initial values dependent and the general transient solutions of the distribution are calculated in term of the characteristics method. Finally, the general initial conditions for the phase oscillators which lie off the poisson submanifold are considered where the symmetric dynamics are governed by the Möbius transformation which are rigorous three dimensional. And therefore one can expect the chaotic behavior of both the symmetric dynamics r2 and the degree of asymmetry r1 occur in the marginal regime. Extensive numerical simulations have been carried out to verify our theoretical analysis. In the following we report our main results both theoretically and numerically. Results Symmetric dynamics. We start by considering the high-order coupled phase oscillators model (2), without loss of generality the range of the coupling strength is restricted to λ > 0 and −∞ < σ < +∞ throughout the paper. The most important characteristic of the current model is the introduction of higher harmonics in the coupling function, and hence the generalized order parameter is needed to symbol the collective behavior of the system [12] which yields rn = Rn e iΘn N 1 X inθj = , e N j=1 (3) for n ∈ integer, where Rn is the magnitude of the complex, n-th order parameter, and Θn is its phase. As named in [12], the amplitude R2 measures the level of cluster synchrony while R1 measures the degree of asymmetry in clustering. Eq. (2) can be rewritten in terms of r2 as θ̇i = ω − λ 2iθi (e − e−2iθi ) + σ · Im(r2 ) , 2i (4) Im represents the imaginary pant. In the thermodynamic limit N → ∞, Eq. (4) is equivalent to the continuity equation as a consequence of the conservation of the number of oscillators, i.e,    ∂ρ λ 2iθ ∂ −2iθ ρ · ω − (e − e + ) + σ · Im(r2 ) = 0. (5) ∂t ∂θ 2i Here ρ(θ, t)dθ gives the fraction of oscillators which lie between θ and θ + dθ at time t with R 2π the appropriate normalization condition 0 ρ(θ, t)dθ = 1, as a result the continuity limit of the generalized order parameter takes the form Z 2π einθ ρ(θ, t)dθ. (6) rn (t) = 0 Additionally, considering the 2π-periodic of θ in the distribution function ρ(θ, t) which allows a Fourier expansion and can be written as ρ(θ, t) = ∞ 1 X rn (t)e−inθ = ρs + ρa . 2π n=−∞ 3 (7) It is obvious that the n-th Fourier coefficient of ρ(θ, t) is just the n-th order parameter Eq. (6). Here ρs is the sum where n is even and is symmetric in the sense that ρs (θ + π, t) = ρs (θ, t), and ρa is the sum where n is odd which is antisymmetric with respect to the translation by π, ρa (θ + π, t) = −ρa (θ, t). Then, substituting the expansion Eq. (7) into the continuity Eq. (5) we obtain a set of two equations   iλ iλ (8) · r2n+3 − r2n−1 + (ω + σIm(r2 )) · r2n+1 , ṙ2n+1 = i(2n + 1) 2 2 for the odd Fourier coefficient and   iλ iλ · r2n+2 − r2n−2 + (ω + σIm(r2 )) · r2n , ṙ2n = i(2n) 2 2 (9) for the even Fourier coefficient. Eq. (8) together with Eq. (9) provide two set of infinite many coupled equations which evolve independently. For instance, the motion of r2 is not only dependent on itself but also governed by r4 and r0 . However, observing the specific form of Eq. (8) and Eq. (9) one notices that Eq. (8) has a trivial invariant manifold solution rn ≡ 0, n ∈ odd, (10) and Eq. (9) has a non-trivial invariant manifold solution rn ≡ r2n , n ∈ even, (11) which is indeed the Ott-Antonsen ansatz [30, 31]. The Ott-Antonsen method yields a special solution for system provided r2 evolves according to a single ordinary differential equation ṙ2 = −λr22 + λ + 2i(ω + σIm(r2 )) · r2 , (12) solution of this kind turn out to form a two-dimensional invariant manifold which is the set of Poisson kernels 1 1 − R22 ρs (2θ, t) = . (13) 2π 1 + R22 − 2R2 cos(Θ2 − 2θ) One central issue in the study is to identify all the possible collective states both steady or nonstationary and reveal various bifurcations as the change of the coupling parameter λ and σ when the initial symmetric distribution has the form of Possion kernels. The Riccati equation Eq. (12) describes the collective symmetric dynamics of Eq. (2) in terms of the second order parameter r2 , and it can be straightforwardly rewritten in cartesian coordinates r2 = x + iy as ẋ = −2(ω + σ · y − λ y)y + λ(1 − x2 − y 2 ) , ẏ = 2(ω + σ · y − λ y)x . (14) (15) In the phase space of the second order parameter, the natural boundary is x2 + y 2 = 1 and a fixed point is determined by the intersection of nullclines ẋ = 0 and ẏ = 0 within the boundary. One re2 calls that the first point is y01 = ω/(λ−σ) and x201 +y01 ≡ 1 which is on the unit circle and defined as the synchrony state, the existence condition for the synchrony state is |λ − σ| ≥ ω. Additionally, the linearization Jacobian matrix of the synchrony state has two eigenvalues δ1 = −2λx01 and 4 δ2 = −2(λ − σ)x01 . The same strategy such as the existence and p the stability conditions can be actually adopted for the second fixed point x02 = 0, y02 = (ω ± ω 2 − λ(λ − 2σ) )/(λ−2σ), which is in the y axis and termed as the splay state [32]. The Jacobian matrix of the second point has two pure imaginary eigenvalues which shows that the splay state is neutral stable to perturbation and is a center in the Ott-Antonsen manifold. Fig. 1 summarizes the results of our analysis of Eq. (14) and Eq. (15). We find that there are four types of regime in the phase diagram, they are the double-center region I, single-center region II, center-synchrony coexistence region III and the synchrony region IV respectively. In addition, various bifurcations and transitions among the states are illustrated when the parameters pass through the boundary curves. There are two kinds of route to synchrony, for example, when we start from II to IV (arrow 1 in the Fig. 1), the transition to synchrony is a first-order phase transition with hysteresis because of the coexistence region III. however, when we turn to the second route (arrow 2 in the Fig. 1) the transition is discontinuity which is absent of hysteresis. The analysis above reveals the low-dimensional collective behavior of symmetric dynamics where the bifurcations and transition among various states are presented. However, all the results are within the framework of Ott-Antonsen invariant manifold including the linear perturbation of the fixed points. In the following, we conduct α through linear stability analysis of the symmetric dynamics in terms of the self-consistent theory [33] where all the boundary curves could be obtained analytically in an alternative way. In particular, we show that when the perturbation is not within the invariant manifold the eigenvalues of the steady states keep the same form above and the splay state is neutral stable in all the directions. A convenient way of solving for the symmetric dynamics is to make the change of variable ϕ = 2θ, which yields a new dynamical equation N 2σ X ϕ̇i = 2ω − 2λ sin ϕi + sin ϕj , N j=1 i = 1, · · · , N. (16) Eq. (16) has the form of linear Josephson junction arrays and it is obviously that the distribution of ϕ is equivalent to the ρs . The synchrony state is a spatially homogeneous fixed point of Eq. (16) defined as ϕ̇i = 0, for all the i and ϕ1 = ϕ2 = · · · = ϕN = ϕ0 which is the simplest attractor. Therefore, this phase leads to a solution of equation ω − (λ − σ) sin ϕ0 = 0. (17) In this case sin ϕ0 = ω/(λ − σ), the Jacobian matrix of Eq. (16) is a circulant matrix [23, 33] which has two kinds of eigenvalues, the first one is δ1 = −2(λ − σ) cos ϕ0 , (18) that corresponds to a spatial homogeneous fluctuation and the other eigenvalues δ2 = −2λ cos ϕ0 , (19) which have (N − 1)-fold degeneracy that correspond to inhomogeneous fluctuations. Incidentally, one notes that x01 ≡ cos ϕ0 , y01 ≡ sin ϕ0 , this result is consistent with a basic fact that the synchrony state is contained in the Ott-Antonsen manifold. Numerical experiment shows that the basin of attraction of the synchrony state in IV is globally. 5 Another scenario is the stationary distribution where the phase ϕ is smoothly distributed over [0, 2π] and from Eq. (16) one observes that the general form of such a stationary distribution is ρs (ϕ) = C , ω̄ − 2λ sin ϕ (20) where ω̄ is the effective frequency which could be determined self-consistently from the equation Z 2π C · sin ϕ , (21) dϕ ω̄ = 2ω + 2σ ω̄ − 2λ sin ϕ 0 √ and C is the normalization constant C = ± ω̄ 2 − 4λ2 /2π when ω̄ > 0, C is a positive one and vice verse. Substituting the expression of C into Eq. (21) and after some calculation. We obtain the solution of ω̄ i h p 2 ω(λ − σ) ± σ 2 (ω 2 − λ2 + 2λσ) . (22) ω̄ = λ − 2σ It should be emphasized that the stationary distribution ρs (ϕ) corresponds to the splay state in the invariant manifold above, because the distribution Eq. (20) has the form of Possion kernels and Z 2π C · sin ϕ hsin ϕi = dϕ = y02 , (23) ω̄ − 2λ sin ϕ 0 and hcos ϕi = Z 2π 0 C · cos ϕ dϕ = 0 ≡ x02 . ω̄ − 2λ sin ϕ (24) The perturbation of the continuum equation for the stationary symmetric distribution ρs (ϕ) is ∂ (δρ) =L̂ δρ(ϕ, t) ∂t Z ∂ ∂ρs 2π ′ =− [(ω̄ − 2λ sin ϕ) δρ(ϕ, t)] + dϕ 2σ sin ϕ′ δρ(ϕ′ , t), ∂ϕ ∂ϕ 0 (25) and the eigenvalue equation Eq. (25) is convenient to treat in the function space ∞ X δρ(ϕ, t) = an (t)e2πinG(ϕ) , ρs (ϕ) n=−∞ (26) where G(ϕ) is the basis function G(ϕ) = Z ϕ ρs (ϕ′ ) dϕ′ , (27) 0 and an (t) are the expansion coefficients. From the stability analysis of the stationary distribution (all the details are included in the supplementary material) we find that in the regime I, II and III, all the infinite many eigenvalues of the operator L are pure imaginary which implies that the stationary distribution is neutrally stable in all the directions and it is not only confined to OttAntonsen manifold. The regime where the stationary distribution is marginally stable is termed as the marginal regime and in this regime there is no particular attractor that the system converges 6 to [33], as a result, the highly non-generic property is associated with the time-reversal symmetry that the Eq. (16) exhibits. When we start in the Ott-Antonsen manifold the trajectory of r2 is two-dimensional closed periodic curve (the insertion of Fig. 1), however, the situation differs substantially when we start with a more general initial condition in the marginal regime as we see in the following part. The steady and transient dynamics. The analysis above investigate the symmetric dynamics by using two kinds of ways, one recalls that when the parameters are in regime IV, the symmetric dynamics converge to steady state, and the two cluster synchrony states emerge sin 2θ0 = ω/(λ − σ), cos 2θ0 > 0 accordingly. Thus, the complete steady state distribution of oscillators is 1 1 ρ0 (θ) = ( + c)δ(θ − θ0 ) + ( − c)δ(θ − θ0 − π), 2 2 1 |c| < , 2 (28) at this time ρ0 (θ) is comprised of two delta functions denoting two clusters of oscillators at θ0 = 0.5 arcsin ω/(λ − σ) and θ0 + π. Hence, the phase oscillators settle to one of the two stable equilibria while the unstable equilibria π/2 − θ0 and 3π/2 − θ0 serve as the boundaries for the basin of attraction. Therefore the degree of asymmetry |r1 | is |r1 | = |2c(cos θ0 + i sin θ0 )| = 2c , (29) which depends on the free parameter c and could be determined approximatively from initial condition, note that 1/2 − c is just the fraction of oscillators in the locked phase θ0 + π, namely 1 c= + 2 Z 3π −θ0 2 π −θ0 2 ρ(θ, t0 ) dθ + ε, (30) ρ(θ, t0 ) is the initial density, ε is the error caused by the Arnold diffusion. Theoretically, when the initial phase is in the one dimensional invariant manifold θ1 = θ2 = · · · = θN the evolution of phase θ(t) for all the oscillators can never pass through the two saddle points π/2−θ0 and 3π/2−θ0 which means ǫ = 0. However, for the more general initial conditions some oscillators which are in the neighborhood of two unstable equilibrium states can bypassing the saddle points (Fig. 2(b) illustrated schematically the mechanism in the low dimensional phase space where N = 3) and therefore ε is non-ignorable. Fig. 2(a) presents the numerical simulation when we choose θ0 = π/6 and ρ(θ, t0 ) = (1+cos θ)/2π, it is clear that the initial phases in the gray area eventually bypassing the saddle points (the pink hollow circle in the horizontal axis) and this interesting phenomenon implies that those oscillators tend to choose a relatively near equilibrium state to settle in the high dimensional phase space while this is forbidden in the one dimensional invariant manifold. When the symmetric dynamics r2 gets to steady state, then the |r1 | dynamics reaches steady state quickly. However, in a large marginal regime of the phase diagram Fig. 1, the dynamics of r2 can never converge to an attractor, the motion of r2 is time-dependent. To capture the dynamics of r1 we can solve the partial differential equation (PDE) Eq. (5) ∂ρ ∂ρ + (ω − λ sin 2θ + σ · Im(r2 )) = λ cos 2θ · ρ , ∂t ∂θ (31) in terms of the characteristics method, when we start in the characteristic curve θ(t, t0 ), the distribution function will be ρ(θ, t) ≡ ρ(θ(t, t0 ), t) and the PDE becomes the ODEs, the characteristic 7 equations are dρ = 2λ cos 2θ · ρ , dt θ̇ = ω − λ sin 2θ + σ · Im(r2 ), (32) (33) when the symmetric dynamics are in the Ott-Antonsen invariant manifold, the motion of r2 is governed by the equation Eq. (12). Generally, the Eq. (32) and Eq. (33) are difficult to solve analytically while for some particular situation such as the orbit of r2 is near to the center point, the amplitude of r2 is small enough that Im(r2 ) approximates a constant, then the expression for the characteristic curves θ(t, t0 ) starting at the initial phase θ0 yields ( ) p p λ + (ω + σIm(r2 ))2 − λ2 · tan[ (ω + σIm(r2 ))2 − λ2 (θ0 + t)] θ(t, t0 ) = arctan (34) ω + σ · Im(r2 ) and the distribution along the characteristic equations is Z t  ′ ′ ρ(t) =ρ0 exp 2λ cos 2θ(t , t0 ) dt t0 "Z # θ(t) ′ dt =ρ0 exp 2λ cos 2θ(t′ , t0 ) · · dθ dθ θ0 = (35) ρ0 (ω − λ sin 2θ0 + σ · Im(r20 )) , ω − λ sin 2θ(t, t0 ) + σ · Im(r2 ) where ρ0 is the initial value of the distribution function, r20 is the initial value of r2 . Theoretically, the general form of density function ρ(θ, t) could be determined by substituting the inverse solution θ0 (θ(t)) Eq. (34) into Eq. (35), and the generalized order parameter could be calculated through the integral Eq. (6) while the difficult is due to the multivalue of anti-trigonometric function Eq. (34). Hence it is convenient to calculate the first-order parameter r1 via the integral Z π ∂θ(t, t0 ) ρ(t)eiθ(t,t0 ) · r1 (t) = · dθ0 , (36) ∂θ0 −π From Fig. 3(a) we find that |r1 (t)| oscillates regularly with a period T =p π (ω + σ · Im(r2 ))2 − λ2 . (37) when the amplitude of r2 is considerable large, the time dependent of r1 is irregular Fig. 3(b). The discussion of the symmetric dynamics Eq. (16) exhibits a two-dimensional Ott-Antonsen manifold proving that the initial symmetric distribution takes the form of Possion kernels Eq. (13). In fact, the significant works point out that the governing equations for the form of Eq. (16) are generated by the action of the Möbius group [34, 35] eiϕj (t) = eiψ eiφj + α , 1 + α∗ eiψ · eiφj 8 (38) α is a complex variable, ψ real, and φj motion constant. The group action partition the N -dimensional state space into three-dimensional invariant manifold and the three parameters Re(α), Im(α), ψ are governed by the following equation ˙ Re(α) = −2(ω + σ · Im(r2 ) − λ · Im(α)) · Im(α) + λ(1 − |α|2 ), ˙ Im(α) = 2(ω + σ · Im(r2 ) − λ · Im(α)) · Re(α), ψ̇ = 2(ω + σ · Im(r2 ) − λ · Im(α)), (39) (40) (41) when the motion constant takes a general distribution, the parameter r2 can be written in terms of α and ψ as [34] ∞ X 2 r2 = α + (|α| − 1) (−1)n c∗n einψ (α∗ )n−1 , (42) n=1 cn is the n-th Fourier expansion coefficient of the distribution of motion constants. For the simple case when the distribution is uniform on [0, 2π] , cn ≡ 0 for all the n , r2 ≡ α, and α decouples from ψ, this implies that the three-dimensional states phase has a two-dimensional invariant submanifold which is indeed the Ott-Antonsen manifold. However in the more typical case that α and ψ are interdependent the three-dimensional equation can exhibit non-general dynamical behavior in the marginal regime. In Fig. 4(a) we use Poincare section at ψ( mod 2π) = 0 to sort out the structure of state space, where the parameters are chosen in regime II of Fig. 1 (λ = 1.5 and σ = 2.0), and the motion constant has a distribution (1 + sin φ)/2π. In the Poincare section, quasiperiodic trajectories appear as closed curves or island chains, periodic trajectories appear as fixed points or periodp points of integer period, and chaotic trajectories fill the remaining regions of the unit disk, the picture of the phase portraits is reminiscent of Hamiltonian chaos and the appearance of the ”quasi-Hamiltonian” properties reflecting the time reversibility symmetry under the transformation t → −t, ψ → −ψ, x → −x, the system is invariant. When we choose the initial value in the chaos regime α(0) = 0.5 + i0.5, ψ(0) = 0.0, the three Lyapunov exponents are λ1 = −0.266 , λ2 = 0 , λ3 = 0.227 respectively. As a result the order parameter r2 and r1 is also chaotic which is initial values sensitive. Fig. 4(b) and (c) depicture the evolution of R1 (t) with two adjacent parameter value where δα(0) = 0.001, δψ(0) = 0 and it is clear that the bias of order parameter with neighboring parameters (illustration in the Fig. 4(b)) will be significant in the long time and both the characteristic curve and numerical simulation are consistent with each other well. Discussion To summarize, we investigated the collective dynamics of globally coupled identical phase oscillator when second harmonics dominate the coupling and solutions can be decomposed into symmetric and antisymmetric part independently. Theoretically, Ott-Antonsen method, linear stability analysis, and characteristic method have been carried out to obtain insights. Together with the numerical simulations, our study presented the following main results. First, we obtain the low-dimensional description of the symmetric dynamics and various regimes are predicted in the phase diagram, including the double-center, the single-center, the center-synchrony coexistence, and the synchrony regime. Second, all the steady states and the boundary curves has been obtained analytically both in terms of the Ott-Antonsen anatz and linear stability analysis. Third, we proved that in the marginal regime the stationary symmetry distribution is only neutrally stable where all the infinitely many eigenvalues are pure imaginary. Finally, the characteristic method 9 has been adopted to obtain the transient dynamics r1 which evolves strongly depends on the initial values. Furthermore, for the general case of initial condition the symmetry dynamics are governed by the Möbius transformation and numerical experiment suggests that three-dimensional invariant manifolds contain neutrally stable chaos. This work provided a complete framework to deal with the high-order coupling phase oscillators model, and the obtained results will enhance our understandings of the dynamical properties of more harmonics coupling phase oscillator system. [1] Pikovsky, A., Rosenblum, M. & Kurths, J. Synchronization: a Universal Concept in Nonlinear Sciences. pp. 279–296 (Cambridge University Press, Cambridge, England, 2001). [2] Buck, J. Synchronous rhythmic flashing of fireflies. II. The Quarterly Review of Biology 63, 265–289 (1988). [3] Georges, B., Grollier, J., Cros, V. & Fert, A. Impact of the electrical connection of spin transfer nanooscillators on their synchronization: an analytical study. Appl. Phys. Lett. 92, 232504 (2008). [4] Kiss, I. Z., Zhai, Y. & Hudson, J. L. Emerging Coherence in a Population of Chemical Oscillators. Science 296, 1676–1678 (2002). [5] Eckhardt, B., Ott, E., Strogatz, S. H., Abrams, D. M. & McRobie, A. Modeling walker synchronization on the Millennium Bridge. Phys. Rev. E 75, 021110 (2007). [6] Néda, Z., Ravasz, E., Vicsek, T., Brechet, Y. & Barabási, A. L. Physics of the rhythmic applause. Phys. Rev. E 61, 6987–6992 (2000). [7] Kuramoto, Y. Chemical Oscillations, Waves and Turbulence. pp. 75–76 (Springer, Berlin, 1984). [8] Strogatz, S. H. From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators. Physica D 143, 1–20 (2000). [9] Acebrón, J. A., Bonilla, L. L., Vicente, C. J. P., Ritort, F. & Spigler, R. The Kuramoto model: A simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77, 137–185 (2005). [10] Arenas, A., Diaz-Guilera, A., Kurths, J., Moreno, Y. & Zhou. C. Synchronization in complex networks. Phys. Rep. 469, 93–153 (2008). [11] Rodrigues, F. A., Peron, T. K. DM., Ji, P. & Kurths, J. The Kuramoto model in complex networks. Phys. Rep. 610, 1–98 (2016). [12] Skardal, P. S., Ott, E. & Restrepo, J. G. Cluster synchrony in systems of coupled phase oscillators with higher-order coupling. Phys. Rev. E 84, 036208 (2011). [13] Komarov, M. & Pikovsky, A. Multiplicity of Singular Synchronous States in the Kuramoto Model of Coupled Oscillators. Phys. Rev. Lett. 111, 204101 (2013). [14] KomarovM. & Pikovsky, A. The Kuramoto model of coupled oscillators with a bi-harmonic coupling function. Physica D 289, 18–31 (2014). [15] Li, K., Ma, S., Li, H. & Yang, J. Transition to synchronization in a Kuramoto model with the first- and second-orderinteraction terms. Phys. Rev. E 89, 032917 (2014). [16] Czolczynski, K., Perlikowski, P., Stefanski, A. & Kapitaniak, T. Synchronization of the self-excited pendula suspendedon the vertically displacing beam. Commun. Nonlinear Sci. Numer. Simul. 18, 386– 400 (2013). [17] Zhang, J., Yuan, Z. & Zhou, T. Synchronization and clustering of synthetic genetic networks: A role for cis-regulatory modules. Phys. Rev. E 79, 041903 (2009). [18] Kiss, I. Z., Zhai, Y. & Hudson, J. L. Predicting Mutual Entrainment of Oscillators with Experiment- 10 [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] Based Phase Models. Phys. Rev. Lett. 94, 248301 (2005). Goldobin, E., Koelle, D., Kleiner, R. & Mints, R. G. Josephson Junction with a Magnetic-Field Tunable Ground State. Phys. Rev. Lett. 107, 227001 (2011). Goldobin, E., Kleiner, R., Koelle, D. & Mints, R. G. Phase Retrapping in a Pointlike ϕ Josephson Junction: The Butterfly Effect. Phys. Rev. Lett. 111, 057004 (2013). Gómez-Gardeñes, J., Gómez, S., Arenas, A. & Moreno, Y. Explosive synchronization transitions in scale-free networks. Phys. Rev. Lett. 106, 128701 (2011). Zou, Y., Pereira, T., Small, M., Liu, Z. & Kurths, J. Basin of Attraction Determines Hysteresis in Explosive Synchronization. Phys. Rev. Lett. 112, 114102 (2014). Xu, C. Gao, J. Sun, Y. Huang, X & Zheng, Z. Explosive or Continuous: Incoherent state determines the route to synchronization. Sci. Rep. 5, 12039 (2015). Coutinho, B. C., Goltsev, A. V., Dorogovtsev, S. N. & Mendes, J. F. F. Kuramoto model with frequency-degree correlations on complex networks. Phys. Rev. E 87, 032106 (2013). Kazanovich, Y. & Borisyuk, R. Synchronization in Oscillator Systems with a Central Element and Phase Shifts. Progress of Theoretical Physics 110, 1047–1057 (2003). Burylko, O., Kazanovich, Y. & Borisyuk, R. Bifurcations in phase oscillator networks with a central element. Physica D 241, 1072–1089 (2011). Kazanovich, Y., Burylko, O. & Borisyuk, R. Competition for synchronization in a phase oscillator system. Physica D 261, 114–124 (2013) Vlasov, V., Zou, Y. & Pereira, T. Explosive synchronization is discontinuous. Phys. Rev. E 92, 012904 (2015). Vlasov, V., Pikovsky, A. & Macau, E. E. N. Star-type oscillatory networks with generic Kuramoto-type coupling: A model for ”Japanese drums synchrony”. Chaos 25, 123120 (2015). Ott, E. & Antonsen, T. M. Low dimensional behavior of large systems of globally coupled oscillators. Chaos 18, 037113 (2008). Ott, E. & Antonsen, T. M. Long time evolution of phase oscillator systems. Chaos 19, 023117 (2009). Watanabe, S. & Strogatz, S. H. Constants of motion for superconducting Josephson arrays. Physica D 74, 197–253 (1994). Golomb, D., Hansel, D., Shraiman, B. & Sompolinsky, H. Clustering in globally coupled phase oscillators. Phys. Rev. A 45, 3516–3530 (1991). Marvel, S. A., Mirollo, R. E. & Strogatz, S. H. Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action. Chaos 19, 043104 (2009). Marvel, S. A. & Strogatz, S. H. Invariant submanifold for series arrays of Josephson junctions. Chaos 19, 013132 (2009). Acknowledgements This work is partially supported by the NSFC grants Nos. 11075016, 11475022, 11135001, and the Scientific Research Funds of Huaqiao University. 11 Author contributions C.X., Y.T.S., J.G., T.Q, S.G.G and Z.G.Z. designed the research; C.X., Y.T.S. and T.Q performed numerical simulations and theoretical analysis; C.X., S.G.G and Z.G.Z. wrote the paper. All authors reviewed and approved the manuscript. Additional information Competing financial interests: The authors declare no competing financial interests. Correspondence and requests for materials should be (zgzheng@bnu.edu.cn), or S.G.G. (guanshuguang@hotmail.com). 12 addressed to Z.G.Z. FIG. 1: Phase diagram of the symmetric dynamics in the parameter space σ and λ . I is double center regime, II is single center regime, III is center-synchrony coexistence regime, IV is synchrony regime, respectively. √ The boundary curves are λ = σ − ω from I to II , λ = σ + ω from II to III or IV , λ = σ + ω 2 + σ 2 from III to IV, in the numerical simulations we set ω = 1. FIG. 2: (a) The initial phase θi vs the final phase θf when we choose the initial phase distribution ρ(θ, t0 ) = (1 + cos θ)/2π and the number of oscillator is N = 10000 in the simulation, when the initial phase is in the gray area these oscillators will bypassing the saddle points (the pink hollow circle in the horizontal axis) which leads to ε 6= 0. (b) The diagram of the Arnold diffusion mechanism in the same parameter with N = 3 in terms of phase space of θ1 and θ2 , for some particular trajectory such as the red line the oscillator will bypassing the saddle point in the invariant manifold. 13 FIG. 3: Transient dynamics of R1 and R2 . (a) The amplitude of R2 is small, R1 oscillates regularly with a period T , σ = 1.0, λ = 1.0, Re(r20 ) = 0, and Im(r20 ) = 0.427. (b) The amplitude of R2 is large, R1 oscillates irregularly, σ = 1.0, λ = 1.4, Re(r20 ) = 0, and Im(r20 ) = −0.5. In the numerical simulation we choose the initial phase distribution ρ(θ, t0 ) = ρs (2θ, t0 )(1 + 0.5 cos(θ)) and N = 10000, the line is calculated by the characteristic theory and the circle is the numerical simulation. FIG. 4: (a) Poincare section of Eq. (39) ∼ Eq. (41) at ψ( mod 2π) = 0, quasiperiodic trajectories appear as closed curves or island chains, periodic trajectories appear as fixed points or period-p points of integer period, and chaotic trajectories fill the remaining region of the unit disk, where ω = 1.0, λ = 1.5, σ = 2.0. (b) the R1 vs t in the chaos regime with α(0) = 0.5 + i 0.5 and ψ(0) = 0.0. (c) the R1 vs t with adjacent parameter value δα(0) = 0.001, δψ(0) = 0. The illustration is the difference of R1 vs t where we can see that the bias of the order parameter with neighboring parameters will be significant in the long time. In the simulation we choose N = 100000, the line is determined in terms of the the characteristic theory and the circle and triangle are calculated by the numerical simulation. 14