Order parameter analysis of synchronization transitions on star networks

TL;DR

This paper introduces an order parameter approach to analyze synchronization transitions in star networks of coupled oscillators, revealing diverse collective states and bifurcations using nonlinear dynamics tools.

Contribution

It develops a low-dimensional dynamical framework for star network synchronization using the Watanabe-Strogatz, Ott-Antonsen, and ensemble order parameter methods.

Findings

Different solutions correspond to various collective states.

Bifurcations indicate transitions among these states.

Analysis tools include time reversibility and stability analysis.

Abstract

Collective behaviors of populations of coupled oscillators have attracted much attention in recent years. In this paper, an order parameter approach is proposed to study the low-dimensional dynam- ical 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 diverse col- lective states, and different bifurcations reveal various transitions among these collective states. The properties of various transitions are revealed in the star-network model by using tools of nonlinear dynamics such as time reversibility analysis and linear stability analysis.

Tables

Table 1. Table 1: The stability conditions of the four fixed points.

Fixed point	Stability condition
$(x_1,y_1)$
	$\lambda >0,\alpha \in (-\pi /2,{\alpha }_0^{-})$
$(x_2,y_2)$	$\lambda >{\lambda }_{sc}^{+},\alpha \in ({\alpha }_0^{+},\pi /2)$
$(x_3,y_3)$	$\lambda >{\lambda }_{sc}^{-},\alpha \in ({\alpha }_0^{-},0)$
$(x_4,y_4)$	always unstable

Equations

$${\dot{\phi }}_j=f{e}^{i{\phi }_j}+g+\bar{f}{e}^{-i{\phi }_j},j=1,\cdots ,N,$$

$$\frac{\partial \rho }{\partial t}+\frac{\partial (\rho \nu )}{\partial \varphi }=0,$$

$$\nu (\varphi ,t)=f{e}^{i\phi }+g+\bar{f}{e}^{-i\phi }.$$

$$\rho (\varphi ,t)=\frac{1}{2\pi }\{1+\sum _{n=1}^{\infty }(\bar{z}(t{)}^n{e}^{in\varphi }+z(t{)}^n{e}^{-in\varphi })\}$$

$$\rho (\varphi )=\frac{1}{2\pi }\frac{1-r^2}{1-2r\cos (\varphi -\Phi )+r^2},$$

$$z=r{e}^{i\Phi }.$$

$$\dot{z}=i(f{z}^2+gz+\bar{f}).$$

$$<e^{i\varphi }>={\int }_0^{2\pi }{e}^{i\varphi }\rho (\varphi ,t)d\varphi .$$

$$z=<e^{i\varphi }>=r{e}^{i\Phi }$$

$$e^{i{\phi }_j(t)}=M_t(e^{i{\theta }_j})$$

$$M_t(w)=\frac{e^{i\psi }w+\eta }{1+\bar{\eta }{e}^{i\psi }w},$$

$$w_j=e^{i{\theta }_j}.$$

$$\dot{\eta }$$

$$\dot{\psi }$$

$$z(t)=\frac{1}{N}\sum _{j=1}^N\frac{e^{i\psi }{e}^{i{\theta }_j}+\eta (t)}{1+\bar{\eta }(t)e^{i\psi }{e}^{i{\theta }_j}},$$

$${\theta }_j=2\pi \frac{j-1}{N},1\le j\le N,$$

$$z(t)=\eta (t)(1+I),$$

$$z(t)\approx \eta (t),N\gg 1.$$

$${\dot{\theta }}_h$$

$${\dot{\theta }}_j$$

$${\dot{\phi }}_i=\Delta \omega -\lambda \sum _{j=1}^K\sin ({\phi }_j+\alpha )-\lambda \sin ({\phi }_i-\alpha ),$$

$$z(t)\equiv r(t)e^{i\Phi (t)}=\frac{1}{K}\sum _{j=1}^K{e}^{i({\phi }_j)}.$$

$${\dot{\phi }}_j=f{e}^{i{\phi }_j}+g+\bar{f}{e}^{-i{\phi }_j},j=1,\cdots ,K,$$

$$\dot{z}=-\frac{\lambda }{2}{e}^{-i\alpha }{z}^2+i(\Delta \omega -\lambda Kr\sin (\Phi +\alpha ))z+\frac{\lambda }{2}{e}^{i\alpha },$$

$$\dot{x}$$

$$+\lambda (K-1)\sin \alpha xy-\Delta \omega y+\frac{\lambda }{2}\cos \alpha ,$$

$$\dot{y}$$

$$-\lambda (K+1)\cos \alpha xy+\Delta \omega x+\frac{\lambda }{2}\sin \alpha .$$

$$x_{1,2}=\frac{-\sin \alpha \Delta \omega \pm A\sin \alpha }{\lambda (2K\cos 2\alpha +1)},$$

$$y_{1,2}=-\frac{-\cos \alpha \Delta \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(\frac{\sin 2\alpha }{2}B-\sin 2{\alpha }^2)}{\lambda \sin \alpha (K^2+2\cos (2\alpha )K+1)},$$

$$y_{3,4}=\frac{-\Delta \omega (-\cos \alpha \pm \sin \alpha B-K\cos \alpha )}{\lambda (K^2+2\cos (2\alpha )K+1)},$$

$$A=\sqrt{-2K{\lambda }^2\cos 2\alpha -{\lambda }^2+\Delta {\omega }^2},$$

$$B=\sqrt{{\lambda }^2+K^2{\lambda }^2+2K{\lambda }^2\cos 2\alpha -\Delta {\omega }^2}.$$

$$\lambda \le {\lambda }_1=\frac{\Delta \omega }{\sqrt{2K\cos 2\alpha +1}},$$

Full text

††preprint: 2016 Conference paper

I Introduction

Understanding the intrinsic microscopic mechanism embedded in collective macroscopic behaviors of populations of coupled units on heterogenous networks has become a focus in a variety of fields, such as the biological neurons circadian rhythm, chemical reacting cells, and even society systems kuramoto1984chemical ; acebron2005kuramoto ; strogatz2000kuramoto ; pikovsky2002synchronization ; dorogovtsev2008critical ; arenas2008synchronization ; zheng1998phase ; zheng2000pre . Numerous different emerging macroscopic states/phases have been revealed, and various non-equilibrium transitions among these states have been observed and studied on heterogenous networks Paley2005 ; Silber1993 ; Strogatz1993 ; lu2014 ; gomez2011explosive ; omel2012nonuniversal ; topaj2002reversibility ; zhou2015 ; zhang2013explosive ; zheng2004book ; yaonan-review .

The transitions among 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 indicate the competitions of miscellaneous attractors and their corresponding basins of attraction in 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  ott2008low 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 an exact three-dimensional dynamics  ws1993 ; ws1994 , but this scheme cannot be extended to general finite systems. The mechanism of the validity of the OA approach was recently studied, and the ensemble order parameter approach is proposed, which extends the OA approach to more general cases such as a finite-number of oscillators and more general coupling forms  gaojian2016eop .

Abrupt or explosive transition from incoherent state to synchronization may occur on networks if the frequencies of oscillators on nodes are positively correlated to the node’s degrees  gomez2011explosive , which has been observed numerically on scale-free networks and experimentally in electronic circuits hu2014exact ; leyva2012explosive . The first-order transition can be changed and more ways of transitions can be observed by adjusting the phase shift among oscillators xucan2015star . 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 zhang2013explosive ; leyva2012explosive ; li2013reexamination ; peron2012explosive ; ji2013cluster ; zhang2014 ; leyva2013explosive ; fop1 .

It is valuable to analytically understand the transitions among various synchrony states on heterogeneous networks. The star topology is the simplest while the key topology in describing the heterogeneity property of complex networks such as the scale-free networks  Bergner2013remote ; Oleksandr2012Bifurcations ; Theesar2013 ; Vladimir2015 . In this paper, we study the collective states and the abundant transitions among these states on a star network by considering the effect of the phase shift among coupled oscillators  xucan2015star ; xucan2016a ; gaojian2016eop ; xucan2016b ; huangxia2016 . The dynamics of star networks of oscillators is analytically studied by building the equations of motion of the order parameter for networks with a finite size, which accomplishes a great reduction from 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 topaj2002reversibility 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.

II The Ott-Antonsen ansatz and the Watanabe-Strogatz approach

We first illustrate the Ott-Antonsen ansatz ott2008low briefly by analyzing the following class of identical oscillators governed by the equations of motion

[TABLE]

where $f$ is a smooth, complex-valued 2$\pi$-periodic function of the phases $\varphi_{1},\cdots,\varphi_{N}$ and the overbar denotes complex conjugate, $\it{g}$ is a real valued function since $\dot{\varphi}_{j}$ is real. In the limit $N\longrightarrow\infty$, by introducing the distribution of phases of oscillators, the evolution of the system (1) is given by the continuity equation

[TABLE]

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

[TABLE]

Suppose $\rho$ is of the form

[TABLE]

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

[TABLE]

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

[TABLE]

The ansatz (4) defines a submanifold in the infinite-dimensional space of the density function $\rho$. 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.  ott2008low is the invariance of the Poisson submanifold, i. e., if the initial phase density is a Poisson kernel, it remains a Poisson kernel for all the time. This can be verfied 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 $e^{in\phi}$ are simultaneously satisfied if and only if $z(t)$ evolves according to

[TABLE]

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 $\phi$ regarded as points $e^{i\phi}$ on the unit circle:

[TABLE]

By substituting Eq. (4) into Eq. (8) one may find that

[TABLE]

for all states on the Poisson submanifold. Then the meaning of $z$ is clear that it represents the order parameter of the system, $r$ is the modulus of it and $\Phi$ is the mean phase of it. However, whether the governing equation Eq. (7) can be used for system with finite size can not be implied from Ott-Antonsen ansatz.

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 ws1993 ; ws1994 , and the governing equations of the system can also be generated by the Möbius group action charles1995 ; marvel2009identical . The class of identical oscillators still governed by the equations of motion Eq. (1), then the oscillators’ phases $\varphi_{j}(t)$ evolve according to the action of the Möbius group on the complex unit cycle

[TABLE]

for $j=1,\ldots,N$, where $M_{t}$ is a one-parameter family of Möbius transformations and $\theta_{j}$ is a constant angle. By parameterizing the one-parameter family of Möbius transformations as

[TABLE]

where $|\eta(t)|<1$ and $\psi(t)\in{R}$, and let

[TABLE]

One then obtains

[TABLE]

With these new variables, one could rewrite the order parameter as

[TABLE]

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

[TABLE]

with which, the order parameter (14) reads

[TABLE]

where $I=(1-|\eta(t)|^{-2})/(1\pm(e^{i\psi}\bar{\eta}(t))^{-N})$, ”$-$” for the case with even N and ”$+$” for the case with odd N. One can verify that for large N, $I\ll 1$, the order parameter could be approximated as

[TABLE]

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

III The Sakaguchi-Kuramoto model on star networks: The order parameter equation

We start 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,...,K$) and the degree of the hub is $k_{h}=K$. Suppose that the natural frequencies of the oscillators are proportional to their degrees, the equations of motion for the hub and leaf nodes read

[TABLE]

where $\theta_{h},\theta_{j}$ and $\omega_{h},\omega$ are instantaneous phases and natural frequencies of the hub and leaf nodes respectively. $\lambda$ is the coupling strength. $K$ is the number of leaf nodes connected with this hub and $\alpha$ is the phase shift. The effect of phase shift among coupled oscillators has been extensively investigated in recent years, while this has been seldom discussed in star networks  sakaguchi1986soluble ; omel2012nonuniversal . Abundant collective dynamics appear in the global coupled model in the presence of a finite phase shift omel2012nonuniversal , where synchrony can decay or incoherence can regain its stability with increasing coupling and multistability between partially synchronized and/or the incoherent state can appear in the globally coupled network.

The coupled phase oscillator system with $\alpha=0$ of the star network was originally to study the characteristics of the explosive synchronization gomez2011explosive , however the process of the synchronization may be influenced with the introduction of phase shift omel2012nonuniversal . By introducing the phase differences between the hub and leaves $\varphi_{j}=\theta_{h}-\theta_{j}$, the phase dynamics on star networks can be transformed to the following phase difference dynamics on an all-connected network,

[TABLE]

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

[TABLE]

It is worth noting that 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 $\varphi_{j}(t)=\varphi(t)$, all the leaf nodes are synchronous to each other while they are asynchronous to the hub oscillator.

It is instructive to rewrite Eq. (19) as

[TABLE]

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

For finite $K$, due 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  gaojian2016eop and the the Watanabe-Strogatz transformation ws1993 ; ws1994 . If the initial phases of oscillators are chosen as (15), we have the dynamical equation for the order parameter $z(t)$ as

[TABLE]

which is just 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.

IV Collective Dynamics of Stationary States

In the following we start discussing the collective dynamics of the star network in terms of the Eq. (22). By setting $z=x+iy$, we can describe the order parameter dynamics in the $x-y$ plane as

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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

[TABLE]

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}\leq 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\times 2$ Jacobian matrix $J$ of the fixed points with elements

[TABLE]

The eigenvalues of the Jacobian matrix are

[TABLE]

The stability conditions of the four fixed points are summarized in Table 1, where The parameters in the table are $\hat{\lambda}_{c}^{f}=\footnotesize{{\Delta\omega}/{\sqrt{2K\cos 2\alpha+1}}}$, $\lambda_{sc}^{+}=\footnotesize{{-\Delta\omega}/({K\cos 2\alpha+1}})$, $\lambda_{sc}^{-}=\footnotesize{{\Delta\omega}/({K\cos 2\alpha+1}})$, $\alpha_{0}^{-}=-\arccos(-1/K)/2$, $\alpha_{0}^{+}=\arccos(-1/K)/2$.

The fixed points of order parameter equation are related to the collective states of coupled phase oscillators. 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 hub and leaf nodes are the same and keep constant as

[TABLE]

implying the global synchronization of leaves and the hub in a star network. The above stability analysis indicates that 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 easily studied.

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 coupled oscillators is bounded by $|z|\leq 1$. If $|z|<1$, the related collective state is called the splay state (SPS) Strogatz1993 ; marvel2009invariant , the phase differences between the hub and leaf nodes satisfy a function relation as

[TABLE]

with $T$ the period of $\varphi(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 move synchronously with a constant time shift.

V Collective Dynamics of Time-dependent States

Long-term solutions of the order parameter equation contains not only the states given by fixed points but also time-dependent states corresponding to periodic solutions. There are two periodic regimes, i.e., the regime $0<\alpha<\pi/2$ and $\lambda<\lambda_{ec}=\lambda_{2}$, and the critical line with $\alpha=0,\pm\pi/2$:

1. The in-phase state

When $\lambda<\lambda_{ec}$ and $0<\alpha<\pi/2$, 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=re^{i\Phi}$:

[TABLE]

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.,

[TABLE]

This state corresponds to the phases of all leaf nodes evolve synchronously, while they are not synchronous to the hub, as shown in Fig. 1(c). In this case, the model Eq. (18) 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 $0<\alpha<\pi/2$ and unstable for $-\pi/2<\alpha<0$.

2. The neutral state

The dynamics at $\alpha=0,\pm\pi/2$ when $\lambda<\lambda_{ec}$ are spacial 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})\neq 0$ in Eq. (30). In this case there are a large class of states in the critical cases with the order parameter $r$ are determined by initial values of $(x,y)$, where long-term behaviors of $z$ depend crucially on initial phases, as shown in Fig. 1(d). We call this state the neutral state(NS) topaj2002reversibility , and the corresponding fixed point is neutrally stable, where all Lyapunov exponents of the fixed point are zero when $\lambda<\lambda_{ec}$, as shown in Fig. 2(a). Because the Kuramoto system is dissipative florian2012disspetive , 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 $\alpha=0$, Eq. (23) can be simplified to

[TABLE]

The fixed points are determined by setting $\dot{x}=0$ and $\dot{y}=0$. When $\lambda<\lambda_{ec}$ 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 the neutral state with $\alpha=0$ are zero. It is worthy to note that 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 topaj2002reversibility . This symmetry endows the system many interesting properties.

Note that, 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 to an attractor, the backward trajectory will evolve to the symmetric repeller of the system. Then the attractor and the repeller of the system emerge in pairs. When the trajectory crosses the invariant set more than once, the forward and backward trajectory will coincide with each other, forming the periodic solution for the system, which is called the reversible trajectory topaj2002reversibility . For any reversible trajectory, the Lyapunov exponents have the sign-symmetry form and the volume of phase space in the vicinity of it are conserved in average as we discovered in numerical simulations.

For our order parameter plane of the system, it is bounded by the unit circle with the invariant set as $x=0,y>0$, the attractor and the repeller emerge in the same time, implying that if the plane only exist one fixed point, it is neither the attractor nor the repeller, i. e., the only fixed point is neutrally stable. Suppose that there only exist one neutrally stable fixed point in the plane, 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 $\alpha=0$ with the region $\lambda<\lambda_{ec}$ as shown in Fig. 2(b). For $\alpha=0$ and $\lambda>\lambda_{ec}$, there is a coexisting region for the synchronous state and the neutral state as the critical cases for the coexistence region for the neutral state and the synchronous state.

All the above possible collective states are summarized in the parameter space ($\alpha$,$\lambda$) as a phase diagram Fig. 3 with the boundaries we get analytically above both from existence and stability conditions. In Fig. 3, four regions of the phase shift $\alpha$ can be identified. For the first region $-\pi/2<\alpha<\alpha_{0}^{-}$, the splay state exists and is stable for any $\lambda$. With the increase of coupling strength $\lambda$, unstable synchronous state exists above the threshold $\lambda>\lambda_{ec}$. For the second region $\alpha_{0}^{-}<\alpha<0$, the splay state exists and is stable within $0<\lambda<\hat{\lambda}_{c}^{f}$, and the synchronous state exist with $\lambda>\lambda_{ec}$ but unstable unless $\lambda>\lambda_{sc}^{-}$. Obviously there exist a co-existing region for the splay state and the synchronous state in the region $-\alpha_{0}^{-}<\alpha<0$, within the coupling interval $\lambda_{sc}^{-}<\lambda<\hat{\lambda}_{c}^{f}$. In the third region where $0<\alpha<\alpha_{0}^{+}$, the splay state is always unstable, the stable synchronous state emerges as the coupling strength $\lambda>\lambda_{ec}$. For the forth region $\alpha_{0}^{+}<\alpha<\pi/2$, the splay state always exist but only stable when $\lambda>\lambda_{sc}^{+}$, and the synchronous state only exist and stable in the region $\lambda_{ec}<\lambda<\lambda_{sc}^{+}$. 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 $0<\alpha<\pi/2$, within the coupling range $0<\lambda<\lambda_{ec}$. The variety of states in the phase diagram leads to various transitions among them.

VI Scenarios of Synchronization Transitions

The phase diagram shown in Fig. 3 presents a great variety of transitions among the different collective dynamical states. One can also find 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.

1. Synchronization transition from the neutral state

We first investigate the synchronization process from the neutral state to the synchronous state for $\alpha=0$. The synchronization process when $\alpha=0$ is discontinuous known as the explosive synchronization which has attracted much attention recently gomez2011explosive . It has been shown that with changing the coupling strength $\lambda$, this kind of transition is abrupt, and there is a hysteretic behavior at the onset of synchronization, and $\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 initial states as shown in Fig. 4(a). The upper limit of $\lambda_{c}^{f}$ is denoted by $\hat{\lambda}_{c}^{f}$. As $\lambda>\hat{\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

[TABLE]

The analytical curve and the simulation results are 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 with A an attractor, C a repeller and B, D neurally stable. Any orbits crossing the nullcline A-B-C will eventually fall to A, and others will hold the property as periodic orbits. It is clear that the stable fixed point A corresponds to the synchronous state. And the basin for the neutral state can be calculated approximated by the circle which has its center in point $D$ and radius as the length of line $B-D$. As $\lambda$ increases, points D and B close to each other and eventually collide at a critical coupling, as shown in Fig. 4(d), and the synchronous state becomes globally attractive.

2. Synchronization transition from the splay state

The synchronization process from the splay state to the synchronous state for $\alpha_{0}^{-}<\alpha<0$ 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

[TABLE]

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 Fig. 5(c,d). Fig. 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 coupling $\lambda$ 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 corresponds to the synchronous state a global attractor.

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 synchronous state are shown in Fig. 6(c,d). As shown in Fig. 6(c), the in-phase state is a limit cycle in the order parameter plane. As $\lambda$ 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).

4. Scenario of desynchronization

In the region $\alpha_{0}^{+}<\alpha<\frac{\pi}{2}$ of the phase diagram 3, one may find 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 the desynchronization, and it is a continuous transition as shown in Fig. 7(a). The order parameter $r$ decrease rapidly at the threshold and effective frequencies of hub and leaf nodes are divided at the same coupling $\lambda$. It is easy to know the route of the de-synchronization is from the synchronous state to the splay state from the view of the phase diagram. The threshold of the de-synchronization is $\lambda_{sc}^{+}=\footnotesize{{-\Delta\omega}/({K\cos 2\alpha+1}})$ as shown in Fig. 7(b), the simulation results are consistent with it obviously.

The dynamical manifestations of the continuous transition from the synchronous state to the splay state are shown in Fig. 7(c,d). Fig. 7(c) exhibits the stable synchronous state of the system when $\lambda<\lambda_{sc}^{+}$, all the orbits in the phase of the order parameter will evolve to the fixed point A eventually. As $\lambda$ increases and larger than the critical coupling $\lambda_{sc}^{+}$, the two nullclines will intersect in four fixed points as shown in Fig. 7(d), the point A loses its stability and a new stable fixed point B which corresponds to the splay state appears. The process from the synchronous state to the splay state is finished by this bifurcation continuously.

VII Conclusion

To summarize, in this paper we study the dynamics of coupled oscillators on a star network with the Sakaguchi-Kuramoto model by resorting to 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 star network can also be approximately described from the Ott-Antonsen ansatz, which is originated from the high symmetry of the topology. By reducing from a high-dimensional phase space to a much lower-dimensional order parameter space without additional approximation, one is able to grasp analytically the essential dynamical mechanism of different scenarios of synchronization. Different solutions of the order parameter equation corresponds to the various collective states of coupled oscillators, and different bifurcations reveal various transitions among those collective states. The process of those transitions are revealed in the plane of order parameter and the critical coupling strengths of them are obtained analytical which are verified by the simulation results.

This work is partially supported by the National Natural Science Foundation of China (Grant No. 11075016 and 11475022) and the Scientific Research Funds of Huaqiao University.