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.

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 learning36,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.

The system we analyze consists of a complete network composed of two populations of different heterogeneous Kuramoto oscillators intra- and interconnected with common coevolutive coupling strengths, as graphically depicted in Fig. 1 and mathematically expressed in the following set of coupled ordinary differential equations:
(1)
(2)
where θ i σ = θ i σ ( t ) ∈ [ 0 , 2 π ) and ω i σ ∈ R represent the phase and natural frequency of the ith oscillator, while N σ and ρ σ denote the number of elements and the modulus of the complex Kuramoto order parameter (3) in the σ-layer (⁠ σ = 1 , 2⁠), respectively. In (1), the term ζ σ σ ′ = ζ σ σ ′ ( t ) ∈ R represents the intracoupling when σ = σ ′ and the intercoupling when σ ≠ σ ′⁠. Regarding (2), we shall consider that only a single of such couplings, but either of them, is coevolutive. Recall that the local Kuramoto order parameter is defined as
(3)
and it describes the centroid of all the phasors e ı ϕ σ⁠, growing larger as a higher level of synchronization is reached in the 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.
FIG. 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 σ ( t ). Although (1) is a complete graph, we present only certain nodes’ connections for simplicity.
View largeDownload slide

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 σ ( t )⁠. Although (1) is a complete graph, we present only certain nodes’ connections for simplicity.

FIG. 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 σ ( t ). Although (1) is a complete graph, we present only certain nodes’ connections for simplicity.
View largeDownload slide

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 σ ( t )⁠. Although (1) is a complete graph, we present only certain nodes’ connections for simplicity.

Close modal
Remark 1
The natural frequency ω i σ of each oscillator is selected from a Cauchy–Lorentz distribution as
(4)
with ω ^ σ and Δ σ being the central frequency and width parameter of each distribution for the corresponding population σ = 1 , … , M⁠. Moreover, the oscillators’ natural frequency distributions (4) may provide sufficient heterogeneity for the network to exhibit chimera states.

Furthermore, the parameter β σ σ ′ represent the phase-lag existing through elements between ( σ = σ ′ ) and across ( σ ≠ σ ′ ) populations, correspondingly. 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 (4). Later, in Remark 7, 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, both for the adaptive inter- and the intracoupling case. Finally, the coevolutive coupling strength ζ σ σ ′ 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 

Remark 2

Using our notation, there are two intracoupling strengths, namely, ζ σ σ⁠, i.e., ζ 11 or ζ 22⁠, and one intercoupling strength ζ σ σ ′⁠, i.e., ζ 12 = ζ 21⁠. Moreover, in order to avoid further confusion, we denote the intracouplings as k σ = ζ σ σ⁠, and the intercoupling as μ = ζ σ σ ′⁠, for σ , σ ′ = 1 , 2⁠.

Thus, in accordance to the results of Ott and Antonsen,41 for two populations without phase-lag, i.e., β σ σ ′ = 0⁠, a mean-field representation of (1) is given by
(5)
Furthermore, considering our two-population problem and by introducing the additional variables ψ = ϕ 2 − ϕ 1⁠, and Ω = ω 2 − ω 1⁠, the mean-field (5) is expressed as
(6)
Remark 3

Observe that, since the selected probability distribution functions are of unimodal Cauchy–Lorentz type as (4), the modulus of the local complex Kuramoto order parameters of (1) directly corresponds to the variables ρ σ in (6).

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.

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.

Proposition III.1
(Invariant and attracting synchronization set)

The set { ( ρ 1 , ρ 2 , ψ , μ ) ∈ K : ρ 2 = 1 } in (6), ( K = [ 0 , 1 ] 2 × [ 0 , 2 π ) × R )⁠, is invariant and attracting if and only if Δ 2 = 0⁠, for k 2 large enough. Moreover, there exists an invariant and attracting set { ( ρ 1 , ρ 2 , ψ , μ ) ∈ K : ρ 2 ( Δ 2 ) = 1 − a Δ 2 + O ( Δ 2 2 ) } with a > 0⁠, and Δ 2 positive and sufficiently small when k 2 > 0 is adequately large.

Proof.
Consider the function
(7)
associated with the flow equation of ρ 2 in (6). Thus, the set { ( ρ 1 , ρ 2 , ψ , μ ) ∈ K : ρ 2 = 1 } is invariant if and only if Δ 2 = 0⁠. Moreover, by taking the derivative of (7) with respect to ρ 2 and evaluating at ( ρ 2 , Δ 2 ) = ( 1 , 0 ) yields
(8)
which is strictly negative in K for k 2 adequately large. Therefore, from the implicit function theorem there exists a set { ( ρ 1 , ρ 2 , ψ , μ ) ∈ K : ρ 2 = q ( Δ 2 ) } such that h 2 ( q ( Δ 2 ) , Δ 2 ) = 0⁠, with q ( 0 ) = 1⁠, q ∈ C 1⁠. In particular, we obtain an expansion of ρ 2 as a regular perturbation on Δ 2 in the form ρ 2 ( Δ 2 ) = 1 − a Δ 2 + O ( Δ 2 2 ) with the invariant flow of (7) given now as
(9)
which is satisfied for a = ( k 2 + μ ρ 1 cos ⁡ ψ ) − 1⁠, positive whenever | k 2 | ≥ | μ | along the interval ρ 1 ∈ ( 0 , 1 )⁠. Moreover, the derivative of (7) evaluated at ρ 2 ( Δ 2 ) = 1 − a Δ 2 + O ( Δ 2 2 )⁠, given by
(10)
approaches (8) in the limit when Δ 2 tends to zero, rendering the set { ( ρ 1 , ρ 2 , ψ , μ ) ∈ K : ρ 2 ( Δ 2 ) = 1 − a Δ 2 + O ( Δ 2 2 ) } attracting with Δ 2 sufficiently small for k 2 sufficiently large.

Therefore, from Proposition III.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 coevolutionary 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 Subsection III A, we study the adaptive intercoupling problem in (6) while considering the intracouplings as constant parameters.

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 (6) is
(11)
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 (11) we obtain
μ cos ⁡ ψ = ρ 1 ( 2 Δ 1 1 − ρ 1 2 − k 1 ) , μ sin ⁡ ψ = − 2 Ω ( ρ 1 3 ρ 1 2 + 1 ) .
Hence, the critical manifold consists of two branches defined by
(12)
with L = [ 0 , 1 ) × [ 0 , 2 π ) × R⁠. Moreover, the Jacobian of (11) evaluated at (12) is
(13)

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.

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

Proposition III.2
(Intercoupling normal hyperbolicity and stability (⁠ | Ω | ≪ 1⁠)

The critical manifold (12) of (11), with Ω = 0⁠, is normally hyperbolic in the interval ρ 1 ∈ [ 0 , 1 ) if and only if k 1 < 2 Δ 1⁠. Moreover, if the critical manifold is normally hyperbolic, then it is attracting for every ρ 1 ∈ [ 0 , 1 )⁠. Finally, for | Ω | ≪ 1⁠, the normal hyperbolicity and stability of the Ω = 0 case are conserved.

Proof. Notice that, when Ω = 0⁠, the Jacobian (13) reduces to

(14)
Thus, the eigenvalues of (14) are given as
λ 1 = − Δ 1 ( 1 + ρ 1 2 1 − ρ 1 2 ) + 1 2 k 1 ( 1 − ρ 1 2 ) , λ 2 = − 1 2 ( 1 + 3 ρ 1 2 1 − ρ 1 2 ) ( 2 Δ 1 − k 1 ( 1 − ρ 1 2 ) ) ,
which are purely real for all ρ 1 ∈ [ 0 , 1 )⁠. Consequently, for normal hyperbolicity, we require that
λ 1 ≠ 0 ⇒ 1 > ( 1 − ρ 1 2 ) 2 1 + ρ 1 2 ≠ 2 Δ 1 k 1 , λ 2 ≠ 0 ⇒ 1 > ( 1 − ρ 1 2 ) ≠ 2 Δ 1 k 1 .
Therefore, the critical manifold (12) is normally hyperbolic in the interval ρ 1 ∈ [ 0 , 1 ) if and only if k 1 < 2 Δ 1⁠. Moreover, for stability notice that if k 1 < 2 Δ 1 then λ 1 , 2 < 0 for every ρ 1 ∈ [ 0 , 1 )⁠. Finally, if now we consider | Ω | positive and sufficiently small, the normal hyperbolicity and attractiveness of the Ω = 0 case are persistent as the eigenvalues λ 1 , 2 of (13) depend differentiably on Ω⁠.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 ρ 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 (6), under Proposition III.1, for k 1 adaptive.

Remark 4
It is important to highlight that the identification of two branches in (12), namely, C 0 + and C 0 −⁠, is done for convenience of notation. Nevertheless, by definition (see Definition A.1 in the  Appendix), every point in these branches conforms the critical manifold, regardless of its location in a particular segment, as shown in Figs. 2 and 4 for the coevolutionary inter- and intracoupling, respectively.
FIG. 2.
Locally exponentially stable intercoupling chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray), 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 μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 ). Parameters: Δ 1 = 0.6, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 2.0, γ μ = 1.0, η μ = 2.5, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 300.
View largeDownload slide

Locally exponentially stable intercoupling chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray), 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 μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 )⁠. Parameters: Δ 1 = 0.6⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1 + 0.01⁠, k 1 = 0.9⁠, k 2 = 2.0⁠, γ μ = 1.0⁠, η μ = 2.5⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 300⁠.

FIG. 2.
Locally exponentially stable intercoupling chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray), 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 μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 ). Parameters: Δ 1 = 0.6, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 2.0, γ μ = 1.0, η μ = 2.5, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 300.
View largeDownload slide

Locally exponentially stable intercoupling chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray), 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 μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 )⁠. Parameters: Δ 1 = 0.6⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1 + 0.01⁠, k 1 = 0.9⁠, k 2 = 2.0⁠, γ μ = 1.0⁠, η μ = 2.5⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 300⁠.

Close modal
As a comparison to the adaptive intercoupling case of Sec. III A, we now analyze the effect of one coevolutive intracoupling. By letting k 1 adaptive and with the result in Proposition III.1, now (6) reduces to
(15)
with the same variable interpretations as before and k 1 as the coevolutive coupling strength. Thus, the critical manifold of (15) is
(16)
Remark 5

Notice that (16) can have complex elements, producing disconnected branches of C 0 ±⁠. 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)

Proposition III.3
(Connected critical manifold)

The critical set (16) of (15) is completely connected if and only if | μ | ≥ | Ω | / 3⁠.

Proof.
In order to produce only real components, we require that
( μ ρ 1 ) 2 − ( 2 Ω 3 ρ 1 2 + 1 ) 2 ≥ 0 ,
which can be rewritten as
(17)
By calculating the derivative of the right-hand side of (17) with respect to ρ 1⁠, its maximum occurs at ρ 1 = 1 / 3⁠, for which the right-hand side of (17) reaches a value of 1 / 2 3⁠. Hence, the critical manifold (16) is completely connected if and only if | μ | ≥ | Ω | / 3⁠.
Corollary III.1

When Ω = 0⁠, the critical manifold (16) is completely connected for every μ ∈ R⁠.

Next, we obtain the Jacobian of (15) with respect to the fast variables, which by evaluation at (16) yields
(18)
with h ( ρ 1 ) := ( μ ρ 1 ) 2 − ( 2 Ω 3 ρ 1 2 + 1 ) 2⁠, non-negative under Proposition III.3. Thus, the following result summarizes the conditions for the normal hyperbolicity and stability of (16).

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

Proposition III.4
(Intracoupling normal hyperbolicity and stability ( | Ω | ≪ 1 )⁠)

The critical manifold (16) of (15) has one normally hyperbolic and attracting branch ( C 0 − )⁠, and one unstable branch ( C 0 + ) along the interval ρ 1 ∈ ( 0 , 1 ) for every μ ≠ 0⁠, with Ω = 0⁠. Moreover, for | Ω | ≪ 1⁠, the normal hyperbolicity and stability conditions of the Ω = 0 case are preserved.

Proof.
Observe that for Ω = 0 the Jacobian (18) reduces to
J f | C 0 ± = [ − 2 Δ 1 ( ρ 1 2 1 − ρ 1 2 ) ± 1 2 ( 1 − ρ 1 2 ρ 1 ) | μ | 0 0 ± 1 2 ( 3 ρ 1 2 + 1 ρ 1 ) | μ | ] ,
with eigenvalues in the form
λ 1 ± = − 2 Δ 1 ( ρ 1 2 1 − ρ 1 2 ) ± 1 2 ( 1 − ρ 1 2 ρ 1 ) | μ | , λ 2 ± = ± 1 2 ( 3 ρ 1 2 + 1 ρ 1 ) | μ | .
Hence, the branch of (16) 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 (18) depend differentiably on Ω⁠.42 
Remark 6

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)

Theorem III.1
(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 ± 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 of 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,44 for ζ σ σ ′⁠, with σ , σ ′ = 1 , 2⁠, denoting the type of coevolutive coupling considered, either inter- ( ζ σ σ ′ = ζ 12 = ζ 21 = μ )⁠, or intracoupling ( ζ 11 = k 1 , ζ 22 = k 2 )⁠. 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 (6) 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 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)

Theorem III.2
(Stable chimera states on a network)

Under normally hyperbolic conditions of the critical manifold C 0 ±⁠, 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 ε 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 )⁠, 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 ε 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.

In Sec. IV, 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 (6) when considering the inter- and intracoupling adaptive scenarios discussed in this section.

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,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 

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 , μ ) and ( ρ 1 , μ ) planes of (1) and (11), with μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 )⁠, where, for distinction purposes, ρ 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 Δ 1 and thus the critical manifold (gray) is normally hyperbolic and attracting for every ρ 1 ∈ ( 0 , 1 )⁠. Moreover, observe that the behavior of 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 (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 μ ˙ = ε μ cos ⁡ ( 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 ρ 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.

FIG. 3.
Locally exponentially stable intercoupling breathing chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray) and response (red). Lower row: network (left) and mean-field (right) order parameter evolution on time. For the network, second (blue) and first (orange) population synchronization level and the filtered version of the latter (red). Adaptive law μ ˙ = ε μ cos ⁡ ( 0.02 t ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.99 , − 0.5 , 1.1 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 9.0, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 1000.
View largeDownload slide

Locally exponentially stable intercoupling breathing chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray) and response (red). Lower row: network (left) and mean-field (right) order parameter evolution on time. For the network, second (blue) and first (orange) population synchronization level and the filtered version of the latter (red). Adaptive law μ ˙ = ε μ cos ⁡ ( 0.02 t ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.99 , − 0.5 , 1.1 )⁠. Parameters: Δ 1 = 1.0⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1 + 0.01⁠, k 1 = 0.9⁠, k 2 = 9.0⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 1000⁠.

FIG. 3.
Locally exponentially stable intercoupling breathing chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray) and response (red). Lower row: network (left) and mean-field (right) order parameter evolution on time. For the network, second (blue) and first (orange) population synchronization level and the filtered version of the latter (red). Adaptive law μ ˙ = ε μ cos ⁡ ( 0.02 t ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.99 , − 0.5 , 1.1 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 9.0, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 1000.
View largeDownload slide

Locally exponentially stable intercoupling breathing chimera. Upper row: network (left) and mean-field (right) projections to the ( R 1 , μ ) and ( ρ 1 , μ ) planes for (1) and (11), respectively. Attracting critical manifold (gray) and response (red). Lower row: network (left) and mean-field (right) order parameter evolution on time. For the network, second (blue) and first (orange) population synchronization level and the filtered version of the latter (red). Adaptive law μ ˙ = ε μ cos ⁡ ( 0.02 t ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.99 , − 0.5 , 1.1 )⁠. Parameters: Δ 1 = 1.0⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1 + 0.01⁠, k 1 = 0.9⁠, k 2 = 9.0⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 1000⁠.

Close modal

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 ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 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 −⁠, 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 ρ 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 Δ 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.

FIG. 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 (15), respectively. Critical manifold regions: attracting (gray), 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 k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.01, ω 1 = 101 / 20, ω 2 = ω 1, μ = 0.1, k 2 = 4.0, γ 1 = 2.5, η 1 = 0.5, ε 1 = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500.
View largeDownload slide

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 (15), respectively. Critical manifold regions: attracting (gray), 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 k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 )⁠. Parameters: Δ 1 = 1.0⁠, Δ 2 = 0.01⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1⁠, μ = 0.1⁠, k 2 = 4.0⁠, γ 1 = 2.5⁠, η 1 = 0.5⁠, ε 1 = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠.

FIG. 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 (15), respectively. Critical manifold regions: attracting (gray), 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 k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.01, ω 1 = 101 / 20, ω 2 = ω 1, μ = 0.1, k 2 = 4.0, γ 1 = 2.5, η 1 = 0.5, ε 1 = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500.
View largeDownload slide

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 (15), respectively. Critical manifold regions: attracting (gray), 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 k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 )⁠. Parameters: Δ 1 = 1.0⁠, Δ 2 = 0.01⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1⁠, μ = 0.1⁠, k 2 = 4.0⁠, γ 1 = 2.5⁠, η 1 = 0.5⁠, ε 1 = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠.

Close modal
Remark 7
Despite the well-known importance of the phase-lag parameter in determining the stability of a chimera state,19,52–54 the assumption βσσ′=0⁠, for σ,σ′=1,2⁠, 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,43 the resulting mean-field equations represent our reduced system in the limit βσσ′→0⁠. 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 β11=β12=β21=β22=π/6⁠. Similarly, in Fig. 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 C0± in gray until it reaches the intersection between Cε and the slow adaptation nullcline. Notice that the apparent shift between Cε and C0± 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.
FIG. 5.
Generation of a stable chimera state on a network with phase-lag and a coevolutive intercoupling strength. Projection to the ( R 1 , μ ) plane (left) and local order parameters (right) for the coevolutive intercoupling case, respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 ). Parameters: Δ 1 = 0.6, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 2.0, γ μ = 1.0, η μ = 2.5, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = π / 6, and N 1 , 2 = 500.
View largeDownload slide

Generation of a stable chimera state on a network with phase-lag and a coevolutive intercoupling strength. Projection to the (R1,μ) plane (left) and local order parameters (right) for the coevolutive intercoupling case, respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law μ˙=εμ(−μ+γμ−ημρ1) and initial conditions (ρ1,ψ,μ)=(0.9,−0.5,2.0)⁠. Parameters: Δ1=0.6⁠, Δ2=0.1⁠, ω1=101/20⁠, ω2=ω1+0.01⁠, k1=0.9⁠, k2=2.0⁠, γμ=1.0⁠, ημ=2.5⁠, εμ=0.02⁠, β11=β12=β21=β22=π/6⁠, and N1,2=500⁠.

FIG. 5.
Generation of a stable chimera state on a network with phase-lag and a coevolutive intercoupling strength. Projection to the ( R 1 , μ ) plane (left) and local order parameters (right) for the coevolutive intercoupling case, respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , − 0.5 , 2.0 ). Parameters: Δ 1 = 0.6, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1 + 0.01, k 1 = 0.9, k 2 = 2.0, γ μ = 1.0, η μ = 2.5, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = π / 6, and N 1 , 2 = 500.
View largeDownload slide

Generation of a stable chimera state on a network with phase-lag and a coevolutive intercoupling strength. Projection to the (R1,μ) plane (left) and local order parameters (right) for the coevolutive intercoupling case, respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law μ˙=εμ(−μ+γμ−ημρ1) and initial conditions (ρ1,ψ,μ)=(0.9,−0.5,2.0)⁠. Parameters: Δ1=0.6⁠, Δ2=0.1⁠, ω1=101/20⁠, ω2=ω1+0.01⁠, k1=0.9⁠, k2=2.0⁠, γμ=1.0⁠, ημ=2.5⁠, εμ=0.02⁠, β11=β12=β21=β22=π/6⁠, and N1,2=500⁠.

Close modal
FIG. 6.
Generation of a stable chimera state on a network with phase-lag and a coevolutive intracoupling strength. Projection to the ( R 1 , k 1 ) plane (left) and local order parameters (right) for the coevolutive intracoupling case, respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.01, ω 1 = 101 / 20, ω 2 = ω 1, μ = 0.1, k 2 = 4.0, γ 1 = 2.5, η 1 = 0.5, ε 1 = 0.02, β 11 = β 12 = β 21 = β 22 = π / 10, and N 1 , 2 = 1000.
View largeDownload slide

Generation of a stable chimera state on a network with phase-lag and a coevolutive intracoupling strength. Projection to the (R1,k1) plane (left) and local order parameters (right) for the coevolutive intracoupling case, respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law k˙1=ε1(−k1+γ1−η1ρ1) and initial conditions (ρ1,ψ,k1)=(0.9,−0.5,4.0)⁠. Parameters: Δ1=1.0⁠, Δ2=0.01⁠, ω1=101/20⁠, ω2=ω1⁠, μ=0.1⁠, k2=4.0⁠, γ1=2.5⁠, η1=0.5⁠, ε1=0.02⁠, β11=β12=β21=β22=π/10⁠, and N1,2=1000⁠.

FIG. 6.
Generation of a stable chimera state on a network with phase-lag and a coevolutive intracoupling strength. Projection to the ( R 1 , k 1 ) plane (left) and local order parameters (right) for the coevolutive intracoupling case, respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law k ˙ 1 = ε 1 ( − k 1 + γ 1 − η 1 ρ 1 ) and initial conditions ( ρ 1 , ψ , k 1 ) = ( 0.9 , − 0.5 , 4.0 ). Parameters: Δ 1 = 1.0, Δ 2 = 0.01, ω 1 = 101 / 20, ω 2 = ω 1, μ = 0.1, k 2 = 4.0, γ 1 = 2.5, η 1 = 0.5, ε 1 = 0.02, β 11 = β 12 = β 21 = β 22 = π / 10, and N 1 , 2 = 1000.
View largeDownload slide

Generation of a stable chimera state on a network with phase-lag and a coevolutive intracoupling strength. Projection to the (R1,k1) plane (left) and local order parameters (right) for the coevolutive intracoupling case, respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue), slow adaptation nullcline (magenta), first population synchronization level (orange) and its filtered version (red), and second population order parameter (blue). Adaptive law k˙1=ε1(−k1+γ1−η1ρ1) and initial conditions (ρ1,ψ,k1)=(0.9,−0.5,4.0)⁠. Parameters: Δ1=1.0⁠, Δ2=0.01⁠, ω1=101/20⁠, ω2=ω1⁠, μ=0.1⁠, k2=4.0⁠, γ1=2.5⁠, η1=0.5⁠, ε1=0.02⁠, β11=β12=β21=β22=π/10⁠, and N1,2=1000⁠.

Close modal
Remark 8
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 (6) 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 Fig. 7 for the coevolutionary intercoupling case. Therefore, the inclusion of the adaptation in a network remarkably expands the capability of the previously known results of such systems. 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 Fig. 7.
FIG. 7.
Synchronization of a partially incoherent system. Projection to the ( R 1 , μ ) plane (left) and local order parameters R 1 , 2 ( t ) (right), respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), and response (red). First and second population synchronization levels (orange and blue) and their filtered version (red and green). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.4 , 0.4 , 0.25 ). Parameters: Δ 1 = 0.3, Δ 2 = 0.2, ω 1 = 0.0, ω 2 = 0.3, k 1 = 0.5, k 2 = 0.7, γ μ = 0.2, η μ = − 2.0, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500. Notice that, since k 1 < 2 Δ 1, the observed partial synchronization is due to the adaptation, as otherwise the first population would tend to an incoherent state.
View largeDownload slide

Synchronization of a partially incoherent system. Projection to the ( R 1 , μ ) plane (left) and local order parameters R 1 , 2 ( t ) (right), respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), and response (red). First and second population synchronization levels (orange and blue) and their filtered version (red and green). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.4 , 0.4 , 0.25 )⁠. Parameters: Δ 1 = 0.3⁠, Δ 2 = 0.2⁠, ω 1 = 0.0⁠, ω 2 = 0.3⁠, k 1 = 0.5⁠, k 2 = 0.7⁠, γ μ = 0.2⁠, η μ = − 2.0⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠. Notice that, since k 1 < 2 Δ 1⁠, the observed partial synchronization is due to the adaptation, as otherwise the first population would tend to an incoherent state.

FIG. 7.
Synchronization of a partially incoherent system. Projection to the ( R 1 , μ ) plane (left) and local order parameters R 1 , 2 ( t ) (right), respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), and response (red). First and second population synchronization levels (orange and blue) and their filtered version (red and green). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.4 , 0.4 , 0.25 ). Parameters: Δ 1 = 0.3, Δ 2 = 0.2, ω 1 = 0.0, ω 2 = 0.3, k 1 = 0.5, k 2 = 0.7, γ μ = 0.2, η μ = − 2.0, ε μ = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500. Notice that, since k 1 < 2 Δ 1, the observed partial synchronization is due to the adaptation, as otherwise the first population would tend to an incoherent state.
View largeDownload slide

Synchronization of a partially incoherent system. Projection to the ( R 1 , μ ) plane (left) and local order parameters R 1 , 2 ( t ) (right), respectively. Attracting critical manifold (gray), slow adaptation nullcline (magenta), and response (red). First and second population synchronization levels (orange and blue) and their filtered version (red and green). Adaptive law μ ˙ = ε μ ( − μ + γ μ − η μ ρ 1 ) and initial conditions ( ρ 1 , ψ , μ ) = ( 0.4 , 0.4 , 0.25 )⁠. Parameters: Δ 1 = 0.3⁠, Δ 2 = 0.2⁠, ω 1 = 0.0⁠, ω 2 = 0.3⁠, k 1 = 0.5⁠, k 2 = 0.7⁠, γ μ = 0.2⁠, η μ = − 2.0⁠, ε μ = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠. Notice that, since k 1 < 2 Δ 1⁠, the observed partial synchronization is due to the adaptation, as otherwise the first population would tend to an incoherent state.

Close modal

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, μ ( 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.

FIG. 8.
Stable breathing chimera state in a two-coevolutionary network. Local order parameters (left) and adaptive coupling strengths (right), respectively. First and second population synchronization level (orange, blue) and their filtered versions (red, green). Adaptive laws μ ˙ = ε ( − 2 ( k 1 − 0.5 ) ) (black), k ˙ 1 = ε ( μ − 0.5 ) (dashed blue), and initial conditions ( ρ 1 , ψ , μ , k 1 ) = ( 0.9 , 0.1 , 0.5 , 0.1 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1, k 2 = 1.0, ε = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500. In this figure, even if k 1 , μ 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.
View largeDownload slide

Stable breathing chimera state in a two-coevolutionary network. Local order parameters (left) and adaptive coupling strengths (right), respectively. First and second population synchronization level (orange, blue) and their filtered versions (red, green). Adaptive laws μ ˙ = ε ( − 2 ( k 1 − 0.5 ) ) (black), k ˙ 1 = ε ( μ − 0.5 ) (dashed blue), and initial conditions ( ρ 1 , ψ , μ , k 1 ) = ( 0.9 , 0.1 , 0.5 , 0.1 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1⁠, k 2 = 1.0⁠, ε = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠. In this figure, even if k 1 , μ 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.

FIG. 8.
Stable breathing chimera state in a two-coevolutionary network. Local order parameters (left) and adaptive coupling strengths (right), respectively. First and second population synchronization level (orange, blue) and their filtered versions (red, green). Adaptive laws μ ˙ = ε ( − 2 ( k 1 − 0.5 ) ) (black), k ˙ 1 = ε ( μ − 0.5 ) (dashed blue), and initial conditions ( ρ 1 , ψ , μ , k 1 ) = ( 0.9 , 0.1 , 0.5 , 0.1 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 101 / 20, ω 2 = ω 1, k 2 = 1.0, ε = 0.02, β 11 = β 12 = β 21 = β 22 = 0, and N 1 , 2 = 500. In this figure, even if k 1 , μ 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.
View largeDownload slide

Stable breathing chimera state in a two-coevolutionary network. Local order parameters (left) and adaptive coupling strengths (right), respectively. First and second population synchronization level (orange, blue) and their filtered versions (red, green). Adaptive laws μ ˙ = ε ( − 2 ( k 1 − 0.5 ) ) (black), k ˙ 1 = ε ( μ − 0.5 ) (dashed blue), and initial conditions ( ρ 1 , ψ , μ , k 1 ) = ( 0.9 , 0.1 , 0.5 , 0.1 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 101 / 20⁠, ω 2 = ω 1⁠, k 2 = 1.0⁠, ε = 0.02⁠, β 11 = β 12 = β 21 = β 22 = 0⁠, and N 1 , 2 = 500⁠. In this figure, even if k 1 , μ 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.

Close modal
FIG. 9.
Canard cycle in a coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = ε μ ( 0.708 065 169 053 771 149 − ρ 1 ), and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 0.0, ω 2 = 0.4, k 1 = 2.0, k 2 = 3.0, ε μ = 0.02, and β 11 = β 12 = β 21 = β 22 = 0.
View largeDownload slide

Canard cycle in a coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = ε μ ( 0.708 065 169 053 771 149 − ρ 1 )⁠, and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 0.0⁠, ω 2 = 0.4⁠, k 1 = 2.0⁠, k 2 = 3.0⁠, ε μ = 0.02⁠, and β 11 = β 12 = β 21 = β 22 = 0⁠.

FIG. 9.
Canard cycle in a coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = ε μ ( 0.708 065 169 053 771 149 − ρ 1 ), and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 0.0, ω 2 = 0.4, k 1 = 2.0, k 2 = 3.0, ε μ = 0.02, and β 11 = β 12 = β 21 = β 22 = 0.
View largeDownload slide

Canard cycle in a coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = ε μ ( 0.708 065 169 053 771 149 − ρ 1 )⁠, and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 0.0⁠, ω 2 = 0.4⁠, k 1 = 2.0⁠, k 2 = 3.0⁠, ε μ = 0.02⁠, and β 11 = β 12 = β 21 = β 22 = 0⁠.

Close modal
FIG. 10.
Transcritical canard cycle in the coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = − ε μ cos ⁡ ψ and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 0.0, ω 2 = ω 1, k 1 = 2.25, k 2 = 3.0, ε μ = 0.02, and β 11 = β 12 = β 21 = β 22 = 0.
View largeDownload slide

Transcritical canard cycle in the coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = − ε μ cos ⁡ ψ and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 0.0⁠, ω 2 = ω 1⁠, k 1 = 2.25⁠, k 2 = 3.0⁠, ε μ = 0.02⁠, and β 11 = β 12 = β 21 = β 22 = 0⁠.

FIG. 10.
Transcritical canard cycle in the coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = − ε μ cos ⁡ ψ and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 ). Parameters: Δ 1 = 0.5, Δ 2 = 0.1, ω 1 = 0.0, ω 2 = ω 1, k 1 = 2.25, k 2 = 3.0, ε μ = 0.02, and β 11 = β 12 = β 21 = β 22 = 0.
View largeDownload slide

Transcritical canard cycle in the coevolutive intercoupling scenario. Projection to the ( ρ 1 , μ ) plane (left) and local order parameter ρ 1 ( t ) (right) of the mean-field (11), respectively. Critical manifold regions: attracting (gray), saddle (dashed green), repelling (blue dotted), and first population synchronization level (red). Adaptive law μ ˙ = − ε μ cos ⁡ ψ and initial conditions ( ρ 1 , ψ , μ ) = ( 0.9 , 0.1 , 0.33 )⁠. Parameters: Δ 1 = 0.5⁠, Δ 2 = 0.1⁠, ω 1 = 0.0⁠, ω 2 = ω 1⁠, k 1 = 2.25⁠, k 2 = 3.0⁠, ε μ = 0.02⁠, and β 11 = β 12 = β 21 = β 22 = 0⁠.

Close modal

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.

The authors have no conflicts to disclose.

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

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

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

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. Therefore, we begin by studying an all-to-all undirected network composed of different populations of heterogeneous Kuramoto phase oscillators, which can be expressed as
(A1)
where θ i σ and N σ denote the state of the ith 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
(A2)
with ω ^ σ and Δ σ being 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, in our context, β σ σ ′ plays the role of a regular perturbation parameter for β σ σ ′ small. The case β σ σ ′ ≈ π 2 is relevant but is left for future research as it leads to extra complications in the mean-field analysis.
Considering the thermodynamic limit (⁠ N σ → ∞ for each σ⁠, but σ finite), the state of (A1) 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
(A3)
where v σ ( θ σ , ω σ , t ) is the continuous version of the angular velocity of the oscillators in each population expressed as
(A4)
Moreover, the well-known Kuramoto complex order parameter is defined as
(A5)
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.39 Subsequently, it is possible to express the angular velocity (A4) in terms of the order parameter (A5), resulting in
(A6)
where the overbar 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⁠,55 and that in general may depend on the complex order parameter (A5). 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.45 In particular, for our problem consisting on a multilayer network of Kuramoto-type phase oscillators with coevolutive inter- and intracouplings, the functions H σ ( t ) are defined as H σ ( t ) = ∑ σ ′ = 1 M k σ σ ′ z σ ′ ( t )⁠, where the effective forcing depends on the population index σ , σ ′ = 1 , 2 , … , M⁠.
Following Ott and Antonsen,41 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.43 Thus, by considering the Fourier series of the form
(A7)
for the density function of each population and by imposing the ansatz (A7) to the continuum equation (A3), we obtain the following set of differential complex equations:
(A8)
where c.c. stands for the complex conjugate of the expression on the right-hand side of (A8). Moreover, since the sum component in (A8) is not zero, it means that the term within the parentheses vanish identically if the ansatz (A7) is a valid solution to the continuity equation (A3). Hence,
(A9)
with z ¯ σ ( t ) = ∫ − ∞ ∞ α σ ( ω σ , t ) g σ ( ω σ ) d ω σ⁠. Furthermore, consider solutions of (A9) with the initial conditions α σ ( ω σ , 0 ) that satisfy: (i) | α σ ( ω σ , t ) | ≤ 1⁠; (ii) α σ ( ω σ , 0 ) is analytically continuable into the lower half plane Im ( ω σ ) < 0⁠; and (iii) | α σ ( ω σ , t ) | → 0 as Im ( ω σ ) → − ∞⁠. If such conditions are satisfied for α σ ( ω σ , 0 )⁠, then they will continue to be satisfied for α σ ( ω σ , t )⁠.41,56 It is important to mention that the Ott–Antonsen ansatz corresponds to a particular case of the technique developed by Watanabe–Strogatz57,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., Δ σ > 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 
By selecting the probability distribution functions as (A2), 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 (A3) produces M-coupled complex ODEs describing the evolution of the complex order parameters as
z ˙ σ + ( Δ σ − ı ω σ ) z σ + 1 2 ∑ σ ′ = 1 M k σ σ ′ ( z ¯ σ ′ z σ 2 − z σ ′ ) = 0.
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⁠,43 produced by orthogonality 2 M coupled real ODEs in the form
(A10)
In particular, our mechanism for the generation of the stable chimera state relies on geometric properties of (A10) with adaptive coupling strengths under a macroscopic slow adaptation, which induces a time-scale separation between the dynamics on and off the network. Next, Section 2 of the Appendix briefly summarizes the principal concepts of fast–slow systems employed in our work.

2. Fast–slow dynamics

A fast–slow system is a singularly perturbed ordinary differential equation in the form
(A11)
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 (A11) as
(A12)
where the prime denotes the derivative with respect to the fast-time t⁠. Hence, (A11) and (A12) 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 (A11) and (A12). Despite being non-equivalent, both problems share a close connection as expressed in the following definition.
(Critical manifold)
Definition A.1
(Critical manifold)
The critical set is defined as
(A13)

Moreover, if C 0 is a submanifold of R m × 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)
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 (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 theorem64,65 stated as follows.

(Fenichel’s Theorem)
Theorem A.1
(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 ∈ C k < ∞⁠, for 0 < ε ≪ 1⁠, then

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

  2. S ε has Hausdorff distance d H ( V , W ) = O ( ε ) from S 0⁠, as ε → 0⁠.

  3. The flow on S ε converges to the flow generated by the reduced problem, known as the slow flow, as ε → 0⁠.

  4. S ε is C k-smooth.

  5. S ε is normally hyperbolic and has the same stability properties with respect to the fast variables as S 0⁠.

  6. S ε is usually not unique. In regions with a fixed distance from ∂ S ε⁠, all manifolds which satisfies conditions 1–5 above are exponentially close to each other with a Hausdorff distance of order O ( exp ( − C / ε ) ) for some C > 0⁠, C ∈ O ( 1 )⁠, as ε → 0⁠.

Thus, any manifold S ε satisfying conditions 1–5 is called a slow manifold.

Remark 9

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

1.
A.
Pikovsky
,
M.
Rosemblum
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
2001
).
Google Scholar
Crossref
Search ADS
 
2.
S.
Boccaletti
,
J.
Kurths
,
G.
Osipov
,
D. L.
Valladares
, and
C. S.
Zhou
, “
The synchronization of chaotic systems
,”
Phys. Rep.
366
,
1
–
101
(
2002
).
https://doi.org/10.1016/S0370-1573(02)00137-0
Google Scholar
Crossref
Search ADS
 
3.
P. J.
Uhlhaas
 and
W.
Singer
, “
Neural synchrony in brain disorders: Relevance for cognitive dysfunctions and pathophysiology
,”
Neuron
52
,
155
–
168
(
2006
).
https://doi.org/10.1016/j.neuron.2006.09.020
Google Scholar
Crossref
Search ADS
PubMed
 
4.
M.
Shafiei
,
S.
Jafari
,
F.
Parastesh
,
M.
Ozer
,
T.
Kapitaniak
, and
M.
Perc
, “
Time delayed chemical synapses and synchronization in multilayer neuronal networks with ephatic inter layer coupling
,”
Commun. Nonlinear Sci. Numer. Simul.
84
,
105175
(
2020
).
https://doi.org/10.1016/j.cnsns.2020.105175
Google Scholar
Crossref
Search ADS
 
5.
S.
Yamaguchi
,
H.
Isejima
,
T.
Matsuo
,
R.
Okura
,
K.
Yagita
,
M.
Kobayashi
, and
H.
Okamura
, “
Synchronization of cellular clocks in the suprachiasmatic nucleus
,”
Science
302
,
1408
–
1412
(
2003
).
https://doi.org/10.1126/science.1089287
Google Scholar
Crossref
Search ADS
PubMed
 
6.
J.
Husse
,
G.
Eichele
, and
H.
Oster
, “
Synchronization of the mammalian circadian timing system: Light can control peripheral clocks independently of the SCN clock
,”
BioEssays
37
,
1119
–
1128
(
2015
).
https://doi.org/10.1002/bies.201500026
Google Scholar
Crossref
Search ADS
PubMed
 
7.
M.
Rohden
,
A.
Sorge
,
M.
Timme
, and
D.
Witthaut
, “
Self-organized synchronization in decentralized power grids
,”
Phys. Rev. Lett.
109
,
064101
(
2012
).
https://doi.org/10.1103/PhysRevLett.109.064101
Google Scholar
Crossref
Search ADS
PubMed
 
8.
A.
Sajadi
,
R. W.
Kenyon
, and
B.
Hodge
, “
Synchronization in electric power networks with inherent heterogeneity up to 100% inverter-based renewable generation
,”
Nat. Commun.
13
,
2490
(
2022
).
https://doi.org/10.1038/s41467-022-30164-3
Google Scholar
Crossref
Search ADS
PubMed
 
9.
J.
Buck
 and
E.
Buck
, “
Mechanisms of rhythmic synchronous flashing of fireflies
,”
Science
159
,
1319
–
1327
(
1968
).
https://doi.org/10.1126/science.159.3821.1319
Google Scholar
Crossref
Search ADS
PubMed
 
10.
A.
Moiseff
 and
J.
Copeland
, “
Firefly synchrony: A behavioral strategy to minimize visual clutter
,”
Science
329
,
181
(
2010
).
https://doi.org/10.1126/science.1190421
Google Scholar
Crossref
Search ADS
PubMed
 
11.
K.
Wiesenfield
,
P.
Colet
, and
S. H.
Strogatz
, “
Frequency locking in Josephson arrays: Connection with the Kuramoto model
,”
Phys. Rev. E
57
,
1563
(
1998
).
https://doi.org/10.1103/PhysRevE.57.1563
Google Scholar
Crossref
Search ADS
 
12.
A. F.
Taylor
,
M. R.
Tinsley
,
F.
Wang
,
Z.
Huand
, and
K.
Showalter
, “
Dynamical quorum sensing and synchronization in large populations of chemical oscillators
,”
Science
323
,
614
–
617
(
2009
).
https://doi.org/10.1126/science.1166253
Google Scholar
Crossref
Search ADS
PubMed
 
13.
D.
Călugăru
,
J. F.
Totz
,
E. A.
Martens
, and
H.
Engel
, “
First-order synchronization transition in a large population of strongly coupled relaxation oscillators
,”
Sci. Adv.
6
,
eabb2637
(
2020
).
https://doi.org/10.1126/sciadv.abb2637
Google Scholar
Crossref
Search ADS
PubMed
 
14.
J.
Pantaleone
, “
Stability of incoherence in an isotropic gas of oscillating neutrinos
,”
Phys. Rev. D
58
,
073002
(
1998
).
https://doi.org/10.1103/PhysRevD.58.073002
Google Scholar
Crossref
Search ADS
 
15.
C.
Bick
,
M.
Goodfellow
,
C. R.
Laing
, and
E. A.
Martens
, “
Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: A review
,”
J. Math. Neurosci.
10
,
9
(
2020
).
https://doi.org/10.1186/s13408-020-00086-9
Google Scholar
Crossref
Search ADS
PubMed
 
16.
M. J. B.
Hauser
, “
Synchronization of glycolytic activity in yeast cells
,”
Curr. Genet.
68
,
69
–
81
(
2022
).
https://doi.org/10.1007/s00294-021-01214-y
Google Scholar
Crossref
Search ADS
PubMed
 
17.
A. T.
Winfree
, “
Biological rhythms and the behavior of populations of coupled oscillators
,”
J. Theor. Biol.
16
,
15
–
42
(
1967
).
https://doi.org/10.1016/0022-5193(67)90051-3
Google Scholar
Crossref
Search ADS
PubMed
 
18.
Y.
Kuramoto
,
Chemical Oscillations, Waves, and Turbulence
(
Springer
,
Berlin
,
1984
).
Google Scholar
Crossref
Search ADS
 
19.
Y.
Kuramoto
 and
D.
Battogtokh
, “
Coexistence of coherence and incoherence in nonlocally coupled phase oscillators
,”
Nonlinear Phenom. Complex Syst.
5
,
380–385
(
2002
).
Google Scholar
 
20.
O. E.
Omel’chenko
, “
The mathematics behind chimera states
,”
Nonlinearity
31
,
R121
(
2018
).
https://doi.org/10.1088/1361-6544/aaaa07
Google Scholar
Crossref
Search ADS
 
21.
D. M.
Abrams
 and
S. H.
Strogatz
, “
Chimera states for coupled oscillators
,”
Phys. Rev. Lett.
93
,
174102
(
2004
).
https://doi.org/10.1103/PhysRevLett.93.174102
Google Scholar
Crossref
Search ADS
PubMed
 
22.
M. J.
Panaggio
 and
D. M.
Abrams
, “
Chimera states: Coexistence of coherence and incoherence in networks of coupled oscillators
,”
Nonlinearity
28
,
R67
(
2015
).
https://doi.org/10.1088/0951-7715/28/3/R67
Google Scholar
Crossref
Search ADS
 
23.
E.
Schöll
, “
Synchronization patterns and chimera states in complex networks: Interplay of topology and dynamics
,”
Eur. Phys. J. Spec. Top.
225
,
891
–
919
(
2016
).
https://doi.org/10.1140/epjst/e2016-02646-3
Google Scholar
Crossref
Search ADS
 
24.
B. K.
Bera
,
S.
Majhi
,
D.
Ghosh
, and
M.
Perc
, “
Chimera states: Effects of different coupling topologies
,”
Europhys. Lett.
118
,
10001
(
2017
).
https://doi.org/10.1209/0295-5075/118/10001
Google Scholar
Crossref
Search ADS
 
25.
N. C.
Rattenborg
,
C. J.
Amlaner
, and
S. L.
Lima
, “
Behavioral, neurophysiological and evolutionary perspectives on unihemispheric sleep
,”
Neurosci. Biobehav. Rev.
24
,
817
–
842
(
2000
).
https://doi.org/10.1016/S0149-7634(00)00039-7
Google Scholar
Crossref
Search ADS
PubMed
 
26.
J. M.
Davidenko
,
A. V.
Pertsov
,
R.
Salomonsz
,
W.
Baxter
, and
J.
Jalife
, “
Stationary and drifting spiral waves of excitation in isolated cardiac muscle
,”
Nature
355
,
349
–
351
(
1992
).
https://doi.org/10.1038/355349a0
Google Scholar
Crossref
Search ADS
PubMed
 
27.
A. V.
Panfilov
, “
Spiral breakup as a model of ventricular fibrillation
,”
Chaos
8
,
57
–
64
(
1998
).
https://doi.org/10.1063/1.166287
Google Scholar
Crossref
Search ADS
PubMed
 
28.
J. C.
González-Avella
,
M. G.
Cosenza
, and
M.
San Miguel
, “
Localized coherence in two interacting populations of special agents
,”
Phys. A
399
,
24
–
30
(
2014
).
https://doi.org/10.1016/j.physa.2013.12.035
Google Scholar
Crossref
Search ADS
 
29.
H.
Sakaguchi
, “
Instability of synchronized motion in nonlocally coupled neural oscillators
,”
Phys. Rev. E
73
,
031907
(
2006
).
https://doi.org/10.1103/PhysRevE.73.031907
Google Scholar
Crossref
Search ADS
 
30.
S.
Olmi
,
A.
Politi
, and
A.
Torcini
, “
Collective chaos in pulse-coupled neural networks
,”
Europhys. Lett.
92
,
60007
(
2011
).
https://doi.org/10.1209/0295-5075/92/60007
Google Scholar
Crossref
Search ADS
 
31.
I.
Omelchenko
,
O. E.
Omel’chenko
,
P.
Hövel
, and
E.
Schöll
, “
When nonlocal coupling between oscillators becomes stronger: Patched synchrony or multichimera states
,”
Phys. Rev. Lett.
110
,
224101
(
2013
).
https://doi.org/10.1103/PhysRevLett.110.224101
Google Scholar
Crossref
Search ADS
PubMed
 
32.
E. A.
Martens
 and
K.
Klemm
, “
Transitions from trees to cycles in adaptive flow networks
,”
Front. Phys.
5
,
1
–
10
(
2017
).
https://doi.org/10.3389/fphy.2017.00062
Google Scholar
Crossref
Search ADS
 
33.
L.
Zino
,
M.
Ye
, and
M.
Cao
, “A coevolutionary model for actions and opinions in social networks,” in
2020 59th IEEE Conference on Decision and Control (CDC)
(IEEE, 2020), pp. 1110–1115.
https://doi.org/10.1109/CDC42340.2020.9303954
34.
L.
Zino
,
M.
Ye
, and
M.
Cao
, “
A two-layer model for coevolving opinion dynamics and collective decision-making in complex social systems
,”
Chaos
30
,
083107
(
2020
).
https://doi.org/10.1063/5.0004787
Google Scholar
Crossref
Search ADS
PubMed
 
35.
R.
Berner
,
S.
Yanchuk
, and
E.
Schöll
, “
What adaptive neuronal networks teach us about power grids
,”
Phys. Rev. E
103
,
042315
(
2021
).
https://doi.org/10.1103/PhysRevE.103.042315
Google Scholar
Crossref
Search ADS
PubMed
 
36.
W.
Gernster
,
R.
Kempter
,
J. L.
van Hemmen
, and
H.
Wagner
, “
A neuronal learning rule for sub-millisecond temporal coding
,”
Nature
383
,
76
–
78
(
1996
).
https://doi.org/10.1038/383076a0
Google Scholar
Crossref
Search ADS
PubMed
 
37.
W.
Gernster
 and
W. M.
Kistler
, “
Mathematical formulations of Hebbian learning
,”
Biol. Cybern.
87
,
404
–
415
(
2002
).
https://doi.org/10.1007/s00422-002-0353-y
Google Scholar
Crossref
Search ADS
PubMed
 
38.
A.
Goriely
,
E.
Kuhl
, and
C.
Bick
, “
Neuronal oscillations on evolving networks: Dynamics, damage, degradation, decline, dementia, and death
,”
Phys. Rev. Lett.
125
,
128102
(
2020
).
https://doi.org/10.1103/PhysRevLett.125.128102
Google Scholar
Crossref
Search ADS
PubMed
 
39.
P.
So
 and
E.
Barreto
, “
Generating macroscopic chaos in a network of globally coupled phase oscillators
,”
Chaos
21
,
033127
(
2011
).
https://doi.org/10.1063/1.3638441
Google Scholar
Crossref
Search ADS
PubMed
 
40.
M.
Ciszak
,
F.
Marino
,
A.
Torcini
, and
S.
Olmi
, “
Emergent excitability in populations of nonexcitable units
,”
Phys. Rev. E
102
,
050201(R)
(
2020
).
https://doi.org/10.1103/PhysRevE.102.050201
Google Scholar
Crossref
Search ADS
 
41.
E.
Ott
 and
T. M.
Antonsen
, “
Low dimensional behavior of large systems of globally coupled oscillators
,”
Chaos
18
,
037113
(
2008
).
https://doi.org/10.1063/1.2930766
Google Scholar
Crossref
Search ADS
PubMed
 
42.
P. D.
Lax
, “Linear algebra and its applications,” 2nd ed. in Pure and Applied Mathematics: A Wiley-Interscience Series of Texts, Monographs and Tracts (John Wiley & Sons, 2007).
43.
D. M.
Abrams
,
R.
Mirollo
,
S.
Strogatz
, and
D. A.
Wiley
, “
Solvable model for chimera states of coupled oscillators
,”
Phys. Rev. Lett.
101
,
084103
(
2008
).
https://doi.org/10.1103/PhysRevLett.101.084103
Google Scholar
Crossref
Search ADS
PubMed
 
44.
P.
Kokotović
,
H.
Khalil
, and
J.
O’Reilly
,
Singular Perturbation Methods in Control: Analysis and Design
(
Society for Industrial and Applied Mathematics
,
1999
).
Google Scholar
Crossref
Search ADS
 
45.
E.
Ott
 and
T. M.
Antonsen
, “
Long time evolution of phase oscillator systems
,”
Chaos
19
,
023117
(
2009
).
https://doi.org/10.1063/1.3136851
Google Scholar
Crossref
Search ADS
PubMed
 
46.
E.
Ott
,
B. R.
Hunt
, and
T. M.
Antonsen
, “
Comment on ‘Long time evolution of phase oscillator systems’
,”
Chaos
21
,
025112
(
2011
).
https://doi.org/10.1063/1.3574931
Google Scholar
Crossref
Search ADS
PubMed
 
47.
A.
Savitzky
 and
M. J. E.
Golay
, “
Smoothing and differentiation of data by simplified least squares procedures
,”
Anal. Chem.
36
,
1627
–
1639
(
1964
).
https://doi.org/10.1021/ac60214a047
Google Scholar
Crossref
Search ADS
 
48.
R. W.
Schafer
, “
What is a Savitzky-Golay filter? [Lecture notes]
,”
IEEE Signal Process. Mag.
28
,
111
–
117
(
2011
).
https://doi.org/10.1109/MSP.2011.941097
Google Scholar
Crossref
Search ADS
 
49.
T.
Böhle
,
C.
Kuehn
, and
M.
Thalhammer
, “
On the reliable and efficient numerical integration of the Kuramoto model and related dynamical systems on graphs
,”
Int. J. Comput. Math.
99
,
31
–
57
(
2022
).
https://doi.org/10.1080/00207160.2021.1952997
Google Scholar
Crossref
Search ADS
 
50.
O. E.
Omel’chenko
, “
Mathematical framework for breathing chimera states
,”
J. Nonlinear Sci.
32
,
22
(
2022
).
https://doi.org/10.1007/s00332-021-09779-1
Google Scholar
Crossref
Search ADS
 
51.
C. R.
Laing
, “
Chimera states in heterogeneous networks
,”
Chaos
19
,
013113
(
2009
).
https://doi.org/10.1063/1.3068353
Google Scholar
Crossref
Search ADS
PubMed
 
52.
N.
Nakagawa
 and
Y.
Kuramoto
, “
Collective chaos on a population of globally coupled oscillators
,”
Prog. Theor. Phys.
89
,
313
–
323
(
1993
).
https://doi.org/10.1143/ptp/89.2.313
Google Scholar
Crossref
Search ADS
 
53.
N.
Nakagawa
 and
Y.
Kuramoto
, “
From collective oscillations to collective chaos in a globally coupled oscillator system
,”
Phys. D
75
,
74
–
80
(
1994
).
https://doi.org/10.1016/0167-2789(94)90275-5
Google Scholar
Crossref
Search ADS
 
54.
M.
Bolotov
,
L.
Smirnov
,
G.
Osipov
, and
A.
Pikovsky
, “
Simple and complex chimera states in a nonlinearly coupled oscillatory medium
,”
Chaos
28
,
045101
(
2018
).
https://doi.org/10.1063/1.5011678
Google Scholar
Crossref
Search ADS
PubMed
 
55.
S. A.
Marvel
,
R. E.
Mirollo
, and
S. H.
Strogatz
, “
Identical phase oscillators with global sinusoidal coupling evolve by Möbius group action
,”
Chaos
19
,
043104
(
2009
).
https://doi.org/10.1063/1.3247089
Google Scholar
Crossref
Search ADS
PubMed
 
56.
E. A.
Martens
,
E.
Barreto
,
S. H.
Strogatz
,
E.
Ott
,
P.
So
, and
T. M.
Antonsen
, “
Exact results for the Kuramoto model with a bimodal frequency distribution
,”
Phys. Rev. E
79
,
026204
(
2009
).
https://doi.org/10.1103/PhysRevE.79.026204
Google Scholar
Crossref
Search ADS
 
57.
S.
Watanabe
 and
S. H.
Strogatz
, “
Integrability of a globally coupled oscillator array
,”
Phys. Rev. Lett.
70
,
2391
–
2394
(
1993
).
https://doi.org/10.1103/PhysRevLett.70.2391
Google Scholar
Crossref
Search ADS
PubMed
 
58.
S.
Watanabe
 and
S. H.
Strogatz
, “
Constant of motion for superconducting Josephson arrays
,”
Phys. D
74
,
197
–
253
(
1994
).
https://doi.org/10.1016/0167-2789(94)90196-1
Google Scholar
Crossref
Search ADS
 
59.
A.
Pikovsky
 and
M.
Rosenblum
, “
Dynamics of heterogeneous oscillator ensembles in terms of collective variables
,”
Phys. D
240
,
872
–
881
(
2011
).
https://doi.org/10.1016/j.physd.2011.01.002
Google Scholar
Crossref
Search ADS
 
60.
S.
Lee
 and
K.
Krischer
, “
Attracting poisson chimeras in two-population networks
,”
Chaos
31
,
113101
(
2021
).
https://doi.org/10.1063/5.0065710
Google Scholar
Crossref
Search ADS
PubMed
 
61.
S.
Jalan
 and
A.
Suman
, “Multiple first-order transitions in simplicial complexes on multilayer systems,”
Phys. Rev. E
106
,
044304
(
2022
).
https://doi.org/10.1103/PhysRevE.106.044304
Google Scholar
Crossref
Search ADS
PubMed
 
62.
M.
Desroches
,
A.
Guillamon
,
E.
Ponce
,
R.
Prohens
,
S.
Rodrigues
, and
A.
Teruel
, “
Canards, folded nodes, and mixed-mode oscillations in piecewise-linear slow-fast systems
,”
SIAM Rev.
58
,
653
–
691
(
2016
).
https://doi.org/10.1137/15M1014528
Google Scholar
Crossref
Search ADS
 
63.
H.
Jardón-Kojakhmetov
 and
C.
Kuehn
, “A survey on the blow-up method for fast-slow systems,” in Mexican Mathematicians in the World: Trends and Recent Contributions, edited by F. Galaz-García, C. González-Tokman, and J. C. Pardo Millán (AMS Press, 2021), pp. 115–160.
64.
N.
Fenichel
, “
Geometric singular perturbation theory for ordinary differential equations
,”
J. Differ. Equ.
31
,
53
–
98
(
1979
).
https://doi.org/10.1016/0022-0396(79)90152-9
Google Scholar
Crossref
Search ADS
 
65.
C. K. R.
T. Jones
, “Geometric singular perturbation theory,” in 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.
66.
C.
Kuehn
,
Multiple Time Scale Dynamics
(
Springer
,
2015
).
Google Scholar
Crossref
Search ADS
 
© 2023 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).
Published open access through an agreement with University of Groningen