Stable Chimera States: A Geometric Singular Perturbation Approach

Over the past decades chimera states have attracted considerable attention given their unexpected symmetry-breaking spatio-temporal nature, simultaneously exhibiting synchronous and incoherent behaviours 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 by two populations of heterogeneous Kuramoto phase oscillators with coevolutive coupling strengths. Moreover, we employ Geometric Singular Perturbation Theory (GSPT) with 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 behaviour 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 co-evolutionary 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 co-evolutionary intracoupling strength and demonstrate the production of stable chimera states which depend on the associated network. Finally, relaxation oscillations and canard cycles, also related to breathing chimeras, are numerically produced under identified conditions due to the geometry of our system.

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 by two populations of heterogeneous Kuramoto phase oscillators with coevolutive coupling strengths.Moreover, we employ Geometric Singular Perturbation Theory (GSPT) with 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 behaviour 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 co-evolutionary 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 co-evolutionary intracoupling strength and demonstrate the production of stable chimera states which depend on the associated network.Finally, relaxation oscillations and canard cycles, also related to breathing chimeras, are numerically produced under identified conditions due to the geometry of our system.

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 Josepshon 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 of almost completely synchronized elements whilst the other 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][23][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 coexist with the other hemispheric 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], and neurology 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 by 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 to design the desired synchronization pattern in the mean-field, which is then observed in the network under the 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 ideas of future work.
The rest of the paper is structured as follows.In section 2 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 section 3, 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 section 4 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 dynamics designed in the mean-field for the coupling strengths of the network.Finally, in section 5 we give our conclusions and discuss future research opportunities as the consequence of the outcomes presented.In particular, we emphasize 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.Additionally, in appendix A 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.

Coevolutionary multilayer Kuramoto network
The system we analyze consists of a complete network composed by two populations of different heterogeneous Kuramoto oscillators intra and interconnected with common coevolutive coupling strengths, as graphically depicted in figure 1 and mathematically expressed in the following set of coupled ordinary differential equations where θ σ i and ω σ i represent the phase and natural frequency of the i-th oscillator, chosen from a unimodal Cauchy-Lorentz distribution as (19), while N σ denotes the number of elements in the σ-layer (σ = 1, 2).Moreover, the parameters β σσ and β σσ ′ represent the phase-lag existing through elements between and across populations, respectively.For the purpose of our work, we consider networks without a phase-lag, i.e. β σσ ′ = 0, for σ, σ ′ = 1, 2, as the necessary level of heterogeneity in the ensemble is already induced by the natural frequencies.Later, in remark 5, we present a detailed explanation on why the assumption β σσ ′ = 0 is considered and its effect on our results.Nevertheless, we display some examples for which β σσ ′ is greater than zero but sufficiently small, as it can be treated as a perturbation parameter in our analysis.Finally, in this work we concentrate only on one coevolutive coupling strength at the time, thus we denote an arbitraty coupling strength with the variable ζ σσ ′ , such that the intracoupling strengths are expressed as while the intercoupling strength is given by (2).Additionally, the coevolutive coupling strength ζ σσ ′ is a slowly adaptive variable which only depend on macroscopic quantities of (1), simplifying the reduction method by considering the couplings as constant parameters in the fast time-scale and avoiding modifications on the mean-field technique employed [39].Thus, in accordance to (27), when Figure 1: Graphic representation of the multilayer coevolutive network (1).The black and blue components represent elements corresponding to the first and second populations, respectively.Moreover, θ σ i denotes the phase of each element in the network, with i = 1, . . ., N σ and σ = 1, 2. The adaptive intercoupling strength (red) is expressed as µ(t) and the intracoupling variable for each population is captured by k σ .Although (1) is a complete graph, we present only certain nodes' connections for simplicity.
Furthermore, considering our two-population problem and by introducing the additional variables ψ = φ 2 − φ 1 , and Ω = ω 2 − ω 1 , the mean-field (3) is expressed as Remark 1. Observe that, since the selected probability distribution functions are of unimodal Cauchy-Lorentz type, the real part of the local Kuramoto order parameter R σ of (1) directly corresponds to the variables ρ σ in (3).
In the following section we analyze the behaviour of (4) under the presence of different types of adaptive coupling strengths, only dependent on macroscopic quantites 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 the use of fast-slow theory without alterations on the reduction method employed.

Analysis
To begin with and in order to facilitate computations, we reduce the dimension of the fast-time mean field (4) by demonstrating that a near synchronization state is attracting when certain parameter conditions are fulfilled.The aforementioned idea is synthesized in the next proposition.
attracting if and only if ∆ 2 = 0, for k 2 large enough.Moreover, there exists an invariant and with a > 0, and ∆ 2 positive and sufficiently small when k 2 > 0 is adequately large.
Therefore, from proposition 3.1 and by choosing k 2 sufficiently large and ∆ 2 sufficiently small, in what follows we simply consider ρ 2 = 1.Moreover, we analyze two different problems, namely, the co-evolutionary intercoupling, i.e., µ 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 the next subsection we study the adaptive intercoupling problem in (4) while considering the intracouplings as constant parameters.

Co-evolutionary Intercoupling
As a first approach, the coevolving intercoupling strength µ effects are analyzed.Therefore, by fixing ρ 2 = 1 and k 1 , k 2 ∈ R, with k 2 sufficiently large, as parameters in both populations of (1), the mean-field ( 4) is where f (ρ 1 , ψ, µ, t) is some adaptive law which only depends on macroscopic quantities of (1), with ε µ as the time-scale separation parameter.Consequently, for the critical manifold of ( 9) we obtain Hence, the critical manifold consists of two branches defined by with L = [0, 1) × [0, 2π) × R.Moreover, the Jacobian of (9) evaluated at (10) is 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., Ω = 0 or Ω = 0, respectively.We start by analyzing the former case in the following proposition.
Hence, proposition 3.2 states that as long as the natural frequencies of the first population are sufficiently sparse, every point in the critical manifold ( 10) is normally hyperbolic and attracting in the interval ρ 1 ∈ (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 (4), under proposition 3.1, for k 1 adaptive.
Remark 2. It is important to highlight that the identification of two branches in (10), namely C + 0 and C − 0 , is done for convenience of notation.Nevertheless, by definition A.1, every point in these branches conforms the critical manifold, regardless its location in a particular segment, as shown in figures 2 and 4 for the coevolutionary inter and intracoupling, respectively.

Coevolutionary Intracoupling
As a comparison to the adaptive intercoupling case of the previous section, we now analyze the effect of one coevolutive intracoupling.By letting k 1 adaptive and with the result in proposition 3.1, now (4) reduces to with the same variable interpretations as before and k 1 as the coevolutive coupling strength.
Thus, the critical manifold of ( 13) is Remark 3. Notice that ( 14) can have complex elements, producing disconnected branches of C ± 0 .
In fact, by definition A.1 the critical manifold is a geometric object of real components.Therefore, we restrict our attention to the case for which the discriminant in ( 14) is non-negative.Proof.In order to produce only real components we require that which can be rewritten as By calculating the derivative of the right-hand side of ( 15) with respect to ρ 1 , its maximum occurs at ρ 1 = 1/ √ 3, for which the function reaches a value of 1/2 √ 3. Hence, the critical manifold ( 14) is completely connected if and only if |µ| ≥ |Ω|/ √ 3.
Next, we obtain the Jacobian of ( 13) with respect to the fast variables, which by evaluation at (14) yields with h(ρ 1 ) := , non-negative under proposition 3.3.Thus, the following result summarize the conditions for the normal hyperbolicity and stability of (14).Proof.Observe that for Ω = 0 the Jacobian ( 16) reduces to with eigenvalues in the form Hence, the branch of ( 14) corresponding to λ − 1,2 , namely C − 0 , is normally hyperbolic and attracting along the interval ρ 1 ∈ (0, 1) for every µ = 0. On the other hand, the branch C + 0 is either of saddle-type or repelling.Finally, for |Ω| ≪ 1 the normal hyperbolicity and stability conditions of the Ω = 0 case persist since the eigenvalues λ ± 1,2 of ( 16) depend differentiably on Ω [40].
Remark 4. The implementation of proposition 3.1 in system (4) sets the problem into a convenient scenario for chimera-like behaviors, in which one population is synchronized whilst the other can produce different patterns, as similar procedures have been applied in the non-adaptive case [41].
Nevertheless, if such a result is neglected, the normal hyperbolicity and stability properties of the critical manifolds ( 10) and ( 14), for both the coevolutive inter and intracoupling cases, (9) and ( 13) 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 (10) and ( 14) 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 ( 9) and ( 13).
Theorem 3.1 (Stable chimera states on the mean field).Given an adequate macroscopic adaptive law in the coevolving intercoupling (9) or intracoupling (13) 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 are satisfied, for ε positive and sufficiently small.
Proof.Once we have determined the normally hyperbolic and attracting conditions of the critical manifold C ± 0 , and recalling that this geometric object is composed by exponentially stable equilibria of the fast problem, we can construct a suitable slow macroscopic adaptive law with an exponentially stable equilibrium point at (ρ 1 , ψ, ζ σσ ′ ) = (ρ * 1 , ψ * , ζ σσ ′ * ), such that if this point is part of C ± 0 , then it is exponentially stable in the perturbed system [42], for ζ σσ ′ , with σ, σ ′ = 1, 2, denoting the type of coevolutive coupling considered, either inter ( . Moreover, by Fenichel's theory there exists a slow manifold for which the dynamics converge to the slow flow as ε → 0, and by regular perturbation results the dynamics of the unperturbed system persist under weak perturbations, so there exists a point in a neighborhood of (ρ * 1 , ψ * , ζ σσ ′ * ) which is an exponentially stable equilibrium of the perturbed system (4) for sufficiently small ε.Finally, recalling that we consider a near synchronization second population, i.e., ρ 2 → 1, 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 which allows us to identify the dynamics occurring in the mean-field system with long time behaviours at the network level.
Theorem 3.2 (Stable chimera states on a network).Under normally hyperbolic conditions of the critical manifold C ± 0 , both for the adaptive intercoupling (9) or intracoupling (13) systems, dynamical behaviours 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 ε sufficiently small and for N large enough.
Proof.It is known that, as long as the probability distributions for the natural frequencies of the oscillators have a finite width, i.e. ∆ σ > 0 (σ = 1, 2, . . ., M ), the attractors for the order parameter dynamics obtained in the reduced Ott-Antonsen manifold (3) are the only attractors of the network (1) [43,44].Moreover, theorem 3.1 demonstrates the preservation of equilibrium points on the perturbed system for ε small enough.Therefore, dynamical behaviours observed in the mean-field reduction, including stable chimera states and breathing chimera states, are present in the network when considering non-homogeneous populations.
In the following section, we show the successful generation of stable chimera states by following our results, with numerical simulations for both the network (1) and the mean-field (4) when considering the inter and intracoupling adaptive scenarios discussed in this section.

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 → ∞, 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 [45], as this least-squares smoothing method reduces noise while preserving the shape and height of waveform peaks [46], 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 [47].
Moreover, observe that the behavior of the order parameters R 1 and ρ 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., ρ 2 ≈ 1 and R 2 ≈ 1, 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 (9), as stated by theorem 3.2.
Similarly, in figure 3 we present a persistent breathing chimera state produced through the same mechanism as before in both ( 1) and ( 9), with μ = ε µ cos (0.02t).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 the order parameters R 1 and ρ 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 [48].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.

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 (13).To begin with, a linear feedback rule k1 = ε 1 (−k 1 + γ 1 − η 1 ρ 1 ) is considered and in figure 4 we present the obtained results.As explained earlier, the critical manifold ( 14) of (13) has one normally hyperbolic and attracting branch, namely C − 0 , while the other branch is at least of saddle type.In figure 4, the stable, saddle and unstable segments of the critical manifold (14) 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 ρ 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 (13), and by means of theorem 3.2, such long term dynamical features are reproduced in the network (1).Finally, it is important to highlight the chimera-like pattern produced in figure 4, as it is  stable for a coupling beyond the critical value, i.e. k 1 > 2∆ 1 , which without any adaptation would usually lead the network to the synchronous attractor, similar to the results obtained for weakly heterogeneous oscillators [49], but now the coevolutive mechanism allows for larger heterogeneity in the network.
Remark 5. Despite the well-known importance of the phase-lag parameter in determining the stability of the chimera state [19,[50][51][52], the assumption β σσ ′ = 0 is done in order to simplify the analysis, as the derived results hold for β σσ ′ 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 [41], the resulting mean-field equations represent our reduced system in the limit 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 figure 5, we present an example of a chimera state in a coevolutionary intercoupling scenario with a phase-lag Similarly, in figure 6, we show a stable chimera state in the coevolutionary intracoupling case for which β 11 = β 12 = β 21 = β 22 = π/10.Naturally, the trajectory follows a manifold C ε sufficiently close to the critical manifold C ± 0 in grey, until it reaches the crossing between C ε and the slow adaptation nullcline.Notice that the apparent shift between C ε and C ± 0 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.
Remark 6.It is important to emphasize that, although there is not an evident and spontaneous symmetry breaking, different synchronization patterns, including the chimera-like behaviors, can be generated for both the inter and intracoupling scenarios in the network (1) and the mean-field (4) by the inclusion of a proper macroscopic adaptive law with an adequate set of parameters.For instance, it is possible to produce a highly synchronized network even when the given parameters, without adaptation, would produce an incoherent system, as depicted in figure 7 for the coevolutionary intercoupling case.Therefore, the inclusion of the adaptation in Figure 4: Locally exponentially stable intracoupling chimera.Upper row: network (left) and mean-field (right) projections to the (R 1 , k 1 ) and (ρ 1 , k 1 ) planes for ( 1) and ( 13), respectively.Critical manifold regions: attracting (grey), saddle (dashed green), repelling (blue dotted), slow adaptation nullcline (magenta) and response (red).Lower row: order parameter evolution in the network (left) and mean-field (right).For the network, second (blue) and first (orange) population synchronization level and the filtered version of the latter (red).Adaptive law k1 = ε 1 (−k 1 + γ 1 − η 1 ρ 1 ) and initial conditions (ρ 1 , ψ, k 1 ) = (0.9, −0.5, 4.0).Parameters:  Finally, we would like to highlight that setting the initial condition of the level of synchrony R 1 (t) in the network is a difficult task as it is determined by the average of the pseudo-randomly generated initial phases of the oscillators, which is why the initial conditions of R 1 (t) and ρ 1 (t) are slightly different in figure 7.

Discussion and conclusions
The term chimera state gather a broad collection of highly interesting 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 strenghts.Particularly, we use the Ott-Antonsen ansatz to derive a mean-field representation of the general network, allowing the use of GSPT results to obtain a deeper insight of 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 attracting.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 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 behaviours 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.Firstly, it is possible to produce oscillating synchronization patterns without a periodically forced adaptive law.For instance, in figure 8 a stable breathing chimera is presented for a two-population system with two-coevolving coupling strengths, µ(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 certain qualitative resemblance.
Following this idea, in figures 9 and 10 two canard cycles for the two-population mean field with coevolving intercoupling scenario are shown.The technicalities and rigorous analysis are part of future work, but here we want to present some numerical evidence.In figure 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 figure 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 behaviours 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 network (1) and the mean-field (9) dynamics 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., µ are both adaptive, thanks to the geometric understanding of the inter and intracoupling critical manifolds, we ensure that the corresponding critical manifold in this case is also normally hyperbolic.This geometric insight allows us to simply propose a periodic adaptation law that leads to this breathing chimera.

A 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.

A.1 Ott-Antonsen reduction
In this work, we analyze a globally connected network formed by different types of heterogeneous Kuramoto oscillators with adaptive coupling strengths.We assume that the coevolutionary laws only depend on macroscopic quantities of the network, which allows us to treat the couplings as fixed parameters for the discussions in this section.Further details are given in section 2. Therefore, we begin by studying an all-to-all network composed by different populations of heterogeneous Kuramoto phase oscillators, which can be expressed in the form where θ σ i and N σ denote the state of the i-th oscillator and the number of elements in the σ-population, respectively.Additionally, k σσ ′ and β σσ ′ are the coupling strength and phase-lag between two layers, with σ = 1, 2, . . ., M .Moreover, ω σ i is the natural frequency of each oscillator in its respective community, selected from M unimodal Cauchy-Lorentz distributions as with ωσ and ∆ σ the central frequency and width parameter of each distribution.Notice that, the heterogeneity in the oscillators' natural frequency distributions itself can lead to the presence of chimera states.Therefore, in what follows we restrict our analysis to the case β σσ ′ = 0.In fact, we will later see that, in our context, β σσ ′ plays the role of a regular perturbation parameter for β σσ ′ small.The case β σσ ′ ≈ π 2 is relevant, but is left for future research because it leads to extra complications in the mean field analysis.
Considering the thermodynamic limit (N σ → ∞ for each σ, but σ finite), the state of (18) can be described by the continuous probability distribution functions f σ (θ σ , ω σ , t), which measure the fraction of oscillators with phases in the interval [θ σ , θ σ + dθ σ ] and natural frequencies between [ω σ , ω σ + dω σ ] at time t.Moreover, since the number of oscillators is preserved for every population at all time, each density function satisfies the continuity equation where v σ (θ σ , ω σ , t) is the continuous version of the angular velocity of the oscillators in each population, expressed as Moreover, the well-known Kuramoto order parameter is defined as In geometric terms, the complex order parameter describes the centroid of all the phasors e ıθ σ and it grows larger as a higher synchronization level is reached in the network [53].Subsequently, it is possible to express the angular velocity (21) in terms of the order parameter (22), resulting in where the over bar denotes the complex conjugate, and the functions H σσ ′ (t) are any smooth, complex-valued 2π−periodic functions of the phases θ σ i , i = 1, . . ., N σ , σ = 1, 2, . . ., M [54], and that in general may depend on the complex order parameter (22).Moreover, as the exact time dependence of H σσ ′ (t) does not affect the derivation of the Ott-Antonsen results, it suffices to consider H σσ ′ (t) as some generic time dependent functions, regardless on how this dependence is determined [43].In particular, for our problem consisting on a multilayer network of Kuramoto-type phase oscillators (1) with coevolutive inter and intracouplings as (2), the functions H σσ ′ (t) are defined as , where the effective forcing depends on the population index σ, σ ′ = 1, 2.
Following Ott and Antonsen [55], we focus our attention on a particular class of density functions f σ (θ σ , ω σ , 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 [41].Thus, by considering Fourier series of the form for the density function of each population, and by imposing the ansatz (24) to the continuum equation ( 20) we obtain the following set of complex equations where c.c. stands for the complex conjugate of the expression in the right-hand side of (25).
By selecting the probability distribution functions as (19), it is possible to calculate zσ by the contour integration in the negative half of the complex plane yielding zσ = α σ (ω σ − ı∆ σ , t) [61].
Substituting this expression in (20) produces M -coupled complex ODEs describing the evolution of the complex order parameter as żσ + (∆ σ − ıω σ ) z σ + 1 2 Finally, with z σ = ρ σ e −ıφσ , for which the negative sign is included in the definition of φ σ in order for the Poisson kernel to converge to δ(θ − φ σ ) and not to δ(θ + φ σ ) as ρ σ → 1 [41], produces by orthogonality 2M coupled real ODEs in the form In particular, our mechanism for the generation of stable chimera state relies on geometric properties of (27) with adaptive coupling strengths under a macroscopic slow adaptation, which induces a time-scale separation between the dynamics on and off the network.Thus, the next subsection briefly summarizes the principal concepts of fast-slow systems employed later on.

A.2 Fast-slow dynamics
A fast-slow system is a singularly perturbed ordinary differential equation in the form ε ẋ = f (x, y, ε), ẏ = g(x, y, ε), where the over-dot represents the derivative with respect to the slow-time τ , x ∈ R m and y ∈ R n denote the fast and slow variables respectively, and 0 < ε ≪ 1 describes the time-scale separation.
Moreover, we consider f : R m × R n × R → R m and g : R m × R n × R → R n of class C k , for k sufficiently large.Equivalently, by defining the fast time t := τ /ε, we can rewrite (28) as x ′ = f (x, y, ε), y ′ = εg(x, y, ε). ( where the prime denotes the derivative with respect to the fast time t.Hence, ( 28) and ( 29) are known as the slow and fast formulations, and in the singular limit (ε = 0) two different problems arise, namely the reduced and the layer problem, respectively from ( 28) and (29).Despite of being non-equivalent both problems share a close connection as expressed in the following definition.
Definition A.1 (Critical manifold).The critical set is defined as Moreover, if C 0 is a submanifold of R m × R n , it is referred to as the critical manifold.Notice that points of (30) are in a direct correspondence to the equilibria set of the fast flow, generated by (29).Additionally, a relevant property that C 0 can have is normal hyperbolicity, given as follows.
Definition A.2 (Normal hyperbolicity).A subset S ⊂ C 0 is normally hyperbolic if the matrix (D x f ) (p, 0) has no eigenvalues with zero real part for all p ∈ 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 ∈ 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 ∈ S. Depending on the normal hyperbolicity condition of (30) different analysis techniques can be employed.For instance, non-hyperbolic points, related to dynamic features such as relaxation oscillations and canards [62], can be studied through the blow-up method [63].On the other hand, for normally hyperbolic cases we have Fenichel's