# A new class of chimeras in locally coupled oscillators with small-amplitude,
high-frequency asynchrony and large-amplitude, low-frequency synchrony

Tasso J. Kaper<EMAIL_ADDRESS>Department of Mathematics and Statistics,
Boston University, Boston, MA 02215, USA  Theodore Vo<EMAIL_ADDRESS>School of Mathematics, Monash University, Clayton, Victoria 3800, Australia

###### Abstract

Chimeras are surprising yet important states in which domains of decoherent
(asynchronous) and coherent (synchronous) oscillations co-exist. In this
article, we report on the discovery of a new class of chimeras, called mixed-
amplitude chimera states, in which the structures, amplitudes, and frequencies
of the oscillations differ substantially in the decoherent and coherent
regions. These mixed-amplitude chimeras exhibit domains of decoherent small-
amplitude oscillations (phase waves) coexisting with domains of stable and
coherent large-amplitude or mixed-mode oscillations. They are observed in a
prototypical bistable partial differential equation with spatially homogeneous
kinetics and purely local, isotropic diffusion. New bifurcations are
identified in which the mixed-amplitude chimeras emerge from, or are
annihilated in, common large-amplitude solutions. Also, key singularities,
folded nodes and folded saddles, arising commonly in multi-scale, bistable
systems play important roles, and these have not previously been studied in
systems with chimeras. The discovery of these mixed-amplitude chimeras is an
important advance for understanding some processes in neuroscience, pattern
formation, and physics which involve both small-amplitude and large-amplitude
oscillations. It may also be of use for understanding some aspects of EEG
recordings from animals that exhibit unihemispheric slow-wave sleep.

> A striking phenomenon in pattern formation occurs when groups of identical
> oscillators exhibit domains of decoherent oscillations coexisting with
> complementary domains of coherent dynamics. These states are known as
> chimera states, named after animals from Greek mythology, which have body
> parts from different known animals. However, unlike the mythological beasts,
> chimera states arise ubiquitously in many areas of physics, chemistry,
> neuroscience, and engineering. This article reports on the discovery of a
> completely new class of chimeras. The structures, amplitudes, and
> frequencies of the oscillations differ substantially in the decoherent and
> coherent domains. Underlying these chimeras are folded singularities, which
> are well-known to be generating mechanisms for complex oscillatory dynamics.
> The discovery of this new class of chimeras was made on a proto-typical
> bistable partial differential equation. These new chimeras have the
> potential to make useful applications, including to unihemispheric slow wave
> sleep, during which the sleeping cerebral hemisphere exhibits synchronous
> high-amplitude, low-frequency rhythms characteristic of sleep, while the
> other hemisphere displays the decoherent low-amplitude, high-frequency
> rhythms characteristic of wakefulness.

##  I Introduction

Discovered in rings of identical phase oscillators KB2002 ; AS2004 , chimera
states are fascinating patterns which exhibit co-existing domains of
decoherent (asynchronous) and coherent (synchronous) oscillations. They occur
broadly in models and experiments in physics, neuroscience, chemistry, and
biology, including in systems of coupled oscillators, integro-differential
equations, and partial differential equations OA2008 ; BPR2010 ; L2010 ; L2015
; OMHS2011 ; TNS2012 ; HMRHOS2012 ; MTFH2013 ; NTS2013 ; SZHK2014 ; YPR2014 ;
ZKS2014 ; OPHSH2015 ; PA2015 ; CCFGR2016 ; KHSKK2016 ; LD2016 ; NRM2017 ;
O2018 ; TRTSE2018 ; KMMG2018 ; ABE2019 ; SA2021 ; Zhang2021 ; SS2014 ; KHK2018
; OMT2008 . Coupling of the identical systems may be non-local, global, local,
or via delay or noise. Reviews and classifications are given in Refs. PA2015,
; KHSKK2016, ; O2018, .

In this article, we report on a new class of chimeras –mixed-amplitude
chimeras– which exhibit regions of decoherent small-amplitude oscillations
(SAOs) and complementary regions of coherent large-amplitude relaxation
oscillations (LAOs) or coherent mixed-mode oscillations (MMOs). (The latter
are combinations of LAOs and SAOs K1995 ; RW2007 ; RW2007 ; BKR2008 ;
DGKKOW2012 ; HKWS2011 .) The decoherent and coherent dynamics differ
qualitatively and quantitatively, with amplitudes separated by at least one
order of magnitude and frequencies widely separated. The SAOs are phase waves
of high frequency, whereas the LAOs and MMOs jump between distinct states and
are of low frequency.

The existence of mixed-amplitude chimeras is surprising, since the LAOs and
MMOs are robust states in bistable systems. Also, as shown here, in parameter
regimes adjacent to the chimera regime large-amplitude states invade the SAO
regions, pushing out the decoherent SAOs there. Hence, the discovery that
decoherent SAOs can coexist with coherent LAOs and MMOs is of substantial
scientific interest.

Along with the minimal chimeras observed in a pair of symmetrically-coupled
Langyel-Epstein oscillators ABE2019 , these mixed-amplitude chimeras are the
first observations of chimeras in which the structures, amplitudes, and
frequencies of the oscillations differ by orders of magnitude in the
decoherent and coherent regions.

In addition, the new chimeras will be important for applications. Many
bistable systems robustly exhibit SAOs, LAOs, and MMOs. See for example Refs.
KS2009, ; W1997, ; vS2003, ; N2002, ; RG2021, ; AFS2004, ; F1960, ; NAY1962, ;
HH1952, ; KS2009, ; M2006, ; AP1996, ; BKR2008, ; DGKKOW2012, ; HKWS2011, ;
RW2007, ; RW2007, ; K1995, ; LRE1990, ; EM2008, ; DBSK2007, ; Wang2010, .
Among these, models of various brain rhythms, cardiac rhythms, chemical waves,
cell cycle transitions, and cell signaling with positive-feedback loops EM2008
; DBSK2007 ; Wang2010 ; AFS2004 ; RG2021 ; LRE1990 , share structural features
of the model here. In particular, it is shown here that a pair of folded
singularities, known as folded nodes and folded saddles, which arise in
systems exhibiting MMOs, are important mechanisms responsible for creating the
mixed-amplitude chimeras.

Also, a prominent application for which chimeras have been suggested as
suitable models is unihemispheric slow-wave sleep (USWS) AS2004 ; Motter2010 ;
Ramlow2019 ; Haugland2015 ; Ma2010 . In certain aquatic mammals and birds, one
cerebral hemisphere sleeps and the other remains awake Rattenborg2000 ;
Rattenborg2006 . This is illustrated by EEG recordings from bottlenose
dolphins during USWS Rattenborg2000 , in which one hemisphere exhibits
synchronous high-amplitude, low-frequency oscillations characteristic of sleep
SMS1993 , whilst the other shows desynchronized low-amplitude, high-frequency
oscillations characteristic of wakefulness. We find here that the main
properties of mixed-amplitude chimeras, with co-existing regions of coherent
LAOs and MMOs and regions of decoherent SAOs are similar, at least in broad
strokes. Of course, this is still far from establishing a possible role for
mixed-amplitude chimeras in USWS. Instead, the existence of mixed-amplitude
chimeras at least demonstrates that chimeras are relevant when the amplitudes
and frequencies of the coherent and decoherent oscillations are widely
different, something that has not been accounted for by the existing examples
and theory of chimeras. Moreover, the mechanisms responsible for creating the
mixed-amplitude chimeras (namely, the folded singularities of bistable
systens) can be robust to small heterogeneities.

This article is organized as follows. In Section II, we study the prototypical
bistable PDE with spatially-homogeneous reaction terms and purely local,
isotropic diffusion, on which the mixed-amplitude chimeras were discovered.
Also, we present the main dynamical properties and bifurcations of these
chimeras. In Section III, we analyze the folded singularities which are the
mechanisms responsible for allowing the identical oscillators to exhibit both
the decoherent SAOs and the coherent MMOs, i.e., the mechanisms responsible
for the creation of the mixed-amplitude chimeras. In Section IV, the regions
of SAO dynamics are examined in more detail, and it is reported that the
recently-discovered phenomenon of delayed Hopf bifurcations plays an important
role. In Section V, we study the bifurcation of these chimeras to a new type
of pattern which consists of co-existing regions of coherent MMOs and coherent
SAOs. Conclusions and discussion are presented in Section VI. Finally, Section
VII contains some open questions.

##  II Mixed-amplitude chimeras in a prototypical bistable PDE

The new mixed-amplitude chimera states are presented for a canonical multi-
scale bistable PDE: the forced van der Pol equation with local diffusion
KS2009 ; vdP1920 ; VKKR1987 ; SE1988 ; KG1995 ; FLB2002 ; GHW2003 ; SW2004 ;
HB2012 ; VW2015 ; BDGKKV2016 ,

$\begin{split}u_{t}&=v-u^{2}-\tfrac{1}{3}u^{3}+D_{u}u_{xx}\\\
v_{t}&=\operatorname{\varepsilon}(a-u+b\cos\theta)+D_{v}v_{xx}\\\
\theta_{t}&=\operatorname{\varepsilon}\omega.\end{split}$ (1)

The activator (or voltage) $u$ is fast, with bistable nullcline form-fvdPol .
The inhibitor (or recovery) $v$ is slow, and $0<\operatorname{\varepsilon}\ll
1$ measures the scale separation in the reaction terms, as is common in
bistable systems. The excitation threshold $a$, forcing amplitude $b$,
frequency $\omega$, and diffusivities $D_{u}$ and $D_{v}$ are spatially
homogeneous, so that the kinetics and diffusion are identical and purely local
at each point $x\in[0,1]$. The boundary conditions are zero flux. Also,
$\omega=1$ unless stated otherwise so that the forcing period is
$T=\tfrac{2\pi}{\operatorname{\varepsilon}}$. The numerical methods and
typical choices of initial data are described in the Appendix.

Figure 1:  Mixed-amplitude chimeras of (1). Coherent MMOs (red regions,
followed by navy blue shadows) occur at random points and random times. All of
the other regions exhibit decoherent SAOs (medium blue). Here, $a=-0.001$,
$b=0.005$, $\operatorname{\varepsilon}=0.01$. (a)
$(D_{u},D_{v})=(0.45,4)\cdot(\Delta x)^{2}$, (b)
$(D_{u},D_{v})=(1.1,11)\cdot(\Delta x)^{2}$. Contour lines at $u=\pm
2.5\sqrt{\operatorname{\varepsilon}}$.

The mixed-amplitude chimeras of (1) exhibit coexisting regions of decoherent
SAOs and coherent MMOs. The SAOs are classical phase waves about the
depolarized state ($u=0$). Diffusion decoheres them. The MMOs consist of
coherent LAOs between the hyperpolarized state (near $u=-3$) and depolarized
state, alternating with coherent SAOs. The MMOs occur either at random points
and times (Fig. 1) or at fixed locations and time-periodically (Fig. 2),
depending on $(D_{u},D_{v})$.

Figure 2:  Mixed-amplitude chimeras of breathing, turbulent type. (a) A fixed
central band of coherent $2T$-periodic MMOs, and decoherent SAOs in the
complementary regions. $(D_{u},D_{v})=(1.4,9)\cdot(\Delta x)^{2}$. Kinetics
parameters as in Fig. 1. (b) Coherent $5T$-periodic MMOs, and decoherent SAOs
centrally. $(D_{u},D_{v})=(1.7,9)\cdot(\Delta x)^{2}$ and $b=0.002$.

These mixed-amplitude chimeras are breathing, turbulent chimeras in the
classification scheme presented in Ref. KHSKK2016, . The spatial correlation
$g_{0}(t)=\int_{0}^{\delta}g(|\hat{D}|,t)\,d|\hat{D}|$ was measured. Here, $g$
is the normalized probability density of $|\hat{D}|$, where $\hat{D}=(\Delta
x)^{2}D$ is the scaled discrete Laplacian, $\delta=0.01D_{m}$ is the threshold
below which the profile is regarded as coherent, and $D_{m}$ is the absolute
maximal curvature. For (1), $g_{0}(t)$ exhibits a pair of fast jumps (one up
and one down) for each spatially coherent MMO. See Fig. 3(a). Moreover, these
pairs of fast jumps occur at random times for the chimeras in Fig. 1 and time-
periodically for those in Fig. 2 (see Fig. 3(a)). Moreover, between jumps,
$g_{0}(t)$ fluctuates randomly, due to the decoherent SAOs, remaining well
inside the interval $(0,1)$, Hence, the solution is a breathing, turbulent
chimera.

The temporal correlation $h_{0}=\sqrt{\int_{\gamma}^{1}h(|\rho|)\,d|\rho|}$ of
Ref. KHSKK2016,  was also measured. Here, $\rho_{ij}$ is the correlation
coefficient for the time series at $x_{i}$ and $x_{j}$, and $h$ is the
probability distribution function of $|\rho|$. We find that, over large times,
$h_{0}$ converges to zero when the MMO clusters occur randomly. This is
consistent with the observations that the positions of the coherent clusters
change in time. In contrast, for the chimeras in which the coherent MMO
clusters are fixed in space, we find that $h_{0}$ is a non-zero constant
inside $(0,1)$. See Fig. 3(a).

Figure 3:  (a) Coherence measures $g_{0}(t)$ (blue) and $h_{0}$ (red) of the
chimera in Fig. 2(a). Time series of $u$ in the (b) coherent MMO region
($x=0.45$), (c) decoherent SAO domain ($x=0.1,0.21$). The amplitudes and
frequencies differ significantly.

The mixed-amplitude chimeras of (1) are observed robustly where $D_{u}$ is
sufficiently less than $D_{v}$; see Fig. 4. The typical chimeras (marked by
red circles) have coherent MMOs which occur at random points and at random
times, with decoherent SAOs in all other regions, just like the chimeras shown
in Fig. 1. There are also wedges in the $D_{u}-D_{v}$ parameter space in which
the chimeras have coherent MMOs that are time-periodic and spatially fixed
(blue circles), examples of which are given by the chimeras shown in Fig. 2.
Moreover, in all of these chimeras, the decoherent SAOs occupy an increasing
portion of $[0,1]$ as $D_{u}$ decreases, within this regime.

Figure 4: Mixed-amplitude chimeras exist for $D_{u}<D_{v}$
($D_{i}={\widetilde{D_{i}}}(\Delta x)^{2},i=u,v$). (a)  Breathing turbulent
chimeras with randomly-occurring coherent MMOs.  Breathing chimeras with time-
periodic, spatially-fixed coherent MMOs. $\times$ Sharp interface solutions.
$\bigtriangleup$ Spatially uniform, time-periodic MMOs. $\Box$ Trigger waves
exhibiting time-periodic annihilation and nucleation of MMOs. Kinetics as in
Fig. 1 (b) In the strip
${\widetilde{D_{v}}}=(1/\sqrt{\operatorname{\varepsilon}}){\widetilde{D_{u}}}-1=10{\widetilde{D_{u}}}-\alpha$,
for $\alpha\in[0.95,1.05]$, many chimeras have time-periodic MMOs. For each
set of parameters, (1) was integrated for a total time of $750T$.

At the right edge of the mixed-amplitude chimeras regime in Fig. 4, the
chimeras bifurcate to spatially-uniform, time-periodic MMOs (Fig. 5(a)) and to
trigger waves (Fig. 5(b)) BDE1995 . Here, $D_{u}$ is no longer small enough
relative to $D_{v}$. The SAO regions are displaced (in part or entirely) by
MMOs and trigger waves. The trigger wave patterns, such as in Fig. 5(b), can
exhibit annihilation points where pairs of traveling waves of MMOs moving in
opposite directions collide and eliminate each other in the interior of the
spatial domain. This annihilation event is typically followed some time later
by a nucleation point at the same spatial location where a pair a trigger
waves emerge and radiate away from the nucleation point in opposite
directions.

Further, we observe that the presence of the spatially coherent MMOs and
trigger waves in this parameter regime immediately adjacent to the regime of
mixed-amplitude chimeras makes the existence of these chimeras all the more
striking. In mixed-amplitude chimeras, the decoherent SAOs resist the invasion
of the stable LAO/MMOs.

At the left edge of the mixed-amplitude chimera regime in Fig. 4, the chimeras
bifurcate from sharp-interface solutions N2002 . As shown in Fig. 5(c), these
sharp interface solutions have steep jumps at fixed locations connecting the
depolarized and hyperpolarized states. Then, on the spatial intervals between
the jumps, these solutions exhibit time-periodic SAOs.

Figure 5: (a) A representative spatially homogeneous, $T$-periodic MMO state
of (1) of the type denoted by $\bigtriangleup$ in Fig. 4. Here,
${\widetilde{D_{u}}}=1$ and ${\widetilde{D_{v}}}=5$. (b) A representative
trigger wave exhibiting time-periodic annihilation and nucleation of MMOs, of
the type denoted by $\Box$ in Fig. 4. Here, ${\widetilde{D_{u}}}=1.7$ and
${\widetilde{D_{v}}}=7$. (c) A representative sharp-interface solution of the
type denoted by $\times$ in Fig. 4. Here, ${\widetilde{D_{u}}}=0.15$ and
${\widetilde{D_{v}}}=3$.

Overall, as shown in Fig. 4, it is also a novelty of mixed-amplitude chimeras
that they bifurcate from large-amplitude solutions, which are ubiquitous in
bistable systems, and not from homogeneous or drift states.

Mixed-amplitude chimeras also exhibit multi-stability and symmetry-breaking.
For example, the chimera state in Fig. 2(a) coexists with three other stable
states of (1) for the same parameters. There is a co-existent spatially-
homogeneous SAO state. Also, there are two other mixed-amplitude chimeras, one
of which resembles the reverse image, with coherent MMOs along the boundaries
and a central strip of decoherent SAOs, and the other of which has time-
periodic coherent MMOs both in the central region, as the chimera in Fig.
2(a), and near the boundaries, as well as decoherent SAOs everywhere else.
Furthermore, symmetry-breaking occurs through bifurcations to chimeras in
which the widths of the MMO regions vary time-periodically.

##  III Folded singularities as a mechanism for mixed-amplitude chimeras

A distinct mechanism, structurally different from that in known chimeras, is
at the heart of the mixed-amplitude chimeras. There is a pair of commonly-
occurring singularities, known as folded nodes (FNs) and folded saddles (FSs),
in the kinetics of bistable PDEs. These singularities arise naturally in the
van der Pol oscillator (1) GHW2003 ; VW2015 ; BDGKKV2016 , and in
neuroscience, physics, chemistry, and electrical engineering models BBE1991 ;
RW2007 ; HKWS2011 ; TTVWB2011 ; DGKKOW2012 ; WMR2013 ; KAA2014 ; MKP2014 ;
RRW2015 ; ABE2019 ; BBU2020 ; SW2004 ; W2005 ; RW2007 ; W2005 ; WMR2013 ;
KAA2014 ; WMR2013 ; BBU2020 .

The FNs and FSs lie on a fold curve between attracting and repelling sheets of
the slow manifold. For (1), they are located at the points $\theta_{N,S}$
where $a+b\cos(\theta_{N,S})=0$, for all $0<|a|<b$, as may be calculated from
the ODE consisting of the reaction terms in (1), see Refs. BDGKKV2016, ; FNFS,
. See Fig. 6.

In the decoherent regions, the time traces of the SAOs are observed to pass
through neighbourhoods of FNs and FSs in alternation, Fig. 6(a,b). Near a FN,
the SAOs are centred on its weak canard SW2004 ; W2005 . As the time traces
approach the next FS, the SAOs are centred about its faux canard WMR2013 . The
weak and faux canards (approximated by the green curves) lie ${\cal
O}(b)$-close to the fold curve, as expected from the ODE VW2015 . In the
passage from an FN to an FS and then through a neighbourhood of the FS, the
amplitude grows until there is an excursion (near the repelling branch of the
slow manifold), into the funnel region (gray) of the next FN. Hence, the time
trace enters the neighbourhood of the next FN, since the funnels are basins of
attraction for the FNs RW2007 , and there it has further SAOs about the weak
canard of that FN.

Figure 6:  Time traces (black) of the chimera in Fig. 2, at $x=0.23$ in (a)
and (b), and at $x=0.41$ in (c) and (d). In the decoherent SAO regions, the
time traces oscillate about the weak canards of the FNs (blue dots) and the
faux canards of the FSs (red dots), near the attracting (upper) sheet and
repelling (middle) sheet of the slow manifold. In the coherent MMO regions,
the time traces also have LAOs down to the lower branch of the attracting slow
manifold.

In the coherent region(s), the time traces are seen to be MMOs induced by the
canards of the folded singularities RW2007 ; TTVWB2011 ; HKWS2011 ; DGKKOW2012
(Fig. 6(c,d)). The SAO segments of these MMOs are also about the FN weak
canard. However, here, after passage near a FN, there is a large-amplitude
excursion toward the lower attracting sheet of the slow manifold.
Subsequently, the traces crawl along this sheet until they reach the lower
fold curve, where they jump back to the upper sheet and into the funnel of the
next FN, and the cycle repeats. (The funnel is enclosed by the true canard
(magenta) of the FS WMR2013 ; WMR2013 ; BBU2020 , the strong canard (cyan) of
the FN W2005 ; RW2007 , and the fold curve.) These MMOs remain spatially
coherent in the presence of diffusion, due to the LAOs.

Figure 7:  Mixed-amplitude chimeras of (1) are robust in the goldilocks zone
(gold) of the $(\omega,a)$ plane. It sits in the middle of the MMO region,
which is enclosed by $a=a_{c}(\operatorname{\varepsilon})\pm
be^{-\varepsilon\omega^{2}/2}$, i.e., by the saddle-node bifurcations of the
FN and FS canards BDGKKV2016 . For the classical van der Pol equation,
$a_{c}(\operatorname{\varepsilon})=-\frac{\operatorname{\varepsilon}}{8}-\frac{3\operatorname{\varepsilon}^{2}}{32}+\mathcal{O}(\operatorname{\varepsilon}^{3})$
BCDD1980 .

Overall, we find that there is a robust zone –which we label as a goldilocks
zone– in which the mixed-amplitude chimera states exist in (1). See Fig. 7. In
this goldilocks zone, the distances between the FNs and FSs are simultaneously
small and large enough. The distances are small enough that, on wide
intervals, the attractor of (1) can pass near the FNs and FSs in alternation,
with the amplitude remaining small and diffusion causing the decoherence. They
are also large enough that, in complementary MMO regions, the return
mechanisms for the LAO segments are robust.

The goldilocks zone shown in Fig. 7 is representative of the zones in the
$(\omega,a)$ plane for a range of values of $b$. Hence, mixed-amplitude
chimeras are robust in (1).

##  IV Delayed Hopf bifurcations underlie the decoherent SAOs

Having presented the novel mixed-amplitude chimeras, the parameter regime in
which they exist, the bifurcations in which they are created or annihilated,
and the folded nodes and folded saddles which make it possible for these
identical oscillators to exhibit both SAO and LAO/MMO dynamics, we next turn
to a fundamental question about the SAO dynamics. Namely, what is the
mechanism by which the solutions in the decoherent regions can remain close to
the fold curve long enough so that they can pass through the neighbourhoods of
both the FN and the FS, hence keeping the oscillation amplitudes small? This
question is important, since it is expected that solutions which pass through
a neighbourhood of a FN generally are repelled away from the fold curve and
rapidly approach the (hyperpolarized) attracting branch of the slow manifold,
resulting in large-amplitudes. It turns out, as we show here, that the SAOs in
the decoherent regions can exist due to delayed passage through Hopf
bifurcations (DHB). That is, DHB makes it possible for large measures of
solutions to stay near the fold curve (and hence maintain small amplitude) so
that they not only pass through a neighbourhood of the FN, but also stay close
enough to pass through a neighbourhood of the FS and into the funnel of the
next FN, so that the cycle can repeat.

To show this, we study (1) with parameters chosen so that the FN and the FS
are close to $\theta=0\ {\rm mod}\ 2\pi$. This includes the parameters in the
goldilocks zone.

First, we show that the local dynamics of (1) may be studied using two coupled
CGL-type equations with a slowly-varying parameter. In particular, we rescale
$u=\sqrt{\operatorname{\varepsilon}}u_{2}$,
$v={\operatorname{\varepsilon}}v_{2}$,
$\theta={\operatorname{\varepsilon}}^{\frac{1}{4}}\theta_{2}$,
$\tau=\sqrt{\operatorname{\varepsilon}}t$,
$a+b=\sqrt{\operatorname{\varepsilon}}\alpha_{2}$, and
$D_{i}=\sqrt{\operatorname{\varepsilon}}d_{i}$ ($i=u,v$), so (1) becomes

$\begin{split}{u_{2}}_{\tau}&={v_{2}}-u_{2}^{2}-\tfrac{1}{3}\sqrt{\operatorname{\varepsilon}}u_{2}^{3}+d_{u}{u_{2}}_{xx}\\\
{v_{2}}_{\tau}&=-u_{2}+\alpha_{2}-\tfrac{1}{2}b\theta_{2}^{2}+{\cal
O}(\sqrt{\operatorname{\varepsilon}}\theta_{2}^{4})+d_{v}{v_{2}}_{xx}\\\
{\theta_{2}}_{\tau}&=\operatorname{\varepsilon}^{\frac{1}{4}}\omega.\end{split}$

This rescaling focuses on the dynamics in the neighborhood of a FN-FS pair.
This local model consists of a fast-slow system, with $(u_{2},v_{2})$ as the
fast variables and $\theta_{2}$ as the slow variable.

We rectify the critical manifold to the $\theta_{2}$-axis using
${\widetilde{u}_{2}}=u_{2}+\mu$ and ${\widetilde{v}_{2}}=v_{2}-\mu^{2}$, where
$\mu=-\alpha_{2}+\frac{b}{2}\theta_{2}^{2}$. Also, we switch to the complex
coordinates $z={\widetilde{u}_{2}}+i{\widetilde{v}}_{2}$ and
${\bar{z}}={\widetilde{u}_{2}}-i{\widetilde{v}}_{2}$, which are more natural
variables for studying the slow passage through a Hopf bifurcation. Hence, (1)
is locally equivalent to

$\begin{split}\begin{bmatrix}z\\\
\overline{z}\end{bmatrix}_{\tau}&=\begin{bmatrix}\mu-i&\mu\\\
\mu&\mu+i\end{bmatrix}\begin{bmatrix}z\\\
\overline{z}\end{bmatrix}-Q\begin{bmatrix}1\\\
1\end{bmatrix}+\begin{bmatrix}d_{+}&d_{-}\\\
d_{-}&d_{+}\end{bmatrix}\begin{bmatrix}z\\\ \overline{z}\end{bmatrix}_{xx}\\\
{\mu}_{\tau}&=\operatorname{\varepsilon}^{\frac{1}{4}}\omega\sqrt{2b(\mu+\alpha_{2})},\end{split}$

where $Q=\tfrac{1}{4}(z+{\bar{z}})^{2}$ and $d_{\pm}=\tfrac{1}{2}(d_{u}\pm
d_{v})$. The origin is an attracting spiral of the $(z,{\bar{z}})$ system for
each $\mu<0$ and a repellor for each $\mu>0$.

Due to the slow growth in $\mu$, this PDE exhibits a dynamic Hopf bifurcation
at $\mu=0$. In particular, solutions initialised with $\mu<0$ do not
immediately diverge from the equilibrium when $\mu$ crosses zero and becomes
positive. Instead, they remain near the repelling state for long times beyond
the dynamic Hopf point (data not shown), in a delayed Hopf bifurcation (DHB).
The time of the delay, past $\mu=0$, is $\mathcal{O}(1)$ in the slow time, so
that solutions spend long times near the repelling sheet of the slow manifold
in the robust SAO regions. Moreover, these times are long enough so that the
solutions can pass near a pair of FN and FS and reach the basin of attraction
(funnel) of the next FN, keeping the oscillation amplitude small.

Remark. For this system of coupled CGL-type PDEs, the observed DHB is
consistent with the recent discovery of DHB in nonlinear PDEs KV2018 ; GKV2020
, including the Hodgkin-Huxley PDE, Brusselator model, and CGL equation. See
also Ref. ADVW2020, .

Remark. The analysis here of DHB in these PDEs was inspired by the ODE
analysis in Ref. VW2015, . There, singularities known as Folded Saddle Nodes
of Type 1 (FSN I), which arise in bifurcations of FNs and FSs, are studied in
ODEs. The analysis in Ref. VW2015,  shows that delayed Hopf bifurcation plays
a central role in ODEs near FSN I bifurcations. Hence, it was natural here to
examine the regime in which the FN and FS of (1) are close together, and to
look for DHB also in the PDE (1).

##  V Bifurcation to coherent-coherent states

The above analyses of the folded singularities (Sec. III) and of the scaling
leading to the coupled CGL-type model for the local dynamics (Sec. IV) also
apply for larger $b$ (relative to $\operatorname{\varepsilon}$), outside the
main chimera regime. With stronger forcing (larger $b$), the attractors also
consist of regions with coherent MMOs and complementary regions of SAOs.
However, we find that with larger $b$ the pure-SAO states are coherent in the
presence of diffusion, rather than decoherent as in the mixed-amplitude
chimeras. These coherent-coherent patterns are multi-mode attractors in the
sense of Ref. VBK2020, .

Examples of these multi-mode attractors, which feature coherent MMOs and
coherent SAOs are shown in Fig. 8. Some of these have multiple sharp
interfaces separating the regions, and others do not have any sharp
interfaces. For example, in Fig. 8(b), one finds an attractor of (1) with
three distinct types of regions: a region of coherent MMOs, a region in which
the time traces are near the hyperpolarized state (near $u=-3$), and a region
of coherent SAOs about the depolarized state. With the knowledge of the FNs
and FSs, one can construct rich patterns.

Furthermore, in these coherent-coherent states that exist for larger $b$, the
time traces of the attractor also exhibit maximal spatiotemporal canards in
the transition intervals between the two distinct states (as for example near
the bulbs in Fig. 8(c)). These maximal spatiotemporal canards mediate the
transitions between regions of distinct oscillatory behaviour, similar to
those in Ref. VBK2020, .

Figure 8: A sample of coherent-coherent multi-mode attractors for various
ratios of $\operatorname{\varepsilon}$ to $b$. (a)
$\operatorname{\varepsilon}=b=0.005$ with $(D_{u},D_{v})=(0.47,4)\cdot(\Delta
x)^{2}$. (b) $\operatorname{\varepsilon}=b=0.01$ with
$(D_{u},D_{v})=(1.71,9)\cdot(\Delta x)^{2}$. (c)
$\operatorname{\varepsilon}=\tfrac{1}{2}b=0.005$ with
$(D_{u},D_{v})=(1.40625,11.25)\cdot(\Delta x)^{2}$. (d)
$\operatorname{\varepsilon}=\tfrac{1}{5}b=0.005$ with
$(D_{u},D_{v})=(0.325,8)\cdot(\Delta x)^{2}$. Finally, for each fixed pair of
$\operatorname{\varepsilon}$ and $b$, the kinetics parameter $a$ is chosen to
lie in the goldilocks zone: (a) $a=-0.001$, (b) $a=-0.001$, (c)
$a=-\tfrac{\operatorname{\varepsilon}}{8}+\tfrac{1}{10}be^{-\operatorname{\varepsilon}\omega^{2}/2}\approx
0.000372503$ and (d)
$a=-\tfrac{\operatorname{\varepsilon}}{8}+\tfrac{1}{20}be^{-\operatorname{\varepsilon}\omega^{2}/2}\approx
0.000621879$.

##  VI Conclusions and discussion

In this article, we presented the novel class of chimeras labeled as mixed-
amplitude chimeras, which consist of co-existing decoherent SAOs and coherent
LAOs or MMOs. This constitutes the first observation of chimeras in PDEs in
which the structures, amplitudes, and frequencies of the oscillations are
widely separated in the two distinct types of regions. With the exception of
Ref. ABE2019, , these chimeras differ from the known chimeras and theory,
where the amplitudes and frequencies in the coherent and decoherent regions
are comparable.

These mixed-amplitude chimeras are surprising, since the coherent LAOs and
MMOs have much larger amplitudes and since they are structurally stable.
Moreover, in parameter regimes immediately adjacent to the chimera regime, the
LAOs and MMOs are observed to invade the decoherent regions and wipe out the
SAOs. Indeed, new bifurcations were also discovered in which mixed-amplitude
chimeras emerge from spatially-uniform MMOs, trigger waves, and sharp-
interface solutions, all of which are ubiquitous in bistable PDEs.

Initial theory has been identified for the mixed-amplitude chimeras. It was
shown here that they are generated by folded node and folded saddle
singularities, common in bistable systems. Moreover, it was discovered that
there is a robust goldilocks zone in which the mixed-amplitude chimeras exist,
where the inter-singularity distances are simultaneously small and large
enough so the decoherent SAOs and coherent MMOs can co-exist. Also, it was
shown that Delayed Hopf Bifurcations in PDEs enables the existence of the SAOs
in the decoherent regions.

The mixed-amplitude chimeras may have important applications in bistable
systems, which are often used to model biological, chemical, and physical
rhythms.

The mixed-amplitude chimeras also bear some of the hallmarks of EEG recordings
taken from animals that exhibit USWS Rattenborg2000 . Namely, the regions of
coherent MMOs consist of large-amplitude, low-frequency activity as observed
in the sleeping hemisphere, and the regions of decoherent SAOs consist of
small-amplitude, high-frequency activity similar to the awake hemisphere.

The substantial contrasts in amplitudes and frequencies observed here for the
mixed-amplitude chimeras of (1) are similar in some respects to the contrasts
reported recently in minimal chimeras in Ref. ABE2019, . There, a two-cell
model consisting of a pair of symmetrically-coupled Lengyel-Epstein ODEs is
studied. The parameters are tuned near a canard point of the L-E system, with
one oscillator exhibiting periodic LAOs and the other aperiodic SAOs.
Moreover, the existence of new types of phase relations between coupled cells
was established, distinct from the classical in-phase and anti-phase
attractors in coupled fast-slow systems. Using analysis similar to that
presented in this article, it can be shown that the underlying singularity in
Ref. ABE2019,  is a (type III) folded saddle-node, and the parameters are in
its goldilocks zone.

The mixed-amplitude chimeras introduced here are distinct from the Multi-Mode
Attractors presented in Ref. VBK2020, . There, the MMAs were introduced in
systems with spatially-dependent parameters. The spatially-dependent applied
currents (or spatially-dependent maximal ion-channel conductivities) give rise
to attractors that exhibit different modes of oscillations in different
spatial regions. Hence, in the systems with MMAs, the oscillators are not all
identical. In contrast, here, the reactions are all spatially homogeneous and
the diffusion is purely local and isotropic, so that the oscillators at each
point are identical. See also the discussion in the introductory section 1 of
Ref. VBK2020, .

##  VII Some open questions

The discovery of mixed-amplitude chimeras raises a number of questions. First,
it is of interest to investigate the mechanism(s) by which some of the mixed-
amplitude chimeras can have coherent MMOs that occur at fixed locations in
space and periodically in time. These are indicated by the blue circles in the
$D_{u}-D_{v}$ parameter plane in Fig. 4. As discussed above, these are
distinct from the typical chimeras (red circles), in which the coherent MMOs
occur at random points and at random times. There may be some type of
resonance responsible for creating these wedges of blue circles within the sea
of red circles.

Second, we ask if some of the known examples of chimeras in other systems can
be continued to parameter regimes in which the structures, amplitudes, and/or
frequencies of the coherent oscillations are substantially different from
those of the decoherent SAOs, as is the case for the mixed-amplitude chimeras
here. The known examples and theory of chimera states in coupled oscillators,
integro-differential systems, and PDEs are based on the amplitude-phase method
or mean-field approach KB2002 ; AS2004 ; PA2015 ; O2018 ; SS2014 , symmetry
and broken symmetry analyses KHK2018 ; ZKS2014 , and the theory Zhang2021 ;
OMT2008 of desynchronization-induced stabilization of coherent oscillations,
among other important results. These methods have been used primarily to find
chimeras in which the amplitudes and frequencies of the oscillations in the
decoherent and coherent regions are comparable, i.e., where they do not differ
by orders of magnitude. For example, the analysis used in the phase-amplitude
method is premised on the fact that the oscillations in the chimera states and
the states from which bifurcate have similar amplitudes and frequencies.
Hence, it would of interest to explore whether any of these may be continued
in parameter space so that the coherent regions exhibit LAOs and MMOs.

In the other direction, one may ask wheher the new mixed-amplitude chimeras
can be continued in parameter space to a regime in which the coherent
oscillations are near homogeneous states (and hence purely of SAO type),
instead of being MMOs or LAOs.

Third, there is a natural question about what the minimal dimension is for a
system to exhibit mixed-amplitude chimeras. We recall that, in ODEs, three is
the minimum number of variables necessary for MMOs. This is the case for the
forced van der Pol reaction terms in (1). It is also the case for the Lengyel-
Epstein ODEs studied in Ref. ABE2019, , where an individual L-E model has two
dependent variables, and the systems studied there consist of two
symmetrically-coupled L-E oscillators (i.e., four variables), or of a single
L-E oscillator with some type of external forcing (i.e., three variables).
Hence, the results of Ref. ABE2019,  and this article suggest that three
variables in the kinetics is expected to be the minimum number needed for
chimera states to exhibit co-existing regions of coherent MMOs and decoherent
SAOs, such as the mixed-amplitude chimeras reported here.

###### Acknowledgements.

The authors thank Karen Corbett for useful conversations and support with
MASSIVE. We also thank Irv Epstein, Ryan Goh, Christian Kuehn, Jonathan
Touboul, Gene Wayne, and Martin Wechselberger for useful comments and
questions. The research of T.V. was partially supported by NSF-DMS 1853342.
The research of T.K. was partially supported by NSF-DMS 1616064. This work was
supported by the MASSIVE HPC facility (www.massive.org.au).

##  VIII Appendix

The numerical simulations were performed using the method of lines, with
adaptive time-stepping suited to stiff problems, and a second-order spatial
discretization. The results were confirmed independently using four other
methods: balanced symmetric Strang operator splitting, with (second-order)
centered finite differences in space and fourth-order Runge-Kutta in time; the
Crank-Nicholson method, which is second-order accurate in both space and time;
a compact finite difference method which is fourth-order accurate in space;
and, a Chebyshev grid method that is spectral in space. All simulations were
run at least for the order of the diffusive time scale, based on the smaller
diffusivity. Unless stated otherwise, $N_{x}=750$ spatial subintervals were
used ($\Delta x=\frac{1}{750}$). Also, convergence was checked with smaller
values of $\Delta x$.

The main method used to generate initial data $(u_{0}(x),v_{0}(x))$ on $[0,1]$
is to compute a limit cycle attractor $(\Gamma_{u}(t),\Gamma_{v}(t))$ (period
$T$) of the kinetics ODE, and then to set
$(u_{0}(x),v_{0}(x))=\left(\Gamma_{u}\left(\frac{t}{nT}\right),\Gamma_{v}\left(\frac{t}{nT}\right)\right)$.
Typical initial data consists of $n$ copies spread out in space, with
$\theta(t=0)=0$. For example, $n=6$ in Fig. 4. Eight other types of initial
data were also used to explore more fully the multi-stable landscape of
solutions, and to find chimera states.

## References

  * (1)
  * (2)
  * (3) Y. Kuramoto and D. Battogtokh, Nonlin. Phen. Complex Sys. 5, 380 (2002). Y. Kuramoto, in Nonlin. Dyn. Chaos: Where Do We Go from Here?, S.J. Hogan, A.R. Champneys, B. Krauskopf, M. di Bernardo, H.M. Osinga, and M.E. Homer (eds.), IOP, Bristol, UK, 209 (2003). Y. Kuramoto and S.I. Shima, Prog. Theor. Phys. Suppl. 150, 115 (2003). S.I. Shima and Y. Kuramoto, Phys. Rev. E 69, 036213 (2004).
  * (4) D.M. Abrams and S.H. Strogatz, Phys. Rev. Lett. 93, 174102 (2004). D.M. Abrams, R.E. Mirollo, S.H. Strogatz, D.A. Wiley, Phys. Rev. Lett. 101, 084103 (2008).
  * (5) E. Ott and T.M. Antonsen, Chaos 18, 037113 (2008). E. Ott and T.M. Antonsen, Chaos 19, 023117 (2009).
  * (6) O.E. Omel’chenko, Y.L. Maistrenko, and P.A. Tass, Phys. Rev. Lett. 100, 044105 (2008).
  * (7) G. Bordyugov, A. Pikovsky, and M. Rosenblum, Phys. Rev. E 82, 035205(R) (2010).
  * (8) C. Laing, Phys. Rev. E 81, 066221 (2010).
  * (9) C. Laing, Phys. Rev. E 92, 050904(R) (2015). C. Laing, Phys. Rev. E 100, 042211 (2019).
  * (10) I. Omelchenko, Y. Maistrenko, P. Hövel, and E. Schöll, Phys. Rev. Lett. 106, 234102 (2011). I. Omelchenko, O. E. Omel’chenko, P. Hövel, and E. Schöll, Phys. Rev. Lett. 110, 224101 (2013).
  * (11) M.R. Tinsley, S. Nkomo, and K. Showalter, Nature Phys. 8, 662 (2012).
  * (12) A.M. Hagerstrom, T. E. Murphy, R. Roy, I. Omelchenko, P. Hövel, and E. Schöll, Nature Phys. 8, 658 (2012).
  * (13) E.A. Martens, S. Thutupalli, A. Fourriere, and O. Hallatschek, PNAS 110, 10563 (2013).
  * (14) S. Nkomo, M.R. Tinsley, and K. Showalter, Phys. Rev. Lett. 110, 244102 (2013).
  * (15) K. Schönleber, C. Zensen, A. Heinrich, and K. Krischer, New J. Phys. 16, 063024 (2014). L. Schmidt, K. Schönleber, K. Krischer, and V. Garcia-Morales, Chaos 24, 013102 (2014). L. Schmidt and K. Krischer, Chaos 25, 064401 (2015).
  * (16) G.C. Sethia and A. Sen, Phys. Rev. Lett. 112, 144101 (2014).
  * (17) A. Yeldesbay, A. Pikovsky, and M. Rosenblum, Phys. Rev. Lett. 112, 144103 (2014).
  * (18) A. Zakharova, M. Kapeller, and E. Schöll, Phys. Rev. Lett. 112, 154101 (2014).
  * (19) I. Omelchenko, A. Provata, J. Hizanidis, E. Schöll, and P. Hövel, Phys. Rev. E 91, 022917b (2015).
  * (20) M.J. Panaggio and D.M. Abrams, Nonlinearity 28, R67 (2015).
  * (21) M.G. Clerc, S. Coulibaly, M.A. Ferre, M.A. Garcia-Nustes, and R.G. Rojas, Phys. Rev. E 93, 052204 (2018).
  * (22) F. P. Kemeth, S. W. Haugland, L. Schmidt, I. G. Kevrekidis, and K. Krischer, Chaos 26, 094815 (2016).
  * (23) B.-W. Li and H. Dierckx, Phys. Rev. E 93, 020202(R) (2016). B.-W. Li, Y. He, L.-D. Li, L. Yang, X. Wang, Arxiv 2012.00983, nlin, Dec 2020.
  * (24) Z.G. Nicolau, H. Riecke, and A.E. Motter, Phys. Rev. Lett. 119, 244101 (2017).
  * (25) F.P. Kemeth, S.W. Haugland, and K. Krischer, Phys. Rev. Lett. 120, 214101 (2018).
  * (26) O.E. Omel’chenko, Nonlinearity 31, R121 (2018).
  * (27) J.F. Totz, J. Rode, M.R. Tinsley, K. Showalter, and H. Engel, Nat. Phys. 14, 282 (2018). J. Rode, J.F. Totz, E. Fengler, and H. Engel, Frontiers Appl. Math. Stat. 5, 31 (2019). J.F. Totz, M.R. Tinsley, H. Engel, and K. Showalter, Nat. Sci. Rep. 10, 7821 (2020).
  * (28) S. Kundu, S. Majhi, P. Muruganandam, D. Ghosh, Eur. Phys. J. Spec. Topics 227, 983–993 (2018).
  * (29) N.M. Awal, D. Bullara, and I.R. Epstein, Chaos 29, 013131 (2019). N.M. Awal and I.R. Epstein, Phys. Rev. E 101, 042222 (2020).
  * (30) I.A. Shepelev and V.S. Anischenko, Comm. Nonlin. Sci., Num. Sim. 93, 105513 (2021).
  * (31) Y. Zhang and A.E. Motter, Phys. Rev. Lett. 126, 094101 (2021).
  * (32) M. Brons, T.J. Kaper, and H.G. Rotstein, Chaos 18, 015101 (2008).
  * (33) M. Desroches, J. Guckenheimer, B. Krauskopf, C. Kuehn, H.M. Osinga, and M. Wechselberger SIAM Rev. 54, 211 (2012).
  * (34) E. Harvey, V. Kirk, M. Wechselberger, and J. Sneyd, J. Nonlin. Sci., 21 639–683 (2011).
  * (35) J.E. Rubin and M. Wechselberger, Bio. Cybernetics 97, 5 (2007). M. Brøns, M. Krupa, and M. Wechselberger, Fields Inst. Commun. 49, 39-63 (2006).
  * (36) M. Koper, Physica D 80, 72 (1995).
  * (37) D. Walgraef, Spatio-Temporal Pattern Formation, Springer, New York (1997).
  * (38) W. van Saarloos, Phys. Rep., 386, 29 (2003).
  * (39) Y. Nishiura, Far-from Equilibrium Dynamics, Amer. Math. Soc., Transl. Math. Monographs, 209 (2002).
  * (40) J. Keener and J. Sneyd, Mathematical Physiology I: Cellular Physiology, 2nd ed., Springer, New York (2009).
  * (41) A. Hodgkin and A. Huxley, J. Physiol. 117, 500 (1952).
  * (42) J. Moehlis, J. Math. Bio. 52, 141 (2006).
  * (43) R. FitzHugh, J. Gen. Physiol. 43, 867 (1960).
  * (44) J. Nagumo, S. Arimoto, and S. Yoshizawa, Proc. IRE 50, 2061 (1962).
  * (45) R.R. Aliev and A.V. Panfilov, Chaos Soliton Fract 7, 293 (1996).
  * (46) D. Angeli, J.E. Ferrell, and E.D. Sontag, PNAS, 101, 1822 (2004).
  * (47) J. Rombouts and L. Gelens, PLOS Comp. Bio., 17, 1008231 (2021).
  * (48) I. Lengyel, G. Rabai, and I.R. Epstein, J. Am. Chem. Soc. 112, 9104 (1990). I. Lengyel and I.R. Epstein, Science 251, 650 (1991).
  * (49) I. Erchova and D.J. McGonigle, Chaos 18, 015115 (2008).
  * (50) C.G. Diniz Behn, E.N. Brown, T.E. Scammell, and N. Kopell, J. Neurophysiol. 97, 3828-3840 (2007).
  * (51) X.-J. Wang, Physiol. Rev. 90, 1195-1268 (2010).
  * (52) A.E. Motter, Nat. Phys. 6, 164–165 (2010).
  * (53) L. Ramlow, J. Sawicki, A. Zakharova, J. Hlinka, J.C. Claussen, and E. Schöll, EPL 126, 50007 (2019).
  * (54) S.W. Haugland, L. Schmidt, and K. Krischer, Sci. Rep.-UK 91, 40006 (2010).
  * (55) R. Ma, J. Wang, and Z. Liu, EPL 91, 40006 (2010).
  * (56) N.C. Rattenborg, C.J.Amlaner, and S.L. Lima, Neurosci. Behav. Rev. 24, 817–842 (2000).
  * (57) N.C. Rattenborg, Naturwissensch. 93, 413–425 (2006).
  * (58) M. Steriade, D.A. McCormick, and T.J. Sejnowski, Science 262, 679-685 (1993).
  * (59) B. van der Pol, Radio Review 1, 701 (1920). B. van der Pol, Lond., Edinb., Dublin Phil. Mag. & J. Sci., Series 7, 3, 65–80 (1927).
  * (60) A.B. Vasil’eva, S.A. Kashchenko, Yu.S. Kolesov, and N.Kh. Rozov, Math. USSR Sbor. 58, 491 (1987).
  * (61) K.R. Schneider and V.M. Evtukhov, Differentsial’nye Urav. 24, 1027 (1988) (Eng. transl.: Diff. Eq. 24, 683)
  * (62) D. Kaplan and L. Glass, Understanding Nonlinear Dynamics, Springer, New York (1995).
  * (63) M.L. Fachinetti, E. de Langre, and F. Biolley Comp. Rend. Acad. Sci., Paris, Ser. II 330 (2002).
  * (64) J. Guckenheimer, K. Hoffman, and W. Weckesser, SIAM J. Appl. Dyn. Sys. 2, 1 (2003); K. Bold, C. Edwards, J. Guckenheimer, S. Guharay, K. Hoffman, J. Hubbard, R. Oliva, and W. Weckesser, ibid. 2, 570 (2003).
  * (65) T. Vo and M. Wechselberger, SIAM J. Math. An. 47, 3235 (2015).
  * (66) J. Burke, M. Desroches, A. Granados, T.J. Kaper, M. Krupa, and T. Vo, J. Nonlin. Sci. 26, 405 (2016).
  * (67) P. Szmolyan and M. Wechselberger, J. Diff. Eq. 200, 69 (2004).
  * (68) X. Han and Q. Bi, Nonlinear Dynamics 68, 275 (2012).
  * (69) The kinetics in (1) are equivalent to the classical $x_{t}=y-\left(\frac{1}{3}x^{3}-x\right)$, $y_{t}=\operatorname{\varepsilon}(\tilde{a}-x+b\cos(\theta))$ under the linear shift $u=x-1$, $v=y+\frac{2}{3}$, and $a=\tilde{a}-1$.
  * (70) F. Buchholtz, M. Dolnik, and I.R. Epstein, J. Phys. Chem. 99, 15093 (1995).
  * (71) W. Teka, J. Tabak, T. Vo, M. Wechselberger, and R. Bertram, J. Math. Neuro. 1, 12 (2011).
  * (72) M. Wechselberger, J. Mitry, and J. Rinzel, in Nonautonomous Dynamical Systems in the Life Sciences, LNM 2012, Springer, 89 (2013). J. Mitry and M. Wechselberger, SIAM J. Appl. Dyn. Sys. 16, 546 (2017).
  * (73) M. Krupa, B. Ambrosio, and M.A. Aziz-Alaouio, Nonlinearity 27, 1555 (2014).
  * (74) K.-L. Roberts, J.E. Rubin, and M. Wechselberger, SIAM J. Appl. Dyn. Sys. 14, 1808 (2015).
  * (75) E. Bossolini, M. Brøns, and K. Uldall-Kristiansen, SIAM Rev. 62, 869 (2020).
  * (76) M. Brøns and K. Bar-Eli, J. Phys. Chem 95, 8706 (1991).
  * (77) P. de Maesschalck, E. Kutafina, and N. Popovic, J. Dyn. Diff. Eq. 26, 955 (2014).
  * (78) M. Wechselberger, SIAM J. Appl. Dyn. Sys. 4, 101 (2005). M. Krupa and M. Wechselberger, J. Diff. Eq. 248, 2841 (2010).
  * (79) The local minimum of the critical manifold $\\{(u,v,\theta)|v=u^{2}+\frac{1}{3}u^{3}\\}$, for $\operatorname{\varepsilon}=0$. On this manifold, the reduced flow is $u(u+2)\dot{u}=a-u+b\cos(\theta)$ and $\dot{\theta}=\omega$, as derived by taking $\frac{d}{dt}$ of the algebraic equation defining the manifold. The regular FN dynamics occur for $b^{2}<a^{2}+\tfrac{1}{64\omega^{2}}$.
  * (80) E. Benoit, J. L. Callot, F. Diener, and M. Diener, Collectanea Mathematica 31, 3 (1980). W. Eckhaus, in Lecture Notes in Math. 985, Springer-Verlag, 449 (1983). A.K. Zvonkin and M.A. Shubin, Russ. Math. Surv. 39, 69 (1984). M. Brøns, in Proc. 9th AIMS Intl. Conf., Disc. Cont. Dyn. Sys., Suppl. 77-83 (2013).
  * (81) T.J. Kaper and T. Vo, Chaos 28, 091103 (2018).
  * (82) R. Goh, T.J. Kaper, and T. Vo, ArXiv DS:2012.10048 (2020).
  * (83) D. Avitabile, M. Desroches, R. Veltz, and M. Wechselberger, SIAM J. Math. An. 52, 5703 (2020).
  * (84) T. Vo, R. Bertram, and T. J. Kaper, Physica D 411, 132544 (2020).


11institutetext: Sberbank, 117997 Moscow, Russia, 22institutetext: National
University of Science and Technology “MISIS”, 119991 Moscow, Russia, 22email:
<EMAIL_ADDRESS>

# Citation network applications in a scientific co-authorship recommender
system

Vladislav Tishin 11 2 2 Artyom Sosedka 11 2 2 Peter Ibragimov 22
Vadim Porvatov 11 2 2

###### Abstract

The problem of co-authors selection in the area of scientific collaborations
might be a daunting one. In this paper, we propose a new pipeline that
effectively utilizes citation data in the link prediction task on the co-
authorship network. In particular, we explore the capabilities of a
recommender system based on data aggregation strategies on different graphs.
Since graph neural networks proved their efficiency on a wide range of tasks
related to recommendation systems, we leverage them as a relevant method for
the forecasting of potential collaborations in the scientific community.

###### Keywords:

Graph machine learning Neural networks Recommender systems Social graphs

##  1 Introduction

Since the advent of scientific communities, there has been a high demand in
the area of collaboration recommendations. Due to the complex nature of
interconnections between researchers, this domain has not reached a successful
automation for a long time.

According to the underlying graph structure of collaboration networks, we
propose to use recently emerged graph neural networks (GNN) to efficiently
predict research cooperation between scientists. This branch of machine
learning has readily proved its outstanding performance in a wide range of
areas related to recommender systems [15]. Such algorithms as Node2Vec [4],
Attri2Vec [16], and GraphSAGE [6] can be trained to capture structural
features of co-authorship network. Embeddings produced by these methods can be
effectively applied to the different forecasting tasks including prediction of
network connections as well [8].

Graph neural networks allow us not only to boost performance in
straightforward link prediction task on co-authorship network, but to improve
the quality on such a forecasting challenge via aggregation of additional
information from the citation graph. Future development of the discussed
pipeline can lead to the simplification of the collaboration assessment
process for the R$\&$D team management.

##  2 Related work

Recommender systems for scientific communities have a long history of
development. Early approaches in this area [13, 1] were based on deterministic
network information, which lacked the ability to represent complex features of
graph data.

Learning algorithms in the area of link prediction resolved a variety of
issues related to capture of graph intricacies. The explicit examples of such
models are local random walk [8] and local naive Bayes [9].

Various metrics such as content similarity LDAcosin [3] were also applied in
this field. Another popular method [11] leverages linear regression on feature
vectors of nodes and set of graph measures [10]. However, implementation of a
learning feature extractor instead of deterministic metrics would
significantly boost the performance of such pipelines. Applying this idea to
the collaboration network structures, usage of graph neural networks becomes
native [17].

Despite the promising results achieved by graph neural networks, different
techniques could be implemented to further increase their efficiency. Those
methods include the alteration of graph topology [14] and the usage of task-
independent techniques like Node2Vec in order to create initial node
representations which serve as better inputs for GNNs [5].

##  3 Data

We use classic HEP-TH dataset [7] as the subject for further development and
extension. It consists of citation and co-authorship graphs obtained from the
arXiv papers published between January 1993 and April 2003 in the area of high
energy physics theory. Unfortunately, there is no connection between these two
parts of the dataset (authors’ IDs were not provided in the citation network)
which makes its initial second part worthless for our purposes.

To our best knowledge, the potential of the HEP-TH citation network was never
explicitly revealed. Due to the presence of highly useful paper metadata (such
as author lists or abstracts), the range of the dataset usage can be
significantly extended. In order to complete our current research, we perform
processing of the citation graph metadata aiming to restore the corresponding
co-authorship network. To preserve homogeneous nature of the reconstructing
graph, we discarded from the citation graph all anonymous papers.

Along with the information about authors and abstracts, HEP-TH involves a
journal reference field with the publisher output information. The presence of
such data provides us an opportunity to extract ISSNs of indexed publications
and further parse scientific metrics (quartile, h-index, and impact factor)
from the ”SCImago Journal & Country Rank” website.

##  4 Methods

Explored architecture consists of two subsequently applied graph neural
networks. The first model generates vector representations of the publications
according to their annotations and the structure of the citation network. As
its input we leverage the directed citation graph $G(V,A,X)$, where
$V=\\{v_{1},v_{2},\dots,v_{n}\\}$ corresponds to the set of graph nodes
(articles), $A:n\times n\rightarrow\\{0,1\\}$ is the adjacency matrix (each
edge encodes citation between two papers), and $X:n\times
d\rightarrow\textbf{R}$ denotes the matrix of node features (vectorized
abstracts via pre-trained FastText [2]).

We perform a set of computational experiments using previously addressed
unsupervised methods GraphSAGE, Node2Vec, and Attri2Vec. In the following, we
briefly discuss each of them to clarify their usage as the feature extractors.

GraphSAGE. This method aggregates information about the set of neighbor nodes
$N(v)$ in order to produce embedding of node $v$

$h^{l+1}_{v}=\sigma(W^{l}\cdot\operatorname{CONCAT}(h^{l}_{v},h^{l+1}_{N(v)})),$
(1)

where $l+1$ is the current number of a convolutional layer, $h^{0}_{v}$ =
$x_{v}$, $W^{l}$ is the matrix of learning parameters, and $h^{l+1}_{N(v)}$
can be extracted by different aggregation functions like Max Pooling or Mean
aggregator.

Node2Vec. This algorithm generates sequences of nodes via second-order random
walks and utilizes them as the input data for a skip-gram model. The skip-gram
generates pairs from input and context nodes in order to cast them to the
feedforward neural network. Its weights can be used as the desired node
embeddings as a result of the following function optimization:

$\max_{f}\sum_{u\in V}\left[-\log Z_{u}+\sum_{n_{i}\in
N(u)}f\left(n_{i}\right)\cdot f(u)\right],$ (2)

where $Z_{u}=\sum_{v\in V}\exp(f(u)\cdot f(v))$ is per-node partition function
and f($\cdot$) corresponds to the mapping function.

Attri2Vec. The last considered model uses the image $f(x_{i})$ of node $v_{i}$
with feature vector $x_{i}$ in the new attribute subspace to predict its
context nodes. In this method, the task is to solve the joint optimization
problem

$\min_{W^{in},W^{\text{out
}}}-\sum_{i=1}^{|V|}\sum_{j=1}^{|V|}n\left(v_{i},v_{j}\right)\log\frac{\exp\left(f\left(v_{i}\right)\cdot
w_{j}^{\text{out }}\right)}{\sum_{k=1}^{|V|}\exp\left(f\left(v_{i}\right)\cdot
w_{k}^{\text{out }}\right)},$ (3)

where $n(v_{i},v_{j})$ is the number of times that $v_{j}$ occurs in $v_{i}$
context within t-window size in the generated set of random walks,
$W_{\text{in}}$ is the weight matrix from the input layer to hidden layer and
$W_{\text{out}}$ is the weight matrix from the hidden layer to the output
layer.

Link prediction step follows after the publications embeddings generating. For
this task, we consider the co-authorship graph
$\hat{G}(\hat{V},\hat{A},\hat{X})$, where
$\hat{V}=\\{\hat{v}_{1},\hat{v}_{2},\dots,\hat{v}_{m}\\}$ is the set of graph
vertexes (authors), $\hat{A}:m\times m\rightarrow\\{0,1\\}$ is the adjacency
matrix (each edge encodes collaboration between the two authors), and
$\hat{X}:m\times k\rightarrow\textbf{R}$ denotes the matrix of node features
(one-hot encoded research interests of authors). In order to supply the
predictive model by additional data about publications of the authors, we
extend each element $\hat{x}_{i}$ of the matrix $\hat{X}$ to
$\hat{x}^{\prime}_{i}$ as follows:

$\hat{x}^{\prime}_{i}=\operatorname{CONCAT}(x_{i},\sum_{e\in e_{i}}e),$ (4)

where $e_{i}$ denotes the set of publications embeddings of ${\rm i}^{th}$
author.

After the concatenation, the extended graph was translated as an input to the
two-layer GraphSAGE with link classifier. It constructs the embedding of the
potential links applying a binary operator to the pair of node embeddings (we
consider L1, L2, Hadamard operator, average, and inner product [12]). Finally,
these link embeddings are passed through the dense classification layer to
obtain probabilities of links existence in the network. The whole pipeline is
illustrated in Figure 1.

Figure 1: Example configuration of the described architecture based on
GraphSAGE convolutions.

The last model is trained by minimizing the binary cross-entropy loss

$L=-\frac{1}{\text{ N }}\sum_{i=1}^{\text{N
}}y_{i}\cdot\log\hat{y}_{i}+\left(1-y_{i}\right)\cdot\log\left(1-\hat{y}_{i}\right),$
(5)

where $N$ is the output size, $y_{i}$ is the true link labels, $\hat{y}_{i}$
is the predicted link existence probabilities.

##  5 Results

We performed training and evaluation on the PC with 1 GPU Tesla V100 and 96 GB
of RAM. The weights of the neural networks were updated by Adam optimizer. We
divided graph edges into train, validation and test samples in ratio 3:1:2.

We evaluated the set of models with various link embedding operators,
representation learning models, and aggregation functions. As the main quality
measurements, we chose binary accuracy, AUC-ROC, and F1-score. Received
embeddings from the models were tested on the supervised GraphSAGE setup
applied to the link prediction task.

Article embedding | Author embedding | LP op. | Accuracy | AUC-ROC | F1-score
---|---|---|---|---|---
– | GraphSAGE (Mean) | L2 | 0.8793 | 0.9442 | 0.8817
FastText | GraphSAGE (Mean) | Had | 0.8844 | 0.9486 | 0.8828
GraphSAGE (Mean) | GraphSAGE (Mean) | Had | 0.8895 | 0.9568 | 0.8911
GraphSAGE (Mean) | GraphSAGE (Mean) | L2 | 0.8928 | 0.9531 | 0.8885
GraphSAGE (Mean) | GraphSAGE (MaxPool) | L1 | 0.8638 | 0.9617 | 0.8489
Table 1: Results of different models on test sample

We conducted experiments with more than 80 different configurations and
represented key results in Table 1. The first two baselines were selected as
the best among the approaches using either vectors of authors interests
without citation network information or just embeddings of the abstracts.
Comparison of these simpler architectures with full models which includes two
graph neural networks could be interpreted as an ablation study.

As it is shown in the table, the proposed aggregation algorithm positively
influences the quality of scientific collaboration forecasting. The model
without any citation data significantly suffers from the lack of expressive
input features as well as the pipeline including only the abstracts of the
papers. Obtained result allow us to report that proposed aggregation strategy
efficiently utilizes citation graph properties.

##  6 Conclusion and Outlook

In the present paper, we briefly introduced the two-stage pipeline for the
collaboration prediction task. Performed computational experiments reveal the
perspective of citation data utilization in sense of co-authorship network
extension. The embeddings generated by GNNs effectively capture the network
properties including its topology and vectorized abstracts represented as
features of the corresponding graph nodes.

Our main contributions in the present work are the following:

  1. 1.

We perform extraction of the co-authorship graph from the corresponding HEP-TH
citation network.

  2. 2.

Presence of useful metadata allows us to parse the scientific significance
measures of the publications (e.g., impact factor).

  3. 3.

We aggregate structural data from the citation graph and apply it to the co-
authorship network in order to evaluate its influence on the link prediction
quality.

Along with the future improvements of link prediction methods in the area of
scientific collaboration, we intend to explore qualitative and quantitative
assessment approaches of emerged links. As the probabilistic estimation of
collaborations is not sufficient, it is important to extend it by less
abstract metrics. In the following, we are going to leverage the average
impact factor and the total number of publications for this task.

## Acknowledgements

We acknowledge fruitful discussions with Natalia Semenova.

## References

  * [1] Alinani, K., Wang, G., Alinani, A., Narejo, D.H.: Who should be my co-author? recommender system to suggest a list of collaborators. In: 2017 IEEE International Symposium on Parallel and Distributed Processing with Applications and 2017 IEEE International Conference on Ubiquitous Computing and Communications (ISPA/IUCC). pp. 1427–1433. IEEE (2017)
  * [2] Bojanowski, P., Grave, E., Joulin, A., Mikolov, T.: Enriching word vectors with subword information. Transactions of the Association for Computational Linguistics 5, 135–146 (2017)
  * [3] Chuan, P.M., Ali, M., Khang, T.D., Dey, N., et al.: Link prediction in co-authorship networks based on hybrid content similarity metric. Applied Intelligence 48(8), 2470–2486 (2018)
  * [4] Grover, A., Leskovec, J.: node2vec: Scalable feature learning for networks. In: Proceedings of the 22nd ACM SIGKDD international conference on Knowledge discovery and data mining. pp. 855–864 (2016)
  * [5] Gupta, C., Jain, Y., De, A., Chakrabarti, S.: Integrating transductive and inductive embeddings improves link prediction accuracy. arXiv preprint arXiv:2108.10108 (2021)
  * [6] Hamilton, W.L., Ying, R., Leskovec, J.: Inductive representation learning on large graphs. In: Proceedings of the 31st International Conference on Neural Information Processing Systems. pp. 1025–1035 (2017)
  * [7] Leskovec, J., Kleinberg, J., Faloutsos, C.: Graph evolution: Densification and shrinking diameters. ACM Trans Knowledge Discov Data 1 (04 2006)
  * [8] Liu, W., Lü, L.: Link prediction based on local random walk. EPL (Europhysics Letters) 89(5), 58007 (2010)
  * [9] Liu, Z., Zhang, Q.M., Lü, L., Zhou, T.: Link prediction in complex networks: A local naïve bayes model. EPL (Europhysics Letters) 96(4), 48007 (2011)
  * [10] Makarov, I., Bulanov, O., Gerasimova, O., Meshcheryakova, N., Karpov, I., Zhukov, L.E.: Scientific matchmaker: collaborator recommender system. In: International Conference on Analysis of Images, Social Networks and Texts. pp. 404–410. Springer (2017)
  * [11] Makarov, I., Bulanov, O., Zhukov, L.E.: Co-author recommender system. In: international conference on network analysis. pp. 251–257. Springer (2016)
  * [12] Makarov, I., Gerasimova, O., Sulimov, P., Korovina, K., Zhukov, L.E.: Joint node-edge network embedding for link prediction. In: van der Aalst, W.M.P., Batagelj, V., Glavaš, G., Ignatov, D.I., Khachay, M., Kuznetsov, S.O., Koltsova, O., Lomazova, I.A., Loukachevitch, N., Napoli, A., Panchenko, A., Pardalos, P.M., Pelillo, M., Savchenko, A.V. (eds.) Analysis of Images, Social Networks and Texts. pp. 20–31. Springer International Publishing, Cham (2018)
  * [13] Sie, R.L., Drachsler, H., Bitter-Rijpkema, M., Sloep, P.: To whom and why should i connect? co-author recommendation based on powerful and similar peers. International Journal of Technology Enhanced Learning 4(1-2), 121–137 (2012)
  * [14] Singh, A., Huang, Q., Huang, S.L., Bhalerao, O., He, H., Lim, S.N., Benson, A.R.: Edge proposal sets for link prediction. arXiv preprint arXiv:2106.15810 (2021)
  * [15] Ying, R., He, R., Chen, K., Eksombatchai, P., Hamilton, W.L., Leskovec, J.: Graph convolutional neural networks for web-scale recommender systems. In: Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining. p. 974–983. KDD ’18, Association for Computing Machinery, New York, NY, USA (2018). https://doi.org/10.1145/3219819.3219890
  * [16] Zhang, D., Yin, J., Zhu, X., Zhang, C.: Attributed network embedding via subspace discovery. Data Mining and Knowledge Discovery 33(6), 1953–1980 (Nov 2019)
  * [17] Zhang, M., Chen, Y.: Link prediction based on graph neural networks. Advances in Neural Information Processing Systems 31, 5165–5175 (2018)