The Sakaguchi-Kuramoto model in presence of asymmetric interactions that break phase-shift symmetry

The fundamental paradigm of the Kuramoto model Kuramoto:1984 ; Strogatz:2000 ; Acebron:2005 ; Gupta:2014 ; Gupta:2018 has been widely employed over the years in studying collective behavior of interacting stable limit-cycle oscillators. The model has been particularly successful in explaining spontaneous collective synchronization, a phenomenon exhibited by large ensembles of coupled oscillators and encountered across disciplines, e.g., in physics, biology, chemistry, ecology, electrical engineering, neuroscience, and sociology Pikovsky:2001 . Examples of collective synchrony include synchronized firing of cardiac pacemaker cells Peskin:1975 , synchronous emission of light pulses by groups of fireflies Buck:1988 , chirping of crickets Walker:1969 , synchronization in ensembles of electrochemical oscillators Kiss:2002 , synchronization in human cerebral connectome Schmidt:2015 , and synchronous clapping of audience Neda:2000 . The Kuramoto model comprises $N$ globally-coupled stable limit-cycle oscillators with distributed natural frequencies interacting symmetrically with one another through the sine of their phase differences. Denoting by $\theta_{j}\in[-\pi,\pi]$ the phase of the $j$-th oscillator, $j=1,2,\ldots,N$, the dynamics of the model is described by a set of $N$ coupled first-order nonlinear differential equations of the form

$$\frac{{\rm d}\theta_{j}}{{\rm d}t}=\omega_{j}+\frac{K}{N}\sum_{k=1}^{N}\sin(\theta_{k}-\theta_{j}),$$

(1)

where $K\geq 0$ is the coupling constant, and $\omega_{j}$ is the natural frequency of the $j$-th oscillator. The frequencies $\{\omega_{j}\}$ denote a set of quenched-disordered random variables sampled independently from a distribution $g(\omega)$ usually taken to be unimodal, that is, symmetric about the mean $\omega_{0}$ and decreasing monotonically and continuously to zero with increasing $|\omega-\omega_{0}|$. In the dynamics (1) the interaction is symmetric: the effect on the $j$-th oscillator due to the $k$-th one being proportional to $\sin(\theta_{k}-\theta_{j})$ is equal in magnitude (but opposite in sign) to the one on the $k$-th oscillator due to the $j$-th one.

In contrast to the dynamics (1), interactions between oscillators may more generally be asymmetric. Considering symmetric interaction in the dynamics is only an approximation that may simplify theoretical analysis, but which may fail to capture important phenomena occurring in real systems. For example, asymmetric interaction leads to novel features such as families of travelling wave states Iatsenko:2013 ; Petkoski:2013 , glassy states and super-relaxation Iatsenko:2014 , and so forth, and has been invoked to discuss coupled circadian neurons Gu:2016 , dynamic interactions Yang:2020 ; Sakaguchi:1988 , etc. A generalization of the Kuramoto model (1) that accounts for asymmetric interaction is the so-called Sakaguchi-Kuramoto model, with the dynamics described by the equation of motion Sakaguchi:1986

$$\frac{{\rm d}\theta_{j}}{{\rm d}t}=\omega_{j}+\frac{K}{N}\sum_{k=1}^{N}\sin(\theta_{k}-\theta_{j}+\alpha),$$

(2)

where $0\leq\alpha<\pi/2$ is the so-called phase lag parameter. It is now easily seen that owing to the presence of an $\alpha\neq 0$, the influence on the $j$-th oscillator due to the $k$-th one is not any more equal in magnitude to the influence on the $k$-th oscillator due to the $j$-th one. As has been the case with the Kuramoto model, the model (2) and its variants have been successfully invoked to study a variety of dynamical scenarios, including, to name a few, disordered Josephson series array Wiesenfeld:1996 , time-delayed interactions Yeung:1999 , hierarchical populations of coupled oscillators Pikovsky:2008 , chaotic transients Wolfrum:2011 , dynamics of pulse-coupled oscillators Pazo:2014 , etc. We note in passing that asymmetric interaction arises in the Kuramoto model with spatially heterogeneous time-delays in the interaction, see Ref. Skardal:2018 . Let us note that both the dynamics (1) and (2) satisfy phase-shift symmetry, whereby rotating every phase by an arbitrary angle same for all leaves the dynamics invariant. In particular, one may implement the transformation $\theta_{j}(t)\to\theta_{j}(t)+\omega_{0}t~{}\forall~{}j,t$, which is tantamount to viewing the dynamics in a frame rotating uniformly with frequency $\omega_{0}$ with respect to an inertial frame, e.g., the laboratory frame.

The paper is organized as follows. In the following section, we describe the model of study, Eq. (3), and summarize our results obtained later in the paper. In Section III, we present our Ott Antonsen ansatz-based analysis of the dynamics (3) and the existence and stability of all the different possible states at long times. In Section IV, we discuss our analytical vis-à-vis numerical results. Finally, in Section V, we draw our conclusions.

In this work, we study a generalization of the Sakaguchi-Kuramoto dynamics (2) by including an interaction term that explicitly breaks the phase-shift symmetry of the dynamics. To this end, we consider the following set of $N$ coupled nonlinear differential equations:

$$\frac{{\rm d}\theta_{j}}{{\rm d}t}=\omega_{j}+\frac{1}{N}\bigg{[}\epsilon_{1}\sum_{k=1}^{N}\sin(\theta_{k}-\theta_{j}+\alpha)+\epsilon_{2}\sum_{k=1}^{N}\sin(\theta_{k}+\theta_{j}+\alpha)\bigg{]},$$

(3)

where the real parameters $\epsilon_{1,2}$ denote the coupling constants. Equation (3) may be obtained as a result of phase reduction analysis of a model of all-to-all-coupled Stuart-Landau oscillators with a symmetry-breaking form of interaction  cj23 ; cj3 .

Equation (3) may also be written as

$$\frac{d\theta_{i}}{dt}=\omega_{i}+\sum_{k=1}^{2}Q_{k}(\theta_{i})\frac{f_{k}}{N}\sum_{j=1}^{N}P_{k}(\theta_{j}),$$

(4)

where we have $f_{1}=\epsilon_{1}+\epsilon_{2}$, $f_{2}=\epsilon_{2}-\epsilon_{1}$, $P_{1}(\theta)=\cos(\theta)$, $P_{2}(\theta)=\sin(\theta)$, $Q_{1}(\theta)=\sin(\theta+\alpha)$ and $Q_{2}(\theta)=\cos(\theta+\alpha)$. Equation (4) resembles the Winfree model winfree , so we may consider our model of study as a variant of a Winfree-type system. The Winfree model played a seminal role in the field of collective synchrony, inspiring the Kuramoto model as well as promoting recent advances in theoretical neuroscience wf-3 ; wf-4 . In spite of the simplifying assumptions of the Winfree model, namely, that of uniform all-to-all weak coupling, analytical solutions have been found only recently using the Ott-Antonsen ansatz wf-5 ; Pazo:2014 , see also wf-6 .

Note that on putting $\epsilon_{2}=0$ in Eq. (3), one recovers the dynamics of the Sakaguchi-Kuramoto model (2) on identifying $\epsilon_{1}$ with the parameter $K\geq 0$ in the latter, and so we take $\epsilon_{1}$ to be positive. Here, we consider $\epsilon_{2}$ to satisfy $\epsilon_{2}\geq 0$. The dynamics (3) has phase-shift symmetry only with the choice $\epsilon_{2}=0$, and so the $\epsilon_{2}$-term acts a phase-shift-symmetry-breaking interaction (With $\epsilon_{2}\neq 0$, the only case of phase-shift symmetry is with respect to the transformation $\theta_{j}\to\theta_{j}+\pi~{}\forall~{}j$.). In contrast to the Sakaguchi-Kuramoto model, the dynamics (3) is not invariant with respect to the transformation $\theta_{j}\to\theta_{j}+\omega_{0}t~{}\forall~{}j,t$ owing to the presence of the $\epsilon_{2}$-term. Consequently, the mean $\omega_{0}$ of the frequency distribution $g(\omega)$ would have a significant influence on the dynamics (3), and cannot be gotten rid of by viewing the dynamics in a frame rotating uniformly with frequency $\omega_{0}$ with respect to the laboratory frame, as is possible with the Sakaguchi-Kuramoto model. We note in passing that the so-called active rotators may also be described by phase oscillators. Active rotators are well-established paradigms for excitable systems, and possess a term that breaks the phase-shift symmetry  Sakaguchi:1988 ; AR2 ; AR3 .

We will consider in this work a unimodal $g(\omega)$. Specifically, we will consider a Lorentzian $g(\omega)$:

$$g(\omega)=\frac{\gamma}{\pi((\omega-\omega_{0})^{2}+\gamma^{2})};~{}~{}\gamma>0.$$

(5)

Analysis of synchronization in the context of the Kuramoto class of models involves introducing a complex-valued order parameter given by Kuramoto:1984

$$z(t)=r(t)e^{\mathrm{i}\psi(t)}\equiv\frac{1}{N}\sum_{j=1}^{N}e^{\mathrm{i}\theta_{j}(t)}.$$

(6)

The quantity $r(t);~{}0\leq r(t)\leq 1$, measures the amount of synchrony present in the system at time $t$, while $\psi(t)\in[-\pi,\pi]$ gives the average phase. Considering the limit $N\to\infty$, both the dynamics (1) and (2) allow for a stationary state to exist at long times. By stationary state is meant a state with time independent $z$. In such a state, the system may be found in a synchronized or an incoherent state. The two states are characterized by the time-independent value that $r(t)$ assumes at long times, namely, a zero value in the incoherent state and a non-zero value in the synchronized state.

Considering the limit $N\to\infty$, this work aims at a detailed characterization of the long-time ($t\to\infty$) limit of the dynamics (3) and understanding of what new features with respect to the Sakaguchi-Kuramoto model are brought in by the introduction of the phase-shift-symmetry-breaking $\epsilon_{2}$-term. An asymmetric phase-shift-symmetry-breaking interaction, such as the one being considered by us, is an hitherto unexplored theme in the context of the Sakaguchi-Kuramoto model, and raises many queries: Does the introduction of the $\epsilon_{2}$-term suffice to allow for the dynamics to have a stationary state, and if so, what is the nature of the stationary state? What is the range of parameter values for which one has a synchronized stationary state? Most importantly, what is the complete bifurcation diagram in the ($\epsilon_{1},\epsilon_{2}$)-plane?

Recently, the dynamics of the Kuramoto model in presence of phase-shift-symmetry-breaking symmetric interaction was shown to unveil a rather rich bifurcation diagram with existence of both stationary and non-stationary synchronized states Chandru:2020 ; the model studied therein corresponds to the dynamics (3) with $\alpha$ set to zero. It is therefore pertinent that we summarize right at the outset what new features does the dynamics (3) offer with respect to (i) the Sakaguchi-Kuramoto model (i.e., the case $\epsilon_{2}=0$ in Eq. (3)), and (ii) the Kuramoto model with an additional phase-symmetry-breaking symmetric interaction (i.e., the case $\alpha=0$ in Eq. (3)). We focus on unimodal frequency distributions. As is well known, the Sakaguchi-Kuramoto model shows either an incoherent ($r(t)$ vanishes as $t\to\infty$) or a synchronized ($r(t)$ as $t\to\infty$ has a nonzero value, and moreover, $r(t)$ does not oscillate as a function of time) state Sakaguchi:1986 . These states are observed only in a co-rotating frame rotating uniformly with frequency $\omega_{0}$ with respect to the laboratory frame. As one tunes the strength of $\epsilon_{1}$, one has typically a continuous synchronization transition between an incoherent and a synchronized state, though for certain specific unimodal frequency distributions and carefully chosen phase lag $\alpha$, one may observe more than one non-trivial synchronization transitions Oleh-SK-1 ; Oleh-SK-2 . Now, coming to the model (ii), it has been observed that in addition to the incoherent and the synchronized state, the dynamics may also exhibit what we called a standing wave state, namely, a state in which $r(t)$ at long times oscillates as a function of time, yielding a non-zero time-independent time average. Note that for the model (ii), all the three mentioned states are observed in the laboratory frame itself. Moreover, the bifurcation diagram in the $(\epsilon_{1},\epsilon_{2})$-plane exhibits a continuous transition between the incoherent and the standing wave state, but instead a first-order transition between the incoherent and the synchronized state, and between the standing wave and the synchronized state. The latter fact implies regions of metastability or coexistence between the incoherent and the synchronized state, and between the standing wave and the synchronized state. In the light of the aforementioned facts, we now summarize the relevant features of the bifurcation diagram in the $(\epsilon_{1},\epsilon_{2})$-plane for the model (3) that we report in this work. The model exhibits in the laboratory frame the incoherent (IC), the synchronized and the standing wave state; for better highlighting of the differences between the last two states, we call them the oscillation death (OD) state and the oscillatory synchronized (OS) state, respectively. However, in contrast to the model (ii), we now have multistability/coexistence between all the different pair of states, that is, there are IC-OD, OS-OD, IC-OS coexistence regions, and in addition, a very peculiar and striking three-state-coexistence region between IC-OS-OD. We believe that this three-state coexistence is new to the literature on Kuramoto and Sakaguchi-Kuramoto models. The highlight of our work is thus the revelation that asymmetric interactions lead to a very rich bifurcation diagram when compared to the original Sakaguchi-Kuramoto model and the Kuramoto model with an additional phase-symmetry-breaking symmetric interaction (i.e., the case $\alpha=0$ in Eq. (3)). In particular, with respect to models (i) and (ii), all the possible transitions in the model (3) are first-order, and, moreover, a new three-state-coexistence region is observed; we may conclude therefore that introducing a phase-shift-symmetry-breaking interaction in the Sakaguchi-Kuramoto model drastically and non trivially modifies the bifurcation diagram with respect to the Sakaguchi-Kuramoto model and with respect to the model studied in Ref. Chandru:2020 .

The rest of the paper is devoted to a derivation of the aforementioned results. For Lorentzian $g(\omega)$, Eq. (5), we use exact analytical results derived by applying the so-called Ott-Antonsen ansatz, combined with numerical integration of the dynamics (3) for large $N$, to support the key result of our work, the bifurcation diagram of Fig. 3. The celebrated Ott-Antonsen ansatz is a remarkable method of analysis introduced to study coupled oscillator systems, which allows to rewrite in the limit $N\to\infty$ the dynamics of coupled networks of phase oscillators in terms of a few collective variables Ott:2008 ; Ott:2009 .

We now analyse the dynamics (3) in the limit $N\to\infty$ by employing the Ott-Antonsen ansatz. In this limit, the dynamics may be characterized by defining a single-oscillator distribution function $f(\theta,\omega,t)$, which is such that $f(\theta,\omega,t)$ is the probability density of finding an oscillator with natural frequency $\omega$ and phase $\theta$ at time $t$. The distribution is $2\pi$-periodic in $\theta$ and is moreover normalized as

$$\int_{0}^{2\pi}{\rm d}\theta~{}f(\theta,\omega,t)=g(\omega)~{}\forall~{}\omega,t.$$

(7)

The evolution of $f$ is given by the continuity equation describing the conservation of number of oscillators of a given frequency under the dynamics (3):

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

(8)

where $v(\theta,\omega,t)\equiv{\rm d}\theta/{\rm d}t$ is the angular velocity at position $\theta$ at time $t$. From Eq. (3), we get

	$\displaystyle v(\theta,\omega,t)=\omega+$	$\displaystyle\frac{\epsilon_{1}}{2\mathrm{i}}[ze^{-\mathrm{i}(\theta-\alpha)}-z^{\star}e^{\mathrm{i}(\theta-\alpha)}]$
		$\displaystyle-\frac{\epsilon_{2}}{2\mathrm{i}}[ze^{\mathrm{i}(\theta+\alpha)}-z^{\star}e^{-\mathrm{i}(\theta+\alpha)}],$		(9)

where $z=z(t)$ is obtained as the $N\to\infty$-limit generalization of Eq. (6):

$$z=\int{\rm d}\omega\int{\rm d}\theta~{}f(\theta,\omega,t)e^{\mathrm{i}\theta}.$$

(10)

In their seminal papers Ott:2008 ; Ott:2009 , Ott and Antonsen pointed out that all the attractors of the Kuramoto model and its many variants and for the case of the Lorentzian $g(\omega)$, Eq. (5), may be found by using the ansatz

$$f(\theta,\omega,t)=\frac{g(\omega)}{2\pi}\left[1+\sum_{n=1}^{\infty}\left(f_{n}(\omega,t)e^{\mathrm{i}n\theta}+{\rm c.c.}\right)\right].$$

(11)

Here, c.c. stands for a term obtained by complex conjugation of the first term within the brackets occurring in the sum, and

$$f_{n}(\omega,t)=\left[a(\omega,t)\right]^{n},$$

(12)

with arbitrary $a(\omega,t)$ satisfying $|a(\omega,t)|<1$, together with the requirements that $a(\omega,t)$ may be analytically continued to the whole of the complex-$\omega$ plane, it has no singularities in the lower-half complex-$\omega$ plane, and $|a(\omega,t)|\to 0$ as ${\rm Im}(\omega)\to-\infty$.

Using the ansatz (11) in Eq. (8), we straightforwardly get

$$\frac{\partial a}{\partial t}+\mathrm{i}\omega a+\frac{\epsilon_{1}}{2}(za^{2}e^{\mathrm{i}\alpha}-z^{\star}e^{-\mathrm{i}\alpha})+\frac{\epsilon_{2}}{2}(ze^{\mathrm{i}\alpha}-z^{\star}a^{2}e^{-\mathrm{i}\alpha})=0,$$

(13)

where we have

$$z^{\star}=\int_{-\infty}^{\infty}a(t,\omega)g(\omega)d\omega=a((\omega_{0}-{\rm i}\gamma),t),$$

(14)

obtained by using in Eq. (10) the ansatz (11) and Eq. (5) for $g(\omega)$, then converting the integral therein into a contour integral, and finally evaluating the contour integral. Equation (13) may then be rewritten as

	$\displaystyle\frac{\partial z}{\partial t}-{\rm i}(\omega_{0}+\mathrm{i}$	$\displaystyle\gamma)z+\frac{\epsilon_{1}}{2}\big{[}|z|^{2}ze^{-\mathrm{i}\alpha}-ze^{\mathrm{i}\alpha}\big{]}$
		$\displaystyle+\frac{\epsilon_{2}}{2}\big{[}z^{\star}e^{-\mathrm{i}\alpha}-z^{3}e^{\mathrm{i}\alpha}\big{]}=0.$		(15)

The above equation expressed in terms of the quantities $r$ and $\psi$, Eq. (6), gives the following two coupled equations:

		$\displaystyle\frac{{\rm d}r}{{\rm d}t}=-\gamma r-\frac{r}{2}(r^{2}-1)(\epsilon_{1}\cos(\alpha)-{\epsilon_{2}}\cos(2\psi+\alpha)),$
		$\displaystyle\frac{{\rm d}\psi}{{\rm d}t}=\omega_{0}+\frac{1}{2}(r^{2}+1)(\epsilon_{1}\sin(\alpha)+{\epsilon_{2}}\sin(2\psi+\alpha)).$

The above equations constitute the Ott-Antonsen-ansatz-reduced order parameter dynamics corresponding to the dynamics (3) in the limit $N\to\infty$ and for the case of the Lorentzian $g(\omega)$, Eq. (5). Let us remark that these equations are invariant under the transformation $(r,\psi)\to(r,\psi+\pi)$, which may be traced back to the fact that the original dynamics (3) is invariant under the transformation $\theta_{j}\to\theta_{j}+\pi~{}\forall~{}j$.

Note that for $\epsilon_{2}=0$ and $\alpha=0$, when one has the Kuramoto model, the two equations in (LABEL:eq:r-dynamics) are decoupled, and the equations then allow for only uniform rotation of $\psi$ with frequency $\omega_{0}$. The case of our interest, namely, $\epsilon_{2}\neq 0$ and $\alpha\neq 0$, is analysed in the subsection that follows. It would prove convenient for the discussion presented therein to define the time average of $r(t)$ in the long-time ($t\to\infty$) limit as

$$R\equiv\lim_{t\to\infty}\frac{1}{\tau}\int_{t}^{t+\tau}{\rm d}t^{\prime}~{}r(t^{\prime}).$$

(17)

In the preceding section, we discussed the existence of the IC, the OS and the OD state. Throughout this work, while discussing numerical integration results for the dynamics (3), unless stated otherwise, we consider the number of oscillators to be $N=10^{5}$ and the natural frequency distribution to be given by the Lorentzian (5) with $\gamma=0.3$ and $\omega_{0}=0.6$. Numerical integration is performed by employing a standard fourth-order Runge-Kutta integration scheme with integration step size ${\rm d}t=0.01$.

In Fig. 1, we have plotted the temporal behavior of the IC, the OS and the OD, for $\alpha=1.3$. The data are obtained from numerical integration of the dynamics (3). At fixed $\epsilon_{2}=2.0$, one has the OD at $\epsilon_{1}=0.5$, the IC at $\epsilon_{1}=1.5$, and the OS at $\epsilon_{1}=2.5$. As expected, we see that the three states are distinguished by the different long-time temporal behavior of $r(t)$. Namely, for the IC, $r(t)$ takes the value zero in the $t\to\infty$ limit. In the OS, $r(t)$ as $t\to\infty$ oscillates in time. For the OD, $r(t)$ as $t\to\infty$ has a non-zero time-independent constant value. Figure 2 shows as a function of $\epsilon_{1}$ the mean-ensemble frequency $f$ and the mean-field frequency $\Omega$ at long times. For details of computation of these quantities, we refer the reader to Ref. Chandru:2020 .

In order to demonstrate the influence of $\alpha$ on the stability of the different states, we now discuss the bifurcation diagram in the ($\epsilon_{1},\epsilon_{2}$)-plane. Figure 3 shows the stable states in the $(\epsilon_{1},\epsilon_{2})$-plane for the model (3) with $\alpha=1.3$. Here the states represent either the IC or the OS or the OD state. The regions representing stable existence of the IC, the OS and the OD state have been constructed by analysing the long-time numerical solution of the dynamics (3) for $N=10^{4}$, namely, by asking at given values of $\epsilon_{1}$ and $\epsilon_{2}$ the quantity $r(t)$ at long times represents which one of the three possible behavior depicted in Fig. 1. The stability boundary for the IC (blue dot-dashed line) is given by the Ott-Antonsen-ansatz-based equation, Eq. (LABEL:eq:stability-ISS). The stability boundary for the OD (black dashed line) is obtained from the Ott-Antonsen-ansatz-based equation, Eq. (LABEL:eq:stability-SSS). The procedure is as follows: For given $\gamma$, $\omega_{0}$ and $\alpha$, we vary at fixed $\epsilon_{2}$ the value of $\epsilon_{1}$ from high to low, and use the first equation in (LABEL:eq:stability-SSS) to ascertain the particular value of $\epsilon_{1}$ when for the first time the equation admits a solution for $r$ in the range $0<r\leq 1$; this particular value of $\epsilon_{1}$ gives the stability threshold of the OD at the chosen value of $\epsilon_{2}$. The process is repeated for different values of $\epsilon_{2}$. The stability boundary for the OS (maroon continuous line) is obtained from the XPPAUT analysis discussed in Section III.1.2 for the Ott-Antonsen-ansatz-based dynamics, Eq. (LABEL:eq:r-dynamics). The detailed procedure involves the following: For given $\gamma$, $\omega_{0}$ and $\alpha$, we vary at fixed $\epsilon_{2}$ the value of $\epsilon_{1}$ from low to high, and use the numerical solution of the Eq. (LABEL:eq:r-dynamics) (rewritten in terms of $x=x(t)$ and $y=y(t)$) to ascertain the particular value of $\epsilon_{1}$ when we first encounter at long times a stable limit cycle in the $(x,y)$-plane. We then repeat the process for different values of $\epsilon_{2}$.

Figure 3 shows in particular the presence of multistability or coexistence regions R1, R2, R3 and R4; these represent multistability between IC-OD, between OS-OD, between IC-OS-OD (coexistence of all the states) and between IC-OS, respectively. At a fixed $\epsilon_{1}$ and on tuning $\epsilon_{2}$ (or vice versa), one observes bifurcations as one crosses the different boundaries. When compared with the bifurcation diagram in case of symmetric interactions that either preserve or break phase-shift symmetry Chandru:2020 , we see that presence of asymmetry in the interaction leads to a very rich bifurcation diagram. The dashed lines L1, L2, and L3 in Fig. 3 indicate the values of $\epsilon_{2}$ for which the plots in Figs. 6 and  7, reported later in the paper, are obtained. Figure 3, the bifurcation diagram of model (3), is the key result of our work.

Figure 4 shows further details of the bifurcation diagram, and in particular, corresponding to the Ott-Antonsen-ansatz-based dynamics (LABEL:eq:r-dynamics), the XPPAUT-generated phase portraits in the $(x,y)$-plane for the IC, the OS and the OD state. The OD corresponds to two symmetrically-placed stable nodes at nonzero $x$ and $y$ (red filled circle); the phase portrait, panel (a), shows an unstable fixed point at $x=0,~{}y=0$ (purple open triangle), while that in panel (b) shows in addition to the unstable fixed point also two saddles (blue open circle). The IC represents a fixed point at $x=0,~{}y=0$ that is a stable spiral (red filled triangle), so that trajectories in course of evolution spiral in to the fixed point, see panel (c). The OS corresponds to a stable limit cycle (black filled square), so that trajectories emanating from the unstable fixed point at $x=0,~{}y=0$ (purple open triangle) as well as those from the region outside the limit cycle eventually approach the limit cycle in course of evolution (see panel (d)). In the figure, we show that the transitions between these states are mediated by either a pitchfork (PF, dot-dashed line), or a saddle-node (SN, dotted line), or a saddle-node limit-cycle (SNLC, continuous line), or a Hopf (HB, continuous line containing filled squares), or a homoclinic (HC, double-dot dashed line) bifurcation. The nature of bifurcations has been obtained from the XPPAUT analysis of the Ott-Antonsen-ansatz-based dynamics (LABEL:eq:r-dynamics).

In Fig. 5, we show the XPPAUT-generated phase-space portraits in each of the coexistence regions R1, R2, R3, R4, based on the Ott-Antonsen-ansatz-based dynamics (LABEL:eq:r-dynamics). In (a), corresponding to R1, we see the coexistence of IC (a stable spiral (red filled triangle)) and the OD (two stable nodes (red filled circle)), together with two saddles (blue open circle). In (b), corresponding to R2, the OS (a stable limit cycle (black filled square)) coexists with the OD (two stable nodes (red filled circle)), together with an unstable fixed point (purple open triangle) and two saddles (blue open circle). In (c), corresponding to R3, the OS (a stable limit cycle (black filled square)) coexists with the IC (a stable spiral (red filled triangle)) and the OD (two stable nodes (red filled circle)), and additionally, there are two saddles (blue open circle). In (d), corresponding to R4, the OS (a stable limit cycle (black filled square)) coexists with the IC (a stable spiral (red filled triangle)) (in this case, there is an unstable limit cycle (not shown) separating the two behavior).

Now that we have seen the complete bifurcation diagram and the phase portraits in the three states, the IC, the OS and the OD, and in regions of their coexistence, we turn to a discussion of the behavior of the quantity $R$ (the time-averaged value of $r(t)$ computed at long times, see Eq. (17)) as a function of adiabatically-tuned $\epsilon_{1}$ for representative values of $\epsilon_{2}$, with the aim to see signatures of bifurcation. We first let the system settle to the stationary state at a fixed value of $\epsilon_{1}$, and then tune $\epsilon_{1}$ adiabatically in time from low to high values while recording the value of $R$ in time; this corresponds to forward variation of $\epsilon_{1}$. Subsequently, we tune $\epsilon_{1}$ adiabatically in time from high to low values (backward variation of $\epsilon_{1}$). Adiabatic tuning ensures that the system remains in the stationary state at all times as the value of $\epsilon_{1}$ changes in time. In Fig. 6, we consider two values of $\epsilon_{2}$: Panels (a) and (b) are for $\epsilon_{2}=2.0$, while panels (c) and (d) are for $\epsilon_{2}=3.4$. Here, the lines correspond to results based on numerical integration of the dynamics (3), while symbols correspond to the analysis of the Ott-Antonsen-ansatz-based dynamics, Eq. (LABEL:eq:r-dynamics), and are obtained as follows: For the IC, the Ott-Antonsen analysis predicts that the IC state with $R=0$ exists only for $\epsilon_{1}$ larger than its threshold value given by Eq. (LABEL:eq:stability-ISS). For the OS, the quantity $R$ is obtained from the XPPAUT analysis of Eq. (LABEL:eq:r-dynamics). For the OD, the values of $R$ are obtained from numerical integration of Eq. (LABEL:eq:stability-SSS).

In panel (a), we observe hysteresis, implying occurrence of the coexistence region R1 between the OD and the IC and the coexistence region R4 between the IC and the OS. This is indeed consistent with Fig. 3 which shows that $\epsilon_{2}=2.0$ allows for the existence of R1 and R4 regions. The hysteresis seen in panel (c) implies occurrence of the coexistence region R2 between the OD and the OS, and is again consistent with Fig. 3 in which we see that $\epsilon_{2}=3.4$ allows for the existence of R2 region. The symbols in panels (b) and (d) represent XPPAUT-generated bifurcation plots for $y$ as a function of $\epsilon_{1}$ with fixed $\epsilon_{2}$ and based on the Ott-Antonsen-ansatz-dynamics (LABEL:eq:r-dynamics), and are a further confirmation of multistability and bifurcation behavior.

In panel (b), we see that when the IC becomes unstable (i) with the decrease of $\epsilon_{1}$, one observes a subcritical pitchfork bifurcation (PF) to give rise to the OD: at the bifurcation point, two symmetric unstable branches (blue open circle, represents saddles) and one stable branch corresponding to the IC (red filled triangle) go over to one unstable branch (purple open triangle, represents unstable fixed points); (ii) with the increase of $\epsilon_{1}$, one observes a saddle-node limit cycle bifurcation (SNLC) to give rise to a stable and an unstable limit cycle coexisting with the IC; the IC loses its stability with increase of $\epsilon_{1}$ via Hopf bifurcation (HB), which gives rise to the stable OS. At HB, the stable and the unstable limit cycle born at SNLC bifurcation collide with each other and disappear. One stable branch corresponding to the IC (red filled triangle) goes over at the bifurcation point to a stable branch corresponding to the OS (black filled square: upper and lower branches in the figure correspond to the maximum and the minimum of the corresponding stable limit cycle). Next, when the OD becomes unstable with the increase of $\epsilon_{1}$, one observes a saddle-node bifurcation (SN) to give rise to the IC: at the bifurcation point, one stable branch corresponding to the OD (red filled circle) and one unstable branch (blue open circle, represent saddles) annihilate each other. On the other hand, when the OS becomes unstable with the decrease of $\epsilon_{1}$, one observes a saddle-node limit-cycle bifurcation (SNLC) to give rise to the IC: at the bifurcation point, a stable branch corresponding to the OS (black filled square: upper and lower branches in the figure correspond to the maximum and the minimum of the corresponding stable limit cycle) and an unstable branch (purple open triangle, represents unstable spirals) annihilate each other.

In panel (d), we see that when the OD becomes unstable with the increase of $\epsilon_{1}$, one observes a saddle-node bifurcation to give rise to the OS: at the bifurcation point, a stable branch corresponding to the OD (red filled circle) and an unstable branch (blue open circle, represent saddles) annihilate each other. Similarly, when the OS becomes unstable with the decrease of $\epsilon_{1}$, one observes a homoclinic bifurcation: at the bifurcation point, a stable branch corresponding to the OS (a stable limit cycle) and represented by black filled squares (upper and lower branch correspond respectively to the maximum and the minimum of the limit cycle) collides with an unstable branch (blue open circle, represents saddles).

In order to witness the R3 region, one has to choose the value of $\epsilon_{2}$ carefully as R3 represents a very narrow region, see Fig. 3; we consider $\epsilon_{2}=3.1$ to witness the R3 region, see Fig. 7. Here, panel (a) (respectively, panel (b)) has been obtained by following the same procedure as the one followed in obtaining panels (a) and (c) (respectively, panels (b) and (d)) of Fig. 6. Here, we see the regions R1 and R2 as well. In panel (a), we see during the forward (adiabatic) variation of $\epsilon_{1}$ and starting with the OD (red filled circle) that the state disappears with the emergence of the OS (black filled square), while during the backward variation, the OS destabilizes to give birth to the IC (red filled triangle) (this confirms the existence of the multistability region R2), and eventually back to the OD. Choosing the IC as the initial state shows during the forward variation of $\epsilon_{1}$ a destabilization to give rise to the OS and during the backward variation a destabilization to give rise to the IC and eventually to the OD (this last behavior confirms the existence of the region R1). This establishes R3 as the region in which the OD, the IC and the OS coexist. In panel (b), we see that when the OD becomes unstable with the increase of $\epsilon_{1}$, one observes a saddle-node bifurcation (SN) to the OS: at the bifurcation point, one stable branch corresponding to the OD (red filled circle) and one unstable branch (blue open circle, represent saddles) annihilate each other. When the OS becomes unstable with the decrease of $\epsilon_{1}$, it bifurcates (i) first as a Hopf bifurcation (HB) to the IC, and (ii) then as a saddle-node limit-cycle bifurcation to the IC. The IC so obtained undergoes with further decrease of $\epsilon_{1}$ a pitchfork bifurcation (PF) to the OD.

We now obtain the temporal behaviour of $r(t)$ for all the four multistability regions while starting from different initial states; the data are obtained from numerical integration of the dynamics (3). Figure 8(a) shows that the dynamics initiated with either the IC or the OS remaining stabilized at these states (the inset shows that both these states remain stabilized at long times) when $\epsilon_{1}$ and $\epsilon_{2}$ have values in the R4 region of the bifurcation diagram. Panel (b) shows that the dynamics initialized with either the IC or the OS or the OD remaining stabilized at these states when $\epsilon_{1}$ and $\epsilon_{2}$ have values in the R3 region of the bifurcation diagram. Panels (c) and (d) represent in a similar manner the expected behavior in the multistability regions R1 and R2, respectively.

We have until now analysed the bifurcation scenario in the $(\epsilon_{1},\epsilon_{2})$-plane for $\alpha=1.3$. To illustrate the effect of varying $\alpha$, we verified qualitatively similar bifurcation diagrams for four different values of $\alpha$, namely, $\alpha=0.1,~{}0.5,~{}1.0$ and $1.5$, shown in Fig. 9, panels (a), (b), (c), (d), respectively. The regions of the different states as well as the boundaries are obtained in the same manner as in Fig. 3. The bifurcation diagram 9(a) for $\alpha=0.1$ is very similar to the one for the dynamics (3) in the absence of asymmetry in the interaction, i.e., with $\alpha=0$ Chandru:2020 . Even with a small increase in $\alpha$ value to $\alpha=0.5$, we see the emergence of regions R3 and R4, see panel (b). Moreover, we see that increase in $\alpha$ values leads to growth of the IC region. In Fig. 10, we report for the dynamics (3) the bifurcation diagram in $(\alpha,\epsilon_{1})$-plane for different values of $\epsilon_{2}$.

In this work, we studied a nontrivial generalization of the Sakaguchi-Kuramoto model of spontaneous collective synchronization, by considering an additional asymmetric interaction in the dynamics that breaks the phase-shift symmetry of the model. We reported the emergence of a very rich bifurcation diagram, including the existence of non-stationary synchronized states arising from destruction of stationary synchronized states. The bifurcation diagram shows existence of regions of two-state as well as three-state coexistence arising from asymmetric interaction in our model. Our results are based on numerical integration of the dynamical equations as well as an exact analysis of the dynamics by invoking the so-called Ott-Antonsen ansatz. It would be of immense interest to explore how inclusion of inertial effects Gupta-inertia in our model modifies the bifurcation diagram 3, an issue we are working on at the moment.

V.K.Chandrasekar and Shamik Gupta formulated the problem. M. Manoranjani carried out the simulations and calculations. All the authors discussed the results and drafted the manuscript.

M.M. wishes to thank SASTRA Deemed University for research funds and for extending infrastructure support to carry out this work. S.G. acknowledges support from the Science and Engineering Research Board (SERB), India under SERB-TARE scheme Grant No. TAR/2018/000023, SERB-MATRICS scheme Grant No. MTR/2019/000560, and SERB-CRG scheme Grant No. CRG/2020/000596. He also thanks ICTP – The Abdus Salam International Centre for Theoretical Physics, Trieste, Italy for support under its Regular Associateship scheme. The work of V.K.C. is supported by the SERB-DST-MATRICS Grant No. MTR/2018/000676 and DST-CRG Project under Grant No. CRG/2020/004353.

The data that support the findings of this study are available from the corresponding author upon reasonable request.