Generalization of the Kuramoto model to the Winfree model by a symmetry breaking coupling

TL;DR

This paper generalizes the Kuramoto model by introducing a symmetry-breaking coupling, creating a variant of the Winfree model, and analyzes the resulting phase transitions and bifurcations across different frequency distributions.

Contribution

The authors develop a generalized Kuramoto model with explicit symmetry breaking, analyze its phase diagram, and derive low-dimensional equations capturing its bifurcation structure.

Findings

Symmetry-breaking coupling influences the stability of synchronized states.

Bistable regions depend on frequency distribution asymmetry.

Analytical bifurcation curves match simulation results.

Abstract

We construct a nontrivial generalization of the paradigmatic Kuramoto model by using an additional coupling term that explicitly breaks its rotational symmetry resulting in a variant of the Winfree Model. Consequently, we observe the characteristic features of the phase diagrams of both the Kuramoto model and the Winfree model depending on the degree of the symmetry breaking coupling strength for unimodal frequency distribution. The phase diagrams of both the Kuramoto and the Winfree models resemble each other for symmetric bimodal frequency distribution for a range of the symmetry breaking coupling strength except for region shift and difference in the degree of spread of the macroscopic dynamical states and bistable regions. The dynamical transitions in the bistable states are characterized by an abrupt (first-order) transition in both the forward and reverse traces. For asymmetric bimodal frequency distribution, the onset of bistable regions depends on the degree of the asymmetry. Large degree of the symmetry breaking coupling strength promotes the synchronized stationary state, while a large degree of heterogeneity, proportional to the separation between the two central frequencies, facilitates the spread of the incoherent and standing wave states in the phase diagram for a low strength of the symmetry breaking coupling. We deduce the low-dimensional equations of motion for the complex order parameters using the Ott-Antonsen ansatz for both unimodal and bimodal frequency distributions. We also deduce the Hopf, pitchfork, and saddle-node bifurcation curves from the evolution equations for the complex order parameters mediating the dynamical transitions. Simulation results of the original discrete set of equations of the generalized Kuramoto model agree well with the analytical bifurcation curves.

Equations

$$\dot{{\theta }_i}$$

$$\dot{{\theta }_i}={\omega }_i+Q({\theta }_i)\sum _{j=1}^{N}P({\theta }_j),$$

$$g(\omega )$$

$${\int }_0^{2\pi }\mathrm{d}\theta f(\theta ,\omega ,t)=g(\omega )\forall \omega .$$

$$\frac{\partial f}{\partial t}+\frac{\partial (fv)}{\partial \theta }=0,$$

$$v(\theta ,\omega ,t)=\omega +\frac{\epsilon }{2i}[(z{e}^{-i\theta }-z^{\star }{e}^{i\theta })+q(z{e}^{i\theta }-z^{\star }{e}^{-i\theta })],$$

$$z={\int }_{-\infty }^{\infty }g(\omega ){\int }_0^{2\pi }f(\theta ,\omega ,t)e^{i\theta }d\theta d\omega .$$

$$f(\theta ,\omega ,t)=\frac{g(\omega )}{2\pi }\left[1+\sum _{n=1}^{\infty }\left({\alpha }_n(\omega ,t)e^{in\theta }+c.c.\right)\right],$$

$${\alpha }_n(\omega ,t)={\left[\alpha (\omega ,t)\right]}^n.$$

$$\frac{\partial \alpha }{\partial t}+i\omega \alpha +\frac{{\epsilon }_1}{2}\left[(z{\alpha }^2-z^{\star })+q(z-z^{\star }{\alpha }^2)\right],$$

$$z^{\star }={\int }_{-\infty }^{\infty }\alpha (t,\omega )g(\omega )d\omega .$$

$$z(t)=a^{\star }({\omega }_0-i\gamma ,t).$$

$$\frac{\mathrm{d}r}{\mathrm{d}t}$$

$$\frac{\mathrm{d}\psi }{\mathrm{d}t}$$

$$~{}~{}M\equiv\begin{bmatrix}-\gamma+\frac{\varepsilon}{2}\big{[}1-q\big{]}&-\omega_{0}\\ \\ ~{}\omega_{0}&-\gamma+\frac{\varepsilon}{2}\big{[}1+q\big{]}\\ \end{bmatrix}.$$

$$\frac{\epsilon q}{2}\cos (2\psi )$$

$$\frac{\epsilon q}{2}\sin (2\psi )$$

$$\frac{{\epsilon }^2{q}^2}{4}$$

$$\tan (2\psi )$$

$$z(t)=\frac{1}{2}[z_1(t)+z_2(t)],$$

$$z_{1,2}(t)={\alpha }^{\star }(\pm {\omega }_0-i{\gamma }_{1,2},t).$$

$${\dot{z}}_1=$$

$$+q((z_1+z_2)z_1^2-(z_1^{\star }+z_2^{\star }))],$$

$${\dot{z}}_2=$$

$$+q((z_1+z_2)z_2^2-(z_1^{\star }+z_2^{\star }))].$$

$$\frac{\mathrm{d}{r}_1}{\mathrm{d}t}$$

$$\displaystyle+q((r_{1}^{2}-1)(r_{1}\cos(2\psi_{1})+r_{2}\cos(\psi_{2}+\psi_{1})))\big{]},~{}~{}~{}$$

$$\frac{\mathrm{d}{\psi }_1}{\mathrm{d}t}$$

$$\displaystyle+q(r_{1}\sin(2\psi_{1})+r_{2}\sin(\psi_{2}+\psi_{1}))\big{]}.$$

$$\frac{\mathrm{d}{r}_2}{\mathrm{d}t}$$

$$\displaystyle+q((r_{2}^{2}-1)(r_{1}\cos(\psi_{2}+\psi_{1})+r_{2}\cos(2\psi_{2})))\big{]},~{}~{}~{}~{}$$

$$\frac{d{\psi }_2}{dt}$$

$$\displaystyle-q(r_{1}\sin(\psi_{2}+\psi_{1})+r_{2}\sin(2\psi_{2}))\big{]}.$$

$${\omega }_0^2=\frac{1}{4}(\epsilon a_1-2a_2+\sqrt{{\epsilon }^2{q}^2{a}_1-4\epsilon a_3^2+4a_3^2{a}_1}),$$

$${\omega }_0^2=$$

Full text

1 Introduction

Symmetry (translational or rotational) prevailing in the coupled dynamical networks due to the coupling geometry manifests in a wide variety of natural systems and in their intriguing macroscopic dynamical states tstp2017 . Nevertheless, symmetry breaking couplings are shown to be a source of a plethora of collective dynamical behavior that are inherent to it and are mostly inaccessible with the symmetry preserving couplings. In particular, networks of the paradigmatic Stuart-Landau oscillators with symmetry breaking coupling have been employed to unravel several collective dynamical states that mimic a variety of collective patterns observed in nature and technology. For instance, symmetry breaking coupling facilitates the transition from the homogeneous to an inhomogeneous steady states akev2013 , symmetry breaking interaction has been identified as an essential feature for the genesis of partially coherent inhomogeneous spatial patterns, namely chimera death state zamk2014 ; tb2015 ; kpvkc2015 . Multicluster oscillation death states have been observed in nonlocally coupled Stuart-Landau oscillators with symmetry breaking coupling ismk2015 . Further, the interplay of the nonisochronicity parameter and the symmetry breaking coupling is found to facilitate the onset of different variants of chimera death state such as multichimera death state and periodic chimera death states in nonlocally coupled Stuart-Landau oscillators kpvkc2016 . The effect of the symmetry breaking coupling has also been investigated on the phenomenon of reviving oscillations wei2017 . Recently, the effect of the symmetry breaking mean-field coupling on the phenomenon of the aging transition has also been investigatedConjugate couplings, a symmetry breaking coupling, have also been widely employed in the literature rkrr2007 ; amit2012 ; wei2022 . Note that the pointed out reports are only a tip of an ice-berg and not an exhaustive list of studies that employed symmetry breaking coupling using the network of the Stuart-Landau oscillators.

Despite the substantial investigations on the effect of the symmetry breaking coupling in networks of Stuart-Landau oscillators, there is a lacunae in understanding the nontrivial role of the symmetry breaking coupling in the phase only models, which indeed allows for exact analytical treatment of the macroscopic dynamical states in most cases. In particular, phase models such as Winfree and Kuramoto models, and their variants have been extensively employed in the literature to investigate the emergence of various intriguing collective dynamical states. Interaction among the phase oscillators in the Winfree model is modeled by a phase-dependent pulse function and a sensitive function. The former characterizes the mean-field, while the latter characterizes the response of the individual oscillators to the mean-field w:1 ; w:2 . Winfree model is one of its kind representing a class of pulse-coupled biological oscillators such as flashing of fire files In2 , applauding audience In7 and many more. Interaction among the phase oscillators in the Kuramoto model is modeled by the sine of difference between the phases of the oscillator and has been widely employed to investigate the emergence of spontaneous synchronization in a wide variety of biological, chemical, mechanical and physical, systems Kuramoto:1984 ; Acebron:2005 ; In1 . Examples include cardiac pacemaker In3 , Josephson junction arrays In6 , and power-grids In8 .

A recent study has generalized the Kuramoto model by including an additional interaction term that breaks the rotational symmetry of the dynamics explicitly and unveiled a rich phase diagram with stationary and standing wave phases due to the symmetry breaking interaction In10 . Specifically, the authors have considered unimodal frequency distributions and revealed the emergence of a stationary state, characterized by time independent amplitude and phase of the complex Kuramoto order parameter, facilitated by the symmetry breaking interaction, which is otherwise absent in the original Kuramoto model that allows for the rotational symmetry of the dynamics. Interesting, in this work, we elucidate that the Kuramoto model can be translated into the Winfree model by the introduction of the additional symmetry breaking coupling and consequently, one can obtain the phase diagrams of both these models simply by tuning the symmetry breaking parameter $q$, thereby bridging the dynamics of both the models. Note that the macroscopic dynamical states of the pulse coupled biological oscillators with different sensitive functions, characterizing the phase-response-curves of biological oscillators, are peculiar to the Winfree model and its generalizations, which are far from reach for the Kuramoto model and its variants. In particular, we consider both the unimodal and bimodal frequency distributions to explore the phase diagrams for various values of the symmetry breaking parameter $q$. On the one hand, we observe the typical phase diagram of the Kuramoto model characterized only by incoherent and standing wave states in the absence of the symmetry breaking interaction for the unimodal frequency distribution. On the other hand, we observe the phase diagram with incoherent state, standing wave pattern along with the synchronized stationary state and bistabilities among them, a typical nature of the Winfree model, for $q=1$. For an intermediate and increasing value of $q\in(0,1)$, one can find the onset of the stationary state and eventually the emergence of bistability among these states in the phase diagram, and enlargement of the bistable regions resulting in the phase diagram of the Winfree model.

All three states are also observed in both Kuramoto and Winfree models for symmetric bimodal frequency distributions along with the region of bistability. The degree of the spread of the different macroscopic dynamical states depends on the strength of the symmetry breaking parameter $q$. Interestingly, for asymmetric bimodal frequency distributions, increase in the degree of asymmetry of the frequency distributions favors the onset of bistable regions even for a rather low values of $q$, which otherwise cannot be observed with the symmetric bimodal and unimodal frequency distributions. We arrive at the phase diagrams by numerical simulation of the original equations of motion. We deduce the reduced low-dimensional evolution equations of motion for the order parameter using the Ott-Antonsen ansatz for both unimodal and bimodal frequency distributions. We also deduce the Hopf, pitchfork and saddle-node bifurcation curves from the governing equations of motion for the order parameters, which mediates the dynamical transitions in the phase diagrams. Homoclinic bifurcation curve is obtained from the XPPAUT software.

The plan of the paper is as follows. In Sec. II, we generalize the Kuramoto model by introducing a symmetry breaking coupling and elucidate that the latter bridges the Kuramoto model and the Winfree model. We deduce the reduced low-dimensional evolution equations for the complex order parameters corresponding to the discrete set of generalized Kuramoto model using the Ott-Antonsen ansatz for both unimodal and bimodal frequency distributions in Sec. III. We also deduce Hopf, pitchfork and saddle-node bifurcation curves from the evolution equations for the complex order parameters in Sec. III, mediating the dynamical transitions among the incoherent, standing wave and synchronized stationary states. In Sec. IV, we discuss the observed dynamical states and their transitions in the various phase diagrams. Finally, in Sec. VI, we summarize the results.

2 Model

We consider a nontrivial generalization of the Kuramoto model by including an interaction term that explicitly breaks the rotational symmetry of the dynamics In10 . The phase $\theta_{i}$ is governed by the set of N ordinary differential equations (ODEs),

[TABLE]

for $i=1,\ldots,N$, where $N\gg 1$. Here $\theta_{i}(t)$ is the phase of the $i$th oscillator at time $t$, $\varepsilon\geq 0$ is the coupling strength, and $q$ is the strength of the symmetry breaking coupling. Note that Eq. (1) reduces to the Kuramoto model by setting $q=0$ and on identifying $\varepsilon$ with the parameter $K>0$. Equation (1) can also be viewed as a variant of the celebrated Winfree model w1 ; w2 ; w3 ; w4 when $q=1$. The Winfree model takes the form

[TABLE]

where $P(\theta_{j})$ is the phase dependent pulse function and the functional form of the response function $Q(\theta)$ characterizes the phase-response curves of certain biological oscillators. From Eq. (1), it easy to recognize that $Q(\theta)={-}2\sin(\theta)$ and $P(\theta)=\cos(\theta)$. It is also evident that the symmetry breaking parameter ‘q’ bridges the Kuramoto and the Winfree models. Equation (1) corresponds to the Kuramoto model when $q=0$ and it corresponds to a variant of the Winfree model when $q=1$, as in Eq. (2). We consider the frequencies of the phase-oscillators are distributed both by the unimodal Lorentzian distribution given as

[TABLE]

and bimodal Lorentzian distribution represented as

[TABLE]

Here $\gamma$, $\gamma_{1}$ and $\gamma_{2}$ are the width parameter (half width at half maximum) of the Lorentzian and $\pm\omega_{0}$ are their central frequencies. Note that $\omega_{0}$ corresponds to the degree of detuning in the system, which is proportional to the separation between the two central frequencies. Note that the bimodal distribution $g(\omega_{0})$ is symmetric about zero when $\gamma_{1}=\gamma_{2}$. It is also to be noted that $g(\omega_{0})$ in Eq. (4) is bimodal if and only if the separation between their central frequencies are sufficiently greater than their widths. To be precise, it is required that $\omega_{0}>\gamma_{1,2}/\sqrt{3}$ for the distribution to be a bimodal, otherwise the classical results of the unimodal distribution holds good.

Heterogeneity in the frequency distribution plays a crucial role in the manifestation of a plethora of collective dynamics in a vast variety of natural systems. In particular, coexisting co-rotating and counter-rotating systems characterized by positive and negative frequencies, respectively, are wide spread in nature. For instance, counter-rotating spirals are observed in protoplasm of the Physarum plasmodium ref1 , counter-rotating vortices are inevitable in the atmosphere and ocean ref2 ; ref3 ; ref4 , in magnetohydrodynamics of plasma flow ref8 , Bose-Einstein condensates ref9 ; ref10 , and in other physical systems ref5 ; ref6 ; ref7 . Very recently, the counter-rotating frequency induced dynamical effects were also reported in the coupled Stuart-Landau oscillator with symmetry preserving as well as symmetry breaking couplings ref11 . The coexistence of co-rotating and counter-rotating oscillators was initially identified by Tabor ref12 , which is followed by a series of work employing co-rotating and counter-rotating oscillators. All these physical systems strongly suggest that counter-rotating time-evolving dynamical systems indeed exist in nature and play a pertinent role in the manifestation of their intriguing collective dynamics.

In the following, we will deduce the low-dimensional evolution equations for the complex macroscopic order parameters corresponding to both the unimodal and bimodal frequency distributions using the Ott-Antonsen (OA) ansatz Ott:2008 ; Ott:2009 . Subsequently, we also deduce the various bifurcation curves facilitating the dynamical transitions among the observed dynamical states in the phase diagrams.

3 Low-dimensional evolution equations for the macroscopic order parameters

We now provide an analysis of the dynamics (1), in the limit $N\to\infty$, by invoking the Ott-Antonsen ansatz. In this limit, the dynamics of the discrete set of equations (1) can be captured by the probability distribution function $f(\theta,\omega,t)$, defined such that $f(\theta,\omega,t){\rm d}\theta$ gives the probability of oscillators with phase in the range $[\theta,\theta+{\rm d}\theta]$ at time $t$. The distribution is $2\pi$-periodic in $\theta$ and obeys the normalization

[TABLE]

Since the dynamics (1) conserves the number of oscillators with a given $\omega$, the time evolution of $f$ follows the continuity equation

[TABLE]

where $v(\theta,\omega,t)$ is the angular velocity of the oscillators. From Eq. (1), we have,

[TABLE]

where $z^{\star}$ denotes the complex conjugate of the macroscopic order parameter defined as

[TABLE]

Now, $f(\theta,\omega,t)$ can be expanded in terms of Fourier series of the form

[TABLE]

where, $\alpha_{n}(\omega,t)$ is the $n$th Fourier coefficient, while c.c. denotes complex conjugation of the preceding sum within the brackets. The normalization condition in (5) is satisfied by the presence of the prefactor of $g(\omega)$ in (9). The Ott-Antonsen ansatz consists in assuming Ott:2008 ; Ott:2009

[TABLE]

Now, it is straightforward to obtain

[TABLE]

where,

[TABLE]

3.1 Unimodal Distribution

Substituting the partial fraction expansion of the unimodal frequency distribution $g(\omega)$ (3) in Eq. (12) and evaluating the integral using an appropriate contour integral, one can obtain the order parameter as

[TABLE]

From (11) and (13), one can obtain the evolution equation for the complex order parameter as

The above evolution equation for the complex order parameter $z(t)=r(t)e^{i\psi(t)}$ can be expressed in terms of the evolution equations in $r$ and $\psi$ as

[TABLE]

The above equations govern the reduced low-dimensional order parameter dynamics, which actually corresponds to the dynamics of the original discrete set of equations (1) in the limit $N\to\infty$ for the unimodal Lorentzian distribution function $g(\omega)$ (3). Now, we discuss the various asymptotic macroscopic dynamical states admitted by Eq. (15).

3.1.1 Incoherent (IC) state:

The incoherent (IC) state is characterized by time independent $z$ satisfying $z=z^{\star}=0$ (thus representing a stationary state of the dynamics (15)); correspondingly, one has $r=0$. The linear stability of such a state is determined by linearizing Eq. (LABEL:eq:z-dynamics) around $z=0$. By representing $z=u$ with u≪1,weobtain\begin{equation}\frac{\partial u}{\partial t}+(\gamma-i\omega_{0})u-\frac{\varepsilon}{2}\big{[}(u)-q(u^{\star})\big{]}=0.\end{equation}Decomposingu=u_x + i u_y

[TABLE] Thematrix

Mhasthecharacteristiceigenvalues\begin{equation}\lambda_{1,2}=\frac{-2\gamma+\varepsilon\pm\sqrt{\Delta}}{2},\end{equation}withΔ=(ε^2q^2-4ω_0^2)$.Notethatwehave$λ_1 > λ_2$.Thestabilitythresholdfortheincoherentstateisthenobtainedbyanalysing$λ_1$asafunctionof$*ε$and$*q$,andseekingtheparticularvalueof$ε$atwhich$λ_1$vanishesforagiven$q$.Thestabilitythresholdcanbeobtainedas\begin{aligned} &\varepsilon_{HB}=2\gamma,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\;\;\;\mbox{for}\;\;\Delta\leq 0,\\ &\varepsilon_{PF}=2\sqrt{\frac{\gamma^{2}+\omega_{0}^{2}}{1+q^{2}}}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\mbox{for}\;\;\Delta>0.\end{aligned}$

3.1.2 Synchronized stationary state (SSS):

Now, we explore the possibility of existence of the synchronized stationary state. Requiring that $r$ and $\psi$ have time-independent non-zero values in this case and hence equating the left hand side of equations (15) to zero, we obtain the two coupled equations for the synchronized stationary state as

[TABLE]

With some algebra, one can obtained the following expressions for the stationary $r$ and $\psi$:

[TABLE]

$r$ and $\psi$ can be calculated for a fixed set of parameters by numerically solving the above set of equations, which is then substituted back into the evolution equations for the low-dimensional order parameters to deduce the characteristic equation. The eigenvalues of the characteristic equation is then used to determine the saddle-node bifurcation curve in the suitable two parameter phase.

3.2 Bimodal Distribution

Now, we will deduce the low-dimensional evolution equations corresponding to the macroscopic order parameters for the original discrete set of equations (1) in the limit $N\to\infty$ for the asymmetric bimodal Lorentzian distribution function $g(\omega)$ (4). Expanding the latter using partial fractions and evaluating the integral in Eq. (12) using appropriate contour integral, one can obtained the complex order parameter as

[TABLE]

where

[TABLE]

Substitution it into Eq. (11) yields two coupled complex ordinary differential equations describing the evolution of two suborder parameters as

[TABLE]

The above evolution equations for the complex order parameters $z(t)_{1,2}=r(t)_{1,2}e^{i\psi(t)_{1,2}}$ can be expressed in terms of the evolution equations in $r_{1,2}$ and $\psi_{1,2}$, as

[TABLE]

and

[TABLE]

The above equations constitute the evolution equations for reduced low-dimensional order parameters corresponding to the dynamics (1) in the limit $N\to\infty$ and for the case of the asymmetric bimodal Lorentzian distribution $g(\omega)$ (4). Now, we discuss the various asymptotic macroscopic dynamical states admitted by Eqs. (25) and (26).

3.2.1 Incoherent state

The incoherent state is defined by $r_{1}$=$r_{2}$=0. A linear stability analysis of the fixed point $(z_{1},z_{2})$ = (0, 0) results in the stability condition,

[TABLE]

where, $a_{1}=\gamma_{1}+\gamma_{2},a_{2}=\gamma_{1}^{2}+\gamma_{2}^{2}$ and $a_{3}=\gamma_{1}-\gamma_{2}$. This stability curve actually corresponds to the pitchfork bifurcation curve across which the fixed point $(z_{1},z_{2})$ = (0, 0) (incoherent state) loses its stability leading to the synchronized stationary state. Note that the incoherent state loses it stability through the Hopf bifurcation, which results in the stability condition

[TABLE]

where, $b_{1}=\gamma_{1}+\gamma_{2},b_{2}=\gamma_{1}\gamma_{2}$ and $b_{3}=\gamma_{1}^{2}+\gamma_{2}^{2}$. The above stability curve corresponds to the Hopf bifurcation curve. The boundary of stable incoherent state is therefore enclosed by both the pitchfork bifurcation and Hopf bifurcation curves.

3.2.2 Synchronized stationary state

Deducing the solution for the synchronized stationary state for the asymmetry bimodal distribution may not be possible as $r_{1}~{}\neq~{}r_{2}$ and $\psi_{1}~{}\neq~{}\psi_{2}$. However, for the symmetry bimodal distribution characterized by $r_{1}~{}=~{}r_{2}$ and $\psi_{1}~{}=~{}-\psi_{2}$, one can deduce the equations for $r$ and $\psi$ as in (20) and obtain the saddle-node bifurcation curves as pointed out in Sec. 3.1.2.

4 Numerical Results

In this section, we will proceed to unravel the macroscopic dynamical states admitted by the generalized Kuramoto model (1) with explicit symmetry breaking coupling by constructing appropriate two parameter phase diagrams and classifying the underlying dynamical states from a numerical analysis of the governing equations of the original discrete model. Specifically, we will unravel the rich phase diagrams of the generalized Kuramoto model, using both unimodal and bimodal frequency distributions, for distinct values of the symmetry breaking parameter $q$. The number of oscillators is fixed as N = $10^{4}$, and we use the standard 4th-order Runge-Kutta integration scheme with integration step size h = 0.01 to solve the generalized Kuramoto model (1). Note that one can break the two-parameter phase into several segments and multiple copies of the same code can be simulated simultaneously for different values of the parameters to generate the data, which can then be concatenated to get the complete phase diagrams with a reasonable workstation.

The initial state of the oscillators ($\theta_{i}$’s ) is distributed with uniform random values between -$\pi$ and +$\pi$. We use the time averaged order parameter $R$ defined as

[TABLE]

where $r(t)=|Z|=|N^{-1}\sum_{j=1}^{N}e^{i\theta_{j}}|$. Incoherent state is characterized by $R=r(t)=0$, while the synchronized stationary state is characterized by $R=r(t)=const.$ Standing wave is characterized by the oscillating nature of $r(t)$. In order to distinguish the synchronized stationary state and the standing wave state more clearly, we use the Shinomoto-Kuramoto order parameter gallego ; ssyk1986

[TABLE]

where $\bar{z}$ denoted the long time average. Shinomoto-Kuramoto order parameter takes $\xi=0$ for the incoherent and synchronized stationary states, whereas it takes nonzero value for the standing wave state.

4.0.1 Phase diagrams for the unimodal distribution

We have depicted phase diagrams in the ($\omega_{0}/\gamma,\varepsilon/\gamma$) parameter space for different values of the symmetry breaking parameter $q$ in Fig. 1 in order to understand the effect of the explicit symmetry breaking interaction on the dynamics of Eq. (1) with unimodal frequency distribution. The phase diagram is demarcated into different dynamical regions using the value of the time averaged order parameter $R$ and the Shinomoto-Kuramoto order parameter $\xi$. Incoherent state (IC), synchronized stationary state (SS) and standing wave (SW), along with the bistable regions (dark and light gray shaded regions) are observed in the phase diagram. The parameter space indicated by light gray shaded region corresponds to the bistable regime between the incoherent and the synchronized stationary states, while that indicated by dark gray shaded region corresponds to the bistable regime between the standing wave and the synchronized stationary states,

Only the incoherent and standing wave states are observed in the phase diagram for $q=0$ (see  1(a)), a typical phase diagram of the Kuramoto model with unimodal frequency distribution. The line connected by the filled squares corresponds to the Hopf bifurcation curve, across which there is a transition from the incoherent state to the standing wave state. Note that a finite value of $q$ results in the loss of the rotational symmetry of the dynamics of the Kuramoto oscillators. Even a feeble value of $q$ manifests the synchronized stationary state in a rather large parameter space at the cost of the standing wave state (see Fig. 1(b) for q=0.1). There is a transition from the incoherent state to the standing wave state via the Hopf bifurcation curve $\varepsilon_{HB}$ (indicated by the line connected by filled squares) as a function of $\varepsilon/\gamma$ for $\omega_{0}/\gamma>0.1$. The standing wave state loses its stability via the homoclinic bifurcation (indicated by the dashed-dotted line) as a function of $\varepsilon/\gamma$ resulting in the synchronized stationary state. There is also a transition from the incoherent state to the synchronized stationary state for $\omega_{0}/\gamma\leq 0.1$ as a function of $\varepsilon/\gamma$ via the pitchfork bifurcation curve $\varepsilon_{PF}$ indicated by the solid line.

Further larger values of the symmetry breaking parameter results in the emergence of the bistability between the standing wave and the synchronized stationary states (indicated by dark shaded region) enclosed by the saddle-node bifurcation curve (indicated by dashed line) and the homoclinic bifurcation curve (see Fig. 1(c) for q=0.5). There is also a bistable region between the incoherent state and the synchronized stationary state (indicated by light grey shaded region) enclosed by the saddle-node bifurcation curve and the pitchfork bifurcation curve. For $q=1$, both the bistable regions enlarged in the phase diagram (see Fig. 1(d)), which is a typical phase diagram of the Winfree model with the unimodal frequency distribution. The phase diagrams for $q=0.5$ and $1.0$ have similar dynamics except for the regime shift and enhanced bistabilities in a larger parameter space. Thus, as the value of $q$ is increased from the null value to the unity, one can observe the transition from the phase diagram of the Kuramoto model to that of the Winfree model. Note that the Hopf, saddle-node and pitchfork bifurcation curves are the analytical bifurcation curves, Eqs. (3.1.1), (3.1.1) and (20) respectively, obtained from the low-dimensional evolution equations for the order parameters deduced in Sec. 3.1. Homoclinic bifurcation curve is obtained from the software XPPAUT xpp .

Time averaged order parameter $R$ and the Shinomoto-Kuramoto order parameter $\xi$ are depicted in Fig. 2 as a function of $\varepsilon/\gamma$ for different values of the symmetry breaking parameter $q$ and $\omega_{0}/\gamma$. The forward trace is indicated by the line connected by open circles, while the backward trace is indicated by the line connected by closed circles. There is a smooth (second order) transition from the incoherent to the standing wave states via the Hopf bifurcation $\varepsilon_{HB}$ at $\varepsilon/\gamma=2$ during both forward and reverse traces for $q=0.0$ and $\omega_{0}/\gamma=1$ as depicted in Figs. 2(a) and 2(d). In addition, to the smooth transition from the incoherent state to the standing wave state via the Hopf bifurcation $\varepsilon_{HB}$ at $\varepsilon/\gamma=2$, there is another smooth transition from the standing wave state to the synchronized stationary state via the homoclinic bifurcation $\varepsilon_{Hc}$ at $\varepsilon/\gamma=2.94$ in both the forward and reverse traces as shown in Fig. 2(b) for $q=0.5$ and $\omega_{0}/\gamma=1$. The transition from the standing wave state to the synchronized stationary state is also corroborated by the sharp fall of the Shinomoto-Kuramoto order parameter $\xi$ to the null value (see Fig. 2(e)). In contrast, there is an abrupt (first order) transition from the incoherent state to the synchronized stationary state at $\varepsilon/\gamma=2$ via the pitchfork bifurcation curve $\varepsilon_{PF}$ for $\omega_{0}/\gamma=1$ during the forward trace, whereas there is an abrupt transition from the synchronized stationary state to the incoherent state at $\varepsilon/\gamma=1.8$ via the saddle-node bifurcation $\varepsilon_{SN}$ during the reverse trace (see Fig. 2(c) for $q=1.0$) elucidating the presence of hysteresis and bistability between the incoherent state and the synchronized stationary state. The Shinomoto-Kuramoto order parameter $\xi$ takes the null value, in the entire range of $\varepsilon/\gamma$ in Fig. 2(f) for $q=1.0$, characterizing both the incoherent and the synchronized stationary states.

The observed dynamical states and their transitions are depicted in the $(q,\varepsilon/\gamma)$ parameter space for different $\omega_{0}/\gamma$ in Fig. 3. The bifurcations mediating the dynamical transitions are similar to those observed in Fig. 1. The phase diagram for $\omega_{0}/\gamma=0.4$ is shown in Fig. 3(a). There is a transition from the incoherent state to the standing wave state via the Hopf bifurcation curve for smaller values of the symmetry breaking parameter as a function of $\varepsilon/\gamma$. Larger values of the symmetry breaking parameter favor the synchronized stationary state in the entire range of $\varepsilon/\gamma$. However, in a narrow range of $q\in(0.36,0.46]$ (see Fig. 3(a)), there is a transition from the incoherent state to the standing wave state and then to the synchronized stationary state. There is also a transition from the incoherent state to the synchronized stationary state in the range of $q\in(0.46,0.6)$. Recall that $\omega_{0}$ quantifies the degree of detuning of the frequency distribution. Increase in the heterogeneity of the frequency distribution promotes bistable regions, incoherent and standing wave states, to a large region of the $(q,\varepsilon/\gamma)$ parameter space. For instance, the phase diagram for $\omega_{0}/\gamma=0.6$ is depicted in Fig. 3(b) elucidates the emergence of the bistable regions and enlarged regions of the incoherent and standing wave states as a function of $q$, a manifestation of increased heterogeneity. Further increase in the $\omega_{0}/\gamma$ enlarges the bistable regions, the incoherent and the standing wave states as depicted in Figs. 3(c) and  3(d) for $\omega_{0}/\gamma=1$ and $1.2$, respectively. These results are in agreement with the phase diagrams in Fig. 1 in the $(\omega_{0}/\gamma,\varepsilon/\gamma)$ parameter space for increasing values of the symmetry breaking parameter. Next, we will explore the effect of symmetric and asymmetric bimodal frequency distributions on the phase diagrams in the following.

4.0.2 Phase diagrams for bimodal distribution

In this section, we analyse the phase space dynamics of the generalized Kuramoto model (1) with symmetric bimodal frequency distribution (4) by setting $\gamma=\gamma_{1}=\gamma_{2}$ for increasing values of the strength of the symmetry breaking coupling. We have depicted the phase diagrams in the $(\omega_{0}/\gamma,\varepsilon/\gamma)$ parameter space for different values of the symmetry breaking parameter $q$ in Fig. 4. Note that the phase space dynamics of the Kuramoto model (see Fig.  4(a) for $q=0$) are similar to those of the Winfree model (see Fig.  4(d) for $q=1$) for the symmetric bimodal frequency distribution except for the regime shift. The dynamical states and the bifurcation curves are similar to those in Fig. 1. Increasing the strength of the symmetry breaking coupling favors the synchronized stationary state and the bistable states in a large region of the parameter space as evident from Fig.  4(b) and  4(c) for $q=0.5$ and $q=0.8$, respectively. Note that a large heterogeneity in the frequency distribution favor the incoherent and the standing wave states in a rather large region of the phase diagram for smaller $q$ and $\varepsilon$ (see Fig.  4(a) for $q=0$). Nevertheless, the synchronized stationary state predominates the phase diagram for larger strength of the symmetry breaking coupling and $\varepsilon$ despite the presence of a large heterogeneity in the frequency distribution (see Fig.  4(d) for $q=1$).

Next, we analyze the phase space dynamics of the generalized Kuramoto model (1) with asymmetric bimodal frequency distribution (4) by increasing the strength of the symmetry breaking coupling and the degree of asymmetry between the bimodal frequency distributions. We have depicted the phase diagrams in the $(\omega_{0}/\gamma_{2},\varepsilon/\gamma_{2})$ parameter space for different values of the symmetry breaking parameter $q$ in Fig. 5. Again, the dynamical states and the bifurcation curves are similar to those in Fig. 1. Phase diagram for $q=0.1$ and $\gamma_{1}/\gamma_{2}=0.6$ is depicted in Fig. 5(a). For most values of $\omega_{0}/\gamma_{2}$, there is a transition from the incoherent state to the synchronized stationary state via the standing wave state and there is no bistability for $\gamma_{1}<\gamma_{2}$. However, there is a transition from the incoherent state to the synchronized stationary state in a large range of $\omega_{0}/\gamma\in(0,1)$ and the emergence of bistable states for $\gamma_{1}>\gamma_{2}$ as depicted in Fig. 5(b) for $\gamma_{1}/\gamma_{2}=1.2$. It is evident that bistable states emerge even for low values of the symmetry breaking coupling when $\gamma_{1}>\gamma_{2}$. Note that bistable states emerge even for $\gamma_{1}<\gamma_{2}$ but for a large strength of the symmetry breaking coupling (see Fig. 5(c) for $q=1$ and $\gamma_{1}/\gamma_{2}=0.6$). The spread of the bistable states increases for $q=1$ and $\gamma_{1}/\gamma_{2}=1.2$ as illustrated in Fig. 5(d). Thus, larger $\gamma_{1}/\gamma_{2}$ and $q$ favor the emergence of the bistable states.

Phase diagrams in the $(q,\varepsilon/\gamma_{2})$ parameter space is depicted in Figs. 6(a) and  6(b) for $\gamma_{1}/\gamma_{2}=0.6$ and $1.2$, respectively, and for $\omega_{0}/\gamma_{2}=1$. The dynamical states and the bifurcation curves are similar to those in Fig. 1. There is a transition from the incoherent state to the synchronized stationary state via the standing wave state for small values of $q$ (see Fig.. 6(a)) similar to that in Fig. 5(a). However, for larger values of $q$ multistability between the standing wave and the synchronized stationary state emerges (dark shaded region in the inset) in addition to the above dynamical transition. For $\gamma_{1}>\gamma_{2}$, there a transition from the incoherent state to the standing wave state along with the bistability among them in a rather narrow range of $q\in(0,0.4)$ as a function of $\varepsilon/\gamma_{2}$ as shown in inset of Fig. 6(b). For $q>0.4$, there is a transition from the incoherent state to the synchronized stationary state with the onset of bistability (light grey shaded region) between them. Phase diagrams in the $(q,\omega_{0}/\gamma_{2})$ parameter space is depicted in Figs. 6(c) and  6(d) for $\gamma_{1}/\gamma_{2}=0.6$ and $1.2$, respectively, for $\varepsilon/\gamma_{2}=2.5$. There is a transition from the synchronized stationary state to the standing wave state as a function of $\omega_{0}/\gamma_{2}$ for $\gamma_{1}<\gamma_{2}$ (see Fig. 6(c)) via the homoclinic bifurcation curve. Both the bistable states emerge when $\gamma_{1}>\gamma_{2}$ as shown in Fig. 6(c) for $\gamma_{1}=1.2$.

5 Conclusions

We have considered a nontrivial generalization of the paradigmatic Kuramoto model by using an additional coupling term that explicitly breaks the rotational symmetry of the Kuramoto model. The strength of the symmetry breaking coupling is found to play a key role in the manifestation of the dynamical states and their transitions along with the onset of bistability among the observed dynamical states in the phase diagram. A typical phase diagram of the Kuramoto model is transformed into a typical phase diagram of the Winfree mode for the unit value of the strength of the symmetry breaking coupling thereby bridging the dynamics of both the Kuramoto and Winfree models. Large values of the strength of the symmetry breaking coupling favor the manifestation of bistable regions and synchronized stationary state in a large region of the phase diagram. The dynamical transitions in the bistable region are characterized by an abrupt (first-order) transition in both the forward and reverse traces. Phase diagrams of both the Kuramoto and Winfree models resemble each other for symmetric bimodal frequency distribution except for the regime shifts and the degree of the spread of the dynamical states and bistable regions. Nevertheless, for asymmetric bimodal frequency distribution one cannot observe the bistable states for low values of the strength of the symmetry breaking coupling when $\gamma_{1}<\gamma_{2}$. In contrast, bistable states emerge even for $\gamma_{1}<\gamma_{2}$ for a large strength of the symmetry breaking coupling. Larger $\gamma_{1}/\gamma_{2}$ and larger $q$ favors the emergence of the bistable states in the case of the asymmetric bimodal frequency distribution. A large $\omega_{0}$ and consequently a large degree of heterogeneity facilitates the spread of the incoherent and standing wave states in the phase diagram for a low strength of the symmetry breaking coupling. However, a large $q$ promotes the spread of the synchronized stationary state and bistable regions in the phase diagram despite the degree of heterogeneity in the frequency distribution. We have deduced the low-dimensional evolution equations for the complex order parameters using the Ott-Antonsen ansatz for both unimodal and bimodal frequency distributions. We have also deduced the Hopf, pitchfork, saddle-node bifurcation curves from the low-dimensional evolution equations for the complex order parameters. Homoclinic bifurcation curve is obtained from XPPAUT software. Simulation results, obtained from the original discrete set of equations agrees well with the analytical bifurcation curves. We sincerely believe that our results will shed more light and enhance our current understanding of the effects of symmetry breaking coupling in the phase models and bridges the dynamics of two distinctly different phase models, which are far from reach otherwise.

6 Acknowledgements

The work of V.K.C. is supported by the DST-CRG Project under Grant No. CRG/2020/004353 and DST, New Delhi for computational facilities under the DST-FIST program (SR/FST/PS- 1/2020/135)to the Department of Physics. M.M. thanks the Department of Science and Technology, Government of India, for provid- ing financial support through an INSPIRE Fellowship No. DST/INSPIRE Fellowship/2019/IF190871. S.G. acknowledges support from the Science and Engineering Research Board (SERB), India under SERB-TARE scheme Grant No. TAR/2018/000023 and SERB-MATRICS scheme Grant No. MTR/2019/000560. He also thanks ICTP – The Abdus Salam International Centre for Theoretical Physics, Trieste, Italy for support under its Regular Associateship scheme. DVS is supported by the DST-SERB-CRG Project under Grant No. CRG/2021/000816.

Data Availability Statement: No Data associated in the manuscript. The data sets on the current study are available from the corresponding author on reasonable request.