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

The celebrated Kuramoto model provides an analytically tractable framework to study spontaneous collective synchronization, and comprises globally-coupled limit-cycle oscillators interacting symmetrically with one another. The Sakaguchi-Kuramoto model is a generalization of the basic model that considers the presence of a phase-lag parameter in the interaction, thereby making it asymmetric between oscillator pairs. Here, we consider a further generalization, by adding an interaction that breaks the rotational symmetry of the model. The highlight of our study is the unveiling of a very rich phase diagram comprising both oscillatory and non-oscillatory synchronized phases as well as an incoherent phase: There are regions of two-phase as well as an interesting and hitherto unexplored three-phase coexistence arising from asymmetric interactions in our model.


Introduction and model of study
The fundamental paradigm of the Kuramoto model [1][2][3][4][5] 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 [6]. Examples of collective synchrony include synchronized firing of cardiac pacemaker cells [7], synchronous emission of light pulses by groups of fireflies [8], chirping of crickets [9], synchronization in ensembles of electrochemical oscillators [10], synchronization in human cerebral connectome [11], and synchronous clapping of audience [12]. 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 θ j ∈ [−π, π] the phase of the j-th oscillator, j = 1, 2, . . . , N , the dynamics of the model is described by a set of N coupled first-order nonlinear differential equations of the form where K ≥ 0 is the coupling constant, and ω j is the natural frequency of the j-th oscillator. The frequencies {ω j } denote a set of quenched-disordered random variables sampled independently from a distribution g(ω) usually taken to be unimodal, that is, symmetric about the mean ω 0 and decreasing monotonically and continuously to zero with increasing |ω − ω 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(θ k − θ 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 traveling wave states [13,14], glassy states and super-relaxation [15], and so forth, and has been invoked to discuss coupled circadian neurons [16], dynamic interactions [17], etc. A generalization of the Kuramoto model (1) that accounts for asymmetric interaction is the so-called Sakaguchi-Kuramoto (SK) model, with the dynamics described by the equation of motion [18] where 0 ≤ α < π/2 is the so-called phase lag parameter. It is now easily seen that owing to the presence of an α = 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 [19], time-delayed interactions [20], hierarchical populations of coupled oscillators [21], chaotic transients [22], dynamics of pulse-coupled oscillators [23], etc. Let us note that both the dynamics (1) and (2) satisfy rotational symmetry, whereby rotating every phase by an arbitrary angle same for all leaves the dynamics invariant. In particular, one may implement the transformation θ j (t) → θ j (t)+ω 0 t ∀ j, t, which is tantamount to viewing the dynamics in a frame rotating uniformly with frequency ω 0 with respect to an inertial frame, e.g., the laboratory frame.
In this work, we study a generalization of the SK dynamics (2) by including an interaction term that explicitly breaks the rotational symmetry of the dynamics. To this end, we consider the following set of N coupled nonlinear differential equations: where the real parameters 1,2 denote the coupling constants. Note that on putting 2 = 0 in Eq. (3), one recovers the dynamics of the SK model (2) on identifying 1 with the parameter K ≥ 0 in the latter, and so we take 1 to be positive. Here, we consider 2 to satisfy 2 ≥ 0. The dynamics (3) has rotational symmetry only with the choice 2 = 0, and so the 2 -term acts a rotational-symmetry-breaking interaction (With 2 = 0, the only case of rotational symmetry is with respect to the transformation θ j → θ j + π ∀ j.). In contrast to the SK model, the dynamics (3) is not invariant with respect to the transformation θ j → θ j + ω 0 t ∀ j, t owing to the presence of the 2 -term. Consequently, the mean ω 0 of the frequency distribution g(ω) 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 ω 0 with respect to the laboratory frame, as is possible with the SK model. We will consider in this work a unimodal g(ω). Specifically, we will consider a Lorentzian g(ω): Analysis of synchronization in the context of the Kuramoto class of models involves introducing a complexvalued order parameter given by [1] The quantity r(t); 0 ≤ r(t) ≤ 1, measures the amount of synchrony present in the system at time t, while ψ(t) ∈ [−π, π] gives the average phase. Considering the limit N → ∞, 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 phase. The two phases are characterized by the time-independent value that r(t) assumes at long times, namely, a zero value in the incoherent phase and a non-zero value in the synchronized phase. Considering the limit N → ∞, this work aims at a detailed characterization of the long-time (t → ∞) limit of the dynamics (3) and understanding of what new features with respect to the SK model are brought in by the introduction of the rotational-symmetry-breaking 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 phase diagram/bifurcation diagram in the ( 1 , 2 )-plane?
Recently, the dynamics of the Kuramoto model in presence of rotational-symmetry-breaking symmetric interaction was shown to unveil a rather rich phase diagram with existence of both stationary and non-stationary synchronized phases [24]; the model studied therein corresponds to the dynamics (3) with α 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 SK model (i.e., the case 2 = 0 in Eq. (3), and (ii) the Kuramoto model with an additional rotational-symmetry-breaking symmetric interaction (i.e., the case α = 0 in Eq. (3)). We focus on unimodal frequency distributions. As is well known, the SK model shows either an incoherent (r(t) vanishes as t → ∞) or a synchronized (r(t) as t → ∞ has a nonzero value, and moreover, r(t) does not oscillate as a function of time) state [18]. These states are observed only in a co-rotating frame rotating uniformly with frequency ω 0 with respect to the laboratory frame. As one tunes the strength of 1 , one has typically a continuous synchronization transition between an incoherent and a synchronized phase, though for certain specific unimodal frequency distributions and carefully chosen phase lag α, one may observe more than one nontrivial synchronization transitions [25,26]. 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 and 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/phase diagram in the ( 1 , 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 phase 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 phase diagram in the ( 1 , 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. 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-phase coexistence region between IC-OS-OD. We believe that this three-phase 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 rotational-symmetry-breaking symmetric interaction (i.e., the case α = 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-phase-coexistence region is observed; we may conclude therefore that introducing an asymmetric rotational-symmetry-breaking interaction in the SK model drastically and nontrivially modifies the phase diagram with respect to the SK model and with respect to the model studied in Ref. [24].
The rest of the paper is devoted to a derivation of the aforementioned results. For Lorentzian g(ω), Eq. (4), we use exact analytical results derived by applying the so-called Ott-Antonsen (OA) 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. 2. The celebrated OA ansatz is a remarkable method of analysis introduced to study coupled oscillator systems, which allows to rewrite in the limit N → ∞ the dynamics of coupled networks of phase oscillators in terms of a few collective variables [27,28].
The paper is organized as follows. In Section 2, we present our OA-ansatz-based analysis of the dynamics (3) and the existence and stability of all the different possible phases at long times. In Section 3, we discuss our analytical vis-à-vis numerical results. Finally, in Section 4, we draw our conclusions.

Analysis of the dynamics (3): The Ott-Antonsen (OA) ansatz
We now analyze the dynamics (3) in the limit N → ∞ by employing the OA ansatz. In this limit, the dynamics may be characterized by defining a single-oscillator distribution function f (θ, ω, t), which is such that f (θ, ω, t) is the probability density of finding an oscillator with natural frequency ω and phase θ at time t. The distribution is 2π-periodic in θ and is moreover normalized as 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): where v(θ, ω, t) ≡ dθ/dt is the angular velocity at posi- where z = z(t) is obtained as the N → ∞-limit generalization of Eq. (5): In their seminal papers [27,28], 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(ω), Eq. (4), may be found by using the ansatz Here, c.c. stands for a term obtained by complex conjugation of the first term within the brackets occurring in the sum, and with arbitrary a(ω, t) satisfying |a(ω, t)| < 1, together with the requirements that a(ω, t) may be analytically continued to the whole of the complex-ω plane, it has no singularities in the lower-half complex-ω plane, and |a(ω, t)| → 0 as Im(ω) → −∞.
Using the ansatz (10) in Eq. (7), we straightforwardly get ∂a ∂t +iωa+ where we have obtained by using in Eq. (9) the ansatz (10) and Eq. (4) for g(ω), then converting the integral therein into a contour integral, and finally evaluating the contour integral. Equation (12) may then be rewritten as The above equation expressed in terms of the quantities r and ψ gives the following two coupled equations: The above equations constitute the OA-ansatz-reduced order parameter dynamics corresponding to the dynamics (3) in the limit N → ∞ and for the case of the Lorentzian g(ω), Eq. (4). Let us remark that these equations are invariant under the transformation (r, ψ) → (r, ψ + π), which may be traced back to the fact that the original dynamics (3) is invariant under the transformation θ j → θ j + π ∀ j. Note that for 2 = 0 and α = 0, when one has the Kuramoto model, the two equations in (15) are decoupled, and the equations then allow for only uniform rotation of ψ with frequency ω 0 . The case of our interest, namely, 2 = 0 and α = 0, is analyzed 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 → ∞) limit as

Analysis of the OA-ansatz-based dynamics
Below we discuss the various states obtained in the long-time (t → ∞) limit of the dynamics (14), equivalently, Eq. (15).

Incoherent (IC) state
The incoherent (IC) state is characterized by time independent z satisfying z = z = 0 (thus representing a stationary state of the dynamics (15)); correspondingly, one has r = 0, and hence, R = 0. The linear stability of such a state is determined by linearizing Eq. (14) around z = 0, by writing z as z = u with |u| 1. We obtain The matrix M has eigenvalues ). Note that we have λ 1 > λ 2 . The stability threshold for the IC is then obtained by analyzing λ 1 as a function of 1 and 2 , and asking for a given 2 the particular value of 1 at which λ 1 vanishes. One thus obtains the stability threshold as here,

Oscillatory Synchronized (OS) state
The oscillatory synchronized (OS) state is characterized by z that is time dependent (thus characterizing a non-stationary state of the dynamics (15)). In this state, the order parameter r(t) oscillates as a function of t, but nevertheless yields a non-zero time-independent time average, R = 0. It is thus distinct from the nonoscillatory synchronized state (i.e., an oscillation death (OD) state, see below) that is also possible as a longtime solution of the dynamics (15). For the OS, one has z that is time independent (thus characterizing a stationary state of the dynamics (15)). and correspondingly, both r(t) and R have non-zero time-independent values, but the former does not oscillate as a function of time. Deriving stability conditions for the OS does not prove easy, unlike the IC and the OD. Hence, we analyzed using Eq. (15) the OS stability using the numerical package XPPAUT [29]. The analysis is pursued by expressing Eq. (15) in the Cartesian coordinates instead of the polar coordinates (r, ψ), as the solutions tend to go out of bound in the latter. Our analysis reveals that the OS represents a stable limit cycle in the (x, y)-plane. In Section 3, we will present XPPAUTgenerated phase-space portraits in the (x, y)-plane for not only the IC, the OS and the OD state, but also for regions in parameter space allowing for coexistence of two or more states. Let us remark in passing that owing to the invariance of the dynamics (15) under the transformation (r, ψ) → (r, ψ + π), it follows that the phase portrait would be symmetric under (x, y) ↔ (−x, −y).

Oscillation Death (OD) state
Considering the dynamics (15), we now explore the possibility of existence of the oscillation death (OD) state.
Requiring in the dynamics that r and ψ have timeindependent non-zero values so that the left hand side of the two equations may be set to zero, we obtain for the OD the two coupled equations The above equations give the following solutions for stationary r and ψ: In the next section, we view the above results in the light of those obtained from numerical integration of the dynamics (3). 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 (4) with γ = 0.3 and ω 0 = 0.6. Numerical integration is performed by employing a standard fourth-order Runge-Kutta integration scheme with integration step size dt = 0.01.
In Fig. 1, we have plotted the temporal behavior of the IC, the OS and the OD, for α = 1.3. The data are obtained from numerical integration of the dynamics (3). At fixed 2 = 2.0, one has the OD at 1 = 0.5, the IC at 1 = 1.5, and the OS at 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 → ∞ limit. In the OS, r(t) as t → ∞ oscillates in time. For the OD, r(t) as t → ∞ has a non-zero time-independent constant value.
In order to demonstrate the influence of α on the stability of the different phases, we now discuss the bifurcation diagram in the ( 1 , 2 )-plane. Figure 2 shows the stable phases in the ( 1 , 2 )-plane for the model (3) with α = 1.3. Here the phases 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 analyzing the long-time numerical solution of the dynamics (3) for N = 10 4 , namely, by asking at given values of 1 and 2 the quantity r(t) at long times represents which one of the three possible behavior depicted in Fig. 1. . The detailed procedure involves the following: For given γ, ω 0 and α, we vary at fixed 2 the value of 1 from low to high, and use the numerical solution of the Eq. (15) (rewritten in terms of x = x(t) and y = y(t)) to ascertain the particular value of 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 2 . Figure 2 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 phases) and between IC-OS, respectively. At a fixed 1 and on tuning 2 (or vice versa), one observes phase transitions/bifurcations as one crosses the different phase boundaries. When compared with the bifurcation diagram in case of symmetric interactions that either preserve or break rotation symmetry [24], we see that presence of asymmetry in the interaction leads to a very rich phase diagram. The dashed lines L1, L2, and L3 in Fig. 2 indicate the values of 2 for which the plots in Figs. 5 and 6, reported later in the paper, are obtained. Figure 2, the phase diagram or the bifurcation diagram of model (3), is the key result of our work.  Figure 3 shows further details of the phase diagram, and in particular, corresponding to the OA-ansatz-based dynamics (15), 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 OA-ansatzbased dynamics (15).
In Fig. 4, we show the XPPAUT-generated phasespace portraits in each of the coexistence regions R1, R2, R3, R4, based on the OA-ansatz-based dynamics (15). 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 phase diagram and the phase portraits in the three phases, 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. (16)) as a function of adiabatically-tuned 1 for representative values of 2 , with the aim to see signatures of bifurcation. We first let the system settle to the stationary state at a fixed value of 1 , and then tune 1 adiabatically in time from low to high values while recording the value of R in time; this corresponds to forward variation of 1 . Subsequently, we tune 1 adiabatically in time from high to low values (backward variation of 1 ). Adiabatic tuning ensures that the system remains in the stationary state at all times as the value of 1 changes in time. In Fig. 5, we consider two values of 2 : Panels (a) and (b) are for 2 = 2.0, while panels (c) and (d) are for 2 = 3.4. Here, the lines corre-spond to results based on numerical integration of the dynamics (3), while symbols correspond to the analysis of the OA-ansatz-based dynamics, Eq. (15), and are obtained as follows: -For the IC, the OA analysis predicts that the IC state with R = 0 exists only for 1 larger than its threshold value given by Eq. (20). 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. 2 which shows that 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. 2 in which we see that 2 = 3.4 allows for the existence of R2 region. The symbols in panels (b) and (d) represent XPPAUTgenerated bifurcation plots for y as a function of 1 with fixed 2 and based on the OA-ansatz-dynamics (15), 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 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 1 , one observes a saddlenode 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 sponding 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 In order to witness the R3 region, one has to choose the value of 2 carefully as R3 represents a very narrow region, see Fig. 2; we consider 2 = 3.1 to witness the R3 region, see Fig. 6. 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. 5. Here, we see the regions R1 and R2 as well. In panel (a), we see during the forward (adiabatic) variation of 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 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 0 0  Fig. 6 The plots correspond to model (3) with α = 1.3, and by considering for ω j 's the Lorentzian distribution (4) with γ = 0.3 and ω 0 = 0.6. The panels correspond to the value 2 = 3.1, represented as the straight line L2 in Fig. 2. Panel (a) shows as a function of adiabatically-tuned 1 the quantity R, namely, the time-averaged order parameter in the long-time limit. The symbols in panel (b) represent XPPAUTgenerated bifurcation plots for y as a function of 1 with fixed the OS coexist. In panel (b), we see that when the OD becomes unstable with the increase of 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 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 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 7 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 1 and 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 1 and 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 analyzed the bifurcation scenario in the ( 1 , 2 )-plane for α = 1.3. To illustrate the effect of varying α, we verified qualitatively similar bifurcation diagrams for four different values of α, namely, α = 0.1, 0.5, 1.0 and 1.5, shown in Fig. 8, panels (a), (b), (c), (d), respectively. The regions of the different phases as well as the phase boundaries are obtained in the same manner as in Fig. 2. The bifurcation diagram 8(a) for α = 0.1 is very similar to the one for the dynamics (3) in the absence of asymmetry in the interaction, i.e., with α = 0 [24]. Even with a small increase in α value to α = 0.5, we see the emergence of regions R3 The regions of the different phases as well as the phase boundaries are obtained in the same manner as in Fig. 2. The regions R1, R2, R3, and R4 represent multistability or coexistence between IC-OD, between OS-OD, between IC-OS-OD and between IC-OS, respectively. and R4, see panel (b). Moreover, we see that increase in α values leads to growth of the IC region. In Fig. 9, we report for the dynamics (3) the phase diagram in (α, 1 )-plane for different values of 2 .

Conclusions
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 rotational symmetry of the model. We reported the emergence of a very rich phase diagram, including the existence of non-stationary synchronized states arising from destruction of stationary synchronized states. The phase diagram shows existence of regions of two-phase as well as three-phase 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 [30]  We show here the various stable phases and their boundaries, with bifurcation behavior observed as one crosses the phase boundaries. The regions representing stable existence of the OD (Oscillation Death state), the IC (Incoherent state) and the OS (Oscillatory Synchronized state) have been constructed by analyzing the long-time numerical solution of the dynamics (3) for N = 10 4 , namely, by asking at given values of 1 the quantity r(t) at long times represents which one of the three possible behavior depicted in Fig. 1. The blue dot-dashed line, the black dashed line, and the maroon continuous line are stability boundaries of the IC, the OD and the OS, respectively, and have been obtained by using the OA-ansatz-based dynamics (15) and following the procedure of its analysis detailed in the text. The regions R1, R2, R3, and R4 represent multistability or coexistence between IC-OD, between OS-OD, between IC-OS-OD and between IC-OS, respectively.
in our model modifies the phase diagram 2, an issue we are working on at the moment.