Over the past decades, chimera states have attracted considerable attention given their unexpected symmetry-breaking spatiotemporal nature and simultaneously exhibiting synchronous and incoherent behaviors under specific conditions. Despite relevant precursory results of such unforeseen states for diverse physical and topological configurations, there remain structures and mechanisms yet to be unveiled. In this work, using mean-field techniques, we analyze a multilayer network composed of two populations of heterogeneous Kuramoto phase oscillators with coevolutive coupling strengths. Moreover, we employ the geometric singular perturbation theory through the inclusion of a time-scale separation between the dynamics of the network elements and the adaptive coupling strength connecting them, gaining a better insight into the behavior of the system from a fast–slow dynamics perspective. Consequently, we derive the necessary and sufficient condition to produce stable chimera states when considering a coevolutionary intercoupling strength. Additionally, under the aforementioned constraint and with a suitable adaptive law election, it is possible to generate intriguing patterns, such as persistent breathing chimera states. Thereafter, we analyze the geometric properties of the mean-field system with a coevolutionary intracoupling strength and demonstrate the production of stable chimera states. Next, we give arguments for the presence of such patterns in the associated network under specific conditions. Finally, relaxation oscillations and canard cycles, seemingly related to breathing chimeras, are numerically produced under identified conditions due to the geometry of our system.

Chimera states have been widely studied due to their unanticipated nature as well as their resemblance to particular synchronization patterns observed in a range of mechanical, biological, and neurological activities. In particular, such interesting dynamics have been modeled and analyzed in networks composed of Kuramoto phase oscillators with diverse structures and configurations. However, most of the existing related studies consider static connections between the nodes of the system, limiting the understanding of the mechanisms involved in the evolution of the overall network. Here, we introduce a novel technique for the generation of different synchronization patterns, including chimera states, through the incorporation of a slow adaptation in the coupling strengths, revealing new underlying mechanisms and, thus, leading to a better description of real-world events involving chimera states.

## I. INTRODUCTION

Synchronization is a crucial element in a large variety of natural and practical processes occurring in daily life.^{1,2} Particularly, synchronization arises in complex systems as neural networks,^{3,4} circadian rhythms,^{5,6} power transmission lines,^{7,8} flashing patterns in fireflies,^{9,10} superconducting Josephson junctions,^{11} chemical oscillations,^{12,13} neutrino oscillations,^{14} and extensive biological processes.^{15,16} A remarkable contribution to this topic was introduced by the Japanese physicist Yoshiki Kuramoto who, inspired by the pioneering work of Arthur Winfree,^{17} presented an archetypal model for the study of spontaneous synchronization, which consists of a complete network of coupled oscillators with randomly distributed natural frequencies and a coupling kernel dependent on the sine of their phase difference.^{18}

Although identical oscillators were expected to only display complete synchronous or incoherent dynamics, Kuramoto and Battogtokh found that, under specific initial conditions, a ring of identical coupled oscillators can separate in two distinguishable spatial regions, one consists of almost completely synchronized elements while the other is composed of partially incoherent oscillators.^{19,20} Later, Abrams and Strogatz named this new unexpected pattern as *chimera state* due to their resemblance to the mythological beast composed of dissonant elements.^{21} Since then, chimera states have been reported in a vast diversity of contexts, gathered in extensive review papers,^{22–24} and various physical phenomena have been related to such behavior. For instance, in nature, many animal species engage in unihemispheric slow-wave sleep in which a highly active region of the brain coexists with the other hemisphere performing in a more erratic manner.^{25} Similarly, the simultaneous presence of highly synchronized and completely incoherent patterns have been reported in ventricular fibrillation,^{26,27} social and cultural trends,^{28} as well as in neuronal models such as non-locally coupled Hodgkin–Huxley oscillators,^{29} leaky integrate-and-fire neurons,^{30} as well as in coupled FitzHugh–Nagumo oscillators,^{31} among several other neural network models.^{22}

In this paper, we present a novel mechanism for the generation of chimera states in a complete network composed of different populations of heterogeneous Kuramoto phase oscillators in which the strength of the coupling connecting the nodes adaptively evolves depending on the macroscopic state of the system. Such coevolutionary networks have been observed in several systems like vascular networks,^{32} opinion dynamics,^{33,34} power grids,^{35} and neural processes involved both in learning^{36,37} and the progression of specific neuro-degenerative diseases.^{38}

Particularly, enabled by our model, we employ the Ott–Antonsen ansatz to obtain a reduced order mean-field representation of the networks’ dynamics. Moreover, by introducing a slow adaptive law, we induce a fast–slow structure, allowing for the first time the identification of stable chimera states with the stable equilibria of a fast–slow system. Additionally, by analyzing the geometric properties of the mean-field, characteristics leading to different synchronization patterns are determined in the reduced manifold, which then, by considering the normally hyperbolic case, are preserved in the network when assuming sufficiently heterogeneous populations in our ensemble with a slowly coevolutive coupling strength, given a suitable large time-scale separation between the dynamics off and on the network. Hence, we provide a novel mechanism that allows one to design the desired synchronization pattern in the mean-field, which is then observed in the network under adequate conditions. In particular, we prove that stable chimera states and persistent breathing chimeras are achieved under a sufficiently slow coevolving law with specific parameter arrangements. Finally, other parameter regimes related to the loss of normal hyperbolicity conditions and leading to complex behavior on the reduced system are identified and compared numerically at the network level, generating exciting future work ideas.

The rest of the paper is structured as follows. In Sec. II, we present our research model with all the necessary assumptions, from which we derive an equivalent mean-field based on the Ott–Antonsen technique in order to analyze the geometric properties of the system from a reduced order representation. Later, in Sec. III, we present the main results of our study for two distinct problems, the coevolutionary inter- and intracoupling, respectively. For each scenario, we determine the geometric properties of the related mean-field and give conditions for the effective generation of chimera states when considering adaptive regimes only dependent on macroscopic quantities to preserve the Ott–Antonsen method. Then, utilizing the attractiveness of the Ott–Antonsen manifold, we provide with arguments for the chimera state to be reflected in the finite-size network. Thereafter, in Sec. IV, we present numerical results as evidence of our research by comparing the dynamics of the mean-field and the network for different macroscopic adaptive laws, showing numerically the successful generation of the desired synchronization patterns by the incorporation of a coevolutive dynamic designed in the mean-field for the coupling strengths of the network. Next, in Sec. V, we give our conclusions and discuss future research opportunities as a consequence of the outcomes presented. In particular, we emphasize on the analysis of the non-hyperbolic case for which relaxation oscillations and canard solutions, related to breathing chimera states due to the patterns numerically observed, have been produced in simulations for both the mean-field and the network. Finally, in the Appendix, we introduce the basic notions supporting our results, including an overview of the Ott–Antonsen ansatz as well as a brief summary on fast–slow dynamics for the interested reader.

## II. COEVOLUTIONARY MULTILAYER KURAMOTO NETWORK

^{39}Such a quantity is of crucial relevance for our study as it dictates the dynamics of the coupling strength since it represents a macroscopic term of the network, and it does not affect the mean-field approach that is presented later.

Furthermore, the parameter $ \beta \sigma \sigma \u2032$ represent the phase-lag existing through elements between $(\sigma = \sigma \u2032)$ and across $(\sigma \u2260 \sigma \u2032)$ populations, correspondingly. For the purpose of our work, we consider networks without a phase-lag, i.e., $ \beta \sigma \sigma \u2032=0$, for $\sigma , \sigma \u2032=1,2$, as the necessary level of heterogeneity in the ensemble is already induced by the natural frequencies (4). Later, in Remark 7, we present a detailed explanation on why the assumption $ \beta \sigma \sigma \u2032=0$ is considered and its effect on our results. Nevertheless, we display some examples for which $ \beta \sigma \sigma \u2032$ is greater than zero but sufficiently small as it can be treated as a perturbation parameter in our analysis, both for the adaptive inter- and the intracoupling case. Finally, the coevolutive coupling strength $ \zeta \sigma \sigma \u2032$ represent the slowly adaptive variable considered which only depend on macroscopic quantities of (1), simplifying the reduction method by considering the remaining couplings as constant parameters in the fast time-scale and avoiding modifications on the mean-field technique employed.^{40}

Using our notation, there are two intracoupling strengths, namely, $ \zeta \sigma \sigma $, i.e., $ \zeta 11$ or $ \zeta 22$, and one intercoupling strength $ \zeta \sigma \sigma \u2032$, i.e., $ \zeta 12= \zeta 21$. Moreover, in order to avoid further confusion, we denote the intracouplings as $ k \sigma = \zeta \sigma \sigma $, and the intercoupling as $\mu = \zeta \sigma \sigma \u2032$, for $\sigma , \sigma \u2032=1,2$.

^{41}for two populations without phase-lag, i.e., $ \beta \sigma \sigma \u2032=0$, a mean-field representation of (1) is given by

In Sec. III, we analyze the behavior of (6) under the presence of different types of adaptive coupling strengths, only dependent on macroscopic quantities of (1). In particular, we benefit from the introduction of a time-scale separation between the dynamics on the nodes and that related to the coupling strength, allowing us to use the fast–slow theory without alterations on the reduction method employed.

## III. ANALYSIS

To begin with and in order to facilitate computations, we reduce the dimension of the fast-time mean-field (6) by demonstrating that a near synchronization state is attracting when certain parameter conditions are fulfilled. The aforementioned idea is synthesized in the next proposition.

### (Invariant and attracting synchronization set)

**(Invariant and attracting synchronization set)**

The set $ { ( \rho 1 , \rho 2 , \psi , \mu ) \u2208 K : \rho 2 = 1}$ in (6), $ ( K = [ 0 , 1 ] 2 \xd7 [ 0 , 2 \pi ) \xd7 R )$, is invariant and attracting if and only if $ \Delta 2=0$, for $ k 2$ large enough. Moreover, there exists an invariant and attracting set $ { ( \rho 1 , \rho 2 , \psi , \mu ) \u2208 K : \rho 2 ( \Delta 2 ) = 1 \u2212 a \Delta 2 + O ( \Delta 2 2 )}$ with $a>0$, and $ \Delta 2$ positive and sufficiently small when $ k 2>0$ is adequately large.

Therefore, from Proposition III.1 and by choosing $ k 2$ sufficiently large and $ \Delta 2$ sufficiently small, in what follows, we simply consider $ \rho 2=1$. Moreover, we analyze two different problems, namely, the coevolutionary intercoupling, i.e., $\mu $ adaptive, and the coevolutive intracoupling, i.e., $ k 1$ adaptive. This distinction is made in order to compare the geometries of both adaptation schemes and their isolated effects in the network. Notice that stable chimera states can be generated in both coevolving scenarios with determined geometric differences as we subsequently demonstrate. To begin with, in Subsection III A, we study the adaptive intercoupling problem in (6) while considering the intracouplings as constant parameters.

### A. Coevolutionary intercoupling

Different results are obtained when considering if the distributions from which the natural frequencies are drawn from are centered at the same value or not, i.e., $\Omega =0$ or $\Omega \u22600$, respectively. We start by analyzing the former case in the following proposition.

#### (Intercoupling normal hyperbolicity and stability ( | Ω | ≪ 1)

**(Intercoupling normal hyperbolicity and stability ( $ |\Omega |\u226a1$)**

The critical manifold (12) of (11), with $ \Omega = 0$, is normally hyperbolic in the interval $ \rho 1\u2208[0,1)$ if and only if $ k 1<2 \Delta 1$. Moreover, if the critical manifold is normally hyperbolic, then it is attracting for every $ \rho 1\u2208[0,1)$. Finally, for $ |\Omega |\u226a1$, the normal hyperbolicity and stability of the $\Omega =0$ case are conserved.

*Proof.* Notice that, when $\Omega =0$, the Jacobian (13) reduces to

^{42}□

Hence, Proposition III.2 states that as long as the natural frequencies of the first population are sufficiently sparse, every point in the critical manifold (12) is normally hyperbolic and attracting in the interval $ \rho 1\u2208(0,1)$, when the centers of the natural frequency distributions are sufficiently close. Observe that this conclusion matches with the global attractiveness condition of the Ott–Antonsen manifold discussed in the Introduction of this work.

Next, we study the case when only the intracoupling of one population evolves on time while preserving constant both the other intracoupling and the intercoupling strength. For convenience, we analyze system (6), under Proposition III.1, for $ k 1$ adaptive.

### B. Coevolutionary intracoupling

Notice that (16) can have complex elements, producing disconnected branches of $ C 0 \xb1$. In fact, by definition (see Definition A.1 in the Appendix), the critical manifold is a geometric object of real components. Therefore, we restrict our attention to the case for which the discriminant in (16) is non-negative.

#### (Connected critical manifold)

When $\Omega =0$, the critical manifold (16) is completely connected for every $\mu \u2208 R$.

#### (Intracoupling normal hyperbolicity and stability ( | Ω | ≪ 1 ))

**(Intracoupling normal hyperbolicity and stability $( |\Omega |\u226a1)$)**

The critical manifold (16) of (15) has one normally hyperbolic and attracting branch $ ( C 0 \u2212 )$, and one unstable branch $ ( C 0 + )$ along the interval $ \rho 1\u2208(0,1)$ for every $\mu \u22600$, with $\Omega =0$. Moreover, for $ |\Omega |\u226a1$, the normal hyperbolicity and stability conditions of the $\Omega =0$ case are preserved.

^{42}

The implementation of Proposition III.1 in system (6) sets the problem into a convenient scenario for chimera-like behaviors, in which one population is synchronized while the other can produce different patterns, as similar procedures have been applied in the non-adaptive case.^{43} Nevertheless, if such a result is neglected, the normal hyperbolicity and stability properties of the critical manifolds (12) and (16), for both the coevolutive inter and intracoupling cases, (11) and (15) respectively, will be different, although the dimension of such manifolds is preserved as long as only one slowly adapting component is considered given the fast–slow structure of the problem.

Consequently, since we have already provided with the conditions for normal hyperbolicity and attractiveness of the critical manifolds (12) and (16) of the adaptive inter- and intracoupling systems, now we present the result which enables us to produce different synchronization patterns, including stable chimera states, by considering a slow macroscopic coevolution in (11) and (15).

#### (Stable chimera states on the mean-field)

**(Stable chimera states on the mean-field)**

Given an adequate macroscopic adaptive law in the coevolving intercoupling (11) or intracoupling (15) case, any long-time synchronization pattern, including stable chimera states, is produced in the perturbed mean-field when the normal hyperbolicity and attracting conditions for the critical manifold $ C 0 \xb1$ are satisfied, for $\epsilon $ positive and sufficiently small.

Once we have determined the normally hyperbolic and attracting conditions of the critical manifold $ C 0 \xb1$, and recalling that this geometric object is composed of exponentially stable equilibria of the fast problem, we can construct a suitable slow macroscopic adaptive law with an exponentially stable equilibrium point at $ ( \rho 1 , \psi , \zeta \sigma \sigma \u2032 )= ( \rho 1 \u2217 , \psi \u2217 , \zeta \sigma \sigma \u2032 \u2217 )$, such that if this point is part of $ C 0 \xb1$, then it is exponentially stable in the perturbed system,^{44} for $ \zeta \sigma \sigma \u2032$, with $\sigma , \sigma \u2032=1,2$, denoting the type of coevolutive coupling considered, either inter- $ ( \zeta \sigma \sigma \u2032 = \zeta 12 = \zeta 21 = \mu )$, or intracoupling $ ( \zeta 11 = k 1 , \zeta 22 = k 2 )$. Moreover, by Fenichel’s theory, there exists a slow manifold for which the dynamics converge to the slow flow as $\epsilon \u21920$, and by regular perturbation results, the dynamics of the unperturbed system persist under weak perturbations, so there exists a point in a neighborhood of $ ( \rho 1 \u2217 , \psi \u2217 , \zeta \sigma \sigma \u2032 \u2217 )$, which is an exponentially stable equilibrium of the perturbed system (6) for sufficiently small $\epsilon $. Finally, recalling that we consider a near synchronization second population, i.e., $ \rho 2\u21921$, it is possible to set the equilibrium point at any position along the critical manifold with an adequate adaptive law. Hence, by persistence, any synchronization behavior, including stable chimera states, can be produced in the mean-field of the weakly perturbed system.

Finally, the next result provides the necessary condition that allows us to identify the dynamics occurring in the mean-field system with long time behaviors at the network level.

#### (Stable chimera states on a network)

**(Stable chimera states on a network)**

Under normally hyperbolic conditions of the critical manifold $ C 0 \xb1$, both for the adaptive intercoupling (11) or intracoupling (15) systems, dynamical behaviors observed in the mean-field reduction, such as stable chimera states and sustained breathing chimera states, are preserved in the network when considering heterogeneous populations with a slowly coevolving and macroscopic coupling strength, for $\epsilon $ sufficiently small and for $N$ large enough.

It is known that, as long as the probability distributions for the natural frequencies of the oscillators have a finite width, i.e., $ \Delta \sigma >0$ $(\sigma =1,2)$, the attractors for the order parameter dynamics obtained in the reduced Ott–Antonsen manifold (5) are the only attractors of the network (1).^{45,46} Moreover, Theorem III.1 demonstrates the preservation of equilibrium points on the perturbed system for $\epsilon $ small enough. Therefore, dynamical behaviors observed in the mean-field reduction, including stable chimera states and breathing chimera states, are present in the network when considering non-homogeneous populations.

## IV. NUMERICAL RESULTS

In this section, we present a series of simulations validating the effectiveness of our results as well as some intriguing patterns observed under certain non-hyperbolic conditions. It is important to emphasize though that the Ott–Antonsen ansatz consider the thermodynamic limit, i.e., $ N 1 , 2\u2192\u221e$, so the network and the mean-field outputs tend to increasingly approach as the number of oscillators is augmented. For this reason, the images corresponding to the network results have two components, the actual output (orange) and a smooth version (red) produced by filtering the original trajectory with a Savitzky–Golay filter,^{47} as this least squares smoothing method reduces noise while preserving the shape and height of waveform peaks,^{48} producing a closer approximation to the mean-field response without dramatically increasing the number of oscillators in our network. Additionally, for the numeric results presented in this work, we adopt the observations made by Böhle, Kuehn, and Thalhammer, which, as they claim, makes the simulations faster compared to classical approaches.^{49}

### A. Coevolutionary intercoupling

First, the behavior of both the network (1) and its mean-field counterpart (11) is compared when different intercoupling coevolving rules are employed. In Fig. 2, we present projections to the $ ( R 1 , \mu )$ and $ ( \rho 1 , \mu )$ planes of (1) and (11), with $ \mu \u02d9= \epsilon \mu (\u2212\mu + \gamma \mu \u2212 \eta \mu \rho 1)$, where, for distinction purposes, $ \rho 1$ represents the mean-field variable in (11) and $ R 1$ describes the measured modulus of the complex Kuramoto’s order parameter in the simulated network (1), respectively. Notice that $ k 1<2 \Delta 1$ and thus the critical manifold (gray) is normally hyperbolic and attracting for every $ \rho 1\u2208(0,1)$. Moreover, observe that the behavior of $ R 1$ and $ \rho 1$ (red) converge to the intersection of the critical manifold and the equilibria set of the slow adaptation (magenta), which corresponds to a highly incoherent region of the first layer and with the assumption that the second population is near synchrony, i.e., $ \rho 2\u22481$ and $ R 2\u22481$, we have effectively produced an exponentially stable chimera state, demonstrating the generation of such pattern in the network (1) by geometric results derived from its mean-field reduction (11), as stated by Theorem III.2.

Similarly, in Fig. 3, we present a persistent breathing chimera state produced through the same mechanism as before in both (1) and (11), with $ \mu \u02d9= \epsilon \mu cos\u2061 ( 0.02 t )$. Observe that, since the fast structure of the problem has not been modified, the critical manifold, as well as its normal hyperbolicity and stability conditions, remain the same and notice that $ R 1$ and $ \rho 1$ oscillate between high and low synchrony regimes, which correspond to the so-called breathing chimera state, characterized by sustained periodical variations on the synchronization level.^{50} Finally, notice that the synchronization level of the second population presents small-amplitude oscillations with the same frequency of the first population due to the intracoupling strength $ k 2$ employed, when this value is increased a stronger interaction in the second population is achieved and the oscillatory pattern is only appreciated in the first population. Therefore, the ratio between both intracoupling strengths $ k 1 , 2$ is selected in such a way that the oscillations in the second populations’ level of synchrony exhibits only small amplitude oscillations.

### B. Coevolutionary intracoupling

Analogously to the intercoupling problem, we examine the numeric results obtained by considering system (1) with an adaptive intracoupling strength and its associated mean-field reduction (15). To begin with, a linear feedback rule $ k \u02d9 1= \epsilon 1(\u2212 k 1+ \gamma 1\u2212 \eta 1 \rho 1)$ is considered and in Fig. 4, we present the obtained results. As explained earlier, the critical manifold (16) of (15) has one normally hyperbolic and attracting branch, namely $ C 0 \u2212$, while the other branch is at least of saddle type. In Fig. 4, the stable, saddle, and unstable segments of the critical manifold (16) are represented by the continuous gray, dashed green, and dotted blue lines, respectively. Observe that we have again initialized the first population of oscillators near complete synchronization and after a transient phase, it is clear that a chimera state is produced as the order parameters $ \rho 1$ and $ R 1$ converge to the intersection between the critical manifold and the slow flow set of equilibria. Therefore, as long as an appropriate macroscopic and slow adaptive law is designed, it is possible to produce diverse patterns due to the geometric properties of the mean-field (15), and by means of Theorem III.2, such long term dynamical features are reproduced in the network (1). Finally, it is important to highlight the chimera-like pattern produced in Fig. 4 as it is stable for a coupling beyond the critical value, i.e., $ k 1>2 \Delta 1$, which without any adaptation would usually lead the network to the synchronous attractor, similar to the results obtained for weakly heterogeneous oscillators,^{51} but now the coevolutive mechanism allows for larger heterogeneity in the network.

^{19,52–54}the assumption $\beta \sigma \sigma \u2032=0$, for $\sigma ,\sigma \u2032=1,2$, is done in order to simplify the analysis, as the derived results hold for $\beta \sigma \sigma \u2032$ greater than zero and sufficiently small. In this regard, following the results for the Ott–Antonsen

*ansatz*applied to an homogeneous two-population network with the same intracouplings but different intercoupling between layers,

^{43}the resulting mean-field equations represent our reduced system in the limit $\beta \sigma \sigma \u2032\u21920$. The aforementioned fact is due to the geometric properties of the system and perturbation results, which is why it is possible to produce different synchronization patterns, including chimera states, even when in the presence of the phase-lag. For instance, in Fig. 5, we present an example of a chimera state in a coevolutionary intercoupling scenario with a phase-lag $\beta 11=\beta 12=\beta 21=\beta 22=\pi /6$. Similarly, in Fig. 6, we show a stable chimera state in the coevolutionary intracoupling case for which $\beta 11=\beta 12=\beta 21=\beta 22=\pi /10$. Naturally, the trajectory follows a manifold $C\epsilon $ sufficiently close to the critical manifold $C0\xb1$ in gray until it reaches the intersection between $C\epsilon $ and the slow adaptation nullcline. Notice that the apparent shift between $C\epsilon $ and $C0\xb1$ is due to the phase-lag considered but, as long as such parameter is sufficiently small, the normally hyperbolic structure is preserved. Therefore, the behavior obtained is in the end a stable chimera state, as can be appreciated from the local order parameters of the aforementioned examples.

## V. DISCUSSION AND CONCLUSIONS

The term *chimera state* gathers a broad collection of highly interesting and unexpected synchronization patterns. In this work, we have shown a new mechanism for the generation of stable chimera states in a complete multilayer network of heterogeneous Kuramoto oscillators through the inclusion of slow coevolutive coupling strengths. Particularly, we use the Ott–Antonsen *Ansatz* to derive a mean-field representation of the general network, allowing the use of geometric singular perturbation theory (GSPT) results to obtain a deeper insight into the dynamics occurring in the system. At first, we have analyzed a two-population network with only one coevolving intercoupling strength, and determined the necessary and sufficient conditions for the critical manifold to be normally hyperbolic and attractive. For such a scenario, we have effectively generated stable chimera states as well as persistent breathing chimera states only by adequately modifying the macroscopic adaptive law dictating the slow coevolving dynamics. Similarly, we studied the case of the same network with one coevolutive intracoupling strength and obtained relevant geometric results that enable the production of stable chimera states, as well as different synchronization patterns, by appropriately modifying the adaptation law employed. Additionally, we have shown the effectiveness of our mechanism for the generation of diverse synchronization patterns in a network of coupled phase oscillators with coevolutive coupling strengths when there exists a phase-lag between the elements of the ensemble. The aforementioned fact holds as long as such a parameter is sufficiently small, as it can be regarded as a perturbation in geometric terms for the given problem. Therefore, the presented numeric results validate our mechanism for the generation of diverse patterns in a complete network of Kuramoto phase oscillators with coevolving macroscopic coupling strengths.

Considering the results we have obtained, it is natural and intriguing to wonder about which kind of behaviors are observed in the non-normally hyperbolic scenario for both types of adaptive couplings. Although rigorous mathematical analysis is currently under development, numerical results have shown remarkable behaviors. First, it is possible to produce oscillating synchronization patterns without a periodically forced adaptive law. For instance, in Fig. 8, a stable breathing chimera is presented for a two-population system with two-coevolving coupling strengths, $\mu (t)$ and $ k 1(t)$. As it can be seen, the periodic oscillation is sustained for the given set of parameters and adaptive laws employed. On the other hand, we suspect that canard cycles and breathing chimeras share a connection since they present a certain qualitative resemblance. Following this idea, in Figs. 9 and 10 two canard cycles for the two-population mean-field with coevolving intercoupling scenarios are shown. The technicalities and rigorous analysis are part of future work, but here we want to present some numerical evidence. In Fig. 9, a canard raising from a singular Hopf bifurcation is depicted. For the selected parameters, the critical manifold folds at a Hopf bifurcation from where the trajectories follow the saddle branch of the critical manifold, giving rise to the depicted canard cycle. Similarly, in Fig. 10, a canard cycle with an oscillatory pattern due to a transcritical bifurcation is shown. Hence, the synchronization level in the first population oscillates periodically, which is why we suspect there exists a connection between canard cycles and breathing chimeras. Nevertheless, these dynamical behaviors are extremely sensitive to initial conditions and parameter variations, rendering highly complex the task of observing them in the network. For this reason, here we show only the results obtained for the mean-field and prefer to present the network results once the non-normally hyperbolic cases are better comprehended. Moreover, we emphasize that the mean-field is only an approximation of the full system and some transient dynamics are not completely captured by this reduction. Thus, a closer correspondence between the dynamics of the network (1) and the mean-field, both (11) or (15), is achieved by increasing the population size with an adequate level of heterogeneity in the network, either by means of the natural frequency distributions or a phase-lag between the oscillators.

## ACKNOWLEDGMENTS

This work was funded by the Centre for Data Science and Systems Complexity (DSSC) of the University of Groningen, as part of the PhD scholarship. The authors thank the editor(s) and anonymous reviewers, whose comments helped us to improve our manuscript.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Luis Guillermo Venegas-Pineda:** Conceptualization (equal); Data curation (equal); Formal analysis (lead); Investigation (lead); Methodology (lead); Writing – original draft (lead); Writing – review & editing (equal). **Hildeberto Jardón-Kojakhmetov:** Conceptualization (equal); Data curation (equal); Supervision (equal); Writing – review & editing (equal). **Ming Cao:** Conceptualization (equal); Data curation (equal); Supervision (equal); Writing – review & editing (equal).

## DATA AVAILABILITY

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

### APPENDIX: PRELIMINARIES

In this section, we summarize the principal methods that are used through the development of our main results. First, we present the Ott–Antonsen ansatz, which produces a mean-field equivalent system of the network under particular considerations. Additionally, we briefly introduce the fast–slow theory for singularly perturbed dynamics and the principal results here employed.

#### 1. Ott–Antonsen reduction

^{39}Subsequently, it is possible to express the angular velocity (A4) in terms of the order parameter (A5), resulting in

^{55}and that in general may depend on the complex order parameter (A5). Moreover, as the exact time dependence of $ H \sigma (t)$ does not affect the derivation of the Ott–Antonsen results, it suffices to consider $ H \sigma (t)$ as some generic time dependent functions, regardless on how this dependence is determined.

^{45}In particular, for our problem consisting on a multilayer network of Kuramoto-type phase oscillators with coevolutive inter- and intracouplings, the functions $ H \sigma (t)$ are defined as $ H \sigma (t)= \u2211 \sigma \u2032 = 1 M k \sigma \sigma \u2032 z \sigma \u2032(t)$, where the effective forcing depends on the population index $\sigma , \sigma \u2032=1,2,\u2026,M$.

^{41}we focus our attention on a particular class of density functions $ f \sigma ( \theta \sigma , \omega \sigma , t )$ in the form of Poisson kernels. In fact, Ott and Antonsen discovered that such kernels satisfy the governing equations if an associated low-dimensional system of ordinary differential equations is satisfied.

^{43}Thus, by considering the Fourier series of the form

*ansatz*(A7) to the continuum equation (A3), we obtain the following set of differential complex equations:

*ansatz*(A7) is a valid solution to the continuity equation (A3). Hence,

^{41,56}It is important to mention that the Ott–Antonsen

*ansatz*corresponds to a particular case of the technique developed by Watanabe–Strogatz

^{57,58}which considers all the harmonics of the frequency distribution function (A7). Thus, although the resulting response in the Ott–Antonsen manifold does not entirely capture the transient behavior of the synchronization level in the network, the stationary states of both systems are effectively the same.

^{59}Furthermore, the Ott–Antonsen manifold is known to be globally attractive when the distribution of natural frequencies is non-homogeneous, i.e., $ \Delta \sigma >0$.

^{45,46}Therefore, every long time behavior of the order parameter can be obtained by analyzing the corresponding mean-field even when weak heterogeneities are considered.

^{60}

^{61}Substituting this expression in (A3) produces $M$-coupled complex ODEs describing the evolution of the complex order parameters as

^{43}produced by orthogonality $2M$ coupled real ODEs in the form

#### 2. Fast–slow dynamics

##### (Critical manifold)

**(Critical manifold)**

Moreover, if $ C 0$ is a submanifold of $ R m\xd7 R n$, it is referred to as the *critical manifold*. Notice that points of (A13) are in direct correspondence to the equilibria set of the fast flow, generated by (A12). Additionally, a relevant property that $ C 0$ can have is normal hyperbolicity, given as follows.

##### (Normal hyperbolicity)

**(Normal hyperbolicity)**

A subset $ S\u2282 C 0$ is normally hyperbolic if the matrix $ ( D x f )(p,0)$ has no eigenvalues with zero real part for all $p\u2208 S$. Moreover, a normally hyperbolic subset is attracting (repelling) if all the eigenvalues of $ ( D x f )(p,0)$ have negative (positive) real part for every $p\u2208 S$. Finally, if $ S$ is normally hyperbolic but neither attracting nor repelling it is of the saddle type.

In contrast, $ S$ is non-hyperbolic if the corresponding matrix $ ( D x f )(p,0)$ has at least one eigenvalue with zero real part for any $p\u2208 S$. Depending on the normal hyperbolicity condition of (A13) different analysis techniques can be employed. For instance, non-hyperbolic points, related to dynamic features such as relaxation oscillations and canard trajectories,^{62} can be studied through the blowup method.^{63} On the other hand, for normally hyperbolic cases, we use Fenichel’s theorem^{64,65} stated as follows.

##### (Fenichel’s Theorem)

**(Fenichel’s Theorem)**

Given a compact normally hyperbolic submanifold (possibly with boundary) $ S= S 0$ of the critical manifold $ C 0$ of (A11), and $f,g\u2208 C k < \u221e$, for $0<\epsilon \u226a1$, then

A locally invariant manifold $ S \epsilon $ diffeomorphic to $ S 0$ exists. Local invariance means that trajectories can enter or escape $ S \epsilon $ only through its boundaries.

$ S \epsilon $ has Hausdorff distance $ d H ( V , W )=O ( \epsilon )$ from $ S 0$, as $\epsilon \u21920$.

The flow on $ S \epsilon $ converges to the flow generated by the reduced problem, known as the slow flow, as $\epsilon \u21920$.

$ S \epsilon $ is $ C k$-smooth.

$ S \epsilon $ is normally hyperbolic and has the same stability properties with respect to the fast variables as $ S 0$.

$ S \epsilon $ is usually not unique. In regions with a fixed distance from $\u2202 S \epsilon $, all manifolds which satisfies conditions 1–5 above are exponentially close to each other with a Hausdorff distance of order $O ( exp ( \u2212 C / \epsilon ) )$ for some $C>0$, $C\u2208O(1)$, as $\epsilon \u21920$.

Regularly, a particular manifold which satisfies conditions 1–5 of Theorem A.1 is referred to as *the* slow manifold. Rigorously, this is not correct since $ S \epsilon $ is not unique, but since all possible choices are exponentially close, it is arbitrary to select one of these manifolds for most analytical and numerical considerations. Moreover, normally hyperbolic submanifolds of the slow problem, such as equilibria and periodic orbits, persist in every slow manifold.^{66}

In the past, sometimes Fenichel’s theorem was referred to as *geometric singular perturbation theory* (GSPT). Nevertheless, nowadays, GSPT includes several geometric techniques such as the fast–slow normal form theory and the blowup method, to name a few.

## REFERENCES

*Synchronization: A Universal Concept in Nonlinear Sciences*

*Chemical Oscillations, Waves, and Turbulence*

*2020 59th IEEE Conference on Decision and Control (CDC)*

*Pure and Applied Mathematics: A Wiley-Interscience Series of Texts, Monographs and Tracts*(John Wiley & Sons, 2007).

*Singular Perturbation Methods in Control: Analysis and Design*

*Dynamical Systems: Lectures Given at the 2nd Session of the Centro Internazionale Matematico Estivo (C.I.M.E.) held in Montecatini Terme, Italy, June 13–22, 1994*, edited by R. Johnson (Springer, Berlin, 1995), pp. 44–118.