We study a system of coupled phase oscillators near a saddle-node on invariant circle bifurcation and driven by random intrinsic frequencies. Under the variation of control parameters, the system undergoes a phase transition changing the qualitative properties of collective dynamics. Using Ott–Antonsen reduction and geometric techniques for ordinary differential equations, we identify heteroclinic bifurcation in a family of vector fields on a cylinder, which explains the change in collective dynamics. Specifically, we show that heteroclinic bifurcation separates two topologically distinct families of limit cycles: contractible limit cycles before bifurcation from noncontractibile ones after bifurcation. Both families are stable for the model at hand.

The Kuramoto model (KM) of coupled phase oscillators provides an important paradigm for studying collective dynamics in systems ranging from neuronal networks to swarms of fireflies to power grids. The classical KM features a remarkable phase transition separating stable mixing dynamics from gradual build-up of synchronization. During the latter phase, the oscillators form a cluster whose coherence (measured by the order parameter) remains approximately constant and increases with coupling strength. If uniformly rotating phase oscillators in the KM are replaced by those close to a saddle-node on invariant circle bifurcation, the order parameter does not stay constant anymore. Instead, it undergoes slow-fast oscillations. Furthermore, for larger values of coupling strength, the system undergoes phase transition, which changes the character of oscillations qualitatively. Previous studies based in part on numerical bifurcation techniques revealed a rich bifurcation structure of the modified KM. In this article, we use Ott–Antonsen reduction and qualitative methods for ordinary differential equations to study collective dynamics in the modified model with an emphasis on slow-fast oscillations of the order parameter. We analytically locate Andronov–Hopf and Bogdanov–Takens bifurcations for the system of equations governing the order parameter and identify the relevant normal form for Bogdanov–Takens bifurcations. Furthermore, we relate the phase transition in the modified model to nonlocal bifurcation in one parameter families of vector fields on a cylinder. The results of this work show that the slow-fast oscillations of the order parameter in the modified KM are shaped by the proximity to both Bogdanov–Takens and heteroclinic bifurcations.

## I. INTRODUCTION

The KM plays a special role in theory synchronization. It provides a framework for studying synchronization and other forms of collective dynamics in systems of coupled oscillators with random parameters. Despite its analytical simplicity, studies of the KM revealed and helped to understand some very nontrivial phenomena in collective dynamics, which are relevant to a range of models in physics and biology (e.g., the onset of synchronization and chimera states^{22,23,38}). Motivated by models featuring type I excitability in mathematical biology, we modify the KM by placing the individual oscillators close to a saddle node on an invariant circle bifurcation. Specifically, we consider the following coupled system:

where $\theta i:R+\u2192T:=R/2\pi R$ is the state of oscillator $i$ at time $t$ and $K\u22650$ is coupling strength. Parameter $\omega i$ controls the frequency of oscillator $i$, if $\omega i>0$, or defines the excitation threshold otherwise. We assume that $\omega i$s are sampled from a unimodal probability distribution with density $g(\omega )$. Equation (1.1) fits into the framework of the coupled active rotators model considered by Shinomoto and Kuramoto in Ref. 37. In contrast to our setting, they used identical oscillators (i.e., $\omega i=const$) albeit forced by small noise. Different variants of the active rotators model with and without noise were studied more recently by other authors.^{1,20,21,42,44} We comment on the relation of our findings to the results in these articles below.

To take a full advantage of the Ott–Antonsen Ansatz^{33} below, we restrict to the Lorentzian distribution with density

where $\u03f52$ and $\delta >0$ are the location and scale parameters, respectively. Lorentzian distribution is commonly used in the studies of collective dynamics in the Kuramoto model and related systems because it fits nicely into the Ott–Antonsen Ansatz (see, e.g., Refs. 32 and 34). There has been a concern that due to its special properties (lack of finite moments), models based on Lorentzian distribution may feature nongeneric scenarios.^{26} We checked that in qualitative form, the bifurcation scenarios reported in this article hold for Gaussian distribution and, thus, are relevant to a large class of models based on unimodal probability distributions. If typical realizations of $\omega i$ are $O(1)$ and positive, then the dynamics of (1.1) is practically the same as in the classical KM. Specifically, for small $K$, the oscillators are highly incoherent. Starting with a certain critical value of $K$, the coherence gradually increases.^{38} Thus, in this article, we focus on the regime when typical values of $\omega i\u2032s$ are small, i.e., when the individual oscillators are close to saddle-node bifurcation.

To describe the collective dynamics of (1.1), we need to remind the reader the definition of the Kuramoto’s order parameter,

The modulus, $\rho n=|hn|$, and the argument, $\varphi n=arg\u2061hn$, of the complex order parameter yield the degree of coherence and the position of the center of mass of the population of oscillators on the unit circle, respectively. The combination of these two quantitates as functions of time provides a good description of collective dynamics. In the classical Kuramoto, after some transients, the modulus of the order parameter approaches a steady value [up to $O(n\u22121/2)$ fluctuations] while the argument drifts with approximately constant velocity. In contrast, in the simulations of the modified KM (1.1), we observe very nonuniform slow-fast oscillations in the modulus and the argument of the complex order parameter (Fig. 1).

We now turn to describe the salient features of the collective dynamics of (1.1). There are three main regimes in the system dynamics:

For small values of $K$, the individual oscillators rotate in the counterclockwise direction, but the overall distribution of oscillators on a unit circle remains practically stationary. The distribution has a peak at $\theta =0$ reflecting the fact that the oscillators slow down in a neighborhood of $\theta =0$ and, thus, spend more time there [see Fig. 2(a)]. The order parameter lies on the real axis very close to $1$ and stays approximately constant [see inset in Fig. 2(a)].

There are two critical values of $K$: $0<KAH<KHC$, which will be shown below to correspond to Andronov–Hopf bifurcation and heteroclinic bifurcation, respectively. For $K\u2208(KAH,KHC)$ the order parameter undergoes slow-fast oscillations (see Fig. 1). Importantly, the argument of the order parameter stays between $\u2212\pi 2$ and $\pi 2$ [see the insets in Figs. 2(b) and 2(c)]. This means that the center of mass of the population of oscillators oscillates around $\theta =0$. The modulus of the order parameter varies between values close to $0$ and $1$ spending more time in the region near $1$ (Fig. 1). Two snapshots of the distribution of oscillators on a unit circle are shown in Fig. 2(b) (coherent phase) and Fig. 2(c) (incoherent phase).

Two representative episodes of the system dynamics for larger values of $K$ are shown in Figs. 2(d) and 2(e). As before, the oscillators slow down and accumulate as they approach the origin and accelerate and spread around once they have passed it. However, there is an important distinction from the previous regime. It is quite pronounced in numerical simulations and can also be seen from static snapshots in Fig. 2. Recall that for small values of $K$ ($K<KHC$), the order parameter never leaves the right half plane. For larger $K$ ($K>KHC$), on the other hand, the order parameter makes a full revolution around the origin in one cycle of oscillations [see the insets in Figs. 2(d) and 2(e)]. The change in the qualitative character of oscillations can be seen from the timeseries of $\varphi n$ for two different values of $K$ in Fig. 1(b). Note that for small $K$, $\varphi n$ undergoes small oscillations, while for large values of $K$, the range of $\varphi n$ covers the entire circle. One can compare this to small oscillations vs full swing revolutions of the pendulum. Below, we show that indeed the two regimes are qualitatively (topologically) distinct and are separated by bifurcation.

In contrast to the classical scenario of transition to synchronization in the original KM, where one observes the formation of a single coherent cluster drifting uniformly around the unit circle,^{23,38} in the modified model (1.1), we see pronounced oscillations in the argument of the order parameter [Figs. 2(b)–2(e)]. The amplitude of these oscillations increases until the oscillations are transformed into a rotational motion of the center of mass of the population of oscillators around the unit circle [Figs. 2(d) and 2(e)]. Another salient feature of this transition is the slow-fast character of the oscillations of the order parameter (Fig. 1).

Periodic regimes in macroscopic dynamics in the model of active rotators were described already in Ref. 37. These regimes were found by numerical simulations. Pulsating oscillations of the order parameter similar to those shown in Fig. 1 were reported in subsequent studies of the coupled active rotators with and without noise.^{20,42,44} In Refs. 20 and 42, the authors derived a system of differential equations for the modulus $\rho $ and the argument $\varphi $ of the order parameter (1.3) in the limit as $n\u2192\u221e$. Using the combination of self-consistent analysis and numerical bifurcation techniques, both studies revealed an Andronov–Hopf, Bogdanov–Takens, and homoclinic bifurcations in the parameter regime relevant to transition to synchronization. Similar results were obtained for a coupled system of noisy rotators in Ref. 44 albeit via a different approximation procedure. Previous studies reveal a rich bifurcation structure underlying macroscopic dynamics in the coupled model. In this article, we focus on the origins of slow-fast oscillations of order parameters and, in particular, on transition from the oscillatory to rotational motion of the order parameter. We use the qualitative methods for ordinary differential equations to elucidate the bifurcation structure of the coupled model. (Note that our setting is slightly different from those of models in Refs. 20, 42, and 44.) Specifically, we locate Andronov–Hopf and Bogdanov–Takens bifurcation analytically and identify the relevant normal form for the Bogdanov–Takens bifurcation. Furthermore, we relate the transition to the rotational motion to a heteroclinic bifurcation of a family of vector fields on a cylinder. A related global bifurcation separating oscillations from rotations of the order parameter in the model of noisy active rotators was reported (but not analyzed) in Ref. 44. In Sec. IV, we provide a detailed analysis of the heteroclinic bifurcation for the model at hand.

There is an extensive literature on the KM. Early articles were devoted mostly to synchronization (cf. Refs. 23 and 38–40). The scope of more recent contributions encompasses rigorous mathematical aspects of synchronization,^{6,7,10} complex spatiotemporal patterns,^{31,35} generalizations or adaptations of the classical KM,^{8,9,27} as well as applications to physical and biochemical systems (see, e.g., Refs. 12, 30, and 36). Our work is close in spirit to Refs. 9 and 43 where the modifications of the KM were used to tackle challenging questions about collective dynamics. Specifically, we wanted to understand collective dynamics in populations of type I excitable oscillators with randomly distributed intrinsic frequencies. Coupled networks of this type have been studied in computational neuroscience (see, e.g., Refs. 14, 15, and 29) albeit in settings that differ from that adopted here. In particular, methods for studying weakly coupled networks (cf. Ref. 17) do not apply to systems with random parameters like (1.1). For coupled systems with different coupling types, collective dynamics in closely related models of coupled theta neurons was analyzed in Refs. 19, 28, and 32. We believe that the analysis of the modified KM used in this article complements the previous studies of coupled active rotators in Refs. 20, 42, and 44 and brings new insights into understanding collective dynamics of type I oscillators.

The article is organized as follows. In Sec. II, we use the Ott–Antonsen Ansatz^{33} to derive a system of two ordinary differential equations, which captures the long time dynamics of the coupled system. In Sec. III, we analyze the reduced system. The analysis uses the unfolding of Bogdanov–Takens bifurcation among other qualitative techniques for ordinary differential equations. Furthermore, we identify heteroclinic bifurcation, which explains phase transition in collective dynamics. The heteroclinic bifurcation is analyzed in Sec. IV. In Sec. V, we relate the analysis of the reduced system to the collective dynamics of (1.1). We conclude with a brief discussion of the main results in Sec. VI.

## II. THE MEAN FIELD LIMIT

In the large $n$ limit, the dynamics of (1.1) is described by the following Vlasov equation (cf. Refs. 11 and 39)

where $f(t,\theta ,\omega )d\theta d\omega $ is the probability of $(\theta ,\omega )\u2208(\theta ,\theta +d\theta )\xd7(\omega ,\omega +d\omega )$ at time $t\u2208R+$. The velocity field is defined by

where

is the continuum order parameter.

From (1.1) with $K=0$, one can determine the stationary distribution of oscillators in phase space,

where $j(\omega )$ is subject to the following equation:

Previous studies of the KM suggest that the Ott–Antonsen Ansatz correctly describes the long time asymptotic behavior of solutions of this model. It has been used to study chimera states^{31} among many other spatiotemporal patterns in the KM and related models.^{20,27,32,34,42} There is a convincing albeit not completely rigorous mathematical argument justifying the use of this Ansatz for studying the long time behavior of solutions of the Kuramoto model. The analysis in the remainder of this article relies on the validity of the Ott–Antonsen Ansatz.

Multiplying both sides by $2$ and rescaling time in (2.5) yields

Equation (2.6) is simpler than the Vlasov equation (2.1), but it is still an integro-partial differential equation. For the Lorentzian density $g(\omega )$ given in (1.2), one can further reduce (2.6) to the following ordinary differential equation (cf. Ref. 33):

where by abuse of notation, we continue to denote $\alpha (t,\u03f52\u2212i\delta )$ by $\alpha (t)$. Furthermore, the order parameter $h$ can now be expressed through $\alpha $,

Equation (2.7) can be written as an ordinary differential equation on $R2$,

where $\alpha =x+iy$. Since we are interested in the macroscopic dynamics of (1.1), we want to track the changes the qualitative behavior of the modulus and the argument of the order parameter. To this end, we rewrite (2.7) in polar coordinates $\alpha =\rho ei\varphi $,

The polar coordinate transformation blows up the origin in the $x$–$y$ plane to a $\rho =0$ circle. This singular change preserves a $1$–$1$ correspondence between the trajectories of (2.9) and (2.10) lying in $R2/{0}$ and $R+\xd7T$, respectively (cf. Ref. 4, Sec. 2.8). The two systems are topologically equivalent in these domains. The right-hand side of (2.10) has a singularity at $\rho =0$. To resolve this singularity, we multiply both sides of the equations in (2.10) by $\rho $ and rescale time to obtain

Systems (2.10) and (2.11) are topologically equivalent on $R+\xd7T$. However, the latter defines a smooth vector field over $R\xd7T$. It can be studied by standard methods of the qualitative theory for differential equations. In particular, this allows us to resolve the singularity of (2.10) at $\rho =0$, which is important for understanding the macroscopic dynamics of the coupled system.

## III. QUALITATIVE ANALYSIS OF THE PLANAR SYSTEM

In this section, we analyze (2.11) when $K,\delta ,$ and $\u03f5$ are nonnegative and small.

There are two fixed points of (2.11),

Additionally, when $\delta =\u03f5=0$, we can explicitly calculate a third fixed point at $(1,0)$. We begin with the local analysis around $(1,0)$. This is followed by linearization about the other two fixed points. After that, we combine this information and identify a nonlocal bifurcation of the vector field responsible for the transformation of collective dynamics in (1.1).

### A. The Bogdanov–Takens singularity

Let $\rho =1\u2212r$ and rewrite (2.11) in the neighborhood of $(1,0)$ keeping terms up to and including order $2$ to obtain

With $\u03f5=\delta =K=0$, (3.2) fits into the normal form of Bogdanov–Takens bifurcation,^{2,5,41}

Rewriting (3.3) in polar coordinates $r=acos\u2061\psi ,\varphi =asin\u2061\psi $, we obtain

Next, we keep $\delta =K=0$ and let $0<\u03f5\u226a1$. By increasing $\u03f5$ from $0$, we observe that the fixed point at the origin bifurcates into two fixed points: $(0,2\u03f5)$ and $(0,\u22122\u03f5)$. Furthermore, the system has an integral (cf. Ref. 16)

Thus, the trajectories are given by two families of circles,

for $c\u2265\u03f5$ and $c\u2264\u2212\u03f5$. As $c2\u2192\u221e$, they limit onto the $r$-axis [Fig. 3(b)].

Next suppose $\u03f5,\delta ,$ and $K$ are positive and small. To study this case, we use the following scaling:

Treating $\nu $ and $\mu $ as positive and small parameters of the same order, we locate the fixed point of (3.7)$(r\xaf,\varphi \xaf)$ with positive $r\xaf$,

Linearization about $(r\xaf,\varphi \xaf)$ yields

Note

The determinant of $A$ is positive for sufficiently small $\mu \u22650$ [recall that $r\xaf=O(1)$]. For $0\u2264\mu <2\nu $, the fixed point $(r\xaf,\varphi \xaf)$ is a stable focus. It becomes unstable for $\mu >2\nu $. At $\mu =2\nu $, the system undergoes a supercritical Andronov–Hopf bifurcation [see Fig. 3(d)].

### B. The heteroclinic bifurcation

We now turn to the remaining two fixed points: $S1=(0,\u2212\pi /2)$ and $S2=(0,\pi /2)$. Linearization about $S1$ yields

The eigenvalues are $\lambda 1=1$ and $\lambda 2=\u22121$ with the corresponding eigenvectors $v1=(1,0)$ and $v2=(1,1+\u03f52)$.

Similarly, linearization about $S2$ yields

The eigenvalues are $\lambda 1=1$ and $\lambda 2=\u22121$ with corresponding eigenvectors $v1=(1,\u22121\u2212\u03f52)$ and $v2=(1,0)$. Consider $D=[\u22121,1]\xd7T$. In the parameter regime of interest, $D$ is positively invariant. Denote the branches of stable and unstable manifolds of $S1$ and $S2$ lying in $D$ by $Ws(S1,2)$ and $Wu(S1,2)$, respectively.

In the remainder of this section, we assume that $0<\u03f5\u226a1$ is fixed. Both $K$ and $\delta $ are also nonnegative and $O(\u03f5)$. Furthermore, we fix $\delta $ and treat $K$ as a control parameter. The linearization about $S1$ and $S2$ shows that both are saddles. The unstable manifold of $S1$ and the stable manifold of $S2$ lie in ${\rho =0}$ (Fig. 4). The unstable manifold of $S2$ and the stable manifold of $S1$ are tangent to $(1,\u22121\u2212\u03f52)$ and $(1,1+\u03f52)$, respectively. Recall that there is another fixed point in $D$: $O\u2248(1\u2212\u03f5,\delta \u2212K)$. For $K>2\delta $, $O$ is an unstable focus. Denote the stable limit cycle born at Andronov–Hopf bifurcation by $PK\u2212$. Note that $Wu(O)$ and $Wu(S2)$ limit onto $PK\u2212$ [Fig. 4(a)].

As $K$ is increased from the Andronov–Hopf bifurcation, $KAH$, the following transformations of the phase portrait take place. The periodic orbit $PK\u2212$ grows in size. $Wu(S2)$ moves up while $Ws(S1)$ moves down [Fig. 4(a)]. At $K=KHC$, they intersect forming a heteroclinic orbit $\Gamma 0$ connecting $S2$ and $S1$ [Fig. 4(b)]. After heteroclinic bifurcation, the limit cycle born at Andronov–Hopf bifurcation disappears blending into a heteroclinic loop. A new limit cycle, $PK+$, is born at $K=KHC+0$ [Fig. 4(c)]. In contrast to $PK\u2212$, $PK+$ is noncontractible. Thus, the heteroclinic bifurcation produces a limit cycle on each side of the bifurcation, i.e., both at $K=KHC\u22120$ and $K=KHC+0$. The bifurcation diagram in Fig. 5 summarizes this information.

To get a first insight into the heteroclinic bifurcation, we numerically computed a Poincaré map for (2.11). Specifically, for values of $K$ near $KHC$, we set up a cross section, $\Sigma $, intersecting the heteroclinic orbit [see a dashed line in plots (a) and (b) in Fig. 6]. Then, we sampled initial conditions from the appropriate region of $\Sigma $ (around the point of intersection with $\Gamma 0$) and followed the corresponding trajectories until their first return to the chosen region of $\Sigma $ [see trajectories plotted in red of Figs. 6(a) and 6(b)]. Note that after bifurcation, the red trajectory makes a full revolution around the cylinder before returning back to $\Sigma $ [see Fig. 6(b)]. The Poincaré maps computed before and after bifurcation are practically identical [see Figs. 6(c) and 6(d)]. Furthermore, the map shows a strong contraction of the vector field near heteroclinic loops. The periodic orbits corresponding to the fixed points of the map shown in Figs. 6(c) and 6(d) are shown in Figs. 6(e) and 6(f), respectively. The two limit cycles approach heteroclinic loops $\Gamma 0\u22c3\Gamma \u2212\xaf$ and $\Gamma 0\u22c3\Gamma +\xaf$ in Hausdorff distance as $K\u2192KHC\u22120$ and $K\u2192KHC+0$, respectively (see Fig. 7). Both limit cycles are stable. Furthermore, numerically computed Poincaré maps do not show any effect of heteroclinic bifurcation. These counterintuitive observations will be explained in Sec. IV by the analysis of the nonlocal bifurcation taking place in (2.11).

## IV. THE BOA CONSTRICTOR BIFURCATION

The main result of this section is given in the following theorem.

*7)*.

Let $\lambda 1,2s:=\lambda 1,2s(0)<0<\lambda 1,2u:=\lambda 1,2u(0)$ be the eigenvalues of the Jacobian $Df(S1,2,0)$. Recall that $\sigma 1,2=|\lambda 1,2s|/\lambda 1,2u$ is called the saddle number. Furthermore, we assume

$\sigma :=\sigma 1\sigma 2\u22601$; and

$\beta \u2032(0)\u22600$, where $\beta (\alpha )$ is a suitably defined split function

*(*see below*)*.Then, for sufficiently small $|\alpha |$, there exists a small neighborhood ofwhich contains a unique limit cycle $P\beta $ bifurcating from $\Gamma $. Moreover, $P\beta $ is contractible for $\beta >0$ and noncontractible for $\beta <0$. It is stable if $\sigma >1$ and unstable otherwise.$\Gamma :=\Gamma \u2212\u22c3\Gamma 0\u22c3\Gamma +\xaf,$

The proof employs a standard scheme for analyzing global bifurcations, which goes back to the proof of the Andronov–Leontovich Theorem (cf. Ref. 24, Theorem 6.1).

If $\sigma =1$ as in the case of (2.11), the stability of bifurcating orbits is determined by positive coefficients $a$, $a\u2212$, and $a+$. Specifically, $P\beta $ is stable if $a(a\u2212)\sigma 2<1$ for $\beta >0$ and $a(a+)\sigma 2$ for $\beta <0$. In this case, the stability of the orbit is determined by the degree of contraction produced by global maps defined by the flow along the heteroclinic loops rather than by that of the local near-to-saddle maps. While obtaining analytical estimates of the contraction of global maps requires tedious calculations, numerical first return maps in Figs. 6(c) and 6(d) unequivocally demonstrate strong contraction.

## V. COLLECTIVE DYNAMICS

Having analyzed the reduced system (2.11), we now return to the description of collective dynamics of (1.1). The analysis in Secs. III and IV shows that the heteroclinic bifurcation separates two topologically distinct families of the limit cycles of (2.11) $PK$ for $K<KHC$ and $K>KHC$. In either case, $PK$ consists of a segment lying in $Vc={0\u2264\rho <c}$ for some $0<c\u226a1$ and an arc outside $Vc$ (see Fig. 7). Thus, in one cycle of oscillations, the composition of the population of oscillators changes from high incoherence ($\rho \u22480$) to high coherence ($\rho \u22481$).

To estimate the duration of each of these phases and the period of oscillations, note that the incoherent phase is determined by the time the periodic trajectory spends near the saddles, which can be easily estimated from (4.2). The trajectory of (4.2) starting from $(\beta ,1)$ leaves the strip $|\xi |\u22641$ after time

Furthermore, at the time of exit

Thus, the trajectory enters the neighborhood of $S2$ with $\eta =O(|\beta |\sigma 1).$ From this, we can estimate the time it spends in the vicinity of $S2$,

[see Fig. 10(a)].

After going back to the original time, for oscillations in (2.10), we obtain

and similarly

Thus, in the original time, the period of oscillations remains finite [see Fig. 10(b)]. Furthermore, since the points on $Ws(S1)$ hit $S1$ in finite time, the heteroclinic orbit reaches the saddle in finite time too. Recall that $\rho =0$ is the image of the origin in $R2$ under polar coordinate transformation. The origin is not a fixed point of (2.9). Thus, the heteroclinic orbit of (2.11) for $K=KHC$ corresponds to a periodic orbit of (2.9) passing through the origin in the Cartesian coordinates (Fig. 12). Thus, the family of periodic orbits $PK,K>KAH$ of (2.11) corresponds to a family of periodic orbits $ZK$ of (2.9). The period of $ZK$ achieves its maximum at $K=KHC$ [see Fig. 10(b)]. The discrepancy in different locations of the maxima in plots (a) and (b) in Fig. 10 is explained by the fact that we dropped higher order terms in (3.2).

Having described the “incoherent” portion of the periodic orbit lying in a neighborhood of ${\rho =0}$, we now turn to the complimentary portion lying in a neighborhood of $\Gamma 0$ (see Fig. 11). First, note that outside a small neighborhood of ${\rho =0}$, the rescaling of time no longer has a qualitative impact on the system’s dynamics. Furthermore, on $\rho =1$, the vector field of (2.11) is pointed downward [see the $\rho $-equation in (2.11)]. On the other hand, the unstable focus $O$ has not moved too far from ${\rho =1}$, where it emerged at Bogdanov–Takens bifurcation (see Fig. 11). Thus, $\Gamma 0$ has to pass in a region between ${\rho =1}$ and the unstable focus $O$ (Fig. 11). This forces $\Gamma 0$ to pass close to $O$ where the vector field is very weak. Consequently, a significant portion of the period is spent near ${\rho =1}$, which correspond to the “coherent” phase of oscillations. This results in an interesting scenario, for which the separation of the timescales in oscillations of the order parameter just before and after heteroclinic bifurcation (Fig. 1) is not due to the proximity of heteroclinic bifurcation, as one would be tempted to assume, but rather to the proximity to Bogdanov–Takens bifurcation. Thus, both Bogdanov–Takens and heteroclinic bifurcations have an impact on the qualitative features of the collective dynamics.

We conclude with several remarks on the relation between (2.9) and (2.11). Recall that we switched to polar coordinates to be able to track the modulus and the argument of the order parameter, which capture macroscopic dynamics. The polar coordinate transformation is a diffeomorphism of $R2/{(0,0)}$ to $R+\xd7T$, but it is not a bijection at the origin. Consequently, (2.9) and (2.11) are topologically equivalent only when restricted to $R2/{(0,0)}$ and $R+\xd7T$, respectively. In particular, Andronov–Hopf and Bogdanov–Takens bifurcations describing local transformations of the vector field (2.11) in closed domains of $R+\xd7T$ translate automatically to corresponding bifurcations of (2.9). On the other hand, the heteroclinic bifurcation, which we analyzed for (2.11) involves the set $\rho =0$, which lies outside $R+\xd7T$. The heteroclinic orbit connecting $\Gamma 0$ (Fig. 7) corresponds to a periodic orbit of (2.9) passing through the origin (Fig. 12). This not a bifurcation of (2.9) as a vector field on $R2$, but it is a bifurcation on $R2/{(0,0)}$. This is a bifurcation because a contractibe periodic orbit before hitting the origin becomes a noncontractibe one after this event. This is border collision bifurcation. Clearly, this bifurcation is a consequence of our using polar coordinate transformation, which is singular at the origin. Nonetheless, both heteroclinic and border collision bifurcations are relevant in the context of the macroscopic dynamics. While in the Cartesian coordinates, the passing of the periodic orbit of (2.9) through the origin is a regular event, it represents a transition point in the description of the macroscopic dynamics. This is the point where oscillations of the center of mass of the population of oscillators are transformed into rotations. At this point, the amplitude and the period of the oscillations of the order parameter reach their respective maximal values, while the modulus of the order parameter reaches its minimal value. Note that the fast dips in the modulus of the order parameter can get arbitrarily close to $0$ provided that $n$ is large enough [Fig. 7(a)]. At the transition point, the oscillators are most dispersed as they undergo fast transitions between their successive stays near $\rho =1$ [Fig. 7(a)]. Therefore, the use of the polar coordinates in the analysis of (2.9) and the analysis of heteroclinic bifurcation in the transformed system are essential for understanding macroscopic dynamics of the coupled system (1.1).

## VI. DISCUSSION

In the present article, we analyzed a modified KM with individual oscillators in the regime near a saddle-node on an invariant circle bifurcation. The modified model features a new type of collective dynamics with alternating phases of high and low coherence. Furthermore, the order parameter exhibits slow-fast oscillations, which reveal a pronounced separation of timescales in collective dynamics. For the most part of the period, the value of the order parameter is close to 1 corresponding to the coherent phase. The long periods of coherence are punctuated by brief intervals of highly incoherent collective dynamics. The Ott–Antonsen reduction and a careful analysis of the reduced system show that the salient features of the collective dynamics are explained by the model’s proximity to local Bogdanov–Takens bifurcation and nonlocal heteroclinic bifurcation. In contrast to previously studied cases, where one limit cycle or two limit cycles of opposite stability appear in a bifurcation of a homoclinic/heteroclinic contour (cf. Refs. 3, 13, 24, and 25), heteroclinic bifurcation for the system at hand generates limit cycles on each side of bifurcation. These limit cycles are topologically distinct (contractible vs noncontractible) and are either both stable or both unstable.

## ACKNOWLEDGMENTS

This work grew out of A.P.’s Research Co-op at Drexel University. G.S.M. and A.P. were supported in part by the National Science Foundation (NSF) under Grant No. DMS 2009233 (to G.S.M.). M.S.M. was supported by a Support of Scholarly Activities Grant at The College of New Jersey.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Georgi S. Medvedev:** Writing – original draft (equal). **Matthew S. Mizuhara:** Writing – original draft (equal). **Andrew Phillips:** Writing – original draft (equal).

## DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

## REFERENCES

*Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector Fields*, Applied Mathematical Sciences Vol. 42 (Springer-Verlag, New York, 1990).

*Weakly Connected Neural Networks*, Applied Mathematical Sciences Vol. 126 (Springer-Verlag, New York, 1997).

*Nonlocal Bifurcations*, Mathematical Surveys and Monographs Vol. 66 (American Mathematical Society, Providence, RI, 1999).

*International Symposium on Mathematical Problems in Theoretical Physics*(Kyoto University, Kyoto, 1975); Lecture Notes in Physics Vol. 39 (Springer, Berlin, 1975), pp. 420–422.

*Elements of Applied Bifurcation Theory*, 3rd ed., Applied Mathematical Sciences Vol. 112 (Springer-Verlag, New York, 2004).