Dynamics in the Sakaguchi-Kuramoto model with bimodal frequency distribution

In this work, we study the Sakaguchi-Kuramoto model with natural frequency following a bimodal distribution. By using Ott-Antonsen ansatz, we reduce the globally coupled phase oscillators to low dimensional coupled ordinary differential equations. For symmetrical bimodal frequency distribution, we analyze the stabilities of the incoherent state and different partial synchronous states. Different types of bifurcations are identified and the effect of the phase lag on the dynamics is investigated. For asymmetrical bimodal frequency distribution, we observe the revival of the incoherent state, and then the conditions for the revival are specified.


Introduction
Collective behaviors emerged out of a large number of interacting units are common in nature. As one type of collective behavior characterizing the phase coherence in nonidentical units, synchronization is well recognized in various systems such as fireflies flashing in unision [1,2], applauding persons in a large audience [3], pedestrians [4,5], and others [6]. Kuramoto model (KM) is the paradigmatic model in the field of synchronization [7,8]. There are two key simplifications in the original KM, which renders the analytical treatments to be possible. Firstly, each unit is treated as a phase oscillator, which is valid for the weak coupling situation where the amplitude information of each unit is inessential to the collective behaviors. The dynamics of each phase oscillator is solely determined by its natural frequency and in turn the frequencies of all oscillators are drawn from a prescribed frequency distribution function g(ω). Secondly, the coupling between units is assumed to be a global one and takes the form of a sinusoidal function with the same strength K. The coupling strength together with the frequency distribution determine the dynamics of KM.
Previously, KM has been intensively investigated. A variety of its generalizations have been proposed and many interesting phenomena have been observed. The repulsive interaction among oscillators (K < 0) may be introduced to KM. Tsimring et al. [9] found that KM with repulsive interaction fails to synchronize. Hong and Strogatz [10,11] treated the coupling strength as an oscillator's ability to response to the mean field and found π synchronous states and novel time-dependent traveling wave synchronous states in the presence of both repulsive and attractive interaction. Yuan et al. further considered the π synchronous state in the presence of correlation between the conformists/contrarians and the natural frequencies of oscillators [12]. Zhang et al. [13] introduced frequency-weighted coupling to KM and found explosive synchronization and chimera-like states. KM has also been extended to complex networks where network topology can affect the synchronization transition. In [14], the authors assigned the natural frequencies of phase oscillators to be the degrees of the nodes they locate on network and found explosive synchronization transition. Recently, KM with higher-order interaction such as biharmonic interaction has drawn some attentions in which infinitely many stable partial synchronous states and a continuum of abrupt desynchronization transition have been identified [15]. The shear is a crucial nonlinear ingredient for complex behaviors in coupled systems [16]. Time delay was also investigated [17][18][19], and small time delay can be approximated by a phase lag parameter β. Along this line, the the phase lag β is introduced into the coupling function as Ksin(θ j − θ i + β) so that KM is generalized to Sakaguchi-Kuramoto model (SKM) [20] and the synchronous dynamics has been investigated [16,[21][22][23][24].
Actually, the original KM is concise enough to display rich dynamics by taking proper frequency distribution g(ω). It has been theoretically shown that the transition to synchronization occurs at K c = 2/[πg(0)] [25] for even and unimodal g(ω). Above K c , the incoherent state yields to a stationary partial synchronous state. For asymmetrical unimodal g(ω), the partial synchronous states are always time-dependent [26]. When g(ω) becomes a bimodal one, increasing coupling strength always first leads to a standing wave state, in which two synchronous clusters of oscillators oscillate at opposite mean frequencies and, then, to traveling wave states, in which synchronous oscillators rotate at the same frequency [27]. Bimodal frequency distributions in the KM were already investigated at different levels [28][29][30][31][32][33][34], and trimodal frequency distribution were also studied [30]. Under proper parameters, KM with bimodal distribution gives rise to discontinuous transitions cross different dynamical states. Martens et al. [33] studied KM with bimodal natural frequency distribution consisting of two equally weighted Lorentzians, and they derived the system's stability diagram. They found three states depending on the parameters and initial conditions, incoherent state, partial synchronous state, and standing wave synchronous states. They also presented analytical results for the bifurcation boundaries between these states. Omel'chenko and colleagues [35] studied SKM with g(ω) being a superposition of two unimodal frequency distributions with the same mean frequency. They found a nonuniversal synchronization transition in which the incoherent state may be revived at stronger coupling strength after it yields to partial synchronous state at K c . Asymmetry has also been studied recently [26,29,31]. For more complicate frequency distribution such as a trimodal one, KM may display collective chaos through a cascade of period-doubling bifurcations [36].
In this work, we study SKM with bimodal natural frequencies distribution. As a natural extension of Ref. [33], the phase-lag parameter β is introduced into the model. The paper is organized as follows. In section 2, we present the model and reduce the coupled phase oscillators to a lowdimensional coupled ordinary differential equations. In section 3, we first study the synchronous dynamics in the model with symmetrical bimodal frequency distribution with an emphasis on the effects of the phase lag. Different dynamical states are analyzed and different types of bifurcations are identified. Then we consider SKM with asymmetrical bimodal frequency distribution. We study the revival phenomenon of the incoherent state and investigate the dependence of revival of the incoherent state on parameters. Summary is made in the last section.

Materials and methods
We consider N phase oscillators with global coupling and the motion equation follows with θ i the phase of oscillator i and K is the global coupling strength. β is the phase lag parameter resulting in rich interesting dynamical phenomena and the model reduces to the original KM at β = 0. ω i is the natural frequency of oscillator i, which is chosen randomly from a probability distribution g(ω). In this work, we assume that the frequency distribution g(ω) takes the form with p 1 + p 2 = 1 and ω 1 = −ω 2 = ω 0 . The parameters Δ 1,2 measure the heterogeneity of oscillators in their natural frequencies. Generally, both the heterogeneity parameter Δ and the phase lag β have strong effects on the synchronous dynamics. However, these two parameters impact on the collective dynamics in different way. Δ is used to measure the fraction of oscillators to be in synchronization. Large Δ always suggests small fraction of phase oscillators to be in synchronization. In contrast, β measures the phase mismatch between the synchronous oscillators and the mean field. Sufficiently large β pushes synchronous phase oscillators to be in antiphase with the mean field, which downgrades synchronization and tends to destroy the coherence in population. Recent work points out that incoherent state may be revival at proper choice of β [35], which suggests the non-monotonic effects of β on the coherence in population. The synchronous dynamics in the model (1) is measured by the complex order parameter, defined as Z ¼ Re iY ¼ 1 N S j e iy j . |Z| = 0 suggests the incoherent state and, otherwise, a synchronous state. Using the order parameter, Eq (1) is reformulated as To study the dynamics, we consider the thermodynamic limit (N ! 1) where Eq (1) can be written in a continuous formulation in terms of a probability density f(θ, ω, t), defined as the fraction of oscillators with natural frequency between ω and ω + dω and phase between θ and θ + dθ at time t, which satisfies the normalization condition The probability density evolves following the continuity equation The order parameter Z in the continuous formalism is reformulated as Since the probability density is periodic in θ, it can be expanded in Fourier series as with c.c. the complex conjugate of the previous term. Ott and Antonsen proposed an ansatz (OA ansatz) [37] that the coefficients f n (ω, t) obey f n (ω, t) = [α(ω, t)] n . Substituting Eq (4) with the ansatz, we obtain For the natural frequency distribution Eq (2), the order parameter Z becomes where we denote z 1, Then the synchronization in the model (1) is characterized by the sub-order parameters z 1,2 (t). The evolution of z i (i = 1, 2) follows Furthermore, we let a j ¼ z � j ¼ r j e À i� j ðj ¼ 1; 2Þ and introduce ψ = ϕ 1 − ϕ 2 . Then substituting them into Eq (11), we have The presence of the phase lag in the model breaks the symmetry between r 1 and r 2 even when Δ 1 = Δ 2 and p 1 = p 2 . Eq (12) consisting of three coupled ordinary differential equations is equivalent to the model (1, 2) and, therefore, the dynamics of the model (1, 2) may be reflected by r 1 , r 2 , and ψ. To be mentioned, the partial synchronous states in the model (1) [or the reduced model (12)] are always time-dependent, periodic or quasiperiodic, for nonzero β. In the reduced model (12), these time-dependent synchronous states are reduced to equilibria or periodic solutions by considering the model in a rotating frame characterizing the time-dependent ϕ 1 . In the following, we claim a solution to be an equilibrium or periodic one according to its behavior in the reduced model (12).

Symmetric frequency distribution
We first consider the symmetric frequency distribution where p 1 = p 2 = 0.5, and Δ 1 = Δ 2 = Δ. We set the coupling strength K = 4 and investigate the effect of the phase lag β on the model dynamics.
We start with the reduced model Eq (11) and investigate the stability of the incoherent state. The incoherent state is defined by z 1 = z 2 = 0. Supposing that the evolution of perturbations to the incoherent state follows δz 1,2 * e λt and substituting them into Eq (11), we may have l 1;2 ¼ e ib À D 0 � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi e 2ib À o 02 0 p ð13Þ For convenience, we assume Re(λ 1 )>Re(λ 2 ). When Re(λ 1 ) becomes positive, the incoherent state becomes unstable. Beyond the bifurcation, Eq (12) gives rise to two new stable equilibria except for the unstable incoherent state, r 1,2 > 0 in one equilibrium, and r 1,2 < 0 in the other which is unrealistic and should be discarded. Therefore, the incoherent state undergoes a supercritical Pitchfork bifurcation when Re(λ 1 ) crosses zero (we denoted it as PB1). Interestingly, when Re(λ 2 ) crosses zero, it induces another pitchfork bifurcation (denoted as PB2) in which two newborn equilibria are unstable and one of them is unrealistic. The pitchfork bifurcations involving the incoherent state occur at the critical curves described by When β = 0, the critical curves (14) are reduced to a semicircle D 0 ¼ 1 � ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 1 À o 02 0 p for o 0 0 < 1, which is related to pitchfork bifurcation, and a line Δ 0 = 1 for o 0 0 > 1 which is related to Hopf bifurcation [33]. Increasing β from zero, the stability regime of the incoherent state shrinks in the plane of Δ 0 and o 0 0 . Then we consider model dynamics by focusing on Eq (12). The equilibria to Eq (12) represent the partial synchronous states and their stabilities can be analyzed by the linear stability method. For β = 0, the partial synchronous state can be acquired rigorously by setting r 1 = r 2 [33]. However, for partial synchronous state, r 1 = r 2 is always not held as β 6 ¼ 0. the equilibria to Eq (12) are obtained by numerical methods and their stabilities are determined by the eigenvalues of the Jacobian matrices at them. To illustrate, we consider the bifurcation diagrams along three parameter paths by setting β = 0.1 and K = 4. Firstly, we consider the parameter path with o 0 0 from 0.4 to 2 at Δ 0 = 0.4. The bifurcation diagrams are presented in Fig 1(a) where r 1 and r 2 are plotted, respectively. Besides the incoherent state which is always unstable along this parameter path, there are at most four equilibria denoted as FP ð1Þ i (i = 1, 2, 3, 4). The eigenvalues of the corresponding Jacobian matrices at these equilibria are plotted in Fig 1(b) and 1(c). As shown, the equilibria FP ð1Þ 1 is stable until o 0 0 ' 1:62 at which it collides with a saddle FP ð1Þ 2 and gives rise to a limit cycle, a standing wave synchronous state, through a SNIPER bifurcation (saddle node infinite period bifurcation). The equilibrium FP ð1Þ 2 is a saddle with a one-dimensional unstable manifold, which is born at o 0 0 ' 1:40 with another saddle FP 1 3 owning a two-dimensional unstable manifold through a saddle-node bifurcation (denoted as SN2). Shortly after SN2, the unstable FP ð1Þ 3 is turned into a saddle-focus. The equilibrium FP ð1Þ 4 has a pair of complex conjugate eigenvalues whose real parts are positive and is an unstable saddlefocus, which is produced by the pitchfork bifurcation (denoted as PB2) of the incoherence state at around o 0 0 ' 0:8 according to Eq (14). Along this parameter path, there are two stable synchronous states, one is represented by FP 1 1 before o 0 0 ¼ 1:62 and the other is represented by a limit cycle [the solid curves in Fig 1(a)].
Secondly, we consider the parameter path with o 0 0 from 1.15 to 1.22 at Δ 0 = 0.95. The bifurcation diagram is presented in Fig 2(a) and the eigenvalues for all equilibria are presented in Fig 2(c) and 2(e). As shown, we find two stable equilibria (FP ð2Þ 1 ,FP ð2Þ 2 ) and one unstable equilibrium (FP 3 ). The bifurcations at which FP ð2Þ 1 and FP ð2Þ 3 annihilate with each other and at which FP ð2Þ 2 and FP ð2Þ 3 are born in pair belong to the saddle-node bifurcation (one is denoted as SN1 and the other is SN2). Along this path, the bistability between FP ð2Þ 1 and FP ð2Þ 2 exists in a range of o 0 0 . The third parameter path is chosen against Δ 0 at o 0 0 ¼ 1:5, which is presented in Fig 2(b), 2 (d) and 2(f). There are two equilibria, FP ð3Þ 1;2 , and a stable periodic solution. FP ð3Þ 1 is a focus, which changes from an unstable to a stable one by colliding with the limit cycle at Δ 0 ' 0.904 through a Hopf bifurcation (denoted as HB). Furthermore, the stable FP ð3Þ 1 disappears at Δ 0 ' 1.09 by turning the unstable incoherent state to being stable one through a pitchfork bifurcation (PB1). The unstable equilibrium FP ð3Þ 2 is always unstable, which results from a pitchfork bifurcation (PB2) of the unstable incoherent state when the real part of its second eigenvalue Re(λ 2 ) crosses zero [see Eq (13)].
Using the above analysis, the phase diagrams in the plane of Δ 0 and o 0 0 at β = 0 and β = 0.1 are presented in Fig 3(a) and 3(b), respectively. Actually, the results at β = 0 have been thoroughly explored [33] and there is only two minor modifications in Fig 3(a). Firstly, we point states. In top panels, red, blue, wine, and dark green symbols are for partial synchronous states FP ð1Þ 1 , FP ð1Þ 2 , FP ð1Þ 3 , and FP ð1Þ 4 , respectively. Thick black and dark green lines represent the maximum and minimum values of r 1 and r 2 for stable standing wave synchronous states. In middle and bottom panels, from left to right, real and imaginary parts of the eigenvalues λ for partial synchronous states from FP ð1Þ 1 to FP ð1Þ 4 are displayed. Squares, circles, and triangles denote eigenvalues λ 1 , λ 2 , and λ 3 , respectively.
https://doi.org/10.1371/journal.pone.0243196.g001 out that the incoherent state loses its stability through a pitchfork bifurcation at low o 0 0 instead of a transcritical bifurcation claimed in Refs. [33,38], which is similar to Refs. [39]. We find that the bifurcations are similar to those in KM with trimodal frequency distribution [30]. Secondly, we include in the phase diagram one more saddle-node bifurcation (SN2) which involves the birth of a pair of unstable saddles. The saddles arising from SN2 were not reported in Ref. [33] in which the authors concerns with the long-term dynamics at β = 0. Interestingly, we find that one of these two saddles becomes stable as β 6 ¼ 0. To be stressed, unstable solutions have no effects on the long-term dynamics of the model dynamics. However, the existence of unstable solutions greatly shapes the topological structure of the underlying phase space and has strong impacts on the transient dynamics of the model. Moreover, under certain conditions, unstable solutions might become stable with the change of parameter and, then, . Solid (open) data points represent stable (unstable) states. In (a), red, blue, and wine lines are for partial synchronous states FP ð2Þ 1 , FP ð2Þ 2 , and FP ð2Þ 3 , respectively. In (b), wine and dark green lines are for partial synchronous states FP ð3Þ 1 and FP ð3Þ 2 , respectively. Thick black and dark green lines refer to the standing wave synchronous state. In the panels from (c) to (f), squares, circles, and triangles denote the real and imaginary parts of eigenvalues λ 1 , λ 2 , and λ 3 , respectively. The inset of (d) shows that HB occurs at a lower Δ 0 than PB2. Note that the incoherent state changes its stability across the pitchfork bifurcation (PB1).
https://doi.org/10.1371/journal.pone.0243196.g002 take effects on the long-term dynamics of the model. Therefore, in the perspective of stability analysis, the exploration of unstable solutions is still necessary.
In Fig 3(b), there are two pitchfork bifurcations involving the incoherence (PB1,PB2), two saddle-node bifurcations involving partial synchronous states (SN1,SN2), and three bifurcations involving limit cycle synchronous states (Hopf bifurcation, homoclinic bifurcation, and SNIPER). The critical curves relating to these bifurcations divide the parameter plane of Δ/4K and ω 0 /4K into several domains. And the phase diagram in Fig 3(b) shows that FP ð1Þ 1 and FP ð2Þ 1 are the same type of solutions while FP ð1Þ 4 and FP ð3Þ 2 are the same type of solutions. The typical evolutions on the plane of r 1 and r 2 from (or towards) the solutions in these different domains are presented in the insets. Compared with Fig 3(a), there are several unique features in Fig 3(b) to be addressed. At β = 0, the two PBs form a continuous semicircle. However, these two PBs become two separated curves. Furthermore, the Hopf bifurcation underlies the transition between the stable incoherent state and the stable limit cycle at β = 0. However, at nonzero β, the Hopf bifurcation occurs between the stable partial synchronous state and the stable limit cycle. In addition, the Hopf bifurcation stays much close to PB2 of the incoherent state. At β = 0, there exists a domain in which the incoherent state coexists with a partial synchronous state. However, no coexistence between the incoherent state and any partial synchronous states at β = 0.1, as shown in Fig 3(b). Instead, there exists the coexistence between two partial synchronous states in the domain enclosed by two saddle-node bifurcations (SN1 and SN2) and HB. As shown in Fig 3(b), there exists a Takens-Bogdanov bifurcation (denoted as TB) where Hopf bifurcation, homoclinic bifurcation (denoted as HC), and saddle-node bifurcation merge. Interestingly, a pair of stable and unstable synchronous states are born at SN2 above TB while a pair of unstable synchronous states occur at SN2 below TB. In addition, SN1 gradually merges with HC to become SNIPER.
In KM with unimodal frequency distribution, increasing the phase lag β always downgrades the coherence in population and, when β = π/2, the critical coupling strength K for the onset of Line color codes: black and red for two pitchfork bifurcations, PB1 and PB2, respectively; green and pink for two saddle-node bifurcations, SN1 and SN2, respectively; blue for HB (Hopf bifurcation); cyan for HC (homoclinic bifurcation). Acronyms: SNIPER for saddle node infinite period; CP for cusp point of SN1 and SN2; TB for Takens-Bogdanov point. To present the topological structure of the phase space in different phase domains, we plot the phase portraits on the (r 1 , r 2 ) plane in several insets with the parameters chosen from different phase domains. The dashed arrows pointing to insets refer to the phase domain represented by the insets. In each inset, several phase portraits (wiggly lines) are plotted with arrows representing the evolution from or towards the solutions in Eq (12). In these insets, solid (open) dots represent stable (unstable) partial synchronous states, while the dark yellow curves represent stable standing wave partial synchronous state denoted by L. The solutions in the same color in different insets are the same solution. K = 4.  Fig 4(c), complicated structure in the phase diagram appears, for example the bistability between different partial synchronous states, the bistability between the partial synchronous states and the standing wave states, and the existence of two Takens-Bagdanov bifurcations.

Asymmetric frequency distribution
Omel'chenko and colleagues have found an interesting phenomenon in a SKM with g(ω) being a superposition of two unimodal frequency distributions with the same mean frequency where the incoherent state may be revived at stronger coupling strength [35]. Liu and colleagues found the same phenomenon in a SKM with g(ω 0 ) being the superposition of two bimodal frequency distributions [40]. Here we show the revival of the incoherent state for asymmetrical bimodal frequency distribution and provide the conditions for better observing revival of the incoherent state. We consider the stability diagrams of the incoherent state on different parameter planes where the stability of the incoherent state is calculated based on Eq (11). With reference to the process of reaching Eq (13), we may have ffi ffi ffi 2 p ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi c 1 þ ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi To be mentioned, ω 0 may be negative when the peak frequency ω 1 is less than the peak frequency ω 2 . Positive and negative ω 0 may exert different impacts on the model dynamics due to the asymmetrical bimodal frequency distribution. Incoherent states change stability with changing parameters at hopf bifurcation or pitchfork bifurcation [25,26,29,31,35]. We can get the bifurcation curves in Fig 5 from Eq (15), in which the more general conditions are considered analytically. Fig 5(a) shows the results on the plane of K and p 1 . For p 1 = 0 (K c ' 0.019), the incoherent state becomes unstable at sufficient weak coupling strength while strong coupling strength is required for p 1 = 1 (K c ' 0.965), which can be seen from Eq (11). Between these two extreme situation, there exists a domain at around p 2 (0.52, 0.86) in which the revival of the incoherent state appears. The stability diagram on the plane of K and Δ 2 /Δ 1 with fixed Δ 1 shows that the revival of the incoherent state requires sufficiently small Δ 2 /Δ 1 and it becomes the most prominent at Δ 2 /Δ 1 = 0 [see Fig 5(b)]. If we measure the revival phenomenon of the incoherent state by the range of the coupling strength K, Fig 5(c) indicates that the superposition of two unimodal distributions with the same mean frequency is not the best candidate for realization of the revival phenomenon. Weak mismatch between the center frequencies of the two unimodal distribution is the optimal for the revival of the incoherent state. Finally, Fig 5(d) suggests that revival of the incoherent state occurs only in SKM with proper phase lag β. We also plot the critical curves, K = 2Δ 1 /cosβ and K = 2Δ 2 /cosβ, for the incoherent state when p 1 = 1 and p 1 = 0. It is interesting to find that these two curves may approximate part of the boundary of the stable incoherent state, which suggests that the revival of the incoherent state is somehow induced by the competition between these two instability mechanisms. To summarize, the revival of the incoherent state studied here requires some conditions. Firstly, the frequency distribution is composed of two unimodal ones and the sufficiently low ratio of their widths is required for the revival of the incoherent state. Secondly, that the fraction of oscillators with the natural frequency from the fat peak in the population is higher than that from the thin peak is required for the revival of the incoherent state. Thirdly, proper choice of β is required. These conditions are similar to those reported in the previous work [35]. Different from the work [35] where the two unimodal distributions share the same central frequency and the frequency distribution is a symmetrical one, the frequency distribution here is a bimodal one and no symmetry on it is required. The results in Fig 5 suggest that the revival of the incoherent state could be a rather popular phenomenon.

Conclusion
In conclusion, we have investigated the globally coupled Sakaguchi-Kuramoto model with bimodal natural frequency distributions. By using Ott-Antonsen ansatz for dimension reduction, we reduce the coupled phase oscillators to a low dimensional coupled ordinary equations. For symmetrical bimodal frequency distribution, we analyze the linear stabilities of the incoherent state and partial synchronous states and identify different types of bifurcations between different dynamical states. Especially, the impacts of the phase lag β on the model dynamics are studied. For example, nonzero β greatly modifies the topological structure of the phase space and unfolds certain bifurcations degenerated at β = 0. More importantly, β impacts on synchronous dynamics in the population in a non-monotonic way. The bifurcation may be unfolded by nonzero. We also study the revival of the incoherent state for the model with asymmetrical bimodal frequency distributions and the conditions for better observing the phenomenon are proposed.