Low dimensional behavior of explosive synchronization on star graphs

Synchronization phenomena, which describe the emergence of collected behavior of large ensembles of coupled units, are always the focus of the intense research in synthetic and natural systems [1–3]. In 1998, Watts and Strognatz showed that the synchronization dynamics were strongly influenced by the structure of the interaction between the units in their seminal work [4], which was the seed of the modern theory of complex networks [5–7]. Since then, the synchronization processes with non trivial interaction patterns have been widely studied [8]. The main interest of these studies was the phase transitions from the incoherent states to the synchronized states, and it has been proven that these phase transitions are usually second-order [8].

In 2011, Gómez-Gardeñes et al [9] found that an explosive synchronization (ES) transition occurs on a scale-free (SF) network [6] for the Kuramoto model [10, 11] when the frequencies of oscillators positively correlate to the degrees of nodes. ES soon attracted many researchers' interest and become a hot topic in the study of dynamics on complex networks. Subsequent papers have further explored the influences of degree-frequency correlations [12, 13] and network structures [14, 15] on ES. ES has also been extended to some other models, such as Rössler units [16] and the weighted Kuramoto model [17]. These studies have greatly improved our understanding about many aspects of ES. However, as far as we know, the intrinsic mathematical mechanism of the emergence of ES on heterogenous networks has not been well clarified yet. To uncover this puzzle, an analytical framework that can deal with the high dimensional equations of the motion of the oscillators system is needed.

Ott and Antonsen introduced an ansatz to study the low dimensional behavior of global coupled Kuramoto oscillators in [18]. Using this anstaz, they derived a finite set of nonlinear ordinary differential equations for the macroscopic evolution of the system, and obtained a closed form solution for the Kuramoto oscillators with a Lorentzian frequency distribution. This anstaz has proven to be very efficient to solve the Kuramoto model and its variants, such as the Kuramoto model with phase lag [19], time delay [20] and bimodal frequency distribution [21].

In this paper, we use the Ott–Antonsen ansatz to study the time evolution of ES of Kuramoto oscillators. For the sake of mathematical treatment, we reduce the problem studied to the analysis of the star configuration, a special structure that grasps the main property of SF networks, namely, the role of hubs. The main work in this paper consists of two parts. Part one (section 2) discusses the time evolution of the Kuramoto oscillators with strict degree-frequency correlation on star configurations. Part two (section 3) extends the model studied to a general case by considering random perturbations to the frequencies of leaf nodes. We present both numerical simulations and analytical solutions in the two parts. Comparison of the analytical and simulation results shows a good agreement.

Let us consider an ensemble of N  +  1 phase oscillators, interacting in a star network via the coupling

$$\begin{eqnarray}&&{{\dot{\mathop{\phi}}}_{N+1}}=\beta +\lambda \frac{\beta}{N}\underset{j=1}{\overset{N}{\mathop \sum}}\sin \left({{\phi}_{j}}-{{\phi}_{N+1}}\right),\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{{\dot{\mathop{\phi}}}_{j}}=1+\lambda \,\sin \left({{\phi}_{N+1}}-{{\phi}_{j}}\right),1\leqslant j\leqslant N,\end{eqnarray} \tag{ 2 }$$

where ${{\phi}_{N+1}}$ and ${{\phi}_{j}}$ are the phases of the hub and leaf nodes, respectively, and λ is the coupling strength. The natural frequencies are set as 1 for the leaves and $\beta >1$ for the hub. The model we proposed here is slightly different from the original degree-frequency correlated Kuramoto model on star graphs [9]. The rescaled coupling term of hub $\lambda \frac{\beta}{N}$ is to guarantee the convergence of the right side of equation (1) as $N\to \infty$. The Kuramoto order parameter that measures the coherence of the system is defined as

$$\begin{eqnarray}&&r=|\frac{1}{N}\underset{j=1}{\overset{N}{\mathop \sum}}{{\text{e}}^{\text{i}{{\phi}_{j}}}}|.\end{eqnarray} \tag{ 3 }$$

The small value of r means an incoherence state while $r\approx 1$ indicates a highly coherence solution.

Let ${{\theta}_{j}}={{\phi}_{j}}-{{\phi}_{N+1}}$ and then we have

$$\begin{eqnarray}&&\dot{{{\theta_{j}}}}\,=1-\beta -\lambda \,\sin \,{{\theta}_{j}}-\lambda \frac{\beta}{N}\underset{l=1}{\overset{N}{\mathop \sum}}\sin \,{{\theta}_{l}},1\leqslant j\leqslant N.\end{eqnarray} \tag{ 4 }$$

In the limit where $N\to \infty$, we can define a probability density function $\rho \left(\theta ,t\right)$ so that the fraction of oscillators with phases between θ and $\theta +\text{d}\theta$ at time t is given by $\rho \left(\theta ,t\right)\text{d}t$. Let $D=\int_{0}^{2\pi}\rho \,\sin \,\theta \text{d}\theta$, then equation (3) can be written as

$$\begin{eqnarray}&&\dot{{{\theta_{j}}}}\,=1-\beta -\lambda \,\sin \,{{\theta}_{j}}-\lambda \beta D,1\leqslant j\leqslant N.\end{eqnarray} \tag{ 5 }$$

Noticing that each oscillator in equation (5) moves with an angular or drift velocity $v=1-\beta -\lambda \,\sin \,\theta -\lambda \beta D$, we can obtain a continuity equation for ρ:

$$\begin{eqnarray}&&\frac{\partial \rho}{\partial t}+\frac{\partial \left(\rho v\right)}{\partial \theta}=0.\end{eqnarray} \tag{ 6 }$$

Using the Ott–Antonsen ansatz, we expand ρ into a Fourier series

$$\begin{eqnarray}&&\rho \left(\theta ,t\right)=\frac{1}{2\pi}\left(1+\left(\underset{n=1}{\overset{\infty}{\mathop \sum}}\alpha {{(t)}^{n}}{{\text{e}}^{\text{i}n\theta}}+\text{c}\text{.c}\right)\right),\end{eqnarray} \tag{ 7 }$$

where c.c stands for the complex conjugate of the preceding term. Substituting this expansion into equation (6), and noticing $D=\frac{\text{i}\left(\alpha -\bar{\alpha}\right)}{2}$ then we obtain

$$\begin{eqnarray}&&\dot{\mathop{\alpha}}+\frac{\lambda \left(\beta +1\right)}{2}{{\alpha}^{2}}+\text{i}\alpha \left(1-\beta \right)-\frac{\lambda \beta}{2}|\alpha {{|}^{2}}-\frac{\lambda}{2}=0.\end{eqnarray} \tag{ 8 }$$

Noticing that $r=|\int_{0}^{2\pi}\rho {{\text{e}}^{\text{i}\theta}}\text{d}\theta |=|\alpha |$, we can write $\alpha =r{{\text{e}}^{\text{i}\psi}}$. Substituting this into equation (8), multiplying by ${{\text{e}}^{\text{i}\psi}}$, and taking the real and imaginary part of the result, thus we obtain

$$\begin{eqnarray}&&\dot{\mathop{r}}=\frac{\lambda}{2}\left(1-{{r}^{2}}\right)\cos \,\psi ,\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&\dot{\mathop{\psi}}=\beta -1-\left(\frac{2\beta +1}{2}r+\frac{1}{2r}\right)\lambda \,\sin \,\psi .\end{eqnarray} \tag{ 10 }$$

Obviously, equation (9) always has a solution r  =  1. Therefore, a trajectory of the above system, starting with an initial condition satisfying r(0)  <  1 cannot cross the unit circle in the complex α-plane, and we have r(t)  <  1 for all finite time t.

We now consider the solution with the restriction r(0)  <  1. First, we can easily obtain the fixed points of the above equations by letting the right side of equation (9) and equation (10) be equal to zero:

$$\begin{eqnarray}&&{{P}_{1,2}}=\left(\frac{\beta -1\pm \sqrt{{{\left(\beta -1\right)}^{2}}-{{\lambda}^{2}}\left(2\beta +1\right)}}{\lambda \left(2\beta +1\right)},\frac{\pi}{2}\right),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{{P}_{3}}=\left(1,\arcsin \frac{\beta -1}{\lambda \left(\beta +1\right)}\right),\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&{{P}_{4}}=\left(1,\pi -\arcsin \frac{\beta -1}{\lambda \left(\beta +1\right)}\right).\end{eqnarray} \tag{ 13 }$$

Here P_1,2 only exist when $0\leqslant \lambda \leqslant \frac{\beta -1}{\sqrt{2\beta +1}}$ while P_3,4 only exist when $\lambda \geqslant \frac{\beta -1}{\beta +1}$.

To study the long term behavior of the system, we first need to know the stability of the points. To do this, we take the Jacobian matrices of the right side of equations (9) and (10) and compute the corresponding eigenvalues at these fixed points. The results are shown directly in table 1 (this process is quite simple, so we do not list details here). It should be specially explained that for $0\leqslant \lambda \leqslant \frac{\beta -1}{\sqrt{2\beta +1}}$ the stability of P₂ cannot be easily obtained by the linearized method because the corresponding eigenvalues are both purely imaginary. Here, we figure out that the system has a first integral: $\beta -1-\lambda \left(\beta +1\right)r\,\sin \,\psi =C|1-{{r}^{2}}{{|}^{\beta +1}}$. These curves which stand for the solutions of the system are a series of periodic orbits surrounding P₂. Then the fixed point P₂ is actually a neutrally stable center of the system.

Table 1. The stability of fixed points of equations (9) and (10).

Fixed point	Existence region	Types	Stability
P1	$\frac{\beta -1}{\beta +1}\leqslant \lambda <\frac{\beta -1}{\sqrt{2\beta +1}}$	Saddle	Unstable
P2	$0<\lambda <\frac{\beta -1}{\sqrt{2\beta +1}}$	Center	Neutrally stable
P3	$\lambda \geqslant \frac{\beta -1}{\beta +1}$	Sink	Stable
P4	$\lambda \geqslant \frac{\beta -1}{\beta +1}$	Source	Unstable

Table 1 actually shows that the critical values $\frac{\beta -1}{\beta +1}$ and $\frac{\beta -1}{\sqrt{2\beta +1}}$ are the bifurcation values of equations (9) and (10). In fact, these two values correspond to the backward ($\lambda _{c}^{b}$) and forward ($\lambda _{c}^{f}$) critical coupling in the explosive synchronized transition. For $0\leqslant \lambda <\lambda _{c}^{b}$, as we conclude above, the solution is a periodic obit around P₂. The order parameter r evolves with time t at different synchronization levels for different initial conditions. For $\lambda >\lambda _{c}^{f}$, the synchronized state P₃ is the only attractor inside the unit circle r  =  1. The solution of the system with arbitrary initial states will approach the synchronized state when time t tends to infinity. For $\lambda _{c}^{b}\leqslant \lambda <\lambda _{c}^{f}$, the synchronized state represented by P₃ and incoherent state represented by P₂ coexist. The parameter r will approach 1 when initial states are close to P₃ and will oscillate around a certain level while initial states are close to P₂. This result in fact shows the mechanism of an important phenomenon in ES—the hysteresis, which means the forward and backward curves of r do not coincide in one certain coupling region.

Now we perform numerical simulations to check the results described above. First, we compute the time evolutions of r by direction simulation of the original system given in equations (1) and (2), and compare the results to the theoretical solutions solved from the reduced system defined by equations (9) and (10). To do this, we fix $\beta =10$ and choose three typical values of λ, 0.6, 1.2, 2.5, which lie in the coupling intervals $(0,\lambda _{c}^{b})$ (incoherent), $(\lambda _{c}^{b},\lambda _{c}^{f})$ (hysteresis), and $(\lambda _{c}^{f},\infty )$ (synchronized), respectively. The initial states ϕ of the original system are randomly drawn from $\left[0,2\sigma \pi \right]$ where σ is a parameter between 0 and 1. To show the influence of the initial states on order parameter r, we perform numerical simulations with two different initial states given by $\sigma =1$ and $\sigma =0.4$ for each λ. The results are presented in figures 1(a)–(c). The simulations confirm the analysis we decribed above, and show a good agreement with the theoretical predictions of our low dimensional reduced system.

Zoom In Zoom Out Reset image size

Then, we want to verify the theoretical solution of forward critical coupling $\lambda _{c}^{f}$ through numerical simulations. We progressively increase the value of λ and calculate the stationary value of the order parameter r with initial ψ randomly drawn from $\left[0,2\pi \right]$ for various $\beta \in \left[10,100\right]$. If the value of r changes sharply from a small value to a high one at certain coupling λ, then that coupling is regarded as the forward critical coupling $\lambda _{c}^{f}$. We plot the value of $\lambda _{c}^{f}$ obtained from the above simulations as well as the analytical result $\lambda _{c}^{f}=\frac{\beta -1}{\sqrt{2\beta +1}}$ in figure 1(d). As we can see, the two curves fit extremely well.

The perfect agreements between simulations and theoretical predictions confirm the validity of the reduced method used in this section. Following this methodology, we extend our study to a more complex case in the next section.

The SF network with a small average degree can be seen as a collection of star graphs. This kind of network approximately has a tree-like structure (with rare loops). The hubs then dominate most of the links while the links between leaves that belong to different hubs are very few. For simplicity, we see the connections between leaves as small random perturbations to their frequency, and then we consider a modified model that can describe the dynamics of star components of the networks in the following:

$$\begin{eqnarray}&&{{\dot{\mathop{\phi}}}_{N+1}}=\beta +\lambda \frac{\beta}{N}\underset{j=1}{\overset{N}{\mathop \sum}}\sin \left({{\phi}_{j}}-{{\phi}_{N+1}}\right),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{\dot{\mathop{\phi}}}_{j}}=1+{{\omega}_{j}}+\lambda \,\sin \left({{\phi}_{N+1}}-{{\phi}_{j}}\right),1\leqslant j\leqslant N,\end{eqnarray} \tag{ 15 }$$

where ${{\omega}_{j}}$ are randomly drawn from a unimodal and symmetric distribution $g\left(\omega \right)$ with zero mean. For simplicity, we assume that the perturbations ${{\omega}_{j}}$ are distributed according to the Cauchy–Lorentz distribution $g\left(\omega \right)=\frac{\Delta}{\pi \left({{\omega}^{2}}+{{\Delta}^{2}}\right)}$. The other parameters about the modified model are the same as the ones we introduced before.

We implement a similar treatment on the above equations as we did in the previous section. First, let ${{\theta}_{j}}={{\phi}_{j}}-{{\phi}_{N+1}}$ be the phase difference between leaf and hub nodes, and then we get the governing equations:

$$\begin{eqnarray}&&\dot{{{\theta_{j}}}}\,={{\omega}_{j}}+1-\beta -\lambda \,\sin \,{{\theta}_{j}}-\lambda \frac{\beta}{N}\underset{l=1}{\overset{N}{\mathop \sum}}\sin \,{{\theta}_{l}}.1\leqslant j\leqslant N.\end{eqnarray} \tag{ 16 }$$

Considering the limit $N\to \infty$, the states of the oscillators at time t can be described as a probability density function $\rho \left(\omega ,\theta ,t\right)$, where

$$\begin{eqnarray}&&\int_{0}^{2\pi}\rho \left(\omega ,\theta ,t\right)\text{d}\theta =g\left(\omega \right).\end{eqnarray} \tag{ 17 }$$

Let $D\left(\omega \right)=\int_{0}^{2\pi}\rho \left(\omega ,\theta ,t\right)\sin \,\theta \text{d}\theta$, then equation (16) can be written as

$$\begin{eqnarray}&&\dot{{{\theta_{j}}}}\,={{\omega}_{i}}+1-\beta -\lambda \,\sin \,{{\theta}_{j}}-\lambda \beta D\left(\omega \right).1\leqslant j\leqslant N.\end{eqnarray} \tag{ 18 }$$

The the evolution of f is governed by the following equation:

$$\begin{eqnarray}&&\frac{\partial \rho}{\partial t}+\frac{\partial \left(\rho v\right)}{\partial \theta}=0.\end{eqnarray} \tag{ 19 }$$

where $v\left(\omega ,\theta ,t\right)=\omega +1-\beta -\lambda \,\sin \,\theta -\lambda \beta D\left(\omega \right)$. Using the Ott–Antonsen ansatz, we expand ρ into a Fourier Series

$$\begin{eqnarray}&&\rho \left(\omega ,\theta ,t\right)=\frac{g\left(\omega \right)}{2\pi}\left(1+\left(\underset{n=1}{\overset{\infty}{\mathop \sum}}\alpha {{\left(\omega ,t\right)}^{n}}{{\text{e}}^{\text{i}n\theta}}+\text{c}\text{.c}\right)\right),\end{eqnarray} \tag{ 20 }$$

where c.c stands for the complex conjugate of the preceding term. Substituting this expansion into equation (19), and noticing $D=\frac{\text{i}\left(\alpha -\bar{\alpha}\right)}{2}$ then we obtain

$$\begin{eqnarray}&&\dot{\mathop{\alpha}}+\frac{\lambda \left(\beta +1\right)}{2}{{\alpha}^{2}}+\text{i}\left(\omega +1-\beta \right)\alpha -\frac{\lambda \beta}{2}|\alpha {{|}^{2}}-\frac{\lambda}{2}=0.\end{eqnarray} \tag{ 21 }$$

Now we consider the solutions of equation (21) with the following restriction: (i) $|\alpha \left(\omega ,t\right)|<1$; (ii) $\alpha \left(\omega ,t\right)$ can be analytically continued from real ω into the complex ω plane; (iii) $\alpha |\left(\omega ,t\right)|\to 0$ as Im$\left(\omega \right)\to -\infty$. As has been shown in [18], if $\alpha \left(\omega ,0\right)$ satisfies these conditions, they are also satisfied by $\alpha \left(\omega ,t\right)$ for all t  >  0. We define a complex order parameter

$$\begin{eqnarray}&&z(t)=\int_{-\infty}^{+\infty}\int_{0}^{2\pi}\rho \left(\omega ,\theta ,t\right)\text{d}\theta \text{d}\omega \end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&=\int_{-\infty}^{+\infty}\alpha \left(\omega ,t\right)g\left(\omega \right)\text{d}\omega ,\end{eqnarray} \tag{ 23 }$$

where we have $r(t)=|z(t)|$. Expanding $g\left(\omega \right)$ as ${{\left(2\pi i\right)}^{-1}}[{{\left(\omega -\Delta i\right)}^{-1}}-{{\left(\omega +\Delta i\right)}^{-1}}]$, and evaluating equation (22) by doing Cauchy integration along a large semicircle in the lower half ω-plane, we have

$$\begin{eqnarray}&&z(t)=\alpha \left(-\text{i} \Delta ,t\right).\end{eqnarray} \tag{ 24 }$$

Substituting this expression into equation (21), we get an ordinary differential equation governing the evolution of z:

$$\begin{eqnarray}&&\dot{\mathop{z}}=-\left(\Delta +\left(1-\beta \right)i\right)z+\frac{\lambda}{2}[1+\beta |z{{|}^{2}}-\left(\beta +1\right){{z}^{2}}].\end{eqnarray} \tag{ 25 }$$

Let $z=r{{\text{e}}^{\text{i}\psi}}$, then the above equation can be written as two equations of order pamperer r and average phase ψ

$$\begin{eqnarray}&&\dot{\mathop{r}}=-\Delta r+\frac{\lambda}{2}\left(1-{{r}^{2}}\right)\cos \,\psi ,\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&\dot{\mathop{\psi}}=\beta -1-\frac{\lambda}{2}\left[\frac{1}{r}+\left(2\beta +1\right)r\right]\sin \,\psi .\end{eqnarray} \tag{ 27 }$$

This two-dimensional system can be completely solved by using the qualitative theory of differential equations. In the following part of this section, we will study how the time evolution of synchronization parameter r varies with parameter $\Delta$ and λ by analyzing the above reduced equations. First, we let the right side of equations (26) and (27) equal zero and then get equations that the fixed points satisfy

$$\begin{eqnarray}&&\frac{{{\Delta}^{2}}{{r}^{2}}}{{{\left(1-{{r}^{2}}\right)}^{2}}}+\frac{{{\left(\beta -1\right)}^{2}}{{r}^{2}}}{{{\left[\left(2\beta +1\right){{r}^{2}}+1\right]}^{2}}}=\frac{{{\lambda}^{2}}}{4},\end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray}&&\psi =\arccos \frac{2\Delta r}{\lambda \left(1-{{r}^{2}}\right)}.\end{eqnarray} \tag{ 29 }$$

Let u  =  r² and we defined a function of u

$$\begin{eqnarray}&&f(u)=\frac{{{\Delta}^{2}}u}{{{(1-u)}^{2}}}+\frac{{{\left(\beta -1\right)}^{2}}u}{{{\left[\left(2\beta +1\right)u+1\right]}^{2}}},\end{eqnarray} \tag{ 30 }$$

where $u\in \left[0,1\right)$. Obviously, we have f (0)  =  0 and $f(u)\to \infty$ as $u\to 1$. The solutions of equation (28) then correspond to the intersections of curve y  =  f (u) and line $y=\frac{{{\lambda}^{2}}}{4}$. Taking the derivative of f , we get

$$\begin{eqnarray}&&{{f}^{\prime}}(u)=\frac{{{\Delta}^{2}}(1+u)}{{{(1-u)}^{3}}}-\frac{{{\left(\beta -1\right)}^{2}}\left[\left(2\beta +1\right)u-1\right]}{{{\left[\left(2\beta +1\right)u+1\right]}^{3}}}.\end{eqnarray} \tag{ 31 }$$

Let ${{f}^{\prime}}(u)=0$ and then we obtain

$$\begin{eqnarray}&&{{\Delta}^{2}}=\frac{{{\left(\beta -1\right)}^{2}}{{(1-u)}^{3}}\left[\left(2\beta +1\right)u-1\right]}{(1+u){{\left[\left(2\beta +1\right)u+1\right]}^{3}}}=g(u).\end{eqnarray} \tag{ 32 }$$

We have ${{f}^{\prime}}(u)>0$ if $g(u)<{{\Delta}^{2}}$ and ${{f}^{\prime}}(u)<0$ if $g(u)>{{\Delta}^{2}}$. The derivative of g is:

$$\begin{eqnarray}&&{{g}^{\prime}}(u)=\frac{8\left(\beta +1\right){{\left(\beta -1\right)}^{2}}{{(1-u)}^{2}}\left[1-\beta u-\left(2\beta +1\right){{u}^{2}}\right]}{{{(1+u)}^{2}}{{\left[\left(2\beta +1\right)u+1\right]}^{4}}}.\end{eqnarray} \tag{ 33 }$$

${{g}^{\prime}}(u)$ has only one root ${{u}_{c}}=\frac{2}{\beta +\sqrt{{{\beta}^{2}}+8\beta +4}}$ in the interval [0, 1). ${{g}^{\prime}}(u)>0$ if $u\in \left[0,{{u}_{c}}\right)$ and ${{g}^{\prime}}(u)<0$ if $u\in \left({{u}_{c}},1\right)$, meaning that g has maximum value at u_c. Let ${{\Delta}_{c}}=\sqrt{g\left({{u}_{c}}\right)}$, then the problem we discussed can be divided into the following two cases according to the relation of $\Delta$ and ${{\Delta}_{c}}$.

Case(i):$\Delta <{{\Delta}_{c}}$. In this case, ${{f}^{\prime}}(u)=0$ has two roots u_c2 and u_c1 ($0<{{u}_{c2}}<{{u}_{c1}}<1$). Then f (u) increases when $u\in \left[0,{{u}_{c2}}\right]\mathop{\cup}^{}\left[{{u}_{c1}},1\right)$ and decreases while $u\in \left({{u}_{c2}},{{u}_{c1}}\right)$ (see figure 2(a)). We let ${{\lambda}_{c1}}=2\sqrt{f\left({{u}_{c1}}\right)}$ and ${{\lambda}_{c2}}=2\sqrt{f\left({{u}_{c2}}\right)}$. Then for $0\leqslant \lambda \leqslant {{\lambda}_{c1}}$, equation (28) has only one solution u₁, corresponding to the unique fixed points ${{P}_{1}}\left({{r}_{1}},{{\psi}_{1}}\right)$ of the system in this case. For ${{\lambda}_{c1}}<\lambda \leqslant {{\lambda}_{c2}}$, two new solutions u₂ and u₃ corresponding to the fixed points ${{P}_{2}}\left({{r}_{2}},{{\psi}_{2}}\right)$ and ${{P}_{3}}\left({{r}_{3}},{{\psi}_{3}}\right)$ occur. Then there are three fixed points u₁, u₂, u₃ coexisting in this case (here ${{u}_{1}}<{{u}_{2}}<{{u}_{3}}$). For $\lambda >{{\lambda}_{c2}}$, the solutions u₁ and u₂ disappear, and then u₃ (P₃) become the only solution (fixed point). To obtain the stability of these fixed points, we take the Jacobian matrixes of equations (26) and (27) and compute the corresponding eigenvalues at these points. For a fixed point $P=\left(r,\psi \right)$, the Jacobian matrix is

$$\begin{eqnarray}J(P)=\left(\begin{array}{c} \frac{\Delta \left(1+{{r}^{2}}\right)}{{{r}^{2}}-1} & \frac{\left(\beta -1\right)\left({{r}^{3}}-r\right)}{\left(2\beta +1\right){{r}^{2}}+1} \\ \frac{\left(1-\beta \right)\left[\left(2\beta +1\right){{r}^{2}}-1\right]}{\left(2\beta +1\right){{r}^{3}}+r} & \frac{\Delta \left[\left(2\beta +1\right){{r}^{2}}+1\right]}{{{r}^{2}}-1} \end{array}\right).\end{eqnarray} \tag{ 34 }$$

Zoom In Zoom Out Reset image size

The eigenvalues of J(P) satisfy the equation ${{x}^{2}}-\text{tr}\left(J(P)\right)x+\det \left(J(P)\right)=0$. Here $\text{tr}\left(J(P)\right)$ and $\det (J(P)$ represent the trace and determinant of J(P), respectively. Obviously, $\text{tr}\left(J(P)\right)<0$ for all 0  <  r  <  1. With some calculation, we have $\det (J)=\left(1-{{r}^{2}}\right)\left[\left(2\beta +1\right){{r}^{2}}+1\right]{{f}^{\prime}}\left({{r}^{2}}\right)$. For ${{P}_{1}}=\left({{r}_{1}},{{\psi}_{1}}\right)$, we have ${{f}^{\prime}}\left(r_{1}^{2}\right)={{f}^{\prime}}\left({{u}_{1}}\right)>0$ and then $\det \left(J\left({{P}_{1}}\right)\right)>0$. So the eigenvalues of J(P₁) that satisfy the equation ${{x}^{2}}-\text{tr}\left(J\left({{P}_{1}}\right)\right)x+\det \left(J\left({{P}_{1}}\right)\right)=0$ are both negative, showing that P₁ is actually a stable node point of the system. For P₂ and P₃, we have ${{f}^{\prime}}\left(r_{2}^{2}\right)={{f}^{\prime}}\left({{u}_{2}}\right)<0$ and ${{f}^{\prime}}\left(r_{3}^{2}\right)={{f}^{\prime}}\left({{u}_{3}}\right)>0$ thus indicating $\det \left(J\left({{P}_{2}}\right)\right)<0$ and $\det \left(J\left({{P}_{3}}\right)\right)>0$. Then J(P₂) has two opposite-sign eigenvalues and J(P₃) has two negative ones, indicating that the fixed point P₂ is a saddle (unstable) point and P₃ is a stable node point.

The stability of a fixed point decides whether it is an attracted state when t tends to infinity. For $0\leqslant \lambda \leqslant {{\lambda}_{c1}}$, P₁ is the only attractor and then the parameter r with an arbitrary initial value will tend to r₁ as $t\to \infty$. For ${{\lambda}_{c1}}<\lambda <{{\lambda}_{c2}}$, two stable states P₁ and P₃ coexist. The parameter r will tend to r₁ (or r₃) when the initial state is close to P₁ (or P₃). For $\lambda \geqslant {{\lambda}_{c2}}$, the fixed point P₃ becomes the only attractor and then r will tend to r₃ with an arbitrary initial state.

This analysis actually indicates that there is a first-order phase transition of order parameter r from an incoherent state represented by P₁ to a synchronized state represented by P₃, where the values ${{\lambda}_{c1}}$ and ${{\lambda}_{c2}}$ correspond to the backward ($\lambda _{c}^{b}$) and forward ($\lambda _{c}^{f}$) critical coupling. The coupling intervals $\left(0,{{\lambda}_{c1}}\right)$, $\left({{\lambda}_{c1}},{{\lambda}_{c2}}\right)$ and $\left({{\lambda}_{c2}},\infty \right)$ correspond to the incoherent region, hysteresis region and synchronized region, respectively.

Case (ii): $\Delta \geqslant {{\Delta}_{c}}$. In this case, we have ${{f}^{\prime}}(u)>0$ for all $0\leqslant u<1$. Then f monotonically increases in the interval [0, 1) (see figure 2(b)). For a given $\lambda >0$, equation (28) has only one solution, corresponding to the unique fixed point $P\left({{r}_{0}},{{\psi}_{0}}\right)$) of equations (26) and (27). Using the same method presented in case (i), we conclude that P is a stable node point of this system. As P₀ is the only global attractor of the system, parameter r(t) will tend to r₀ as $t\rightarrow\infty$ with an arbitrary initial state. From figure 2(b), we can see that r₀ should gradually increase from 0 to 1 as λ varies from 0 to $\infty$. There is not a nonzero critical point for the transition from the incoherent state to the synchronized state. In other words, we say that the transition is continuous and the critical coupling point ${{\lambda}_{c}}$ is zero.

The above results show that the synchronization phase transitions of the system depend on the width of perturbation distribution $\Delta$. Small perturbations of frequencies ($\Delta <{{\Delta}_{c}}$) do not change the first-order synchronization transition of the original model in section 2, while large perturbations of frequencies weaken the degree-frequency correlation and then lead to a continuous synchronization transition. Besides, as $\Delta \to 0$, the model in this section described by equations (14) and (15) should approximate to the model without perturbations described by equations (1) and (2). We check this conclusion in the following. First, it is easy to check that ${{\lim}_{\Delta \to 0}}{{u}_{c2}}=\frac{1}{2\beta +1}$ and ${{\lim}_{\Delta \to 0}}{{u}_{c1}}=1$. Then we have ${{\lim}_{\Delta \to 0}}{{\lambda}_{c2}}={{\lim}_{\Delta \to 0}}2\sqrt{g\left({{u}_{c2}}\right)}=\frac{\beta -1}{\sqrt{2\beta +1}}$, which is just the forward critical coupling $\lambda _{c}^{f}$ of the model in section 2. For u_c1, we have ${{\Delta}^{2}}\sim A\left(\beta \right){{\left(1-{{u}_{c1}}\right)}^{3}}$ as $\Delta \to 0$, here $A\left(\beta \right)$ is a function of β. Then we have ${{\lim}_{\Delta \to 0}}{{\lambda}_{c1}}={{\lim}_{\Delta \to 0}}2\sqrt{g\left({{u}_{c1}}\right)}=\frac{\beta -1}{\beta +1}$, which is actually the backward critical coupling $\lambda _{c}^{b}$ of the model described by equations (1) and (2). Moreover, the critical value of r as $\Delta \to 0$ is ${{\lim}_{\Delta \to 0}}{{r}_{c1}}=1$ and ${{\lim}_{\Delta \to 0}}{{r}_{c2}}=\frac{1}{\sqrt{2\beta +1}}$, which are just the same as the results of the model in section 2. The agreement of these two models when perturbations tend to zero partly confirms the validity of the method we used.

All of the results described above are obtained by analyzing the reduced ODE system governed by equations (26) and (27), and are therefore subject to the restrictions described therein. It is reasonable to ask if these results agree with the dynamics of the original system given in equations (14) and (15). To check this, a series of direct simulations of equations (14) and (15) using the fourth-order Runge–Kutta numerical integration are performed in the following.

Figure 3 presents the time evolutions of parameter r for different values of $\Delta$ and λ. The blue solid lines stand for the numerical simulations from the original model given in equations (14) and (15) with $N=10\,000$, while the red dashed lines are the analytical solutions obtained from equations (26) and (27). We fix $\beta =10$ and then we have ${{\Delta}_{c}}=1.44$. We choose two typical values of $\Delta$, 0.5 and 2.5 which correspond to case (i) and case (ii), respectively. Two different values of λ are chosen for each $\Delta$ to stand for the weak coupling case and the strong coupling case. As we can see, the simulation results and the analytical solutions fit extremely well in all the four pictures, indicating that the reduced ODE system governed by equations (26) and (27) can quite precisely describe the time evolution of the order parameter r as N is large.

Zoom In Zoom Out Reset image size

Figure 4 shows the phase transition diagrams of the order parameter r. The parameter β here is set to be 10. The upper two figures are the simulation results of r versus coupling λ for $\Delta =0.5$ (figure 4(a)) and $\Delta =1.5$ (figure 4(b)). For comparison, the lower two figures present the analytical solutions of r versus λ solved from equation (28). Figure 4(a) clearly shows a sharp jump of r and a hysteresis region between the forward and backward diagrams, which verify the existence of the first-order transition. The backward and forward processes in figure 4(b) coincide completely, showing that phase transition is actually continuous. Analytical results in figure 4(c) and (d) agree with corresponding simulations in figures 4(a) and (b). The backward and forward critical coupling obtained from the simulation in figure 4(a) is $\lambda _{c}^{b}=1.56$ and $\lambda _{c}^{f}=1.96$. They are quite close to the analytical results presented in figure 4(c), where $\lambda _{c}^{b}=1.558$ and $\lambda _{c}^{f}=1.976$.

Zoom In Zoom Out Reset image size

For $\Delta <{{\Delta}_{c}}$, it has been shown that there is a first-order phase transition which is labeled by the two inconsistent critical couplings $\lambda _{c}^{b}$, $\lambda _{c}^{f}$ and the existence of the hysteresis region. We are interested in how the value of critical coupling and the width of hysteresis depend on $\Delta$. To do this, we solve equation (32) for a fixed β and each $\Delta <{{\Delta}_{c}}$ to get the roots u_c1 and u_c2. Then we have $\lambda _{c}^{b}=2\sqrt{g\left({{u}_{c1}}\right)}$ and $\lambda _{c}^{f}=2\sqrt{g\left({{u}_{c2}}\right)}$. Figure 5(a) presents the values of $\lambda _{c}^{b}$ and $\lambda _{c}^{f}$ for $\beta =10$ when $\Delta$ varies from 0 to ${{\Delta}_{c}}=1.44$. We can see $\lambda _{c}^{f}$ slowly increases with $\Delta$ while $\lambda _{c}^{b}$ increases much faster until, at $\Delta ={{\Delta}_{c}}=1.44$, the two lines overlap. The $\Delta -\lambda$ plane is then divided into three regions—the synchronized region (blue), hysteresis region (yellow) and incoherent region (purple). Figure 5(b) shows the hysteresis width $\lambda _{c}^{f}-\lambda _{c}^{b}$ as a function of $\Delta$. Here we can see it shows a descending trend starting from a certain level at $\Delta =0$ and ending with $\lambda _{c}^{f}-\lambda _{c}^{b}=0$ at $\Delta ={{\Delta}_{c}}$.

Zoom In Zoom Out Reset image size

Finally, we want to show how the parameter $\Delta$ and β codetermine the order of the phase transition of the system. We plot ${{\Delta}_{c}}$ as a function of β from 1 to 100 in figure 6. Then the phase diagram is divided into two parts. In the upper part where $\Delta \geqslant {{\Delta}_{c}}$, the phase transition is continuous while in the lower part where $\Delta <{{\Delta}_{c}}$, transition is first-order (explosive). We can see that the boundary $\Delta ={{\Delta}_{c}}\left(\beta \right)$ looks like a straight line. In fact, when $\beta \gg 1$, ${{u}_{c}}\approx \frac{2}{2\beta +1}$, and then by equation (32) we have ${{\Delta}_{c}}\approx \frac{1}{3\sqrt{3}}\left(\beta -1\right)$. We plot $\Delta =\frac{1}{3\sqrt{3}}\left(\beta -1\right)$ as a dashed line in this figure and we can see that the two lines are quite close to each other.

Zoom In Zoom Out Reset image size

In this paper, we have analyzed ES on star graphs. Using the Ott–Antonsen ansatz, we first established a two-dimensional reduced system that can describe the evolution of the global behavior of the oscillators. We studied the dynamical behavior of this reduced system and showed that the phase transition from incoherence to synchronization is actually discontinuous. We also presented direct simulations of the original model and compared them to analytical results solved from the reduced system, and showed they have a good agreement. Then, we extended our model studied to a more general case where there are random perturbations to the frequencies of leaf nodes. We studied how the type of synchronization transition is decided by the perturbation strength $\Delta$ and the frequency of hub β. The results showed that for small perturbations the synchronization phase transition is still discontinuous while for large perturbations the transition becomes continuous.

Our works have completely solved the ES of Kuramoto oscillators on star graphs and have uncovered the intrinsic mathematical mechanism behind the first-order transition. For the SF network with tree-like structure, the modified model in section 3 can roughly describe the dynamic of the oscillators of star components. Then the results in this section can partly explain the mechanism of ES in this kind of network. Although it is still hard to completely solve the puzzle of ES on general SF networks, our work is still of important value in understanding the explosive phase transitions in large ensembles of coupled units, and moreover, it also points out some feasible directions and provides inspiration for a deeper investigation in the future.

This work is supported by the Major Program of National Natural Science Foundation of China (11290141), NSFC (11201018) and the international cooperation project no. 2010DFR00700.