Networks of coupled dynamical units give rise to collective dynamics such as the synchronization of oscillators or neurons in the brain. The ability of the network to adapt coupling strengths between units in accordance with their activity arises naturally in a variety of contexts, including neural plasticity in the brain, and adds an additional layer of complexity: the dynamics on the nodes influence the dynamics of the network and vice versa. We study a minimal model of Kuramoto phase oscillators including a general adaptive learning rule with three parameters (strength of adaptivity, adaptivity offset, adaptivity shift), mimicking learning paradigms based on spike-time-dependent plasticity. Importantly, the strength of adaptivity allows to tune the system away from the limit of the classical Kuramoto model, corresponding to stationary coupling strengths and no adaptation and, thus, to systematically study the impact of adaptivity on the collective dynamics. We carry out a detailed bifurcation analysis for the minimal model consisting of $N=2$ oscillators. The non-adaptive Kuramoto model exhibits very simple dynamic behavior, drift, or frequency-locking; but once the strength of adaptivity exceeds a critical threshold non-trivial bifurcation structures unravel: A symmetric adaptation rule results in multi-stability and bifurcation scenarios, and an asymmetric adaptation rule generates even more intriguing and rich dynamics, including a period-doubling cascade to chaos as well as oscillations displaying features of both librations and rotations simultaneously. Generally, adaptation improves the synchronizability of the oscillators. Finally, we also numerically investigate a larger system consisting of $N=50$ oscillators and compare the resulting dynamics with the case of $N=2$ oscillators.

## I. INTRODUCTION

The synchronization of coupled oscillators is a fascinating manifestation of self-organization—indeed, self-emergent synchronization is a central process to a spectacular range of natural and technological systems, including the beating of the heart,^{16} flashing fireflies,^{17} pedestrians on a bridge locking their gait,^{18} genetic clocks,^{19} pendulum clocks hanging on a beam,^{20} mechanical oscillators,^{21,22} superconducting Josephson junctions,^{23} chemical oscillations,^{24,25} metabolic oscillations in yeast cells,^{26} life cycles of phytoplankton,^{27} and networks of neurons in the brain.^{3,4,28}

A desirable property in real-world oscillator networks is the presence of synchronization whenever it ensures the proper functioning of the network: in the realm of technology, for instance, the AC between generators and consumers in a power grid needs to stay synchronized to ensure ideal power transmission and ultimately avoid power blackouts;^{5} and wireless networks require synchronized clocks^{29} to ensure safe data transmission; in the realm of biology, the proper functioning of the heart requires that the rhythmic electric activation of cardiac cells stays coordinated.^{16}

While the network interaction between dynamic units may give rise to intriguing collective behaviors such as synchronization, adding the ability to adapt the (coupling) weights on a network according to the dynamics on its nodes leads to co-evolutionary network dynamics.^{30} Indeed—“Intelligence is the ability to adapt to change”^{31}—the adaptive dynamics of co-evolutionary networks may allow to increase their functional robustness.^{32,33} Adaptive co-evolutionary networks appear in a wide range of systems, including the vascular network,^{34,35} the glymphatic network of the brain,^{36,37} osteocyte network formation,^{38} social networks,^{39} and, in particular, in neural networks in the brain where neural plasticity plays an important role for learning,^{15,40} but also for the progression of certain neuro-degenerative diseases.^{41} In the current context, this raises the question if adaptivity in a coupled oscillator network can increase its ability to synchronize.

We devise a special variation of the Kuramoto model with adaptive coupling weights. While previous studies have considered such systems from various perspectives, we here introduce a parameterization that allows to systematically deviate from the traditional Kuramoto system with stationary uniform coupling and systematically study the effects of adaptation and the related bifurcation behavior.

This article is structured as follows. In Sec. II, we explain our adaptive Kuramoto model and how we parameterize it. In Sec. III, we consider $N=2$ oscillators and carry out a detailed bifurcation analysis for several important parameter cases: (Sec. III A) Non-adaptive limit (classical Kuramoto model); (Sec.III B) Symmetric adaptation ( $\alpha =\beta =0$); (Sec. III C) Asymmetric adaptation. In Sec. IV, we numerically investigate larger oscillator systems with $N=50$ and compare the resulting dynamic behavior with the results for $N=2$. Section V concludes with a summary and discussion of our results.

## II. MODEL

### A. General model

#### 1. Coupling function and adaptation rule

*adaptation offset*, $ a 0\u22600$, and the

*adaptivity*$ a 1$ (i.e., the strength of adaptation). The adaptation rule is inspired by models of spike-time-dependent plasticity, a Hebbian synaptic learning rule,

^{42}and can be approximated in phase oscillator networks

^{43}as we do here. The shape of the adaptation rule induces different types of adaptation or learning.

^{9,10,12}In this context, the

*adaptation shift*$\beta $ tunes the type of interaction. Restricting our attention to positive adaptivity, $ a 1>0$, the following adaptation regimes may be distinguished (see also Fig. 1):

Symmetric adaptation with in-phase gain ( $\beta =0$) amplifies (or suppresses) the coupling strength between oscillators that have a phase difference close to 0 (or $\pi $).

Symmetric adaptation with anti-phase gain ( $\beta =\pi $) amplifies (or suppresses) the coupling strength between oscillators that have a phase difference close to $\pi $ (or 0).

Non-symmetric adaptation ( $\beta \u22600,\pi $) amplifies (or suppresses) coupling strengths between oscillators with non-trivial phase difference; e.g., for the odd adaptation function with $\beta =\pi /2$, coupling strengths associated with phase difference close to $\u2212\pi /2$ (or $\pi /2$) are amplified (or suppressed).

Note that the system has the parameter symmetry $(a,\beta )\u21a6(\u2212a,\beta +\pi )$. Thus, if we consider $\beta =0$ for both positive ( $ a > 0$) and negative ( $a<0$) adaptivity, we need not include the case of $\beta =\pi $ in our investigation. In addition to these special cases, we are also interested in investigating mixed-type learning rules such as the case $\beta =\pi /4$, see Sec. III C 3.

Finally, note that symmetric/asymmetric adaptation corresponds to symmetric/asymmetric coupling, i.e., undirected/directed network weights.

#### 2. Reduced/rescaled model and its symmetries

^{44}Furthermore, we may reduce the number of parameters involved by appropriately rescaling parameters and variables,

*(relative) adaptivity*, $a:= a 1/ a 0$, we then obtain the governing equations for the model,

Since $g(\varphi )$ and $ A(\varphi )$ only depend on phase differences, (5) is invariant to shifts in constant phase and frequency, i.e., $ \varphi l(t)\u21a6 \varphi l(t)+\Phi +\Omega t$ with constants $\Phi ,\Omega \u2208 R$.

#### 3. Non-adaptive case (*a* **=** 0)

^{7,32,45}In the continuum limit, $N\u2192\u221e$, when the coupling strength is subcritical ( $\kappa < \kappa c$), the order parameter converges to 0, $ |Z |\u21920$, corresponding to incoherent oscillations; vice versa, for supercritical coupling ( $\kappa > \kappa c$), the order parameter approaches a non-zero value corresponding to partial synchronization, $ |Z |\u2192C>0$. (Note that for finite oscillator systems, the order parameter is subject to pseudo-random fluctuations and it is preferable to consider the time-averaged order parameter, $\u27e8 |Z(t) | \u27e9 t$. In the subcritical regime, fluctuations will prevent the order parameter from precisely attaining 0.) For the continuum limit $N\u2192\u221e$, the threshold for this transition occurs at $ \kappa c=2/(\pi g(0))$,

^{32,33}where $g(\omega )$ is a unimodal frequency distribution with maximum $g(0)$ and width $\sigma $. Note that we normalized the coupling to $1$; thus, the synchronization threshold, $ \sigma c$, is indirectly determined via $1=2/(\pi g(0))$. Vice versa, if $a\u22600$, we allow the network to be adaptive and we move away from the manifold corresponding to the Kuramoto model with stationary coupling. We will also use the order parameter (6) to characterize the adaptive dynamics for $N=50$ in Sec. IV.

### B. Two oscillator model

Due to (8a), we may restrict the parameter range to $\u2212\pi /2<\beta \u2264\pi /2$; due to (8b), we may restrict the parameter range to $\u2212\pi /2<\alpha \u2264\pi /2$; furthermore, (8c) allows to restrict either the range of $\beta $ or $\alpha $ to positive values, and we chose to restrict $0\u2264\beta \u2264\pi /2$. However, observe that all symmetry transformations also affect other parameters and variables; e.g., the symmetry (8b) effectively swaps (near-)phase-locked states where $\varphi $ is close to zero, with (near-)antiphase states where $\varphi $ is close to $\pi $. Thus, if a phase-locked stationary state is stable for $\alpha =0$ with adaptivity $a$, then an antiphasic stationary state is also stable for $\alpha =\pi $ with $\u2212a$. There are several other symmetries, see the Appendix, some of which are reflected in the stability diagrams [Figs. 3–6(a)].

Furthermore, we consider $\epsilon >0$. As our analysis will show, we find interesting nontrivial dynamic behavior (at least) if $\epsilon =0.2$. We use this value throughout the analysis unless specified otherwise.

## III. ANALYSIS

### A. Non-adaptive limit (classical Kuramoto model)

*drifting*oscillators. When the frequency mismatch is smaller, $ |\omega |<1$, we have the two equilibria

### B. Symmetric adaptation (*β* **=** 0, *π*)

*β*

*π*

#### 1. Reduction

Equation (14) admits for two degenerate cases: (i) if $\omega =cos\u2061\alpha =0$, $\varphi $ is constant and, thus, $\kappa (t)\u21921+acos\u2061\varphi $ exponentially fast and (ii) if $cos\u2061\alpha =0$ with $\omega \u22600$, we have $ \varphi ( t ) = \varphi 0 + \omega t$ so that oscillators drift apart and $\kappa (t)$ forever undergoes an oscillation.

#### 2. Stability analysis

We seek to determine bifurcation curves in $(\omega ,a)$-parameter-space. Our system is planar so that saddle-node and Hopf bifurcations are determined via the conditions $det ( J )=0$ and $ tr ( J )=0$ with $det ( J )>0$, respectively. Bogdanov–Takens bifurcations occur at the intersection of Hopf and saddle-node bifurcations. Eliminating $\varphi $ and $\kappa $ corresponding to equilibria in $det ( J )$ or $ tr ( J )$ appears not to be feasible since (18) neither can be solved for $\kappa $ in closed form nor is there a suitable substitution to achieve the elimination. Instead, to determine the desired bifurcation curves, we seek a parameterization of $\omega $ and $a$ in terms of one of the equilibrium variables, $\kappa $.

##### 1. Saddle-node bifurcation curve

##### 2. Hopf bifurcation curve

##### 3. Homoclinic bifurcation of (librational) limit cycles

When a librational limit cycle (born in a Hopf bifurcation) collides with a saddle point, it may eventually get destroyed in a homoclinic bifurcation. Homoclinic bifurcations have been determined numerically as follows. The associated homoclinic bifurcation has been determined by inspection of phase portraits while varying parameters, i.e., for a fixed value of a by determining the value of $\omega $ at which the limit cycle is destroyed in saddle collision. When limit cycles in this planar system are unstable, Eqs. (15) have been numerically integrated backwards in time while continuing $\omega $. Doing so for several values of $a$ allowed to construct the associated homoclinic bifurcation curve shown in Fig. 3.

##### 4. Bogdanov–Takens points

Curves of saddle-node, Hopf, and homoclinic bifurcations intersect in a codimension 2 Bogdanov–Takens bifurcation point (BT). This point is found by seeking solutions that simultaneously satisfy equilibrium conditions together with $det ( J )= tr ( J )=0$.

##### 5. Cusp bifurcation points

Another codimension 2 bifurcation point is the cusp bifurcation (CP) where two saddle-node bifurcations meet. Thus, we find that CP points are located at $( \omega CP, a CP)=(0,\xb11)$.

##### 6. Heteroclinic bifurcation and loss of stability of drift cycle

A similar procedure as in Sec. III B 2 c was used to determine bifurcations of the drift cycle, the rotational limit cycle (limit cycle wrapping around the cylinder). This drift cycle is destroyed in a (global) heteroclinic bifurcation (Het, see Fig. 3). However, for $a<0$, the drift cycle loses appears to lose its stability (D-loss, see Fig. 3) before a saddle collision may occur; instead, trajectories approach a stable spiral.

#### 3. Stability and bifurcation diagrams

We explain the bifurcation scenarios and associated stability boundaries based on our previous analysis. We first observe that the bifurcation and stability diagram (see Figs. 2 and 3, respectively) exhibit certain symmetries regarding the reflection $\omega \u21a6\u2212\omega $. It is easy to check that Eqs. (15), (18), and (19) obey the symmetry (A4), which effectively swaps the labels of the two oscillators and, thus, preserves equilibria, their stabilities, and bifurcations (as seen in Figs. 2 and 3). A further symmetry regarding $a\u21a6\u2212a$ (A5) preserves equilibria, but the stability of the equilibria is not preserved as can be seen by inspecting the Jacobian in (20) [more specifically, saddle-node and cusp bifurcations are preserved, while Hopf bifurcations are not, see Figs. 3, 2(b), and 2(c)]. See also the Appendix for further details and explanations.

The non-adaptive case ( $a=0$) represents the limit of the classical Kuramoto model (see Sec. III A). This case illustrates the perhaps most fundamental bifurcation organizing the structure of the stability diagram on a “global” scale. When intrinsic frequencies are sufficiently similar, $\omega = \omega 1\u2212 \omega 2<1$, oscillators lock their frequency, i.e., synchronization is achieved; vice versa, if intrinsic frequencies are too dissimilar, $\omega >1$, oscillator phases drift apart. Close to $\omega =0$, oscillators are nearly in-phase (compare with Fig. 2 where $a\u22600$); but the phase difference $\varphi $ increases with larger $ |\omega |$ until the locking breaks apart in a saddle-node bifurcation (black curve) and a drifting state appears. This saddle-node bifurcation separates the stability regions for frequency-locked (L, light gray shaded region) and drifting (D, white region) solutions and can be followed throughout the stability diagram in Fig. 3. Note that this drift solution corresponds to a rotation on the cylinder $ T\xd7 R$.

For non-zero adaptivity ( $a\u22600$), the structure of the stability diagram in Fig. 3 is organized around two bifurcation points of codimension two, from where additional bifurcation curves emanate: (i) cusp points (CP, black dot) located at $( \omega CP, a CP)=(0,\xb11)$ and (ii) Bogdanov–Takens points (BT, purple dot). Furthermore, several parameter regions of bistability (dark gray shaded regions, yellow regions) appear near the cusp bifurcation points (dark gray shaded regions) and the saddle-node bifurcations leading to drifting solutions (yellow regions).

##### 1. Positive adaptivity (*a* > 0)

For smaller frequency mismatch (around $ |\omega |<1$) and intermediate adaptivity ( $0<a<1$), frequency-locked solutions (L, light gray shaded region in Fig. 3) with positive coupling strengths exist, stable nodes [black curve in Fig. 2(b)] closer to in-phase and saddles (black dashed curve) closer to anti-phase configurations. For $a\u22651$, two additional saddle-node bifurcation curves [black dots in Fig. 2(a)/black curves in Fig. 3] emerge in a cusp bifurcation (CP, black dot in Fig. 3) located at $( \omega CP, a CP)=(0,1)$. Traversing these saddle-node bifurcation curves gives birth to an additional stable node and a saddle with negative coupling, see Fig. 2(a). Thus, for $a>1$, above the cusp point, between the inner SN bifurcation curves, we have a region of bistability (dark gray shaded region) with two stable locked solutions characterized by near in- and anti-phase configurations, respectively.

For larger frequency mismatch (around $ |\omega |>1$), as we diminish $ |\omega |$ and traverse the saddle-node bifurcation giving rise to the lock state (L), we first observe a region of bistability (D/L, yellow region in Fig. 3) between the lock state (L) and a drift cycle (D). As the frequency mismatch $ |\omega |$ is further diminished, the drift cycle, which already was present in the region D for larger $ |\omega |$ collides with the saddle version of the lock state emanating from the saddle-node bifurcation and is, thus, destroyed in a heteroclinic bifurcation. (labeled Het in Fig. 3; see also Sec. III B 2 c). Note this heteroclinic bifurcation is global and expected to depend on $\u03f5$.

##### 2. Negative adaptivity (*a* < 0)

The symmetry (A5) swaps near in-phase with anti-phase configurations and changes the stability of equilibria. The overall bifurcation structure is, thus, similar to the one observed for $a>0$; however, their structure is more intricate as explained in the following.

First we restrict our attention to smaller frequency mismatch ( $ |\omega |<1$). For $\u22121<a<0$, the saddle-node bifurcations [black dot in Fig. 2(c)/black curve in Fig. 3] now give birth to a saddle and an unstable node. The unstable node turns into an unstable spiral, which later gains stability in a Hopf bifurcation (HB). The resulting stable spiral [red line in Fig. 2(c)] now, however, corresponds to a near in-phase state with positive coupling (but smaller compared to $a>0$). We observe a cusp bifurcation point (CP, black dot in Fig. 3) at $( \omega CP, a CP)=(0,\u22121)$. Below the cusp point ( $a<\u22121$), a Bogdanov–Takens point (BT $ 2$, purple dot in Fig. 3) appears on the associated saddle-node bifurcation [black dot in Fig. 2(d)/black curve in Fig. 3], which gives rise to a subcritical Hopf bifurcation [blue dot in Fig. 2(d)/blue curve in Fig. 3] that stabilizes the second lock state, a stable spiral with small negative coupling and $\varphi $ near $\xb1\pi /4$ [red line in Fig. 2(d)]. The unstable limit cycles (librations) emanating from the Hopf bifurcation get destroyed in a homoclinic bifurcation (purple curve).

For larger frequency mismatch ( $\omega >1$), we observe a region of bistability (D/L, yellow region in Fig. 3) between lock (L) and drift (D) states similar to the one observed for positive adaptivity, $a>0$. However, there are important differences. First, below the Bogdanov–Takens point BT $ 1$, the frequency-locked solution (L) loses stability in a subcritical Hopf bifurcation (blue curve in Fig. 3). More specifically, as shown in Figs. 2(c) and 2(d), the saddle-node bifurcation (black dot) gives birth to a saddle (black dashed curve) and an unstable node, which first becomes an unstable spiral (red dashed curve) and then becomes a stable spiral (red line) in a Hopf bifurcation (blue dot).

Second, for adaptivity values roughly above the BT $ 1$ point, [the precise location of the origin for this boundary remains unclear (see also detection method in Sec. III B 2)] the drift cycle D is destroyed in a heteroclinic bifurcation; this scenario, however, appears to be different for adaptivity $a$ below BT $ 1$: On the boundary of the D/L (labeled D-loss in Fig. 3; see also Sec. III B 2 f) and L regions (light gray shaded region)—before any heteroclinic bifurcation occurs—the drift (limit) cycle appears to lose its stability. Thus, trajectories emanating from the unstable drift end up spiraling into the stable spiral that arises in the Hopf bifurcation related to BT $ 1$.

To summarize the respective differences regarding stability changes occurring for $a<0$ and $a>0$, the presence of Bogdanov–Takens bifurcation points BT $ 1$ and BT $ 2$ diminishes the size of the L locking region and the 2L bistability region (dark gray shaded region in Fig. 3). Furthermore, the size of the bistable D/L region with lock and drift is also effectively diminished.

### C. Asymmetric adaptation

We now allow for arbitrary values of $\beta $; in general, this may lead to asymmetric coupling strengths, $ \kappa 12\u2260 \kappa 21$, corresponding to directed network weights.

#### 1. Stability analysis

Finally, if $\mu =\delta =0$, two eigenvalues are zero and the equilibrium is a Bogdanov–Takens point. Thus, simultaneously solving $\mu =\delta =0$ together with the fixed point conditions (28), given a value of $\epsilon $ we are able to obtain $(\varphi , \kappa 12, \kappa 21,a,\omega )$ at a Bogdanov–Takens point.

#### 2. Dynamics for *β* **=** *π*/2

*β*

*π*

##### 1. Phase-lag *α* = 0 and *α* = *π*

##### 2. Arbitrary phase-lag (*α* ≠ 0, *π*)

It is easily checked that (34a) has the parameter symmetry $(a,\alpha )\u21a6(\u2212a,\u2212\alpha )$. Since we were allowed to limit $\alpha \u2208(\u2212\pi /2,\pi /2]$ in the original governing Eq. (7), we may restrict $\alpha $ further to either $[\u2212\pi /2,0]$ or $[0,\pi /2]$ while observing this symmetry.

To keep the analysis manageable, we consider only $\alpha =\u2212\pi /10$. The resulting stability diagram is shown in Fig. 4(b). Even though this is just a small deviation, the bifurcation landscape differs drastically from the case where $\alpha =0$. As was the case for $\alpha =\beta =0$ (Fig. 3), the stability diagram is symmetric regarding $\omega \u21a6\u2212\omega $ [preserving equilibria and their stabilities and bifurcations, see (A4)]. Moreover, equilibria and the SN and cusp bifurcations are symmetric about $a\u21a6\u2212a$ [preserving equilibria and the SN and cusp bifurcations, see (A5)]. See the Appendix for further details and explanations.

#### 3. Dynamics for *β* **=** *π*/4

*β*

*π*

Next, we consider the case where $\beta =\pi /4$. The stability diagram for $\alpha =0$ (not shown) is qualitatively identical to the one obtained for $\alpha =\beta =0$ (Fig. 3). Therefore, we instead consider small deviations from $\alpha =0$, i.e., $\alpha =\xb1\pi /10$.

##### 1. Phase-lag *α* = *π*/10

First, we consider the stability diagram for $\alpha =\pi /10$, see Fig. 5. There are three distinct Bogdanov–Takens points [BT $ 1$ and BT $ 2$ in the main plot of Fig. 5 and (two) BT $ 3$ in the inset], which organize the bifurcation structure. The Hopf curves (blue curves) emanating from all three Bogdanov–Takens points (BT $ 1$, BT $ 2$, BT $ 3$) are subcritical. These Hopf curves are adjacent to SN curves (black curves) giving birth to saddles and unstable nodes. The unstable nodes turn into unstable spirals, and, when undergoing the Hopf bifurcations, into stable spirals. The homoclinic curves (purple curves) adjacent to the SN curves destroy the unstable (libration) limit cycles created in the subcritical Hopf bifurcations.

##### 2. Phase-lag *α* = −*π*/10

Second, we consider the stability diagram for $\alpha =\u2212\pi /10$ in Fig. 6. The resulting dynamics become far more involved when compared to the previously considered cases. The Hopf curve (blue curve) that emanates from the BT $ 1$ point [purple dot in Fig. 6(a)] contains a Generalized Hopf (GH) point located at ( $ \omega GH$, $ a GH$) [blue dot in Fig. 6(a)], which separates the Hopf curve (blue curve) into a supercritical segment (below GH) and a subcritical segment (above GH).

###### 1. *Regime of subcritical Hopf bifurcations (**a* > *a*_{GH})

*a*>

*a*

_{GH})

The subcritical Hopf curve [blue curve in Fig. 6(a)] produces unstable (libration) limit cycles, which may undergo a saddle-node-of limit cycles (SNLC) bifurcation (black dashed curve) before they die in a homoclinic bifurcation (purple curve). The homoclinic bifurcation curve, which originates in the point BT $ 1$ (purple dot), was found via numerical continuation using MatCont.^{46} Below the GH point (blue dot), there are two branches of SNLC bifurcations which meet in a cusp of SNLCs (CPC, green dot). Thus, between the GH and the CPC points, limit cycles produced in the Hopf bifurcation (blue curve) undergo *two* SNLC bifurcations before dying in the homoclinic bifurcation (purple curve). One of the SNLC branches (black dashed curve) meets the homoclinic bifurcation curve (purple curve) in a point which we call SLH (red dot).

###### 2. *Regime of supercritical Hopf bifurcations (**a* < *a*_{GH})

*a*<

*a*

_{GH})

The supercritical Hopf bifurcations [blue curve in Fig. 6(a)] below the GH (blue dot) produce stable (libration) limit cycles (LB), which can only exist on the same side of the supercritical Hopf curve as the unstable spirals, i.e., in the blue region between the HB curve (blue) and the HC curve (purple) that destroys the limit cycles. However, the stable limit cycle can be destroyed/destabilized before reaching the HC bifurcation (purple curve) in an intriguing bifurcation scenario. We let $a=\u22121.9425$ and follow the green line in the diagram [Fig. 6(a)] while varying $\omega $ and show example trajectories for chosen values of $\omega $ in Fig. 7.

After the supercritical Hopf bifurcation [right blue dot in Fig. 6(b)], i.e., for increasing $\omega $, the stable limit cycle LB [blue in Figs. 6(b), 6(d), and 7(a)] undergoes a Period-Doubling (PD) bifurcation [left green dashed line in Fig. 6(d)]. The resulting stable limit cycle [blue in Figs. 6(d) and 7(b)] has period 2. The original, period-1 cycle has become unstable [blue dashed curve in Fig. 6(d)] in the PD [left green dashed line in Fig. 6(d)]; the associated branch is ultimately destroyed in the homoclinic bifurcation [purple dot in Figs. 6(b) and 6(c)]. Remarkably, note that the libration cycle LB co-exists with a (rotational) drift limit cycle in the region of period-1 and period-2 cycles, see magenta curves in Fig. 6(d) and magenta trajectories in Figs. 7(a) and 7(b). Unlike the previously described D state, these drift cycles consist of both rotations and librations, which is why refer to them as mixed oscillations (MOs). While the precise nature of the bifurcation mechanism giving rise to its birth on the left for smaller $\omega $ remains open, inspection of trajectories for larger $\omega $ on the right strongly suggests that this mixed oscillation is destroyed in a collision with the unstable period-1 version of the libration cycle [blue dashed curve in Fig. 6(d)].

As we increase $\omega $ further, the stable period-2 cycle undergoes a period-doubling cascade (PDs are not explicitly marked to improve readability) resulting in a chaotic attractor with $\varphi $ bounded, i.e., it can (still) be seen as a libration LB [blue in Figs. 6(c), 6(d), and 7(c)]. To characterize this chaotic motion, we numerically estimated the maximum Lyapunov exponent $ \lambda max$, and we found that $ \lambda max>0$ (close to zero) throughout the (non-)chaotic ranges of $\omega $ [Figs. 6(c) and 6(d)]. (Fluctuations in $ \lambda max$ while varying $\omega $ are due to the fact that phase points of the original and the perturbed trajectory can lie anywhere on the chaotic attractor by the simulation time $ \lambda max$ is computed.)

This chaotic attractor can be seen as a chaotic libration, since $\varphi $ does not rotate around the cylinder $ T\xd7 R$ [Fig. 7(c)]. Remarkably, as $\omega $ is further increased, the chaotic attractor at some point also includes oscillations characterized as *rotations*, where $\varphi $ revolves around the cylinder [Fig. 7(d)] (in a fashion that is reminiscent of *phase slips*). More precisely, the trajectory librates a number of times, until it escapes and becomes a rotation with drift character, while, however, remaining chaotic in nature. Presumably, this behavior is akin to a mixed mode oscillation seen characteristic of slow-fast systems, which is a reminiscence from low $\epsilon $ (recall that we use a fairly large $\epsilon =0.2$ here). Eventually, as $\omega $ is increased even further, the cycle ceases to display chaoticity and becomes a periodic oscillation including both librations (of a certain period) and (drift) rotations [Fig. 7(e)], which we therefore refer to as *mixed oscillations* (MOs; not to be confused with mixed-mode oscillations^{47}). Ultimately, as we further increase $\omega $, the librations disappear and the cycle is a purely rotational drift cycle [Fig. 7(h)]. This is remarkable, since—as already noted further above—a stable mixed oscillation cycle was already present for smaller $\omega $ [magenta curves in Figs. 6(c), 6(d), 7(a), and 7(b)], which was destroyed in a bifurcation scenario different from the scenario seen here. Thus, these (rotational/librational) drift cycles just described are indeed distinct, as they seem to emerge from the libration cycles originally born in the supercritical Hopf bifurcation. To summarize, stable oscillations encountered are all either non-drifting libration cycles (LB, periodic or chaotic), purely rotational drift states (D, periodic), or mixed rotational/librational drift oscillations (MO, periodic or chaotic).

Finally, note that states containing both librations and rotations exist, for $a=\u22121.9425$, in the range $\omega \u2208[0.69456,0.78299)$ [Figs. 6(c) and 7(d)–7(g)], where subintervals hosting periodic cycles alternate with subintervals hosting (aperiodic) chaotic motion. This was found via quasi-continuation of the stable oscillatory states in the interval $\omega \u2208[0.693,1.25]$. Interestingly, the respective number of rotations and librations per period of the oscillation is different for each periodic subinterval [see Figs. 7(e)–7(g)]. Finally, based on quasi-continuation of the drift states in the direction of $\omega \u21920$ for several values of $a$, we find that stable purely rotational drift states (D) exist only in the white regions and yellow regions in Fig. 6(a). While the blue region *may* contain all aforementioned types of LB and MO cycles, none of these types emerges everywhere inside of the blue region; for readability, we do not distinguish existence regions for each type individually.

## IV. SIMULATIONS FOR *N* = 50 OSCILLATORS

In Sec. III, we have studied in detail the behavior of the system (5) for $N=2$ oscillators. The question abounds whether some of the behavior observed for only $N=2$ oscillators carries over to a larger version of the system with $N=50$ oscillators. In order to address this question, we restrict our attention to the simplest case with parameters $\alpha =\beta =0,\epsilon =0.2$, (see Figs. 2 and 3 for the results obtained for $N=2$). For two oscillators, we could tune the frequency mismatch, $\omega $; for the many oscillator systems, we instead let the intrinsic frequencies $ \omega l$ follow a normal distribution with zero mean and standard deviation $\sigma $. Furthermore, instead of drawing these frequencies randomly, we constructed them via an equidistant set of points mapped to the interval $[\u22123\sigma ,3\sigma ]$ via the inverse error function. This symmetrization eliminates potential random effects that due to finite size effects might otherwise obscure the dynamic effects we are interested in, such as to render comparison to the case of $N=2$ more immediate.

We study the dynamic asymptotic behavior of this system on a grid in $(\sigma ,a)$ parameter space, the results of which are summarized in Fig. 8(a). We may distinguish four states: the (frequency-)locked state (L) and the antipodal state (AP), two frequency-locked states where all $d \varphi l/dt,d \kappa l m/dt$ tend to zero after some transient; a drift (D) state with all oscillator phases drifting a part; and the partial synchrony state (PS) which can be seen as a mixture of both L and D states, in the sense that the $N$ oscillators are split into a frequency-locked and a drifting subset. We describe these states in the order of appearance as we increase $\sigma $ for $a=1.7$ as shown in Fig. 8.

For the *antipodal/anti-phase state* (AP) [Fig. 8(b)], the phases of oscillators split into two groups inside which phases differ only very little from one another, but in between the groups phases differ by about $\pi $, i.e., the two oscillator groups are in an anti-phase (antipodal) configuration with respect to one another on the phase circle; $d \varphi l/dt,d \kappa l m/dt$ tend to zero asymptotically in time. This special configuration of the phases also impacts the coupling on the network: $ \kappa l m$ that represent couplings between oscillator pairs residing in the same population (distinct populations) assume values close to the maximum possible value $1+a$ (the minimum possible value $1\u2212a$), see black (white) $ \kappa l m$ values in Fig. 8(b). Note that these values for $ \kappa l m$ are heterogeneously distributed near these values as they accommodate for the non-identical phase configuration (the figure compresses the variation due to the chosen gray scale). Note that the antipodal state occurs when the distribution of $ \omega l$ is sufficiently narrow and the adaptivity $a>0$ is sufficiently strong. This state may be interpreted as a corollary to the 2L co-existence observed for $N=2,\alpha =\beta =0,a\u226b0, |\omega |\u226a1$ [see Fig. 2(a)], where the two oscillators either occupy a phase configuration with 0 or $\pi $ difference; here, on the other hand, oscillators form groups of varying size that mutually adhere to either of the two configurations. Note that a similar correspondence applies to the coupling strengths both at the upper and the lower end of the applicable $\kappa $ range.

For the *frequency-locked state* L [Fig. 8(c)], all phases are close to one another. This configuration implies that the $ \kappa l m$ are close to $1+a$ (and bounded by that value from above). Note that individual coupling strengths $ \kappa l m$ for $l$ and $m$ are not identical but heterogeneously distributed: the coupling strengths have adapted to accommodate for the non-identical phase configuration. In Figs. 8(b) and 8(c), the gray scale for the values of $ \kappa l m$ ranges from $1\u2212a$ to $1+a$ so that the distributed values cannot easily be distinguished. Topologically, this essentially corresponds to an all-to-all network. This state is, thus, interpreted as the corollary of the frequency-locked state L for $N=2$, see Fig. 2(b) for small $ |\omega |$.

For the *partially synchronous state* (PS) [in Fig. 8(d)], oscillators split into two groups where one is locked (blue dots for phases), while the other contains drifters (red dots for phases). Correspondingly, coupling strengths $ \kappa l m$ where $l$ and/or $m$ can be associated with drifting oscillators and are, thus, oscillatory [red curves in Fig. 8(d)]. By contrast, the locked oscillators [blue dots in Fig. 8(d)] display vanishingly small dynamic frequencies; thus, the $ \kappa l m$ for which neither $l$ nor $m$ are drifting oscillate close to $1+a$ and with very small amplitude [blue time traces in Fig. 8(d)]. We note that the correspondence to states observed for $N=2$ is less obvious; however, an extension of the D/L bistable region for $N=2$ oscillators to larger systems could give space to accommodate for configurations where a group of oscillators populate a lock state while the remaining oscillators adhere to drifting oscillations.

For the *drifting state* (D), oscillator phases are *incoherent* and drift apart from each other in unlocked motion while the coupling strengths $ \kappa l m$ appear to attain a (quasi-)stationary distribution centered around 1 (not shown). We interpret this state as the corollary to the D state for $N=2$.

Just as for the $N=2$ state, the locked state L occurs for both positive and negative $a$. For $N=50$ and $a<0$ with $\sigma $ small, we observe locked states (corresponding to stable foci for $N=2$) whose $ \kappa l m$ are mostly close to the minimal value $1\u2212a$, in analogy to the case of $N=2$ oscillators [see Fig. 2(c)].

For the $N=2$ system, we found that (co)existence of one (L) or two (2L) stable locked state(s) is guaranteed for $ |a |$ sufficiently large. Interestingly, Fig. 8(a) suggests that there is a stable locked state for $a>0$ sufficiently large; yet, for $a<0$ it appears as if the region of stable locked states is diminished in favor of partially synchronous, and, even drift states (incoherence). Similarly, we observed for $N=2$ and $a<0$ that the locking region is diminished when compared to $a>0$, which is due to a subcritical Hopf bifurcation that destabilizing the frequency-locked branch (albeit not as much as strongly as for $N=50$)—we hypothesize that this effect is more pronounced for larger oscillator numbers (and ensuing bifurcations). It would be interesting to investigate this aspect in future research.

## V. DISCUSSION

The model in Eq. (5) is a Kuramoto–Sakaguchi model with phase-lag $\alpha $ and coupling strengths that adapt according to a learning rule with adaptivity strength $a$ and adaptation shift $\beta $. The parameter $a$ allows to tune away from the limit of the non-adaptive classical Kuramoto model with stationary coupling ( $a=0$) and systematically investigate the effect of adaptivity on the collective dynamics.

Our bifurcation analysis concerns the case of $N=2$ oscillators, for which we—chiefly—may distinguish three paradigms of adaptivity: (i) the non-adaptive limit with stationary coupling ( $a=0$); (ii) symmetric adaptation ( $\beta =0$ or $\beta =\pi $); and (iii) asymmetric adaptation (arbitrary $\beta $). The non-adaptive limit ( $a=0$) trivially reduces to the classical Kuramoto model, where, for frequency locking to occur, the frequency mismatch must be smaller than the coupling strength between oscillators.

Considering the case of symmetric coupling ( $\beta =0$), we first observe that deviating from the non-adaptive limit ( $a=0$) with non-zero adaptivity ( $a\u22600$) leads to an overall larger locking region (L), i.e., a larger frequency mismatch is required to break locking. It is interesting to note that in the context of forced oscillations, the mode-locking region corresponds to an Arnol’d tongue; indeed, the strength of adaptivity could be related to a forcing strength acting on one of the oscillators, at least in the limit of slow adaptation ( $\epsilon \u21920$). Furthermore, smaller frequency mismatch $ |\omega |$ and sufficiently large adaptivity $a$ allow for bistable regions 2L where two frequency locked modes exist, i.e., anti-phase (antipodal) configurations co-exist in addition to in-phase configurations—indeed, considering stationary coupling in, e.g., Eq. (15) effectively results in a bi-harmonic coupling function which may lead to such bi-stability.^{48} For larger $ |\omega |$ and $a\u22600$, a further bistability region D/L appears where locked states co-exist with drift cycles (rotations). While this general picture prevails for both positive and negative adaptivity, the situation is slightly more complicated in the region with $a<0$, where additional bifurcations diminish the width of the various (bi-)stable regions (L, 2L and D/L). Drift cycles correspond to rotational limit cycles around the cylindrical phase space, whereas Hopf bifurcations give rise to librational limit cycles; however, in the case of symmetric coupling ( $\beta =0,\pi $), these cycles remain always unstable.

We found that similar bifurcation scenarios occur for networks with asymmetric coupling for $\beta \u22600$ and $a>0$; but, even more intriguing dynamics are possible for $a<0$ and $\beta =\pi /4,\alpha =\u2212\pi /10$. In short, a generalized Hopf bifurcation (GH) may be present, around which a complex structure of bifurcations is organized [including sub-/supercritical Hopf, saddle-node of limit cycles (SNLC), cusp of cycles (CPC), a homoclinic and SLH bifurcations, similar to a scenario seen for a system of Theta neurons, see Ref. 49]; in particular, for $a<0$ the GH may give rise to a supercritical Hopf bifurcation leading to stable librational limit cycles. These libration cycles may undergo a further period-doubling cascade to chaos. Moreover, remarkably, period-1 and period-2 libration cycles may co-exist with a rotational drift cycle characterized by a non-trivial winding number; this rotational cycle appears to be destroyed in a collision with the unstable period-1 cycle generated by the first period-doubling bifurcation, before the stable libration cycle becomes chaotic. However, the period-doubling cascade of these librations features a surprise. As $\omega $ increases, after the librations have become chaotic, parameter regions exist where the librations are altered into a mixed oscillation, which combines features of both libration and rotation. The winding number of the libration gradually changes, and mixed oscillations alternate between chaotic and periodic, until for large $\omega $ one finds a regular period-1 rotational drift cycle.

One may expect that the analysis of the $N=2$ oscillator system might be able to capture certain aspects of the dynamic behavior seen in larger systems. Our simulations for adaptive networks with symmetric coupling ( $\beta =0$) with $N=50$ revealed that such a correspondence exists, albeit with some limitations. First, considering $a>0$, we find a relatively good correspondence between antipodal/anti-phase (AP), locked (L), and drifting/incoherent states (D) seen for $N=50$ with the bistable locking region (2L), the locking region (L), and drifting region (D) seen for $N=2$; the partially synchronous (PS) state for $N=50$ may be related to the D/L region for $N=2$. Second, increasing strength of adaptivity $ |a |$ leads to a widening of the locking regions for $N=2$; this effect appears to be present for $N=50$ only when $a>0$, but not for $a<0$.

We point out similarities between the system studied here and systems subject to previous research. These systems are contained as special cases in our model (5). The relationships are more obvious when considering the unscaled version of our adaptive model (4): Letting $ a 0=0$ with $\alpha =\beta =0$, we recover the model Eqs. (3) and (4) studied by Seliger *et al.*, Ref. 8. Similarly, if we set $ a 0=0$ with $ \omega 1=\u2026= \omega N=0, a 1=1$, we recover the model by Berner *et al.*^{10} Our adaptation rule generalizes the model in Ref. 8 by inclusion of an adaptation shift $\beta \u22600$ (see Ref. 10); but in contrast to Ref. 10, we introduced a nontrivial adaptation offset $ a 0$ and strength of adaptation $ a 1$. This choice allowed us to systematically deviate from the classical Kuramoto–Sakaguchi model with stationary coupling, and thereby, to address the question of how varying levels of adaptivity impacts the dynamics of the network. Furthermore, unlike Ref. 10, we allow for non-identical frequencies ( $\omega \u22600$ and $\sigma \u22600$, respectively) which breaks symmetries in the system—thus, this enables us to study the width of the locking region. We note that the choice $ a 0=1$ does not allow us to retrieve the model in Ref. 8 since $ a 0\u21920$ constitutes a singular limit.

Nevertheless, some of the basic dynamic behavior, such as the presence of frequency locking or antipodal (antiphasic) cluster states seen in Ref. 10, naturally carry over to our system as long as intrinsic frequencies are kept relatively close to each other. Vice versa, breaking the symmetry via disordered frequencies ( $\omega \u22600$) opens the door for bifurcations resulting in novel states, such as small amplitude librations and the associated period-doubling cascade reported here. A recent study^{50} carries out an analysis of our system with $ a 0=0$ and $ a 0=1$ in the slowly adapting limit ( $\epsilon \u21920$); however, this analysis separates the dynamics of the coupling strengths into a planar problem, and thus, the resulting dynamics cannot display the transition to chaos or the mixed oscillations containing both librations and rotations that we found here. Finally, we note that Kasatkin and Nekorkin^{51} studied the same oscillator system with $N=2$ as in this study, except that frequencies are identical ( $\omega =0$); in this case, these authors also observed chaotic dynamics.

Our simulations for $N=50$ revealed that oscillators may split into two groups, which experience distinct strong or weak coupling, depending on the specific phase configuration that the oscillators assume. This phenomenon is particularly prominent for the states we called AP and PS. We note that several studies addressed the effect of non-uniform, but *time-independent* coupling on oscillator networks. In order to understand the occurrence of chimera states, past research performed bifurcation analysis for $N=4$ oscillators^{52} and for larger $N$, or in the continuum limit ( $N\u2192\u221e$).^{53–57} However, drawing analogies to our model is not straightforward. First, the configurations seen for the AP or PS state in Fig. 8 display groups with different oscillator numbers, but aforementioned studies consider oscillator groups of equal size. Second, the adaptive capability of the systems renders the dynamics intrinsically more complicated. For instance, the adaptive rule in Eq. (5b) may give rise to bi-harmonic coupling function (depending on the parameter choice) and, thus, to bi-stability of configurations, see Fig. 2, a feature that is due to higher harmonic interaction.^{48} [This becomes evident, e.g., by assuming a stationary solution for (15).] Moreover, the adaptivity may stabilize certain non-uniform coupling configuration, while models without adaption simply enforce the non-uniform coupling. This is a fundamental difference—as we have seen, adaptive dynamics generates nontrivial dynamic behavior via bifurcations that are absent in the non-adaptive limit (Sec. III A). There exist other models with bimodal frequency distributions where the oscillator population may effectively split into two groups adhering to distinct dynamic behavior,^{45,58–61} and even models that combine bimodally distributed frequencies with non-uniform coupling strength;^{62} however, how to correlate the observed states with the behavior observed for these systems seems less obvious.

We plan to further study adaptive Kuramoto–Sakaguchi oscillator networks with large oscillator numbers. To this end, a variety of questions are interesting to pursue. For instance, it will be interesting to see if similar complex dynamics, such as the period-doubling cascade seen for $N=2$ with asymmetric coupling ( $\beta \u22600,\pi $; see Sec. III C), persists for larger systems. Moreover, developing a (finite) mean-field theory as well as a continuum limit description for Eq. (5), e.g., using graphons/graphops and other methods,^{33,63–67} would be very useful to study the dynamics of Eq. (5), but remains a formidable challenge. Finally, it would be interesting to study and relate adaptive oscillator networks to more realistic settings such as those that occur in a biological, neuro-scientific, or experimental context.^{21,25,43} These represent formidable open problems for future research.

## ACKNOWLEDGMENTS

We thank M. P. Sørensen, C. Bick, S. Yanchuk, and R. Cestnik for useful discussions and comments throughout this study. We acknowledge the DTU International Graduate School for support via the EU-COFUND project “Synchronization in Co-Evolutionary Network Dynamics (SEND).”

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Benjamin Jüttner:** Conceptualization (equal); Formal analysis (equal); Methodology (equal); Writing – original draft (equal); Writing – review & editing (equal). **Erik A. Martens:** Conceptualization (equal); Formal analysis (equal); Methodology (equal); Supervision (lead); Writing – original draft (lead); Writing – review & editing (equal).

## DATA AVAILABILITY

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

### APPENDIX: SYMMETRIES

In summary, symmetry (A4) preserves equilibria, their stabilities, as well as all bifurcations since we are merely renaming the oscillators; the symmetry (A5) preserves SN and cusp bifurcations, but not Hopf bifurcations. Note that these two symmetries are evident in all stability diagrams shown in this article, i.e., Figs. 3–6(a).

## REFERENCES

*Sync: The Emerging Science of Spontaneous Order*

*Synchronization: A Universal Concept in Nonlinear Sciences*

*Progress in Industrial Mathematics at ECMI 2018*, edited by I. Faragó, F. Izsák, and P. L. Simon (Springer International Publishing, 2019), pp. 147–155.

*Adaptive networks*(Springer, 2009), pp. 231–251.

*Chemical Oscillations, Waves, and Turbulence*