# Stable Chimera States: A Geometric Singular Perturbation Approach.

Luis Guillermo Venegas-Pineda<EMAIL_ADDRESS>Bernoulli Institute
for Mathematics, Computer Science and Artificial Intelligence, University of
Groningen, Nijenborgh 9, Groningen, The Netherlands. Hildeberto Jardón-
Kojakhmetov  Bernoulli Institute for Mathematics, Computer Science and
Artificial Intelligence, University of Groningen, Nijenborgh 9, Groningen, The
Netherlands. Ming Cao  Engineering and Technology Institute Groningen,
University of Groningen, Nijenborgh 9, Groningen, The Netherlands.

## Abstract

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

##  1 Introduction

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

####

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

####

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

####

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 chimeras
are determined. In particular, we prove that stable chimera states and
persistent breathing chimeras are achieved under a sufficiently slow co-
evolving law with specific parameter arrangements. Finally, other parameter
regimes related to the loss of normal hyperbolicity conditions and leading to
complex behavior on the reduced system are identified and compared numerically
at the network level, generating exciting ideas of future work.

####

The rest of the paper is structured as follows. In section 2 we introduce the
basic notions supporting our results, including an overview of the Ott-
Antosnen ansatz as well as a brief summary on fast-slow dynamics. Later, in
section 3 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. In section 4, we present the main results of our
study for two distinct problems, the co-evolutionary inter- and intra-
coupling, respectively. For each scenario, we determine the geometric
properties of the related mean field and give conditions for the effective
generation of chimera states when considering adaptive regimes only dependent
on macroscopic quantities to preserve the Ott-Antonsen method. Then, utilizing
the attractiveness of the Ott-Antonsen manifold, we provide with arguments for
the chimera state to be reflected in the finite-size network. Thereafter, in
section 5 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 coevolution in
the coupling strengths of the network. Finally, in section 6 we give our
conclusions and discuss future research opportunities as the consequence of
the outcomes presented. In particular, we emphasize the analysis of the non-
hyperbolic case for which relaxation oscillations and canard solutions,
related to breathing chimera states due to the patterns numerically observed,
have been produced in simulations for both the mean field and the network.

##  2 Preliminaries

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

###  2.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 co-evolutionary laws only depend on macroscopic quantities
of the network, which allows us to treat the couplings as fixed parameters for
the discussions in this section. Further details are given in section 3.
Therefore, we begin by studying an all-to-all network composed by different
populations of heterogeneous Kuramoto phase oscillators, which can be
expressed in the form

$\dot{\theta}_{i}^{\sigma}=\omega_{i}^{\sigma}+\sum_{\sigma^{\prime}=1}^{M}\frac{k_{\sigma\sigma^{\prime}}}{N_{\sigma^{\prime}}}\sum_{j=1}^{N_{\sigma^{\prime}}}\sin{\left(\theta_{j}^{\sigma^{\prime}}-\theta_{i}^{\sigma}+\beta_{\sigma\sigma^{\prime}}\right)},$
(1)

where $\theta_{i}^{\sigma}$ and $N_{\sigma}$ denote the state of the $i$-th
oscillator and the number of elements in the $\sigma$-population,
respectively. Additionally, $k_{\sigma\sigma^{\prime}}$ and
$\beta_{\sigma\sigma^{\prime}}$ are the coupling strength and phase-lag
between two layers, with $\sigma=1,2,\dots,M$. Moreover, $\omega_{i}^{\sigma}$
is the natural frequency of each oscillator in its respective community,
selected from $M$ unimodal Cauchy-Lorentz distributions as

$g^{\sigma}(\omega_{\sigma})=\frac{1}{\pi}\left[\frac{\Delta_{\sigma}}{(\omega_{\sigma}-\hat{\omega}_{\sigma})^{2}+\Delta_{\sigma}^{2}}\right],$
(2)

with $\hat{\omega}_{\sigma}$ and $\Delta_{\sigma}$ the central frequency and
width parameter of each distribution. Notice that, 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 $\beta_{\sigma\sigma^{\prime}}=0$.

####

Considering the thermodynamic limit ($N_{\sigma}\rightarrow\infty$), the state
of (1) can be described by the continuous probability distribution functions
$f^{\sigma}\left(\theta^{\sigma},\omega_{\sigma},t\right)$, which measure the
fraction of oscillators with phases in the interval
$[\theta^{\sigma},\theta^{\sigma}+\textnormal{d}\theta^{\sigma}]$ and natural
frequencies between
$[\omega_{\sigma},\omega_{\sigma}+\textnormal{d}\omega_{\sigma}]$ at time $t$.
Moreover, since the number of oscillators is preserved for every population at
all time, each density function satisfies the continuity equation

$\frac{\partial f^{\sigma}}{\partial
t}+\frac{\partial}{\partial\theta^{\sigma}}\left(f^{\sigma}v^{\sigma}\right)=0,$
(3)

where $v^{\sigma}\left(\theta^{\sigma},\omega^{\sigma},t\right)$ is the
continuous version of the angular velocity of the oscillators in each
population, expressed as

$v^{\sigma}\left(\theta^{\sigma},\omega_{\sigma},t\right)=\omega_{\sigma}+\frac{1}{2\imath}\left(z_{\sigma}e^{-\imath\theta^{\sigma}}-\bar{z}_{\sigma}e^{\imath\theta^{\sigma}}\right),$

with the over-bar representing the complex conjugate and $z_{\sigma}$ the
Kuramoto order parameter in the form

$z_{\sigma}(t)=\rho_{\sigma}e^{\imath\phi_{\sigma}}=\sum_{\sigma^{\prime}=1}^{M}k_{\sigma\sigma^{\prime}}\int_{0}^{2\pi}e^{\imath\theta^{\sigma^{\prime}}}f^{\sigma^{\prime}}\left(\theta^{\sigma^{\prime}},\omega_{\sigma^{\prime}},t\right)d\theta^{\sigma^{\prime}}.$

In geometric terms, the complex order parameter describes the centroid of all
the phasors $e^{\imath\theta^{\sigma}}$ and it grows larger as a higher
synchronization level is reached in the network [38].

####

Following Ott and Antonsen [39], we focus our attention on a particular class
of density functions
$f^{\sigma}\left(\theta^{\sigma},\omega_{\sigma},t\right)$ 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 [40]. Thus, by considering
Fourier series of the form

$f^{\sigma}\left(\theta^{\sigma},\omega_{\sigma},t\right)=\frac{g^{\sigma}(\omega)}{2\pi}\left[1+\sum_{n=1}^{\infty}\alpha_{\sigma}^{n}(\omega_{\sigma},t)e^{\imath
n\theta^{\sigma}}+\sum_{n=1}^{\infty}\bar{\alpha}^{n}_{\sigma}(\omega_{\sigma},t)e^{-\imath
n\theta^{\sigma}}\right],$ (4)

for the density function of each population, and by imposing the ansatz (4) to
the continuum equation (3) we obtain the next set of differential complex
equations

$0=\sum_{n=1}^{\infty}n\alpha_{\sigma}^{n-1}e^{\imath
n\theta^{\sigma}}\left(\dot{\alpha}_{\sigma}+\imath\omega_{\sigma}\alpha_{\sigma}+\frac{1}{2}\left[\alpha_{\sigma}^{2}z_{\sigma}e^{-\imath\alpha_{\sigma}}-\bar{z}_{\sigma}e^{\imath\alpha_{\sigma}}\right]\right)+\textnormal{c.c.},$
(5)

where c.c. stands for the complex conjugate of the expression in the right-
hand side of (5). Moreover, since the sum component in (5) is not zero, it
means that the term within the parentheses vanish identically if the ansatz
(4) is a valid solution to the continuity equation (3). Hence,

$\dot{\alpha}_{\sigma}+\imath\omega_{\sigma}\alpha_{\sigma}+\frac{1}{2}\left[\alpha_{\sigma}^{2}z_{\sigma}e^{-\imath\alpha_{\sigma}}-\bar{z}_{\sigma}e^{\imath\alpha_{\sigma}}\right]=0,$
(6)

with
$z_{\sigma}=\int_{-\infty}^{\infty}\alpha_{\sigma}(\omega_{\sigma},t)g^{\sigma}(\omega_{\sigma})$.
Furthermore, consider solutions of (6) with the initial conditions
$\alpha_{\sigma}(\omega_{\sigma},0)$ that satisfy: (i)
$|\alpha_{\sigma}(\omega_{\sigma},t)|\leq 1$; (ii)
$\alpha_{\sigma}(\omega_{\sigma},0)$ is analytically continuable into the
lower half plane $\text{Im}(\omega_{\sigma})<0$; and (iii)
$|\alpha_{\sigma}(\omega_{\sigma},t)|\rightarrow 0$ as
$\text{Im}(\omega_{\sigma})\rightarrow-\infty$. If such conditions are
satisfied for $\alpha_{\sigma}(\omega_{\sigma},0)$, then they will continue to
be satisfied for $\alpha_{\sigma}(\omega_{\sigma},t)$ [39, 41]. It is
important to mention that the Ott-Antonsen ansatz corresponds to a particular
case of the technique developed by Watanabe-Strogatz [42, 43] which considers
all the harmonics of the frequency distribution function (4). 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 [44]. Furthermore,
the Ott-Antonsen manifold is known to be globally attracting when the
distribution of natural frequencies is non-homogeneous, i.e.,
$\Delta_{\sigma}>0$ [45, 46]. Therefore, every long time behavior of the order
parameter can be obtained by analyzing the corresponding mean-field even when
weak heterogeneities are considered [47].

####

By selecting the probability distribution functions as (2), it is possible to
calculate $\bar{z}_{\sigma}$ by the contour integration in the negative half
of the complex plane yielding
$\bar{z}_{\sigma}=\alpha_{\sigma}\left(\omega_{\sigma}-\imath\Delta_{\sigma},t\right)$
[48]. Substituting this expression in (3) produces $M$-coupled complex ODEs
describing the evolution of the complex order parameter as

$\dot{z}_{\sigma}+\left(\Delta_{\sigma}-\imath\omega_{\sigma}\right)z_{\sigma}+\frac{1}{2}\sum_{\sigma^{\prime}=1}^{M}k_{\sigma\sigma^{\prime}}\left(\bar{z}_{\sigma^{\prime}}z_{\sigma}^{2}-z_{\sigma^{\prime}}\right)=0.$

Finally, with $z_{\sigma}=\rho_{\sigma}e^{-\imath\phi_{\sigma}}$, for which
the negative sign is included in the definition of $\phi_{\sigma}$ in order
for the Poisson kernel to converge to $\delta(\theta-\phi_{\sigma})$ and not
to $\delta(\theta+\phi_{\sigma})$ as $\rho_{\sigma}\rightarrow 1$ [40],
produces by orthogonality $2M$ coupled real ODEs in the form

$\displaystyle\begin{split}\dot{\rho}_{\sigma}&=-\Delta_{\sigma}\rho_{\sigma}+\frac{1}{2}\left(1-\rho_{\sigma}^{2}\right)\sum_{\sigma^{\prime}=1}^{M}k_{\sigma\sigma^{\prime}}\rho_{\sigma^{\prime}}\cos{\left(\phi_{\sigma^{\prime}}-\phi_{\sigma}\right)},\\\
\dot{\phi}_{\sigma}&=-\omega_{\sigma}+\frac{1}{2}\frac{\rho_{\sigma}^{2}+1}{\rho_{\sigma}}\sum_{\sigma^{\prime}=1}^{M}k_{\sigma\sigma^{\prime}}\rho_{\sigma^{\prime}}\sin{\left(\phi_{\sigma^{\prime}}-\phi_{\sigma}\right)}.\end{split}$
(7)

In particular, our mechanism for the generation of stable chimera state relies
on geometric properties of (7) with adaptive coupling strengths under a
macroscopic slow adaptation, which induces a time-scale separation between the
dynamics on and off the network. Thus, the next section briefly summarizes the
principal concepts of fast-slow systems employed later on.

###  2.2 Fast-slow dynamics

A fast-slow system is a singularly perturbed ordinary differential equation in
the form

$\displaystyle\begin{split}\varepsilon\dot{x}&=f(x,y,\varepsilon),\\\
\enspace\dot{y}&=g(x,y,\varepsilon),\end{split}$ (8)

where the over-dot represents the derivative with respect to the slow-time
$\tau$, $x\in\mathbb{R}^{m}$ and $y\in\mathbb{R}^{n}$ denote the fast and slow
variables respectively, and $0<\varepsilon\ll 1$ describes the time-scale
separation. Moreover, we consider
$f:\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{m}$
and
$g:\mathbb{R}^{m}\times\mathbb{R}^{n}\times\mathbb{R}\rightarrow\mathbb{R}^{n}$
of class $C^{k}$, for $k$ sufficiently large. Equivalently, by defining the
fast time $t\coloneqq\tau/\varepsilon$, we can rewrite (8) as

$\displaystyle\begin{split}x^{\prime}&=f(x,y,\varepsilon),\\\
y^{\prime}&=\varepsilon g(x,y,\varepsilon).\end{split}$ (9)

where the prime denotes the derivative with respect to the fast time $t$.
Hence, (8) and (9) are known as the slow and fast formulations, and in the
singular limit ($\varepsilon=0$) two different problems arise, namely the
reduced and the layer problem, respectively from (8) and (9). Despite of being
non-equivalent both problems share a close connection as expressed in the
following definition.

######  Definition 2.1 (Critical manifold).

The critical set is defined as

$\mathcal{C}_{0}=\left\\{(x,y)\in\mathbb{R}^{m}\times\mathbb{R}^{n}:f(x,y,0)=0\right\\}.$
(10)

Moreover, if $\mathcal{C}_{0}$ is a submanifold of
$\mathbb{R}^{m}\times\mathbb{R}^{n}$, it is referred to as the _critical
manifold_. Notice that points of (10) are in a direct correspondence to the
equilibria set of the fast flow, generated by (9). Additionally, a relevant
property that $\mathcal{C}_{0}$ can have is normal hyperbolicity, given as
follows.

######  Definition 2.2 (Normal hyperbolicity).

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

In contrast, $\mathcal{S}$ is non-hyperbolic if the corresponding matrix
$\left(\textnormal{D}_{x}f\right)(p,0)$ has at least one eigenvalue with zero
real part for any $p\in\mathcal{S}$. Depending on the normal hyperbolicity
condition of (10) different analysis techniques can be employed. For instance,
non-hyperbolic points, related to dynamic features such as relaxation
oscillations and canards [49], can be studied through the blow-up method [50].
On the other hand, for normally hyperbolic cases we have Fenichel’s Theorem
[51, 52], stated as follows.

######  Theorem 2.1 (Fenichel’s Theorem).

Given a compact normally hyperbolic submanifold (possibly with boundary)
$\mathcal{S}=\mathcal{S}_{0}$ of the critical manifold $\mathcal{C}_{0}$ of
(8), and $f,g\in C^{k<\infty}$, for $0<\varepsilon\ll 1$, then:

  1. 1.

There exists a locally invariant manifold $\mathcal{S}_{\varepsilon}$
diffeomorphic to $\mathcal{S}_{0}$. Local invariance means that trajectories
can enter or escape $\mathcal{S}_{\varepsilon}$ only through its boundaries.

  2. 2.

$\mathcal{S}_{\varepsilon}$ has Hausdorff distance
$\mathrm{d_{H}}\left(V,W\right)=O\left(\varepsilon\right)$ from
$\mathcal{S}_{0}$, as $\varepsilon\rightarrow 0$.

  3. 3.

The flow on $\mathcal{S}_{\varepsilon}$ converges to the flow generated by the
reduced problem, known as the slow flow, as $\varepsilon\rightarrow 0$.

  4. 4.

$\mathcal{S}_{\varepsilon}$ is $C^{k}$-smooth.

  5. 5.

$\mathcal{S}_{\varepsilon}$ is normally hyperbolic and has the same stability
properties with respect to the fast variables as $\mathcal{S}_{0}$.

  6. 6.

$\mathcal{S}_{\varepsilon}$ is usually not unique. In regions with a fixed
distance from $\partial\mathcal{S}_{\varepsilon}$, all manifolds which
satisfies conditions 1-5 above are exponentially close to each other with a
Hausdorff distance of order $O\left(\textnormal{exp}(-C/\varepsilon)\right)$
for some $C>0$, $C\in O(1)$, as $\varepsilon\rightarrow 0$.

Thus, any manifold $\mathcal{S}_{\varepsilon}$ satisfying conditions 1-5 is
called a slow manifold.

######  Remark.

Regularly, a particular manifold which satisfies conditions $1-5$ of theorem
2.1 is referred as _the_ slow manifold. Rigorously, this is not correct since
$\mathcal{S}_{\varepsilon}$ 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 [53].

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 fast-slow normal form theory and the blow-up
method, to name a few.

####

Hence, with the previously summarized results we are now ready to present the
research model studied in our work. In the next section, we introduce the
network analyzed and its corresponding mean-field reduction, obtained through
the Ott-Antonsen ansatz with all its necessary considerations. In particular,
we emphasize that our network is composed by heterogeneous oscillators in
order for the Ott-Antonsen manifold to be globally attracting.

##  3 Co-evolutionary multilayer Kuramoto network

The system we analyze consists of a complete network composed by two
populations of different heterogeneous Kuramoto oscillators intra- and inter-
connected with common coevolutive coupling strengths, as graphically depicted
in figure 1 and mathematically expressed in the following set of coupled
ordinary differential equations

$\dot{\theta}_{i}^{\sigma}=\omega_{i}^{\sigma}+\frac{k_{\sigma}}{N_{\sigma}}\sum_{j=1}^{N_{\sigma}}\sin{\left(\theta_{j}^{\sigma}-\theta_{i}^{\sigma}\right)}+\frac{\mu}{N_{\sigma^{\prime}}}\sum_{l=1}^{N_{\sigma^{\prime}}}\sin{\left(\theta_{l}^{\sigma^{\prime}}-\theta_{i}^{\sigma}\right)},$
(11)

where $\theta_{i}^{\sigma}$ and $\omega_{i}^{\sigma}$ represent the phase and
natural frequency of the $i$-th oscillator, chosen from a unimodal Cauchy-
Lorentz distribution as (2), while $N_{\sigma}$ denotes the number of elements
in the $\sigma$-layer ($\sigma=1,2$). Additionally, the coupling strengths
$\mu$, $k_{1}$, and $k_{2}$ are slow adaptive variables which only depend on
macroscopic quantities of (11), simplifying the reduction method by
considering the couplings as constant parameters in the fast time-scale and
avoiding modifications on the mean-field technique employed [54]. Thus, in
accordance to (7) a mean-field representation of (11) is given by

$\displaystyle\begin{split}\dot{\rho}_{\sigma}&=-\Delta_{\sigma}\rho_{\sigma}+\frac{1}{2}\left(1-\rho_{\sigma}^{2}\right)\left(k_{\sigma}\rho_{\sigma}+\mu\rho_{\sigma^{\prime}}\cos{(\phi_{\sigma^{\prime}}-\phi_{\sigma})}\right),\\\
\dot{\phi}_{\sigma}&=-\omega_{\sigma}+\frac{1}{2}\mu\frac{\rho_{\sigma}^{2}+1}{\rho_{\sigma}}\rho_{\sigma^{\prime}}\sin{(\phi_{\sigma^{\prime}}-\phi_{\sigma})}.\end{split}$
(12)

Furthermore, considering our two-population problem and by introducing the
additional variables $\psi=\phi_{2}-\phi_{1}$, and
$\Omega=\omega_{2}-\omega_{1}$, the mean-field (12) is expressed as

$\displaystyle\begin{split}\dot{\rho}_{1}&=-\Delta_{1}\rho_{1}+\frac{1}{2}\left(1-\rho_{1}^{2}\right)\left(k_{1}\rho_{1}+\mu\rho_{2}\cos{\psi}\right),\\\
\dot{\rho}_{2}&=-\Delta_{2}\rho_{2}+\frac{1}{2}\left(1-\rho_{2}^{2}\right)\left(k_{2}\rho_{2}+\mu\rho_{1}\cos{\psi}\right),\\\
\dot{\psi}&=-\Omega-\frac{1}{2}\mu\left(\frac{\rho_{1}^{2}+\rho_{2}^{2}+2\rho_{1}^{2}\rho_{2}^{2}}{\rho_{1}\rho_{2}}\right)\sin{\psi}.\end{split}$
(13)

Figure 1: Graphic representation of the multilayer coevolutive network (11).
The black and blue components represent elements corresponding to the first
and second populations, respectively. Moreover, $\theta_{i}^{\sigma}$ denotes
the phase of each element in the network, with $i=1,\dots,N_{\sigma}$ and
$\sigma=1,2$. The adaptive intercoupling strength (red) is expressed as
$\mu(t)$ and the intracoupling variable for each population is captured by
$k_{\sigma}$. Although (11) is a complete graph, we present only certain
nodes’ connections for simplicity.

######  Remark.

Observe that, since the selected probability distribution functions are of
unimodal Cauchy-Lorentz type, the real part of the local Kuramoto order
parameter $R_{\sigma}$ of (11) directly corresponds to the variables
$\rho_{\sigma}$ in (12).

In the following section we analyze the behaviour of (13) under the presence
of different types of adaptive coupling strengths, only dependent on
macroscopic quantites of (11). In particular, we benefit from the introduction
of a time-scale separation between the dynamics on the nodes and that related
to the coupling strength, allowing us the use of fast-slow theory without
alterations on the reduction method employed.

##  4 Analysis

To begin with and in order to facilitate computations, we reduce the dimension
of the fast-time mean field (13) 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 4.1 (Invariant synchronization set).

The set
$\left\\{\left(\rho_{1},\rho_{2},\psi,\mu\right)\in\mathbb{K}:\rho_{2}=1\right\\}$
in (13), with $\mathbb{K}=[0,1]^{2}\times[0,2\pi)\times\mathbb{R}$, is
invariant and attracting if and only if $\Delta_{2}=0$, for $k_{2}$ large
enough. Furthermore, there exists an invariant and attracting set
$\left\\{\left(\rho_{1},\rho_{2},\psi,\mu\right)\in\mathbb{K}:\rho_{2}(\Delta_{2})=1-a\Delta_{2}+O\left(\Delta_{2}^{2}\right)\right\\}$
with $a>0$, and $\Delta_{2}$ positive and sufficiently small when $k_{2}>0$ is
adequately large.

###### Proof.

Consider the function

$h_{2}(\rho_{2},\Delta_{2})\coloneqq-\Delta_{2}\rho_{2}+\frac{1}{2}\left(1-\rho_{2}^{2}\right)\left(k_{2}\rho_{2}+\mu\rho_{1}\cos{\psi}\right),$
(14)

associated with the flow equation of $\rho_{2}$ in (13). Thus, the set
$\left\\{\left(\rho_{1},\rho_{2},\psi,\mu\right)\in\mathbb{K}:\rho_{2}=1\right\\}$
is invariant if and only if $\Delta_{2}=0$. Moreover, by taking the derivative
of (14) with respect to $\rho_{2}$ and evaluating at
$(\rho_{2},\Delta_{2})=(1,0)$ yields

$\frac{\partial
h_{2}}{\partial\rho_{2}}\bigg{\rvert}_{(1,0)}=-\left(k_{2}+\mu\rho_{1}\cos{\psi}\right),$
(15)

which is strictly negative in $\mathbb{K}$ for $k_{2}$ adequately large.
Therefore, from the implicit function theorem there exists a set
$\\{\left(\rho_{1},\rho_{2},\psi,\mu\right)\in\mathbb{K}:\rho_{2}=q(\Delta_{2})\\}$
such that $h_{2}\left(q(\Delta_{2}),\Delta_{2}\right)=0$, with $q(0)=1$, $q\in
C^{1}$. In particular, we obtain an expansion of $\rho_{2}$ as a regular
perturbation on $\Delta_{2}$ in the form
$\rho_{2}(\Delta_{2})=1-a\Delta_{2}+O\left(\Delta_{2}^{2}\right)$ with the
invariant flow of (14) given now as

$0=\left(-1+ak_{2}+a\mu\rho_{1}\cos{\psi}\right)\Delta_{2}+O\left(\Delta_{2}^{2}\right),$
(16)

which is satisfied for $a=\left(k_{2}+\mu\rho_{1}\cos{\psi}\right)^{-1}$,
positive whenever $|k_{2}|\geq|\mu|$ along the interval $\rho_{1}\in(0,1)$.
Moreover, the derivative of (14) evaluated at
$\rho_{2}\left(\Delta_{2}\right)=1-a\Delta_{2}+O\left(\Delta_{2}^{2}\right)$,
given by

$\frac{\partial
h_{2}}{\partial\rho_{2}}\bigg{\rvert}_{\rho_{2}=1-a\Delta_{2}+O\left(\Delta_{2}^{2}\right)}=-\Delta_{2}+\frac{1}{2}k_{2}\left(a\Delta_{2}(2-3a\Delta_{2})-2\right)-\mu\rho_{1}(1-a\Delta_{2})\cos{\psi},$
(17)

approaches (15) in the limit when $\Delta_{2}$ tends to zero, rendering the
set
$\\{\left(\rho_{1},\rho_{2},\psi,\mu\right)\in\mathbb{K}:\rho_{2}(\Delta_{2})=1-a\Delta_{2}+O(\Delta_{2}^{2})\\}$
attracting with $\Delta_{2}$ sufficiently small for $k_{2}$ sufficiently
large. ∎

Therefore, from proposition 4.1 and by choosing $k_{2}$ sufficiently large and
$\Delta_{2}$ sufficiently small, in what follows we simply consider
$\rho_{2}=1$. Moreover, we analyze two different problems, namely, the co-
evolutionary intercoupling, i.e., $\mu$ adaptive, and the coevolutive
intracoupling, i.e., $k_{1}$ adaptive. This distinction is made in order to
compare the geometries of both adaptation schemes and their isolated effects
in the network. Notice that stable chimera states can be generated in both
coevolving scenarios with determined geometric differences, as we subsequently
demonstrate. To begin with, in the next subsection we study the adaptive
intercoupling problem in (13) while considering the intracouplings as constant
parameters.

###  4.1 Co-evolutionary Intercoupling

As a first approach, the coevolving intercoupling strength $\mu$ effects are
analyzed. Therefore, by fixing $\rho_{2}=1$ and $k_{1},k_{2}\in\mathbb{R}$,
with $k_{2}$ sufficiently large, as a parameters in both populations of (11)
the mean-field (13) is

$\displaystyle\begin{split}\dot{\rho}_{1}&=-\Delta_{1}\rho_{1}+\frac{1}{2}\left(1-\rho_{1}^{2}\right)\left(k_{1}\rho_{1}+\mu\cos{\psi}\right),\\\
\dot{\psi}&=-\Omega-\frac{1}{2}\mu\left(\frac{3\rho_{1}^{2}+1}{\rho_{1}}\right)\sin{\psi},\\\
\dot{\mu}&=\varepsilon_{\mu}f\left(\rho_{1},\psi,\mu,t\right),\end{split}$
(18)

where $f(\rho_{1},\psi,\mu,t)$ is some adaptive law which only depends on
macroscopic quantities of (11), with $\varepsilon_{\mu}$ as the time-scale
separation parameter. Consequently, for the critical manifold of (18) we
obtain

$\displaystyle\begin{split}\mu\cos{\psi}&=\rho_{1}\left(\dfrac{2\Delta_{1}}{1-\rho_{1}^{2}}-k_{1}\right),\\\
\mu\sin{\psi}&=-2\Omega\left(\frac{\rho_{1}}{3\rho_{1}^{2}+1}\right).\end{split}$

Hence, the critical manifold consists of two branches defined by

$\mathcal{C}_{0}^{\pm}=\left\\{\left(\rho_{1},\psi,\mu\right)\in\mathbb{L}:\mu=\pm\rho_{1}\sqrt{\left(\frac{2\Delta_{1}}{1-\rho_{1}^{2}}-k_{1}\right)^{2}+\left(\frac{2\Omega}{3\rho_{1}^{2}+1}\right)^{2}}\right\\},$
(19)

with $\mathbb{L}=\left[0,1\right)\times\left[0,2\pi\right)\times\mathbb{R}$.
Moreover, the Jacobian of (18) evaluated at (19) is

$J_{f}\Big{\rvert}_{\mathcal{C}_{0}^{\pm}}=\begin{bmatrix}-\Delta_{1}\left(\frac{1+\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)+\frac{1}{2}k_{1}\left(1-\rho_{1}^{2}\right)&\Omega\rho_{1}\left(\frac{1-\rho_{1}^{2}}{1+3\rho_{1}^{2}}\right)\\\
\frac{\Omega}{\rho_{1}}\left(\frac{3\rho_{1}^{2}-1}{3\rho_{1}^{2}+1}\right)&-\frac{1}{2}\left(\frac{1+3\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\left(2\Delta_{1}-k_{1}(1-\rho_{1}^{2})\right)\end{bmatrix}.$
(20)

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

######  Proposition 4.2 (Intercoupling normal hyperbolicity and stability
($|\Omega|\ll 1$)).

The critical manifold (19) of (18), with $\Omega=0$, is normally hyperbolic in
the interval $\rho_{1}\in[0,1)$ if and only if $k_{1}<2\Delta_{1}$. Moreover,
if the critical manifold is normally hyperbolic then it is attracting for
every $\rho_{1}\in[0,1)$. Finally, for $|\Omega|\ll 1$, the normal
hyperbolicity and stability of the $\Omega=0$ case are conserved.

###### Proof.

Notice that, when $\Omega=0$, the Jacobian (20) reduces to

$J_{f}\Big{\rvert}_{\mathcal{C}_{0}^{\pm}}=\begin{bmatrix}-\Delta_{1}\left(\frac{1+\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)+\frac{1}{2}k_{1}\left(1-\rho_{1}^{2}\right)&0\\\
0&-\frac{1}{2}\left(\frac{1+3\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\left(2\Delta_{1}-k_{1}(1-\rho_{1}^{2})\right)\end{bmatrix}.$
(21)

Thus, the eigenvalues of (21) are given as

$\displaystyle\lambda_{1}$
$\displaystyle=-\Delta_{1}\left(\frac{1+\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)+\frac{1}{2}k_{1}\left(1-\rho_{1}^{2}\right),$
$\displaystyle\lambda_{2}$
$\displaystyle=-\frac{1}{2}\left(\frac{1+3\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\left(2\Delta_{1}-k_{1}(1-\rho_{1}^{2})\right),$

which are purely real for all $\rho_{1}\in[0,1)$. Consequently, for normal
hyperbolicity we require that

$\displaystyle\lambda_{1}$ $\displaystyle\neq 0\Rightarrow
1>\frac{\left(1-\rho_{1}^{2}\right)^{2}}{1+\rho_{1}^{2}}\neq\frac{2\Delta_{1}}{k_{1}},$
$\displaystyle\lambda_{2}$ $\displaystyle\neq 0\Rightarrow
1>\left(1-\rho_{1}^{2}\right)\neq\frac{2\Delta_{1}}{k_{1}}.$

Therefore, the critical manifold (19) is normally hyperbolic in the interval
$\rho_{1}\in[0,1)$ if and only if $k_{1}<2\Delta_{1}$. Moreover, for stability
notice that if $k_{1}<2\Delta_{1}$ then $\lambda_{1,2}<0$ for every
$\rho_{1}\in[0,1)$. Finally, if now we consider $|\Omega|$ positive and
suffiently small the normal hyperbolicity and attractiveness of the $\Omega=0$
case are persistent as the eigenvalues $\lambda_{1,2}$ of (20) depend
differentiably on $\Omega$ [55]. ∎

Hence, proposition 4.2 states that as long as the natural frequencies of the
first population are sufficiently sparse, every point in the critical manifold
(19) is normally hyperbolic and attracting in the interval $\rho_{1}\in(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 (13), under
proposition 4.1, for $k_{1}$ adaptive.

######  Remark.

It is important to highlight that the identification of two branches in (19),
namely $\mathcal{C}_{0}^{+}$ and $\mathcal{C}_{0}^{-}$, is done for
convenience of notation. Nevertheless, by definition 2.1, every point in these
branches conforms the critical manifold, regardless its location in a
particular segment, as shown in figures 2 and 4 for the co-evolutionary inter-
and intra-coupling, respectively.

###  4.2 Co-evolutionary Intracoupling

As a comparison to the adaptive intercoupling case of the previous section, we
now analyze the effect of one coevolutive intracoupling. By letting $k_{1}$
adaptive and with the result in proposition 4.1, now (13) reduces to

$\displaystyle\begin{split}\dot{\rho}_{1}&=-\Delta_{1}\rho_{1}+\frac{1}{2}\left(1-\rho_{1}^{2}\right)\left(k_{1}\rho_{1}+\mu\cos{\psi}\right),\\\
\dot{\psi}&=-\Omega-\frac{1}{2}\mu\left(\frac{3\rho_{1}^{2}+1}{\rho_{1}}\right)\sin{\psi},\\\
\dot{k}_{1}&=\varepsilon_{1}f\left(\rho_{1},\psi,k_{1},t\right),\end{split}$
(22)

with the same variable interpretations as before and $k_{1}$ as the
coevolutive coupling strength. Thus, the critical manifold of (22) is

$\mathcal{C}_{0}^{\pm}=\left\\{\left(\rho_{1},\psi,k_{1}\right)\in\mathbb{L}:k_{1}=\frac{2\Delta_{1}}{1-\rho_{1}^{2}}\pm\sqrt{\left(\frac{\mu}{\rho_{1}}\right)^{2}-\left(\frac{2\Omega}{3\rho_{1}^{2}+1}\right)^{2}}\right\\}.$
(23)

######  Remark.

Notice that (23) can have complex elements, producing disconnected branches of
$\mathcal{C}_{0}^{\pm}$. In fact, by definition 2.1 the critical manifold is a
geometric object of real components. Therefore, we restrict our attention to
the case for which the discriminant in (23) is non-negative.

######  Proposition 4.3 (Connected critical manifold).

The critical set (23) of (22) is completely connected if and only if
$|\mu|\geq|\Omega|/\sqrt{3}$.

###### Proof.

In order to produce only real components we require that

$\left(\frac{\mu}{\rho_{1}}\right)^{2}-\left(\frac{2\Omega}{3\rho_{1}^{2}+1}\right)^{2}\geq
0,$

which can be rewritten as

$\frac{|\mu|}{2|\Omega|}\geq\frac{\rho_{1}}{3\rho_{1}^{2}+1}.$ (24)

By calculating the derivative of the right-hand side of (24) with respect to
$\rho_{1}$, its maximum occurs at $\rho_{1}=1/\sqrt{3}$, for which the
function reaches a value of $1/2\sqrt{3}$. Hence, the critical manifold (23)
is completely connected if and only if
$\left|\mu\right|\geq|\Omega|/\sqrt{3}$. ∎

######  Corollary 4.1.

When $\Omega=0$, the critical manifold (23) is completely connected for every
$\mu\in\mathbb{R}$.

Next, we obtain the Jacobian of (22) with respect to the fast variables, which
by evaluation at (23) yields

$J_{f}\Big{\rvert}_{\mathcal{C}_{0}^{\pm}}=\begin{bmatrix}-2\Delta_{1}\left(\frac{\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\pm\frac{1}{2}\left(1-\rho_{1}^{2}\right)h\left(\rho_{1}\right)&\Omega\rho_{1}\left(\frac{1-\rho_{1}^{2}}{3\rho_{1}^{2}+1}\right)\\\
\frac{\Omega}{\rho_{1}}\left(\frac{3\rho_{1}^{2}-1}{3\rho_{1}^{2}+1}\right)&\pm\frac{1}{2}(3\rho_{1}^{2}+1)h\left(\rho_{1}\right)\end{bmatrix},$
(25)

with
$h(\rho_{1})\coloneqq\sqrt{\left(\frac{\mu}{\rho_{1}}\right)^{2}-\left(\frac{2\Omega}{3\rho_{1}^{2}+1}\right)^{2}}$,
non-negative under proposition 4.3. Thus, the following result summarize the
conditions for the normal hyperbolicity and stability of (23).

######  Proposition 4.4 (Intracoupling normal hyperbolicity and stability
$(|\Omega|\ll 1)$).

The critical manifold (23) of (22) has one normally hyperbolic and attracting
branch $\left(\mathcal{C}_{0}^{-}\right)$, and one unstable branch
$\left(\mathcal{C}_{0}^{+}\right)$ along the interval $\rho_{1}\in(0,1)$ for
every $\mu\neq 0$, with $\Omega=0$. Moreover, for $|\Omega|\ll 1$, the normal
hyperbolicity and stability conditions of the $\Omega=0$ case are preserved.

###### Proof.

Observe that for $\Omega=0$ the Jacobian (25) reduces to

$J_{f}\Big{\rvert}_{\mathcal{C}_{0}^{\pm}}=\begin{bmatrix}-2\Delta_{1}\left(\frac{\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\pm\frac{1}{2}\left(\frac{1-\rho_{1}^{2}}{\rho_{1}}\right)|\mu|&0\\\
0&\pm\frac{1}{2}\left(\frac{3\rho_{1}^{2}+1}{\rho_{1}}\right)|\mu|\end{bmatrix},$

with eigenvalues in the form

$\displaystyle\lambda_{1}^{\pm}$
$\displaystyle=-2\Delta_{1}\left(\frac{\rho_{1}^{2}}{1-\rho_{1}^{2}}\right)\pm\frac{1}{2}\left(\frac{1-\rho_{1}^{2}}{\rho_{1}}\right)\left|\mu\right|,$
$\displaystyle\lambda_{2}^{\pm}$
$\displaystyle=\pm\frac{1}{2}\left(\frac{3\rho_{1}^{2}+1}{\rho_{1}}\right)\left|\mu\right|.$

Hence, the branch of (23) corresponding to $\lambda_{1,2}^{-}$, namely
$\mathcal{C}_{0}^{-}$, is normally hyperbolic and attracting along the
interval $\rho_{1}\in(0,1)$ for every $\mu\neq 0$. On the other hand, the
branch $\mathcal{C}_{0}^{+}$ is either of saddle-type or repelling. Finally,
for $|\Omega|\ll 1$ the normal hyperbolicity and stability conditions of the
$\Omega=0$ case persist since the eigenvalues $\lambda_{1,2}^{\pm}$ of (25)
depend differentiably on $\Omega$ [55]. ∎

Consequently, since we have already provided with the conditions for normal
hyperbolicity and attractiveness of the critical manifolds (19) and (23) 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 (18) and
(22).

######  Theorem 4.1 (Stable chimera states on the mean field).

Given an adequate macroscopic adaptive law in the coevolving intercoupling
(18) or intracoupling (22) 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
$\mathcal{C}_{0}^{\pm}$ are satisfied, for $\varepsilon$ positive and
sufficiently small.

###### Proof.

Once we have determined the normally hyperbolic and attracting conditions of
the critical manifold $\mathcal{C}_{0}^{\pm}$, and recalling that this
geometric object is composed by exponentially stable equilibria of the fast
problem, we can construct a suitable slow macroscopic adaptive law with an
exponentially stable equilibrium point at
$\left(\rho_{1},\psi,K\right)=\left(\rho_{1}^{*},\psi^{*},K^{*}\right)$, such
that if this point is part of $\mathcal{C}_{0}^{\pm}$, then it is
exponentially stable in the perturbed system [56], with $K$ denoting the type
of coevolutive coupling considered, either inter or intracoupling. Moreover,
by Fenichel’s theory there exists a slow manifold for which the dynamics
converge to the slow flow as $\varepsilon\rightarrow 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
$\left(\rho_{1}^{*},\psi^{*},K^{*}\right)$ which is an exponentially stable
equilibrium of the perturbed system (13) for sufficiently small $\varepsilon$.
Finally, recalling that we consider a near synchronization second population,
i.e., $\rho_{2}\rightarrow 1$, it is possible to set the equilibrium point at
any position along the critical manifold with an adequate adaptive law. Hence,
by persistence any synchronization behavior, including stable chimera states,
can be produced in the mean-field of the weakly perturbed system. ∎

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

######  Theorem 4.2 (Stable chimera states on a network).

Under normally hyperbolic conditions of the critical manifold
$\mathcal{C}_{0}^{\pm}$, both for the adaptive intercoupling (18) or
intracoupling (22) systems, dynamical behaviours observed in the mean-field
reduction, such as stable chimera states and sustained breathing chimera
states, are preserved in the network when considering heterogeneous
populations with a slowly coevolving and macroscopic coupling strength, for
$\varepsilon$ 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. $\Delta_{\sigma}>0$
$(\sigma=1,2,\dots,M)$, the attractors for the order parameter dynamics
obtained in the reduced Ott-Antonsen manifold (12) are the only attractors of
the network (11) [45, 46]. Moreover, theorem 4.1 demonstrates the preservation
of equilibrium points on the perturbed system for $\varepsilon$ small enough.
Therefore, dynamical behaviours observed in the mean-field reduction,
including stable chimera states and breathing chimera states, are present in
the network when considering non-homogeneous populations. ∎

In the following section, we show the successful generation of stable chimera
states by following our results, with numerical simulations for both the
network (11) and the mean-field (13) when considering the inter and
intracoupling adaptive scenarios discussed in this section.

##  5 Numerical results

In this section we present a series of simulations validating the
effectiveness of our results as well as some intriguing patterns observed
under certain non-hyperbolic conditions. It is important to emphasize though
that the Ott-Antonsen ansatz consider the thermodynamic limit, i.e.
$N_{1,2}\rightarrow\infty$, 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 [57], as this least-
squares smoothing method reduces noise while preserving the shape and height
of waveform peaks [58], 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 [59].

###  5.1 Co-evolutionary Intercoupling

First, the behavior of both the network (11) and its mean-field counterpart
(18) are compared when different intercoupling coevolving rules are employed.
In figure 2, we present projections to the $\left(R_{1},\mu\right)$ and
$\left(\rho_{1},\mu\right)$-planes of (11) and (18) with
$\dot{\mu}=\varepsilon_{\mu}(-\mu+\gamma_{\mu}-\eta_{\mu}\rho_{1})$. Notice
that $k_{1}<2\Delta_{1}$ and thus the critical manifold (grey) is normally
hyperbolic and attracting for every $\rho_{1}\in(0,1)$. Moreover, observe that
the behavior of the order parameters $R_{1}$ and $\rho_{1}$ (red) converge to
the intersection of the critical manifold and the equilibria set of the slow
adaptation (magenta), which corresponds to a highly incoherent region of the
first layer and with the assumption that the second population is near
synchrony, i.e. $\rho_{2}\approx 1$ and $R_{2}\approx 1$, we have effectively
produced an exponentially stable chimera state, demonstrating the generation
of such pattern in the network (11) by geometric results derived from its
mean-field reduction (18), as stated by theorem 4.2.

####

Similarly, in figure 3 we present a persistent breathing chimera state
produced through the same mechanism as before in both (11) and (18), with
$\dot{\mu}=\varepsilon_{\mu}\cos{(0.02t)}$. Observe that, since the fast
structure of the problem has not been modified, the critical manifold, as well
as its normal hyperbolicity and stability conditions, remain the same and
notice that the order parameters $R_{1}$ and $\rho_{1}$ oscillate between high
and low synchrony regimes, which correspond to the so called breathing chimera
state, characterized by sustained periodical variations on the synchronization
level [60]. 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.

Figure 2: Locally exponentially stable intercoupling chimera. Upper row:
network (left) and mean-field (right) projections to the
$\left(R_{1},\mu\right)$ and $\left(\rho_{1},\mu\right)$ planes for (11) and
(18), respectively. Attracting critical manifold (grey), 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
$\dot{\mu}=\varepsilon_{\mu}(-\mu+\gamma_{\mu}-\eta_{\mu}\rho_{1})$ and
initial conditions $(\rho_{1},\psi,\mu)=(0.99,-0.5,3.0)$. Parameters:
$\Delta_{1}=1.0$, $\Delta_{2}=0.1$, $a=0.2$, $\omega_{1}=101/20$,
$\omega_{2}=\omega_{1}+0.01$, $k_{1}=0.9$, $k_{2}=9.0$, $\gamma_{\mu}=2.5$,
$\eta_{\mu}=10.0$, $\varepsilon_{\mu}=0.02$, and $N_{1,2}=1000$.

Figure 3: Locally exponentially stable intercoupling breathing chimera. Upper
row: network (left) and mean-field (right) projections to the
$\left(R_{1},\mu\right)$ and $\left(\rho_{1},\mu\right)$ planes for (11) and
(18), respectively. Attracting critical manifold (grey), 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 $\dot{\mu}=\varepsilon_{\mu}\cos{(0.02t)}$ and initial conditions
$(\rho_{1},\psi,\mu)=(0.99,-0.5,1.1)$. Parameters: $\Delta_{1}=1.0$,
$\Delta_{2}=0.1$, $a=0.2$, $\omega_{1}=101/20$, $\omega_{2}=\omega_{1}+0.01$,
$k_{1}=0.9$, $k_{2}=9.0$, $\varepsilon_{\mu}=0.02$, and $N_{1,2}=1000$.

###  5.2 Co-evolutionary Intracoupling

Analogously to the intercoupling problem, we examine the numeric results
obtained by considering system (11) with an adaptive intracoupling strength
and its associated mean-field reduction (22). To begin with, a linear feedback
rule $\dot{k}_{1}=\varepsilon_{1}(-k_{1}+\gamma_{1}-\eta_{1}\rho_{1})$ is
considered and in figure 4 we present the obtained results. As explained
earlier, the critical manifold (23) of (22) has one normally hyperbolic and
attracting branch, namely $\mathcal{C}_{0}^{-}$, while the other branch is at
least of saddle type. In figure 4, the stable, saddle and unstable segments of
the critical manifold (23) are represented by the continuous gray, dashed
green and dotted blue lines, respectively. Observe that we have again
initialized the first population of oscillators near complete synchronization
and after a transient phase, it is clear that a chimera state is produced as
the order parameters $\rho_{1}$ and $R_{1}$ converge to the intersection
between the critical manifold and the slow flow set of equilibria. Therefore,
as long as an appropriate macroscopic and slow adaptive law is designed it is
possible to produce diverse patterns due to the geometric properties of the
mean-field (22), and by means of theorem 4.2, such long term dynamical
features are reproduced in the network (11).

Figure 4: Locally exponentially stable intracoupling chimera. Upper row:
network (left) and mean-field (right) projections to the
$\left(R_{1},k_{1}\right)$ and $\left(\rho_{1},k_{1}\right)$ planes for (11)
and (22), respectively. Critical manifold regions: attracting (grey), saddle
(dashed green), repelling (blue dotted), slow adaptation nullcline (magenta)
and response (red). Lower row: order parameter evolution in the network (left)
and mean-field (right). For the network, second (blue) and first (orange)
population synchronization level and the filtered version of the latter (red).
Adaptive law $\dot{k}_{1}=\varepsilon_{1}(-k_{1}+\gamma_{1}-\eta_{1}\rho_{1})$
and initial conditions $(\rho_{1},\psi,k_{1})=(0.99,-0.5,3.5)$. Parameters:
$\Delta_{1}=1.0$, $\Delta_{2}=0.1$, $a=0.2$, $\omega_{1}=101/20$,
$\omega_{2}=\omega_{1}+0.01$, $\mu=0.3$, $k_{2}=9.0$, $\gamma_{1}=2.5$,
$\eta_{1}=10.0$, $\varepsilon_{1}=0.02$, and $N_{1,2}=1000$.

##  6 Discussion and conclusions

The term chimera state gather a broad collection of highly interesting
unexpected synchronization patterns. In this work, we have shown a new
mechanism for the generation of stable chimera states in a complete multilayer
network of heterogeneous Kuramoto oscillators through the inclusion of slow
coevolutive coupling strenghts. Particularly, we use the Ott-Antonsen ansatz
to derive a mean-field representation of the general network, allowing the use
of GSPT results to obtain a deeper insight of the dynamics occurring in the
system. At first, we have analyzed a two-population network with only one
coevolving intercoupling strength, and determined the necessary and sufficient
conditions for the critical manifold to be normally hyperbolic and attracting.
For such a scenario, we have effectively generated stable chimera states as
well as persistent breathing chimera states only by adequately modifying the
macroscopic adaptive law dictating the slow coevolving dynamics. Similarly, we
studied the case of the same network with one coevolutive intracoupling
strength and obtained relevant geometric results that enable the production of
stable chimera states, as well as different synchronization patterns, by
appropriately modifying the adaptation law employed. Therefore, the presented
numeric results validate our mechanism for the generation of diverse patterns
in a complete network of Kuramoto phase oscillators.

####

Considering the results we have obtained, it is natural to wonder about which
kind of behaviours are observed in the non-normally hyperbolic scenario for
both types of adaptive couplings. Although rigorous mathematical analysis is
currently under development, numerical results have shown remarkable patterns.
For instance, in figure 5 a canard trajectory in the intercoupling case is
produced for both the network (11) and the mean-field (18) by introducing a
periodic adaptation law given as $\dot{\mu}=-\varepsilon_{\mu}\cos{\psi}$,
with the periodic dynamics of the angle $\psi$ as depicted in figure 6. The
deformation of the critical manifold occurs as the coupling strength $k_{1}$
increases beyond the value of $2\Delta_{1}$, giving rise to the required
bifurcations for the existence of these particular solutions. However, canards
are highly sensitive to initial conditions and parameter variations, which is
why observing the precise same behavior in the network is a notably complex
task. 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. Nevertheless, a closer correspondence between the network (11)
and the mean-field (18) dynamics is achieved by increasing the population size
and with an adequate level of heterogeneity in the network, either by means of
the natural frequency distributions or a phase-lag between the oscillators.
Finally, it is noteworthy to mention the oscillatory pattern present in the
synchronization level, leading to the question if there exists indeed a
connection between canard solutions and breathing chimera states. The analysis
of such dynamics and their possible applications are open problems that are
being treated, as well as other ideas such as the consideration of multiple
simultaneous adaptive couplings and different mechanisms for the generation of
more intriguing synchronization patterns.

Figure 5: Conjectured canard cycle. Upper row: network (left) and mean-field
(right) projections to the $\left(R_{1},\mu\right)$ and
$\left(\rho_{1},\mu\right)$ planes for (11) and (18), respectively. Critical
manifold regions: attracting (grey), saddle (dashed green), repelling (blue
dotted), and response (red). Lower row: order parameter evolution in the
network (left) and mean-field (right). For the network, first (orange)
population synchronization level and its filtered version (red). Adaptive law
$\dot{\mu}=-\varepsilon_{\mu}\cos{\psi}$ and initial conditions
$(\rho_{1},\psi,\mu)=(0.99,0.4680,0.5)$. Parameters: $\Delta_{1}=0.1$,
$\Delta_{2}=0.1$, $a=0.2$, $\omega_{1}=101/20$, $\omega_{2}=\omega_{1}$,
$k_{1}=5.0$, $k_{2}=5.0$, $\varepsilon_{\mu}=0.02$, and $N_{1,2}=5000$.

Figure 6: Dynamics of the angle $\psi$ for the conjectured canard cycle.
Network (left) and mean-field (right) variable $\psi(t)$ for (11) and (18),
respectively. Adaptive law $\dot{\mu}=-\varepsilon_{\mu}\cos{\psi}$ and
initial conditions $(\rho_{1},\psi,\mu)=(0.9,0.4680,0.5)$. Parameters:
$\Delta_{1}=0.1$, $\Delta_{2}=0.1$, $a=0.2$, $\omega_{1}=101/20$,
$\omega_{2}=\omega_{1}$, $k_{1}=5.0$, $k_{2}=5.0$, $\varepsilon_{1}=0.02$, and
$N_{1,2}=5000$.

## References

## References

  * [1] A. Pikovsky, M. Rosemblum and J. Kurths  “Synchronization: A universal concept in nonlinear sciences”  Cambridge University Press, 2001
  * [2] S. Boccaletti, J. Kurths, G. Osipov, D. L. Valladares and C. S. Zhou  “The synchronization of chaotic systems”  In _Physics Reports_ 366, 2002, pp. 1–101
  * [3] P. J. Uhlhaas and W. Singer  “Neural Synchrony in Brain Disorders: Relevance for Cognitive Dysfunctions and Pathophysiology”  In _Neuron_ 52, 2006, pp. 155–168
  * [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”  In _Communications in Nonlinear Science and Numerical Simulation_ 84, 2020
  * [5] S. Yamaguchi, H. Isejima, T. Matsuo, R. Okura, K. Yagita, M. Kobayashi and H. Okamura  “Synchronization of Cellular Clocks in the Suprachiasmatic Nucleus”  In _Science_ 302, 2003
  * [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”  In _Bioessays_ 37, 2015
  * [7] M. Rohden, A. Sorge, M. Timme and D. Witthaut  “Self-Organized Synchronization in Decentralized Power Grids”  In _Physical Review Letters_ 109, 2012
  * [8] A. Sajadi, R. W. Kenyon and B. Hodge  “Synchronization in electric power networks with inherent heterogeneity up to 100% inverter-based renewable generation”  In _Nature Communications_ 13, 2022
  * [9] J. Buck and E. Buck  “Mechanisms of Rhythmic Synchronous Flashing of Fireflies”  In _Science_ 159.3821, 1968, pp. 1319–1327
  * [10] A. Moiseff and J. Copeland  “Firefly Synchony: A Behavioral Strategy to Minimize Visual Clutter”  In _Science_ 329, 2010, pp. 181
  * [11] K. Wiesenfield, P. Colet and S. H. Strogatz  “Frequency locking in Josephson arrays: Connection with the Kuramoto model”  In _Physical Review E_ 57.2, 1998
  * [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”  In _Science_ 323.5914, 2009
  * [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”  In _Science Advances_ 6, 2020
  * [14] J. Pantaleone  “Stability of incoherence in an isotropic gas of oscillating neutrinos”  In _Physical Review D_ 58, 1998
  * [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”  In _The Journal of Mathematical Neuroscience_ 12, 2020
  * [16] M. J. B. Hauser  “Synchronization of glycolytic activity in yeast cells”  In _Current Genetics_ 68, 2022, pp. 69–81
  * [17] A. T. Winfree  “Biological rhythms and the behavior of populations of coupled oscillators”  In _Journal of Theoretical Biology_ 16.1, 1967, pp. 15–42
  * [18] Y. Kuramoto  “Chemical Oscillations, Waves, and Turbulence”  Springer, Berlin, 1984
  * [19] Y. Kuramoto and D. Battogtokh  “Coexistence of coherence and incoherence in nonlocally coupled phase oscillators”  In _Nonlinear Phenomena in Complex Systems_ 5.380, 2002
  * [20] D. M. Abrams and S. H. Strogatz  “Chimera states for coupled oscillators”  In _Physical Review Letters_ 93, 2004
  * [21] M. J. Panaggio and D. M. Abrams  “Chimera states: coexistence of coherence and incoherence in networks of coupled oscillators”  In _Nonlinearity_ 28, 2015
  * [22] E. Schöll  “Synchronization patterns and chimera states in complex networks: Interplay of topology and dynamics”  In _The European Physical Journal Special Topics_ 225, 2016, pp. 891–919
  * [23] B. K. Bera, S. Majhi, D. Ghosh and M. Perc  “Chimera states: effects of different coupling topologies”  In _Europhysics Letters_ 118, 2017
  * [24] N. C. Rattenborg, C. J. Amlaner and S. L. Lima  “Behavioral, neurophysiological and evolutionary perspectives on unihemispheric sleep”  In _Neuroscience and Biobehavioral Reviews_ 24, 2000
  * [25] J. M. Davidenko, A. V. Pertsov, R. Salomonsz, W. Baxter and J. Jalife  “Stationary and drifting spiral waves of excitation in isolated cardiac muscle”  In _Nature_ 355, 1992
  * [26] A. V. Panfilov  “Spiral breakup as a model of ventricular fibrillation”  In _Chaos_ 8, 1998, pp. 57–64
  * [27] J. C. González-Avella and M. G. Cosenza and M. San Miguel  “Localized coherence in two interacting populations of special agents”  In _Physica A_ 399, 2014, pp. 24–30
  * [28] H. Sakaguchi  “Instability of synchronized motion in nonlocally coupled neural oscillators”  In _Physical Review E_ 73, 2006
  * [29] S. Olmi, A. Politi and A. Torcini  “Collective chaos in pulse-coupled neural networks”  In _Europhysics Letters_ 92.6, 2011
  * [30] 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”  In _Physical Review Letters_ 110, 2013
  * [31] E. A. Martens and K. Klemm  “Transitions from trees to cycles in adaptive flow networks”  In _Frontiers in Physics_ 5, 2017, pp. 1–10
  * [32] 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)_ , 2020, pp. 1110–1115
  * [33] L. Zino, M. Ye and M. Cao  “A two-layer model for coevolving opinion dynamics and collective decision-making in complex social systems”  In _Chaos_ 30, 2020
  * [34] R. Berner, S. Yanchuk and E. Schöll  “What adaptive neuronal networks teach us about power grids”  In _Physical Review E_ 103, 2021
  * [35] W. Gernster and R. Kempter and J. L. van Hemmen and H.ΩWagner  “A neuronal learning rule for sub-millisecond temporal coding”  In _Nature_ 383, 1996
  * [36] W. Gernster and W. M. Kistler  “Mathematical formulations of Hebbian learning”  In _Biological Cybernetics_ 87, 2002, pp. 404–415
  * [37] A. Goriely, E. Kuhl and C. Bick  “Neuronal oscillations on evolving networks: dynamics, damage, degradation, decline, dementia, and death”  In _Physical Review Letters_ 125, 2020
  * [38] P. So and E. Barreto  “Generating macroscopic chaos in a network of globally coupled phase oscillators”  In _Chaos_ 21, 2011
  * [39] E. Ott and T. M. Antonsen  “Low domensional behavior of large systems of globally coupled oscillators”  In _Chaos_ 18, 2008
  * [40] D. M. Abrams, R. Mirollo, S. Strogatz and D. A. Wiley  “Solvable Model for Chimera States of Coupled Oscillators”  In _Physical Review Letters_ 101, 2008
  * [41] 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”  In _Physical Review E_ 79, 2009
  * [42] S. Watanabe and S. H. Strogatz  “Integrability of a globally coupled oscillator array”  In _Phys. Rev. Lett._ , 1993, pp. 2391–2394
  * [43] S. Watanabe and S. H. Strogatz  “Constant of motion for superconducting Josephson arrays”  In _Physica D_ , 1994, pp. 197–253
  * [44] A. Pikovsky and M. Rosenblum  “Dynamics of heterogeneous oscillator ensembles in terms of collective variables”  In _Physica D_ , 2011, pp. 872–881
  * [45] E. Ott and T. M. Antonsen  “Long time evolution of phase oscillator systems”  In _Chaos_ 19, 2009
  * [46] E. Ott, B. R. Hunt and T. M. Antonsen  “Comment on ”Long time evolution of phase oscillator systems””  In _Chaos_ 21, 2011
  * [47] S. Lee and K. Krischer  “Attracting Poisson chimeras in two-population networks”  In _Chaos_ , 2021
  * [48] S. Jalan and A. Suman  “Multiple first-order transitions in simplicial complexes on multilayer networks”  In _Accepted paper_ , 2022
  * [49] M. Desroches, A. Guillamon, E. Ponce, R. Prohens, S. Rodrigues and A.E. Teruel  “Canards, Folded Nodes, and Mixed-Mode Oscillations in Piecewise-Linear Slow-Fast Systems”  In _SIAM Review_ 58.4, 2016, pp. 653–691
  * [50] H. Jardón-Kojakhmetov and C. Kuehn  “A survey on the blow-up method for fast-slow systems”  American Mathematical Society, 2021, pp. 115–160
  * [51] N. Fenichel  “Geometric Singular Perturbation Theory for Ordinary Differential Equations”  In _Journal of Differential Equations_ 31, 1979, pp. 53–98
  * [52] 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_ Berlin, Heidelberg: Springer Berlin Heidelberg, 1995, pp. 44–118
  * [53] C. Kuehn  “Multiple Time Scale Dynamics”  Springer, 2015
  * [54] M. Ciszak, F. Marino, A. Torcini and S. Olmi  “Emergent excitability in populations of nonexcitable units”  In _Phys. Rev. E._ 102, 2020
  * [55] P. D. Lax  “Linear Algebra and its Applications”  Pure and Applied Mathematics: A Wiley-Interscience Series of Texts, Monographs and Tracts, 2007
  * [56] P. Kokotović, H.K. Khalil and J. O’Reilly  “Singular Perturbation Methods in Control: Analysis and Design”  Society for Industrial and Applied Mathematics, 1999
  * [57] A. Savitzky and M. J. E. Golay  “Smoothing and Differentiation of Data by Simplified Least Squares Procedures”  In _Analytical Chemistry_ , 1964, pp. 1627–1639
  * [58] R. W. Schafer  “What Is a Savitzky-Golay Filter? [Lecture Notes]”  In _IEEE Signal Processing Magazine_ 28.4, 2011, pp. 111–117
  * [59] 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”  In _International Journal of Computer Mathematics_ 99.1, 2022, pp. 31–57
  * [60] O. E. Omel’chenko  “Mathematical Framework for Breathing Chimera States”  In _Journal of Nonlinear Science_ , 2022