Order parameter analysis of synchronization transitions on star networks

Abstract

The collective behaviors of populations of coupled oscillators have attracted significant attention in recent years. In this paper, an order parameter approach is proposed to study the low-dimensional dynamical mechanism of collective synchronizations, by adopting the star-topology of coupled oscillators as a prototype system. The order parameter equation of star-linked phase oscillators can be obtained in terms of the Watanabe–Strogatz transformation, Ott–Antonsen ansatz, and the ensemble order parameter approach. Different solutions of the order parameter equation correspond to the diverse collective states, and different bifurcations reveal various transitions among these collective states. The properties of various transitions in the star-network model are revealed by using tools of nonlinear dynamics such as time reversibility analysis and linear stability analysis.

Graphical abstract

Keywords

Kuramoto model / synchronization / order parameter / Ott–Antonsen ansatz / star network

1 Introduction

Understanding the intrinsic microscopic mechanism embedded in collective macroscopic behaviors of populations of coupled units on heterogeneous networks has become the focus in a variety of fields, such as the biological neurons circadian rhythm, chemical reacting cells, and even society systems [1-8]. Numerous different emerging macroscopic states/phases have been revealed, and various non-equilibrium transitions among these states have been observed and studied on heterogeneous networks [9-19].

The transitions among the different collective states on heterogeneous networks exhibit the typical feature of multistability, i.e., these states may coexist for a group of given parameters and depend on the choice of initial conditions. This interesting behavior is closely related to the first-order phase transition, and multistability in the discontinuous transitions indicates the competition of miscellaneous attractors and their corresponding basins of attraction in the phase space. For a network of coupled oscillators, the microscopic description of the dynamics of oscillators should be made in a high-dimensional phase space, which is very difficult to deal with. The key point in understanding macroscopic transitions is the projection of the dynamics from this high-dimensional space to a much lower-dimensional subspace. This can be executed by introducing appropriate order parameters and building their dynamical equations. Ott and Antonsen [20] proposed an ansatz to project the infinite-dimensional dynamics to a low-dimensional manifold called the Ott–Antonsen (OA) manifold, which has been successfully applied to systems composed of large numbers of oscillators. However, strictly speaking, the OA manifold analysis cannot be applied to finite-oscillator systems. Watanabe and Strogatz introduced the Möbius transformation for finite-size systems with specific symmetries to obtain exact three-dimensional dynamics [21, 22]; however, this scheme cannot be extended to general finite systems. The mechanism of the validity of the OA approach was studied recently, and the ensemble order parameter approach is proposed, which extends the OA approach to more general cases such as finite-number of oscillators and more general coupling forms [23].

Abrupt or explosive transition from the incoherent state to synchronization may occur on networks if the frequencies of the oscillators on the nodes are positively correlated to the node’s degrees [13], which has been observed numerically on scale-free networks and experimentally in electronic circuits [24, 25]. The first-order transition can be changed and more ways of transitions can be observed by adjusting the phase shift among the oscillators [26]. Numerous efforts have been made to understand the mechanism of explosive synchronization from different viewpoints such as the topological structures of networks, the coupling functions among nodes, and so on [17, 25, 27-32].

An analytical understanding of the transitions among the various synchrony states in heterogeneous networks is very valuable. In addition to being the simplest, the star topology is the key topology that can describe the heterogeneity property of complex networks such as scale-free networks [33-36]. In this paper, the collective states and the abundant transitions among these states on a star network are studied by considering the effect of phase shift among coupled oscillators [23, 26, 37-39]. The dynamics of star networks of oscillators are analytically studied by building the equations of motion of the order parameter for networks with a finite size, which accomplishes a great reduction from the microscopic high-dimensional phase dynamics of coupled oscillators to a macroscopic low-dimensional dynamics. Based on the order parameter dynamics, we further reveal numerous transitions among different collective states in this model by using tools of nonlinear dynamics such as time reversibility analysis [15] and linear stability analysis. We found three typical processes of the transitions to the synchronous state, i.e., the transitions from the neutral state, the in-phase state, or the splay state to the synchronous state, and a continuous process of desynchronization and a group of hybrid phase transitions that are discontinuous with no hysteresis.

2 The Ott–Antonsen ansatz and the Watanabe–Strogatz approach

We first illustrate the Ott–Antonsen ansatz [20] briefly by analyzing the following class of identical oscillators governed by the equations of motion

\begin{array}{}\text{(1)}& {\dot{\phi }}_j=f{\mathrm{e}}^{\mathrm{i}{\phi }_j}+g+\overline{f}{\mathrm{e}}^{-\mathrm{i}{\phi }_j},j=1,\mathrm{\cdots },N,\end{array}

where f is a smooth, complex-valued 2π-periodic function of the phases {\phi }_1,\mathrm{\cdots },{\phi }_N and the overbar denotes complex conjugate, g is a real valued function since {\dot{\phi }}_j is real. By introducing the distribution of phases of oscillators in the limit N⟶\mathrm{\infty }, the evolution of the system (1) is given by the continuity equation

\begin{array}{}\text{(2)}& \frac{\partial \rho }{\partial t}+\frac{\partial (\rho \nu )}{\partial \varphi }=0,\end{array}

where \rho (\varphi ,t) is the phase distribution function, and \rho (\varphi ,t)\mathrm{d}\varphi gives the fraction of phases that lie between ϕ and \varphi +\mathrm{d}\varphi at time t. The velocity field is the Eulerian version of Eq. (1),

\begin{array}{}\text{(3)}& \nu (\varphi ,t)=f{\mathrm{e}}^{\mathrm{i}\phi }+g+\overline{f}{\mathrm{e}}^{-\mathrm{i}\phi }.\end{array}

Suppose ρ is of the form

\begin{array}{}\text{(4)}& \rho (\varphi ,t)=\frac{1}{2\pi }\left\{1+\sum _{n=1}^{\mathrm{\infty }}[\overline{z}{(t)}^n{\mathrm{e}}^{\mathrm{i}n\varphi }+z{(t)}^n{\mathrm{e}}^{-\mathrm{i}n\varphi }]\right\}\end{array}

for some unknown function z that is independent of ϕ. Note that Eq. (4) is an algebraic rearrangement of the usual form for the Poisson kernel

\begin{array}{}\text{(5)}& \rho (\varphi )=\frac{1}{2\pi }\frac{1-r^2}{1-2r\cos (\varphi -\Phi )+r^2},\end{array}

where the complex number z can be expressed in the complex plane as

\begin{array}{}\text{(6)}& z=r{\mathrm{e}}^{\mathrm{i}\Phi }.\end{array}

The ansatz (4) defines a submanifold in the infinite-dimensional space of the density function ρ. This Poisson submanifold is two-dimensional and is parameterized by the complex number z, or equivalently, by the polar coordinates r and \Phi. An intriguing point discovered by Ott et al. [20] is the invariance of the Poisson submanifold, i.e., if the initial phase density is a Poisson kernel, it always remains a Poisson kernel. This can be verified by substituting the velocity field (3) and the ansatz (4) into the continuity equation (2). It can be found that the amplitude equations for each harmonic {\mathrm{e}}^{\mathrm{i}n\varphi } are simultaneously satisfied if and only if  z(t) evolves according to

\begin{array}{}\text{(7)}& \dot{z}=\mathrm{i}(f{z}^2+gz+\overline{f}).\end{array}

This equation can be recast in a more physically meaningful form in terms of the complex order parameter defined as the centroid of the phases ϕ regarded as points {\mathrm{e}}^{\mathrm{i}\varphi } on the unit circle:

\begin{array}{}\text{(8)}& ⟨{\mathrm{e}}^{\mathrm{i}\varphi }⟩={\int }_0^{2\pi }{\mathrm{e}}^{\mathrm{i}\varphi }\rho (\varphi ,t)d\varphi .\end{array}

By substituting Eq. (4) into Eq. (8), it is seen that

\begin{array}{}\text{(9)}& z=⟨{\mathrm{e}}^{\mathrm{i}\varphi }⟩=r{\mathrm{e}}^{\mathrm{i}\Phi }\end{array}

for all states on the Poisson submanifold. Thus, it is clear that z represents the order parameter of the system, r is its modulus, and \Phi is its mean phase. However, from Ott–Antonsen ansatz, it cannot be determined whether the governing equation Eq. (7) can be used for a system with finite size.

For a finite number of oscillators N, the original microscopic dynamical state can be reduced to a macroscopic collective state by the Watanabe–Strogatz approach [21, 22], and the governing equations of the system can be generated by the Möbius group action [40, 41]. The class of identical oscillators is still governed by the equations of motion Eq. (1); therefore, the oscillators’ phases {\phi }_j(t) evolve according to the action of the Möbius group on the complex unit cycle

\begin{array}{}\text{(10)}& {\mathrm{e}}^{\mathrm{i}{\phi }_j(t)}=M_t({\mathrm{e}}^{\mathrm{i}{\theta }_j})\end{array}

for j=1,\mathrm{\dots },N, where M_t is a one-parameter family of the Möbius transformations and {\theta }_j is a constant angle. By parameterizing the one-parameter family of Möbius transformations as

\begin{array}{}\text{(11)}& M_t(w)=\frac{{\mathrm{e}}^{\mathrm{i}\psi }w+\eta }{1+\overline{\eta }{\mathrm{e}}^{\mathrm{i}\psi }w},\end{array}

where and \psi (t)\in R, and assuming that

\begin{array}{}\text{(12)}& w_j={\mathrm{e}}^{\mathrm{i}{\theta }_j},\end{array}

we get

\begin{array}{}\text{(13a)}& \dot{\eta }=\mathrm{i}(f{\eta }^2+g\eta +\overline{f}),\end{array}

\begin{array}{}\text{(13b)}& \dot{\psi }=f\eta +g+\overline{f}\overline{\eta }.\end{array}

With these new variables, the order parameter can be rewritten as

\begin{array}{}\text{(14)}& z(t)=\frac{1}{N}\sum _{j=1}^N\frac{{\mathrm{e}}^{\mathrm{i}\psi }{\mathrm{e}}^{\mathrm{i}{\theta }_j}+\eta (t)}{1+\overline{\eta }(t){\mathrm{e}}^{\mathrm{i}\psi }{\mathrm{e}}^{\mathrm{i}{\theta }_j}},\end{array}

Eqs. (13) and (14) can describe the system with arbitrary initial conditions as \eta (0),\psi (0), and N constants {\theta }_j,1\le j\le N. The order parameter (14) could be simplified further by choosing the constants

\begin{array}{}\text{(15)}& {\theta }_j=2\pi \frac{j-1}{N},1\le j\le N,\end{array}

with which, the order parameter (14) reads

\begin{array}{}\text{(16)}& z(t)=\eta (t)(1+I),\end{array}

where I=(1-{|\eta (t)|}^{-2})/[1\pm {({\mathrm{e}}^{\mathrm{i}\psi }\overline{\eta }(t))}^{-N}], “-” for the case with even N and “+” for the case with odd N. For large N, I≪1, the order parameter could be approximated as

\begin{array}{}\text{(17)}& z(t)\approx \eta (t),N\gg 1.\end{array}

Based on the above analysis, the dynamics of the system with finite size can be described by the same equation as the governing equation (7), which is obtained from the Ott–Antonsen ansatz for the system with infinite size. Then, Eq. (7) can be used to explore the low-dimensional collective behaviors of the system with finite size.

3 The Sakaguchi–Kuramoto model on star networks: The order parameter equation

We started with a star network of coupled phase oscillators with nonzero phase shift as our working model. In the star network with one hub and K leaves, the degree of the leaves is k_i=1 (i=1,\mathrm{\dots },K) and the degree of the hub is k_h=K. By assuming that the natural frequencies of the oscillators are proportional to their degrees, the equations of motion for the hub and leaf nodes can be written as

\begin{array}{}\text{(18)}& \begin{array}{l}{\dot{\theta }}_h={\omega }_h+\lambda {\displaystyle \sum _{j=1}^K\sin }({\theta }_j-{\theta }_h-\alpha ),\\ {\dot{\theta }}_j=\omega +\lambda \sin ({\theta }_h-{\theta }_j-\alpha ),1\le j\le K,\end{array}\end{array}

where {\theta }_h,{\theta }_j, and {\omega }_h,\omega are instantaneous phases and natural frequencies of the hub and leaf nodes respectively. λ is the coupling strength. K is the number of leaf nodes connected with this hub and α is the phase shift. The effect of phase shift among the coupled oscillators has been extensively investigated in recent years, while this has been seldom discussed in the case of star networks [14, 42]. Abundant collective dynamics appear in the global coupled model in the presence of a finite phase shift [14], where synchrony can decay, or incoherence can regain its stability with increased coupling and multistability between partially synchronized states, and/or the incoherent state can appear in the globally coupled network.

The coupled phase oscillator system with α=0 of the star network was originally designed to study the characteristics of explosive synchronization [13]; however the process of synchronization may be affected by the introduction of phase shift [14]. By introducing the phase differences between the hub and leaves {\phi }_j={\theta }_h-{\theta }_j, the phase dynamics of star networks can be transformed to the following phase difference dynamics on an all-connected network,

\begin{array}{}\text{(19)}& {\dot{\phi }}_i=\mathrm{\Delta }\omega -\lambda \sum _{j=1}^K\sin ({\phi }_j+\alpha )-\lambda \sin ({\phi }_i-\alpha ),\end{array}

where 1\le i\le K. We further define the order parameter of the all-connected network to describe the degree of synchronization as

\begin{array}{}\text{(20)}& z(t)\equiv r(t){\mathrm{e}}^{\mathrm{i}\Phi (t)}=\frac{1}{K}\sum _{j=1}^K{\mathrm{e}}^{\mathrm{i}{\phi }_j}.\end{array}

The star network becomes globally synchronous if the modulus r(t)=1 and the mean phase \Phi (t)=const. If the modulus r(t)=1 while the mean phase \Phi (t) is periodic, which corresponds to the state with {\phi }_j(t)=\phi (t), all the leaf nodes are synchronous to each other while they are asynchronous to the hub oscillator.

Eq. (19) can be rewritten as

\begin{array}{}\text{(21)}& {\dot{\phi }}_j=f{\mathrm{e}}^{\mathrm{i}{\phi }_j}+g+\overline{f}{\mathrm{e}}^{-\mathrm{i}{\phi }_j},j=1,\mathrm{\cdots },K,\end{array}

where \mathrm{i} denotes the imaginary unit, f=\mathrm{i}{\displaystyle \frac{\lambda }{2}}{\mathrm{e}}^{-\mathrm{i}\alpha }, and g=\mathrm{\Delta }\omega -\lambda Kr\sin (\Phi +\alpha ).

For finite K, owing to the high topological symmetry of the star network, the collective behaviors of the system can be analyzed by deriving the low-dimensional dynamical equations in terms of both the ensemble order parameter approach [23] and the Watanabe–Strogatz transformation [21, 22]. If the initial phases of the oscillators are chosen as (15), we have the dynamical equation for the order parameter z(t) as

\begin{array}{}\text{(22)}& \dot{z}=-\frac{\lambda }{2}{\mathrm{e}}^{-\mathrm{i}\alpha }{z}^2+\mathrm{i}[\mathrm{\Delta }\omega -\lambda Kr\sin (\Phi +\alpha )]z+\frac{\lambda }{2}{\mathrm{e}}^{\mathrm{i}\alpha },\end{array}

which is the OA result (18) in terms of the order parameter. Different solutions of Eq. (22) build correspondences with diverse collective states of the coupled oscillator system.

4 Collective dynamics of stationary states

In this section, we discuss the collective dynamics of the star network in terms of Eq. (22). By setting z=x+\mathrm{i}y, the order parameter dynamics in the x−y plane can be described as

\begin{array}{}\text{(23)}& \begin{array}{l}\begin{array}{ll}\dot{x}=\hfill & \lambda \left(\frac{1}{2}+K\right)\cos \alpha y^2-\frac{\lambda }{2}\cos \alpha x^2\hfill \\ \hfill & +\lambda (K-1)\sin \alpha xy-\text{Δ}\omega y+\frac{\lambda }{2}\cos \alpha ,\hfill \end{array}\hfill \\ \begin{array}{ll}\dot{y}=\hfill & \lambda \left(\frac{1}{2}-K\right)\sin \alpha x^2-\frac{\lambda }{2}\sin \alpha y^2\hfill \\ \hfill & -\lambda (K+1)\cos \alpha xy+\text{Δ}\omega x+\frac{\lambda }{2}\sin \alpha .\hfill \end{array}\hfill \end{array}\end{array}

The steady-state solutions are determined by setting \dot{x}=0 and \dot{y}=0, which results in four fixed points noted by (x_i,y_i) with

\begin{array}{}\text{(24)}& \begin{array}{l}{x}_{1,2}=\frac{-\sin \alpha \text{Δ}\omega \pm A\sin \alpha }{\lambda [2K\cos (2\alpha )+1]},\\ y_{1,2}=-\frac{-\cos \alpha \text{Δ}\omega \pm A\cos \alpha }{\lambda [2K\cos (2\alpha )+1]},\\ x_{3,4}=\frac{\sin \alpha }{\lambda }+\frac{\frac{\sin (2\alpha )}{2}B\pm K\left[\frac{\sin (2\alpha )}{2}B-\sin (2{\alpha }^2)\right]}{\lambda \sin \alpha [K^2+2\cos (2\alpha )K+1]},\\ y_{3,4}=\frac{-\text{Δ}\omega (-\cos \alpha \pm \sin \alpha B-K\cos \alpha )}{\lambda [K^2+2\cos (2\alpha )K+1]},\end{array}\end{array}

where “+” represents the fixed points (x_{1,3},y_{1,3}), “-” represents the fixed points (x_{2,4},y_{2,4}), and

\begin{array}{}\text{(25)}& \begin{array}{l}A=\sqrt{-2K{\lambda }^2\cos (2\alpha )-{\lambda }^2+\text{Δ}{\omega }^2},\\ B=\sqrt{{\lambda }^2+K^2{\lambda }^2+2K{\lambda }^2\cos (2\alpha )-\text{Δ}{\omega }^2}.\end{array}\end{array}

The existence condition for the fixed points are determined by Eq. (25), where -2K{\lambda }^2\cos (2\alpha )-{\lambda }^2+\mathrm{\Delta }{\omega }^2\ge 0, and {\lambda }^2+K^2{\lambda }^2+2K{\lambda }^2\cos (2\alpha )-\mathrm{\Delta }{\omega }^2\ge 0. For the fixed points (x_{1,2},y_{1,2}), the existence condition can be given as

\begin{array}{}\text{(26)}& \lambda \le {\lambda }_1=\frac{\mathrm{\Delta }\omega }{\sqrt{2K\cos (2\alpha )+1}},\end{array}

and for the fixed points (x_{3,4},y_{3,4}), the existence condition is

\begin{array}{}\text{(27)}& \lambda \ge {\lambda }_2=\frac{\mathrm{\Delta }\omega }{\sqrt{K^2+2K\cos (2\alpha )+1}}.\end{array}

For the fixed points (x_{3,4},y_{3,4}), there is an additional natural restriction relation

\begin{array}{}\text{(28)}& x^2+y^2=1,\end{array}

while for the fixed points (x_{1,2},y_{1,2}), x^2+y^2 may be greater or lower than 1. The definition of z implies that only those fixed points satisfying x^2+y^2\le 1 are reasonable.

Linear stability analysis can be applied to the fixed points (x_i,y_i),i=1,2,3,4 by computing the eigenvalues of the 2×2 Jacobian matrix J of the fixed points with elements

\begin{array}{}\text{(29)}& \begin{array}{l}{J}_{11}=-\lambda \cos \alpha x_i+\lambda (K-1)\sin \alpha y_i,\\ J_{12}=\lambda (1+2K)\cos \alpha y_i+\lambda (K-1)\sin \alpha x_i-\text{Δ}\omega ,\\ J_{21}=\lambda (1-2K)\sin \alpha x_i-\lambda (K+1)\cos \alpha y_i+\text{Δ}\omega ,\\ J_{22}=-\lambda \sin \alpha y_i-\lambda (K+1)\cos \alpha x_i.\end{array}\end{array}

The eigenvalues of the Jacobian matrix are

\begin{array}{}\text{(30)}& {\beta }_{1,2}=\frac{J_{11}+J_{22}\pm \sqrt{{(J_{11}+J_{22})}^2-4(J_{11}{J}_{22}-J_{12}{J}_{21})}}{2}.\end{array}

The stability conditions of the four fixed points are summarized in Table 1, where the parameters in the table are {\widehat{\lambda }}_c^f=\mathrm{\Delta }\omega /\sqrt{2K\cos \left(2\alpha \right)+1}, {\lambda }_{sc}^{+}=-\mathrm{\Delta }\omega /\left[K\cos \left(2\alpha \right)+1\right], {\lambda }_{sc}^{-}=\mathrm{\Delta }\omega /\left[K\cos \left(2\alpha \right)+1\right], {\alpha }_0^{-}=-\arccos (-1/K)/2, {\alpha }_0^{+}=\arccos (-1/K)/2.

The fixed points of the order parameter equation are related to the collective states of the coupled phase oscillators. The fixed points (x_{3,4},y_{3,4}) with |z|=1 correspond to the synchronous state (SS) of the system, where all the phase differences between the hub and leaf nodes are constant and have the same value, as

\begin{array}{}\text{(31)}& {\phi }_j(t)=const,1\le j\le K,\end{array}

implying the global synchronization of the hub and leaves in a star network. The above stability analysis indicates that the fixed point (x_4,y_4) corresponds to the unstable synchronous state and the fixed point (x_3,y_3) corresponds to the stable synchronous state, and their stability can be studied easily.

The fixed points (x_{1,2},y_{1,2}) with the modulus |z|=\sqrt{x^2+y^2}>1 is unphysical, because the order parameter z of the coupled oscillators is bounded by |z|\le 1. If , the related collective state is called the splay state (SPS) [11, 43], and the phase differences between the hub and leaf nodes satisfy a function relation as

\begin{array}{}\text{(32)}& {\phi }_j(t)=\phi \left(t+\frac{jT}{K}\right),1\le j\le K\end{array}

with T the period of \phi (t), as shown in Fig. 1(a). This kind of state physically represents the collective state where all the leaf oscillators in the star network are in synchronous motion with a constant time shift.

5 Collective dynamics of time-dependent states

Long-term solutions of the order parameter equation contain not only the states given by fixed points, but also the time-dependent states corresponding to the periodic solutions. There are two periodic regimes: the regime and , and the critical line with \alpha =0,\pm \pi /2:

5.1 The in-phase state

When and , as shown in Fig. 1(b), the largest Lyapunov exponent is zero and the other exponents are negative, implying a stable limit-cycle solution. This solution can be conveniently found by transforming the Eq. (22) to polar coordinates as z=r{\mathrm{e}}^{\mathrm{i}\Phi }:

\begin{array}{}\text{(33)}& \begin{array}{l}\dot{r}=-\frac{\lambda }{2}(r^2-1)\cos (\Phi +\alpha ),\\ \dot{\Phi }=-\frac{\lambda }{2}\left(r+\frac{1}{r}\right)\sin (\Phi -\alpha )+\text{Δ}\omega -\lambda Kr\sin (\Phi +\alpha ).\end{array}\end{array}

There is a limit cycle solution with r=1 and periodic phase \Phi (t), which is called the in-phase state (IPS), where all the phase differences between leaves and the hub are the same and time dependent, i.e.,

\begin{array}{}\text{(34)}& {\phi }_j(t)=\phi (t),1\le j\le K.\end{array}

This state corresponds to the phases of all the leaf nodes evolve synchronously, while they are not synchronous to the hub, as shown in Fig. 1(c). In this case, the model Eq. (3) is reduced to the case with K=1 and the stability of the state can be obtained by Floquet theory for limit cycle as it is stable for and unstable for .

5.2 The neutral state

The dynamics at \alpha =0,\pm \pi /2 when are special cases and correspond to the critical dynamical states of the system, where the related fixed point is found to be neutrally stable where Re({\beta }_{1,2})=0 and Im({\beta }_{1,2})\ne 0 in Eq. (30). In this case, there is a large class of critical states with the order parameter r determined by the initial values of (x,y), and the long-term behavior of z depend crucially on the initial phases, as shown in Fig. 1(d). We call this state the neutral state(NS) [15], and the corresponding fixed point is neutrally stable, where all Lyapunov exponents of the fixed point are zero when , as shown in Fig. 2(a). The Kuramoto system is dissipative [44]; therefore, the existence of a large class of neutral states are counterintuitive, which implies that the phase space of this state contains an integrable Hamiltonian system family of periodic orbits, shown in Fig. 2(b).

To understand the mechanism of these neutral states, we resort to the analysis of the order parameter equations (23). For the case of α=0, Eq. (23) can be simplified to

\begin{array}{}\text{(35)}& \begin{array}{l}\dot{x}=\lambda \left(K+\frac{1}{2}\right)y^2-\frac{\lambda }{2}{x}^2-\text{Δ}\omega y+\frac{\lambda }{2},\\ \dot{y}=-\lambda (K+1)xy+\text{Δ}\omega x.\end{array}\end{array}

The fixed points are determined by setting \dot{x}=0 and \dot{y}=0. When only the fixed point (x_1,y_1) exists inside the unit cycle in the plane as shown in Fig. 2(b). The fixed point is neutrally stable, and all the Lyapunov exponents of the neutral state with α=0 are zero. If we define a time reversal transformation as R:(t,x,y)\mapsto (-t,-x,y), the dynamical equations (35) remain invariant. Hence, they are called the time-reversible dynamical system or the quasi-Hamiltonian system [15]. This symmetry endows the system with many interesting properties.

The time reversal transformation R can be resolved into R=TW with T:t\mapsto -t and W:(x,y)\mapsto (-x,y). Hence, the invariant set for W is the y axis with x=0,y>0. For any trajectory crossing this invariant set, according to the time reversal symmetry, the forward trajectory and the backward trajectory are symmetric. If the forward trajectory evolves into an attractor, the backward trajectory will evolve into the symmetric repeller of the system. Then, the attractor and the repeller of the system will emerge in pairs. When the trajectory crosses the invariant set more than once, the forward and backward trajectories will coincide with each other, forming the periodic solution for the system, which is called the reversible trajectory [15]. For any reversible trajectory, the Lyapunov exponents have the sign-symmetry form and the volume of phase space in its vicinity are conserved in the form of averages as shown in the numerical simulations.

The order parameter plane of our system is bounded by the unit circle with an invariant set as x=0,y>0. Therefore, the attractor and repeller emerge at the same time, implying that if only one fixed point exists in the plane, it is neither the attractor nor the repeller, i.e.; the only fixed point is neutrally stable. If only one neutrally stable fixed point exists in the plane, and the trajectories are vagrant and must cross the invariant set more than once, then those trajectories are closed and periodic. This is what happens when α=0 with the region as shown in Fig. 2(b). For α=0 and \lambda >{\lambda }_{ec}, there is a coexisting region for the synchronous and neutral states as the critical cases in the coexistence region for the neutral and synchronous states.

All the above possible collective states are summarized in the parameter space (α,λ) as a phase diagram in Fig. 3 having the boundaries that were obtained analytically from both the existence and stability conditions. In Fig. 3, four regions of the phase shift α can be identified. For the first region , the splay state exists and is stable for any λ. With the increase of coupling strength λ, the unstable synchronous state exists above the threshold \lambda >{\lambda }_{ec}. For the second region , the splay state exists and is stable within , and the synchronous state exists with \lambda >{\lambda }_{ec}; however, it is unstable unless \lambda >{\lambda }_{sc}^{-}. Obviously there exists a co-existing region for the splay state and the synchronous state in the region , within the coupling interval . In the third region, where , the splay state is always unstable, the stable synchronous state emerges as the coupling strength \lambda >{\lambda }_{ec}. For the fourth region , the splay state always exists, however, it is stable only when \lambda >{\lambda }_{sc}^{+}, and the stable synchronous state only exists in the region . The neutral state exists as a particular case for the phase shift, \alpha =0,\pm \pi /2, and the in-phase state is always stable in the region , within the coupling range . The variety of states in the phase diagram leads to various transitions among them.

6 Scenarios of synchronization transitions

The phase diagram shown in Fig. 3 presents a great variety of transitions among the different collective dynamical states. It is observed that some of the states coexist with each other at the same parameter. These coexisting states may lead to abrupt transitions among them and hysteresis behaviors, while the others lead to continuous transitions.

6.1 Synchronization transition from the neutral state

We first investigated the synchronization process from the neutral state to the synchronous state for α=0. The synchronization process when α=0 is discontinuous, which is known as the explosive synchronization, has attracted significant attention recently [13]. It was shown that by changing the coupling strength λ, this transition is abrupt, and there is a hysteretic behavior at the onset of synchronization. {\lambda }_c^b and {\lambda }_c^f are the backward and forward critical coupling strengths respectively, where {\lambda }_c^b={\lambda }_2 and {\lambda }_c^f depends on the initial states as shown in Fig. 4(a). The upper limit of {\lambda }_c^f is denoted by {\widehat{\lambda }}_c^f. As \lambda >{\widehat{\lambda }}_c^f, the synchronization state is globally attractive. It is difficult to understand this process on the basis of the self consistent method, especially for the hysteresis behavior and coexisting region.

The critical coupling corresponds to the upper limit of {\lambda }_c^f, which can be determined as

\begin{array}{}\text{(36)}& {\widehat{\lambda }}_c^f=\frac{\mathrm{\Delta }\omega }{\sqrt{2K+1}}.\end{array}

From the analytical curve and the simulation results given in Fig. 4(b), it is clear that the results conform with the curve.

In the bistable regime, as shown in Fig. 4(c), the nullclines \dot{x}=0 (the red lines) and \dot{y}=0 (the green lines) have four intersections labeled by A-D, where A is an attractor, C is a repeller, and B and D are neurally stable. Any orbit crossing the nullcline A-B-C will eventually fall to A, and the others will hold the property as periodic orbits. It is clear that the stable fixed point A corresponds to the synchronous state. Moreover, the basin for the neutral state can be calculated using an approximation by the circle, which has its center at point D and the radius equal to the length of the line B-D. As λ increases, points D and B move closer to each other and eventually collide at the critical coupling, as shown in Fig. 4(d), and the synchronous state becomes globally attractive.

6.2 Synchronization transition from the splay state

The synchronization process from the splay state to the synchronous state for is found to be discontinuous. Numerical computations reveal that this kind of transition is abrupt with hysteresis at the onset of synchronization as shown in Fig. 5(a). The abrupt transition implies that there are two critical coupling strengths {\lambda }_c^b and {\lambda }_c^f, where {\lambda }_c^b={\lambda }_{sc}^{-} and {\lambda }_c^f depend on the basin of attraction. The upper limit of {\lambda }_c^f can be determined by analyzing the inverse saddle-node bifurcation as

\begin{array}{}\text{(37)}& \widehat{{\lambda }_c^f}=\frac{\mathrm{\Delta }\omega }{\sqrt{2K\cos (2\alpha )+1}}.\end{array}

As shown in Fig. 5(b), the simulation results are consistent with the analytical curve.

The dynamical manifestations of the discontinuous transition from the splay state to the synchronous state are shown in Figs. 5(c) and (d). Figure 5(c) exhibits the coexistence of the splay state and the synchronous state as the stable fixed points D and A respectively. The basins of attraction of the splay state and the synchronous state are separated by the saddle point B. When the coupling λ increases, as shown in Fig. 5(d), the saddle point B and the attractor D collide and disappear via an inverse saddle-node bifurcation, and this discontinuous transition makes the fixed point A correspond to the synchronous state a global attractor.

6.3 Synchronization transition from the in-phase state

The route of synchronization from the in-phase state to the synchronous state for \alpha >0 is shown in Fig. 6(a). The critical coupling strength of this continuous transition {\lambda }_{ec} is determined by Eq. (27). It can be found from Fig. 6(b) that the simulation results agree well with the analytical curve.

The dynamical manifestations of the transition from the in-phase state to the synchronous state are shown in Figs. 6(c) and (d). As shown in Fig. 6(c), the in-phase state is a limit cycle in the order parameter plane. As λ increases, the stable fixed point A corresponding to the synchronous state emerges on the limit cycle. The transition from the in-phase state to the synchronous state takes place continuously through a saddle-node bifurcation, as shown in Fig. 6(d).

6.4 Scenario of desynchronization

In the region of the phase diagram 3, the synchronous state is unstable when \lambda >{\lambda }_{sc}^{+} and the stable splay state emerges, which is contrary to our conventional belief that the system will always be synchronous if the coupling strength is large enough. The transition is called desynchronization, and it is a continuous transition as shown in Fig. 7(a). The order parameter r decreases rapidly at the threshold and effective frequencies of hub and leaf nodes are divided at the same coupling λ. From the phase diagram, it can be seen that the route of the de-synchronization is from the synchronous state to the splay state. The threshold of de-synchronization is {\lambda }_{sc}^{+}=-\mathrm{\Delta }\omega /\left[K\cos \left(2\alpha \right)+1\right] as shown in Fig. 7(b); this is consistent with the simulation results.

The dynamical manifestations of the continuous transition from the synchronous state to the splay state are shown in Figs. 7(c) and (d). From Fig. 7(c), which illustrates the stable synchronous state of the system, it is seen that when , all the orbits in the phase of the order parameter will eventually evolve to the fixed point A. As λ increases and becomes larger than the critical coupling {\lambda }_{sc}^{+}, the two nullclines will intersect at the four fixed points as shown in Fig. 7(d). The point A will lose its stability and a new stable fixed point B, which corresponds to the splay state, appears. This process of transition from the synchronous state to the splay state is continuously completed by this bifurcation.

7 Conclusion

To summarize, in this paper we studied the dynamics of coupled oscillators on a star network with the Sakaguchi–Kuramoto model using the dynamical order parameter equation that can be obtained in terms of different approaches, e.g., the ensemble order parameter approach and the Watanabe–Strogatz approach. The order parameter equation obtained for the star network can be approximately described from the Ott–Antonsen ansatz, which originates from the high symmetry of the topology. The essential dynamical mechanisms of different scenarios of synchronization can be understood analytically by the reduction from a high-dimensional phase space to a much lower-dimensional order parameter space without any additional approximation. Different solutions of the order parameter equation correspond to the various collective states of the coupled oscillators, and the different bifurcations reveal various transitions among those collective states. These transition processes were revealed in the plane of order parameter and their critical coupling strengths were obtained analytically, and verified by the simulation results.

References

[1] Y. Kuramoto, Chemical Oscillations, Waves and Turbulence, Springer Science and Business Media, 2012

[2] J. A. Acebrón, L. L. Bonilla, C. J. Pérez Vicente, F. Ritort, and R. Spigler, The Kuramoto model: A simple paradigm for synchronization phenomena, Rev. Mod. Phys. 77(1), 137 (2005)

[3] S. H. Strogatz, From Kuramoto to Crawford: Exploring the onset of synchronization in populations of coupled oscillators, Physica D 143(1–4), 1 (2000)

[4] A. Pikovsky, M. Rosenblum, and J. Kurths, Synchronization: A Universal Concept in Nonlinear Sciences, Cambridge: Cambridge University Press, 2001

[5] S. N. Dorogovtsev, A. V. Goltsev, and J. F. F. Mendes, Critical phenomena in complex networks, Rev. Mod. Phys. 80(4), 1275 (2008)

[6] A. Arenas, A. Diaz-Guilera, J. Kurths, Y. Moreno, and C. Zhou, Synchronization in complex networks, Phys. Rep. 469(3), 93 (2008)

[7] Z. Zheng, G. Hu, and B. Hu, Phase slips and phase synchronization of coupled oscillators, Phys. Rev. Lett. 81(24), 5318 (1998)

[8] Z. Zheng, G. Hu, and B. Hu, Collective phase slips and phase synchronizations in coupled oscillator systems, Phys. Rev. E 62(1), 402 (2000)

[9] D. A. Paley, N. E. Leonard, and R. Sepulchre, Oscillator models and collective motion: Splay state stabilization of self-propelled particles, in: Proc. 51st IEEE Conf. Decision Control, pp 3935-3940 (2005)

[10] M. Silber, L. Fabiny, and K. Wiesenfeld, Stability results for in-phase and splay-phase states of solid-state laser arrays, J. Opt. Soc. Am. B 10(6), 1121 (1993)

[11] S. H. Strogatz and R. E. Mirollo, Splay states in globally coupled Josephson arrays: Analytical prediction of Floquet multipliers, Phys. Rev. E 47(1), 220 (1993)

[12] L. Lü, C. Li, W. Wang, Y. Sun, Y. Wang, and A. Sun, Study on spatiotemporal chaos synchronization among complex networks with diverse structures, Nonlinear Dyn. 77(1–2), 145 (2014)

[13] J. Gómez-Gardeñes, S. Gomez, A. Arenas, and Y. Moreno, Explosive synchronization transitions in scale-free networks, Phys. Rev. Lett. 106(12), 128701 (2011)

[14] O. E. Omel’chenko and M. Wolfrum, Nonuniversal transitions to synchrony in the Sakaguchi–Kuramoto model, Phys. Rev. Lett. 109(16), 164101 (2012)

[15] D. Topaj and A. Pikovsky, Reversibility vs. synchronization in oscillator lattices, Physica D 170(2), 118 (2002)

[16] L. Zhou, C. Wang, Y. Lin, and H. He, Combinatorial synchronization of complex multiple networks with unknown parameters, Nonlinear Dyn. 79(1), 307 (2015)

[17] X. Zhang, X. Hu, J. Kurths, and Z. Liu, Explosive synchronization in a general complex network, Phys. Rev. E 88, 010802(R) (2013)

[18] Z. Zheng, Spatiotemporal Dynamics and Collective Behaviors in Coupled Nonlinear Systems, Beijing: Higher Education Press, 2004 (in Chinese)

[19] N. Yao and Z. Zheng, Chimera states in spatiotemporal systems: Theory and applications, Int. J. Mod. Phys. B 30(1), 163002 (2016)

[20] E. Ott and T. M. Antonsen, Low dimensional behavior of large systems of globally coupled oscillators, Chaos 18(3), 037113 (2008)

[21] S. Watanabe and S. H. Strogatz, Integrability of a globally coupled oscillator array, Phys. Rev. Lett. 70(16), 2391 (1993)

[22] S. Watanabe and S. H. Strogatz, Constants of motion for superconducting Josephson arrays, Physica D 74(3–4), 197 (1994)

[23] J. Gao, C. Xu, Y. Sun, and Z. Zheng, Order parameter analysis for low-dimensional behaviors of coupled phase oscillators, Sci. Rep. 6, 30184 (2016)

[24] X. Hu, S. Boccaletti, W. Huang, X. Zhang, Z. Liu, S. Guan, and C. H. Lai, Exact solution for first-order synchronization transition in a generalized Kuramoto model, Sci. Rep. 4, 7262 (2014)

[25] I. Leyva, R. Sevilla-Escoboza, J. M. Buldu, I. Sendina-Nadal, J. Gomez-Gardenes, A. Arenas, Y. Moreno, S. Gomez, R. Jaimes-Reategui, and S. Boccaletti, Explosive first-order transition to synchrony in networked chaotic oscillators, Phys. Rev. Lett. 108(16), 168702 (2012)

[26] C. Xu, J. Gao, Y. Sun, X. Huang, and Z. Zheng, Explosive or continuous: Incoherent state determines the route to synchronization, Sci. Rep. 5, 12039 (2015)

[27] P. Li, K. Zhang, X. Xu, J. Zhang, and M. Small, Reexamination of explosive synchronization in scale-free networks: The effect of disassortativity, Phys. Rev. E 87(4), 042803 (2013)

[28] T. K. Peron and F. A. Rodrigues, Explosive synchronization enhanced by time-delayed coupling, Phys. Rev. E 86(1), 016102 (2012)

[29] P. Ji, T. K. Peron, P. J. Menck, F. A. Rodrigues, and J. Kurths, Cluster explosive synchronization in complex networks, Phys. Rev. Lett. 110(21), 218701 (2013)

[30] L. Zhang, J. Chen, B. Sun, Y. Tang, M. Wang, Y. Li, and S. Xue, Nonlinear dynamic evolution and control in a new scale-free networks modeling, Nonlinear Dyn. 76(2), 1569 (2014)

[31] I. Leyva, A. Navas, I. Sendina-Nadal, J. A. Almendral, J. M. Buldu, M. Zanin, D. Papo, and S. Boccaletti, Explosive transitions to synchronization in networks of phase oscillators, Sci. Rep. 3, 1281 (2013)

[32] C. Wang, A. Pumir, N. B. Garnier, and Z. Liu, Explosive synchronization enhances selectivity: Example of the cochlea, Front. Phys. 12(5), 128901 (2017)

[33] A. Bergner, M. Frasca, G. Sciuto, A. Buscarino, E. J. Ngamga, L. Fortuna, and J. Kurths, Remote synchronization in star networks, Phys. Rev. E 85(2), 026208 (2012)

[34] O. Burylko, Y. Kazanovich, and R. Borisyuk, Bifurcations in phase oscillator networks with a central element, Physica D 241(12), 1072 (2012)

[35] S. J. S. Theesar, M. R. K. Ariffin, and S. Banerjee, Synchronization and a secure communication scheme using optical star network, Opt. Laser Technol. 54, 15 (2013)

[36] V. Vlasov, A. Pikovsky, and E. E. N. Macau, Star-type oscillatory networks with generic Kuramoto-type coupling: A model for Japanese drums synchrony, Chaos 25(12), 123120 (2015)

[37] C. Xu, Y. Sun, J. Gao, T. Qiu, Z. Zheng, and S. Guan, Synchronization of phase oscillators with frequency-weighted coupling, Sci. Rep. 6, 21926 (2016)

[38] C. Xu, H. Xiang, J. Gao, and Z. Zheng, Collective dynamics of identical phase oscillators with high-order coupling, Sci. Rep. 6, 31133 (2016)

[39] X. Huang, J. Gao, Y. Sun, Z. Zheng, and C. Xu, Effects of frustration on explosive synchronization, Front. Phys. 11(6), 110504 (2016)

[40] C. J. Goebel, Comment on “Constants of motion for superconductor arrays”, Physica D 80(1–2), 18 (1995)

[41] S. A. Marvel, R. E. Mirollo, and S. H. Strogatz, Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action, Chaos 19(4), 043104 (2009)

[42] H. Sakaguchi and Y. Kuramoto, A soluble active rotater model showing phase transitions via mutual entertainment, Prog. Theor. Phys. 76(3), 576 (1986)

[43] S. A. Marvel and S. H. Strogatz, Invariant submain-fold for series arrays of Josephson junctions, Chaos 19(1), 013132 (2009)

[44] F. Dorfler and F. Bullo, Exploring synchronization in complex oscillator networks, in: Proc. 51st IEEE Conf. Decision Control, pp 7157-7170 (2012)

Rights & permissions

Higher Education Press and Springer-Verlag Berlin Heidelberg 2017