# Consistency between black hole and mimetic gravity
– Case of $(2+1)$-dimensional gravity –

Shin’ichi Nojiri<EMAIL_ADDRESS>Department of Physics,
Nagoya University, Nagoya 464-8602, Japan
&
Kobayashi-Maskawa Institute for the Origin of Particles and the Universe,
Nagoya University, Nagoya 464-8602, Japan  G. G. L. Nashed<EMAIL_ADDRESS>Centre for Theoretical Physics, The British University, P.O. Box 43, El
Sherouk City, Cairo 11837, Egypt

###### Abstract

We show that the mimetic theory with the constraint
$g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=1$ cannot realize the
black hole geometry with the horizon(s). To overcome such issue, we may change
the mimetic constraint a little bit by
$\omega(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1,$ where
$\omega(\phi)$ is a function of the scalar field $\phi$. As an example, we
consider $(2+1)$-dimensional mimetic gravity with the mimetic potential and
construct black hole (BH) solutions by using this modified constraint. We
study three different classes: In the first class, we assume the Lagrange
multiplier and mimetic potential are vanishing and obtain a BH solution that
fully matches the BH of GR despite the non-triviality of the mimetic field
which ensures the study presented in JCAP 01 (2019) 058. In the second class,
we obtain a BH having constant mimetic potential and a non-trivial form of the
Lagrange multiplier. In the third class, we obtain a new BH solution with non-
vanishing values of the mimetic field, the Lagrange multiplier, and the
mimetic potential. In any case, the solutions correspond to the space-time
with only one horizon but we show that the formalism for the constraint works.

Mimetic gravitational theory, analytic (2+1)-dimension BHs.

###### pacs:

04.50.Kd, 04.25.Nx, 04.40.Nr

##  I Introduction

The study of a black hole (BH) is considered the most ravishing and
astonishing astrophysical topic which appeared as the solution of Einstein’s
general relativity (GR). These topics yield a robust background for
discovering different subjects of physics comprising thermodynamics, quantum
gravity, paramagnetism-ferromagnetism phase transition, superconducting phase
transition, super-fluids, spectroscopy, condensed matter physics, holographic
hypothesis, information theory, and so on. The discovery of the well-known
confirmed correspondence between gravitational field in $N$-dimensional anti-
de Sitter (AdS) space-time and the conformal field theory (CFT), existing on
the boundary of $(N-1)$-dimensional space-time, which is known as
$\mathrm{AdS}_{N}/\mathrm{CFT}_{N-1}$ correspondence or gauge/gravity duality,
makes the study of BH solutions be a hot topic. The understanding of the
puzzle of information paradox that was studied in detail in Hawking  _et al._
(2016) and the shadow of the super-massive BHs whose results given by $M87$
event horizon telescope Akiyama  _et al._ (2019) yield researchers to do more
effort to understand the nature of BHs. The existence of soft gravitons and/or
soft photons, on the BH horizon and their quantum state information which is
stored on a holographic plate at the future boundary of the horizon, are
discussed in Hawking  _et al._ (2016). The entropy and microscopic
construction of BH near the horizon are discussed through the soft hairs in
Afshar  _et al._ (2016); Haco  _et al._ (2018, 2019); Grumiller  _et al._
(2020); Sheykhi (2021).

This paper is concerned with the study of the mimetic gravitational theory
coupled with mimetic potential and the Lagrange multiplier. The mimetic theory
in $(3+1)$-dimensions is a new gravitational theory constructed recently as a
new prescription of the dark matter (DM) issue Chamseddine and Mukhanov (2013,
2017a). It is shown that the mimetic theory can have an additional
longitudinal degree of freedom to the gravitational field although the degree
of freedom is not dynamical. The additional longitudinal degree of freedom is
responsible for the mimetic DM in the case of $(3+1)$-dimensions. Recently, it
is discussed that an amendment of mimetic gravity in $(3+1)$-dimensions can
solve the issue of the cosmological singularities Chamseddine and Mukhanov
(2017a) and the singularity at the center of a BH Chamseddine and Mukhanov
(2017b). Moreover, it is shown that the mimetic gravitational theory in
$(3+1)$-dimensions forecasts that gravitational wave (GW) can spread with a
speed of light, confirming fully consistent with the results of the event,
GW170817, and its optical counterpart Casalino  _et al._ (2018, 2019).
Moreover, it is shown that mimetic theory can investigate the flat rotation
curves of spiral galaxies Vagnozzi (2017); Sheykhi and Grunau (2019).
Recently, a mimetic theory has gained a lot of attention in the domain of
cosmology Chamseddine  _et al._ (2014); Dutta  _et al._ (2018); Abbassi  _et
al._ (2018); Zhong and Sáez-Chillón Gómez (2018); Matsumoto (2016); Nojiri and
Odintsov (2014); Odintsov and Oikonomou (2016a); Nojiri  _et al._ (2016);
Sadeghnezhad and Nozari (2017); Nashed (2021); Gorji  _et al._ (2019, 2018a);
Bouhmadi-Lopez  _et al._ (2017); Gorji  _et al._ (2018b); Firouzjahi  _et al._
(2018); Sebastiani  _et al._ (2017); Chamseddine  _et al._ (2019) and in the
domain of astrophysics Deruelle and Rua (2014); Myrzakulov and Sebastiani
(2015); Myrzakulov  _et al._ (2016); Astashenok  _et al._ (2015); Odintsov and
Oikonomou (2015); Nojiri  _et al._ (2017); Odintsov and Oikonomou (2018);
Oikonomou (2016a); Nashed (2014); Gorji  _et al._ (2020); Nashed  _et al._
(2019); Chen  _et al._ (2018); Nashed and Nojiri (2021a, b, c); Nashed (2018,
2018); Ben Achour  _et al._ (2018); Zheng  _et al._ (2017); Elizalde  _et al._
(2020); Shen  _et al._ (2019); Sheykhi (2020).

In the mimetic gravitational theory, the conformal degree of freedom of
gravitational field can be isolated by introducing a relation between the
auxiliary metric $\bar{g}^{\alpha\beta}$, the physical metric
$g_{\alpha\beta}$, and a mimetic field $\phi$ which is a scalar through the
form,

$\displaystyle
g_{\alpha\beta}=\mp\left(\bar{g}^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi\right)\bar{g}_{\alpha\beta}\,.$
(1)

Eq. (1) has the scale invariance
$\bar{g}_{\mu\nu}\to\mathrm{e}^{\sigma}\bar{g}_{\mu\nu}$ with a parameter
$\sigma$. Equation (1) yields that the mimetic scalar field should give,

$\displaystyle g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi=\mp
1\,.$ (2)

In this paper, we choose the signature of the metric as $(+,-,-,-)$. Therefore
if we choose the negative sign in the r.h.s. of Eq. (2),
$\left(\partial_{\alpha}\phi\right)$ is a space-like vector and if we choose
the positive sign, $\left(\partial_{\alpha}\phi\right)$ is time-like. When the
cosmology is discussed, the positive sign has been often used and the mimetic
matter can be interpreted as dust but when we consider the static objects like
black holes, the negative sign has been used well. Recently, many amendments
of mimetic theory have been constructed like, mimetic gravity with the
Lagrange multipliers Cid and Labrana (2012); Capozziello  _et al._ (2010) and
$f(R)$ mimetic gravity. Additionally, the existence of mimetic potential and
the Lagrange multiplier assisted the possibility for understanding different
topics in the frame of cosmology, Odintsov and Oikonomou (2015). The mimetic
theory of gravity is also extended to $f(R)$ mimetic theory of gravity Nojiri
and Odintsov (2014); Odintsov and Oikonomou (2016a); Oikonomou (2016b); Nashed
(2010); Oikonomou (2016c, a); Nashed (2008); Myrzakulov and Sebastiani (2016);
Odintsov and Oikonomou (2015, 2016b, 2016c); Nojiri  _et al._ (2017); Odintsov
and Oikonomou (2018); Chen  _et al._ (2021) and the Gauss-Bonnet mimetic
gravity Astashenok  _et al._ (2015); Oikonomou (2015); Nashed (2007); Zhong
and Elizalde (2016); Zhong and Sáez-Chillón Gómez (2018). More specific, a
unified formulation of early inflation and late-time acceleration in the
$f(R)$ mimetic gravity was formulated in Nojiri  _et al._ (2016) where
confirmation of the inflationary era is established contrary to the $f(R)$
gravity.

In this paper, we show that the mimetic constraint (2) is not consistent with
the black hole geometry with the horizon(s). In order to solve this problem,
we modify the constraint (2) by introducing a function $\omega(\phi)$ as
follows,

$\displaystyle\omega(\phi)g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi=\mp
1\,,$ (3)

which is locally equivalent to (2) as we will see in the next section but
$\omega(\phi)$ plays the important role when crossing the horizon. We should
also note that by introducing $\omega(\phi)$, we can treat the case that the
time-like vector $\left(\partial_{\alpha}\phi\right)$ changes to the space-
like vector and vice versa in a solution even if we fix the negative or
positive sign in the r.h.s. of Eq. (3). We apply this formalism to
$(2+1)$-dimensional mimetic gravity and show that the formalism works well. We
should also note that by the constraint (3), Eq. (1) can be rewritten as
$g_{\alpha\beta}=\mp\omega(\phi)\left(\bar{g}^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi\right)\bar{g}_{\alpha\beta}$
and it is clear that there remains the scale invariance
$\bar{g}_{\mu\nu}\to\mathrm{e}^{\sigma}\bar{g}_{\mu\nu}$.

The derivation of $(2+1)$-dimensional BH solutions of GR in AdS space-time,
which is known as the BTZ (Banados-Teitelboim-Zanelli) BHs Banados  _et al._
(1992), is considered as the most fascinating achievement in BH physics. The
metric of the BTZ black hole is given by

$\displaystyle
ds^{2}=-f(r)dt^{2}+\frac{1}{f(r)}dr^{2}+r^{2}\left(h(r)\,dt+d\theta\right)^{2}\,,\quad\textrm{where}\quad
f(r)=\Lambda r^{2}-M+\frac{J}{4r^{2}}\,,\textrm{and}\quad
h(r)=-\frac{J}{2r^{2}}\,.$ (4)

Here, $\Lambda$, $M$, and $J$ are the cosmological constant and the black
hole’s mass and angular momentum. The $(2+1)$-dimensional BH solution of the
Einstein GR gravitational theory yields an easy model to explain and discuss
some basic problems as quantum gravity, BH thermodynamics, holographic
superconductors, string and gauge/gravity duality, in the frame of
$\mathrm{AdS}_{3}/\mathrm{CFT}_{2}$ Carlip (1995); Ashtekar  _et al._ (2003);
Sarkar  _et al._ (2006); Witten (1998); Carlip (2005). The coincidence of the
quasinormal modes in the $\mathrm{AdS}_{3}/\mathrm{CFT}_{2}$ with the poles of
the correlation function in the dual CFT gives quantitative proof for
$\mathrm{AdS}_{3}/\mathrm{CFT}_{2}$ Birmingham  _et al._ (2002). Moreover, the
BTZ BHs are important tools to increase our knowledge of gravitational
interaction in lower-dimensional space-times Witten (2007), especially, in the
quantum gravity in $(2+1)$-dimension. The geometries of the spinning
$(2+1)$-dimensional BHs have been discussed in Banados  _et al._ (1993). The
study of BTZ that has an electric charge, mass, and angular momentum has been
analyzed in Martinez  _et al._ (2000); Fernando and Mansouri (1997).
Furthermore, the $(2+1)$-dimensional BHs can supply a strong background to
investigate one-dimensional holographic superconductors Ren (2010); Liu  _et
al._ (2011); Kord Zangeneh  _et al._ (2017); Mohammadi  _et al._ (2018);
Ghotbabadi  _et al._ (2018). There are many studies with different topics on
$(2+1)$-dimensional BH solutions of gravitational field done in the
literatures (see e.g. Panotopoulos and Rincón (2018); Rincón and Panotopoulos
(2018); Clement (1996); Carlip (1995); Nojiri and Odintsov (1998); Emparan
_et al._ (2000); Cadoni and Setare (2008); Parsons and Ross (2009); Hendi
(2012); Sheykhi  _et al._ (2014); Xu (2014); Grumiller and Merbis (2020);
Capozziello  _et al._ (2010); Sebastiani  _et al._ (2017)).

The arrangement of the present study is as follows: In Section II, we show
that the mimetic constraint (2) is inconsistent with the black hole geometry
with the horizon(s) and give a solution for this problem as given in (3). In
Section III, we apply the above formalism to the $(2+1)$-dimensional mimetic
gravity coupled with mimetic potential and the Lagrange multiplier. In this
section, we employ the field equations of mimetic theory to a
$(2+1)$-dimensional space-time that has unequal metric potentials and presents
the non-linear differential equations which have five unknown functions, i.e.,
two of the metric potentials, mimetic scalar field, mimetic potential field,
and the Lagrange multiplier. We solved this system of non-linear differential
equations using three restrictions of the mimetic potential and the Lagrange
multiplier and obtain three different analytic BH solutions. We summarize and
discuss the results obtained in this study in the final Section IV.

##  II Inconsistency between mimetic gravity with black hole geometry

The term “mimetic DM” in $(3+1)$-dimensional space-time is delivered in the
scientific society by Mukhanov and Chamseddine Chamseddine  _et al._ (2014).

In this paper, we consider the black hole space-time but the mimetic
constraint (2) is not consistent with the black hole space-time with
horizon(s). In (2), we choose the minus sign as follows,

$\displaystyle g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1\,.$
(5)

Even if we choose the plus sign in (2), the situation is not changed. We
consider the $D$ dimensional sphericaly symmetric and static space-time with
following line element,

$\displaystyle
ds^{2}=b(r)dt^{2}-\frac{dr^{2}}{b_{1}(r)}-r^{2}d\Omega_{D-2}^{2}\,.$ (6)

Here $d\Omega_{D-2}^{2}$ is the line element of $D-2$ dimensional unit sphere.
If we may also assume $\phi=\phi(r)$, the mimetic constraint (5) has the
following form

$\displaystyle b_{1}\left(\phi^{\prime}\right)^{2}=1\,.$ (7)

The equation has no solution if $b_{1}$ is negative, $b_{1}<0$. In the case of
black hole geometry, $b$ vanishes and changes its signature at the horizon. As
we will see below, to avoid the curvature singularity, $b_{1}$ also must
vanish at the horizon. Therefore even if $b_{1}$ is positive, $b_{1}$ is
negative in general, as in the Schwarzschild black hole or the BTZ black hole,
which conflicts with Eq. (7). This tells that the mimetic theory with the
constraint (5) cannot realize the black hole geometry with the horizon(s).

In order to avoid the above problem, we may change the mimetic constraint in
(5), a little bit as in (3),

$\displaystyle\omega(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1\,.$
(8)

Here $\omega(\phi)$ is a function of $\phi$. If $\omega(\phi)$ is positive, we
may define a new scalar field $\tilde{\phi}$ by $\tilde{\phi}=\int
d\phi\sqrt{\omega(\phi)}$, the constraint (8) is reduced to the form of (5),

$\displaystyle
g^{\rho\sigma}\partial_{\rho}\tilde{\phi}\partial_{\sigma}\tilde{\phi}=-1\,.$
(9)

The signature of $\omega(\phi)$ can be, however, changed in general. If we may
also assume $\phi=\phi(r)$ and the space-time is given by (6), instead of (7),
the constraint (8) has the following form,

$\displaystyle b_{1}\omega(\phi)\left(\phi^{\prime}\right)^{2}=1\,.$ (10)

Then for a solution of $\phi$ where $\omega(\phi)$ is positive when $b_{1}$ is
positive and $\omega(\phi)$ is negative when $b_{1}$ is negative, the
constraint (10) is consistent even inside the horizon. As an example, we
consider a simple case where

$\displaystyle\omega(\phi)=\frac{1}{\phi}\,.$ (11)

Near the horizon, $b_{1}$ should behaves as

$\displaystyle b_{1}(r)\sim b_{0}\left(r-r_{\mathrm{h}}\right)\,.$ (12)

Here $r_{\mathrm{h}}$ is the radius of the horizon and $b_{0}$ is a positive
constant. Then a solution of (10) with (11) is given by

$\displaystyle\phi\sim\frac{r-r_{\mathrm{h}}}{b_{0}}\,.$ (13)

Then the scalar field $\phi$ and therefore $\omega(\phi)$ change the sign at
the horizon and Eq. (10) is consistent even inside the horizon.

In case there are several horizons, the problem might not be solved only by
the choice in (11). A way to solve the problem in this case, we may choose

$\displaystyle\omega(\phi)=b_{1}(\phi)\,.$ (14)

In this case, the solution of (10) is simply given by

$\displaystyle\phi=r\,.$ (15)

Therefore it is clear that the problem is solved by the choice of (14).
However, this choice may look rather artificial because it looks like we have
assumed the solution from the beginning. Anyway, the possibility of the choice
(14) tells that a model is realizing the solution of Eq. (10).

##  III $(2+1)$-dimensional BH solutions in mimetic gravity coupled with
mimetic potential and Lagrange multiplier

In this section, we apply the formalism proposed in the last Section II to
$(2+1)$-dimensional mimetic gravity and show that the formalism works well.
The action of the mimetic gravitational theory that has the Lagrange
multiplier $\lambda(\phi)$ and mimetic potential $V(\phi)$ has the following
form Chamseddine  _et al._ (2014); Nojiri and Odintsov (2014); Gorji  _et al._
(2020),

$\displaystyle
S=\frac{1}{2\kappa^{2}}\int\mathrm{d}x^{3}\sqrt{-g}\left\\{R\left(g_{\mu\nu}\right)-2\Lambda-V(\phi)+\lambda\left(\omega(\phi)g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi+1\right)\right\\}+S_{\mathrm{matt}}\,,$
(16)

where $\kappa^{2}$ is the Einstein gravitational constant which in the
relativistic units equals $\kappa^{2}=8\pi$ and $S_{\mathrm{matt}}$ is the
action of the matter fields. Furthermore, we assume the auxiliary scalar
mimetic field, $\phi$, relies on the radial coordinate $r$. The variation of
the action (16) w.r.t. the metric tensor $g_{\mu\nu}$, yield the field
equations,

$\displaystyle
0=R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}\left(R-2\Lambda\right)+\frac{1}{2}g_{\mu\nu}\left\\{\lambda\left(\omega(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi+1\right)-V(\phi)\right\\}-\lambda\omega(\phi)\partial_{\mu}\phi\partial_{\nu}\phi+8\pi
T_{\mu\nu}\,,$ (17)

with $T_{\mu\nu}$ being the energy-momentum tensor of the matter fluids. In
the present study, we only consider the case $T_{\mu\nu}=0$. Moreover, the
variation of the action w.r.t. the auxiliary mimetic scalar field $\phi$,
gives

$\displaystyle
2\nabla^{\mu}(\lambda\omega(\phi)\partial_{\mu}\phi)-\lambda\omega^{\prime}(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi+V^{\prime}(\phi)=0\,,$
(18)

where the “prime” means the differentiation w.r.t. the auxiliary scalar
$\phi$. Finally, the variation of the action (16) w.r.t. the Lagrange
multiplier $\lambda$, gives

$\displaystyle\omega(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1\,,$
(19)

which coincides with Eq. (8). We should note, the scalar field equation (18)
can be obtained from (17) and (19) when $T_{\mu\nu}=0$ or by using the
conservation law. Therefore we do not use (18) hereafter.

Now, we are ready to apply the field equations of the mimetic theory, (17),
and (18) to a $(2+1)$-dimensional space-time that has two unknown functions
and takes the form in the coordinates $(t,r,\theta)$ as,

$\displaystyle
ds^{2}=b(r)dt^{2}-\frac{dr^{2}}{b_{1}(r)}-r^{2}\,d\theta^{2}\,.$ (20)

Here $b(r)$ and $b_{1}(r)$ are unknown functions of radial coordinate $r$,
that will fixed through the use of the field equations (17), and (18). Using
Eq. (20), we get the Ricci scalar as,

$\displaystyle
R={\frac{b^{\prime}b^{\prime}_{1}br+2b_{1}b^{\prime\prime}br-b^{\prime
2}b_{1}r+2b_{1}b^{\prime}b+2b^{\prime}_{1}b^{2}}{2b^{2}r}}\,.$ (21)

The expression (21) tells that the scalar curvature $R$ might diverge when
$b=0$ because there is $b^{2}$ in the denominator. If $R$ diverges, the radius
where $b$ vanishes does not correspond to the horizon but the singularity. To
check if the radius with $b=0$ corresponds to the horizon or the singularity,
we assume, as in the standard Schwarzschild space-time $b_{1}$ is written in
the form $b_{1}=b_{2}b$, where $b_{2}$ is finite when $b=0$. Then the scalar
curvature $R$ in (21) is rewritten as

$\displaystyle
R=\frac{b^{\prime}b^{\prime}_{2}r+2b_{2}b^{\prime\prime}r+2b^{\prime}_{2}b+4b_{2}b^{\prime}}{2r}\,.$
(22)

Therefore if $b_{2}$ vanishes when $b$ vanishes, the scalar curvature $R$ is
finite and the radius where $b$ vanishes corresponds to the horizon.

Applying the field equations (17) and (18) to space-time (20), after using Eq.
(21), we obtain the following non-linear differential equations:

The $(t,t)$-component of the field equation (17) is

$\displaystyle 0=\frac{b^{\prime}_{1}+Vr-2r\Lambda}{2r}\,,$ (23)

the $(r,r)$-component is

$\displaystyle 0=\frac{b_{1}b^{\prime}-2\lambda bb_{1}r\omega\phi^{\prime
2}-2br\Lambda+V\,br}{2rb}\,,$ (24)

the $(\theta,\theta)$-component is

$\displaystyle
0=\frac{bb^{\prime}_{1}b^{\prime}+2b^{\prime\prime}bb_{1}-b^{\prime
2}b_{1}-4b^{2}\Lambda+2V\,b^{2}}{4b^{2}}\,,$ (25)

where $b\equiv b(r)$, $b_{1}\equiv b_{1}(r)$, $V\equiv V(r)$ 111Due to the
nature of our study in this research, we assume the mimetic potential and the
Lagrange multiplier to depend on the radial coordinate $r$.,
$\lambda\equiv\lambda(r)$, $b^{\prime}=\frac{db}{dr}$,
$b^{\prime}_{1}=\frac{db_{1}}{dr}$, $V^{\prime}=\frac{dV}{dr}$,
$\phi^{\prime}=\frac{d\phi}{dr}$, and $\lambda^{\prime}=\frac{d\lambda}{dr}$.
The above non-linear differential equations, i.e., (23), (24), and (25) are
going to be solved for the following three different classes:

Case I: $V(r)=\lambda(r)=0$.
When $V(r)=\lambda(r)=0$, the analytic solution of the non-linear differential
equations (23), (24), and (25) after resealing the constants of integration,
takes the following form,

$\displaystyle b(r)=b_{1}(r)=\Lambda r^{2}+C=\Lambda r^{2}-M\,,\quad C=-M\,.$
(26)

If we choose $\omega(\phi)$ as in (11), we obtain

$\displaystyle\phi(r)=\frac{1}{2}\left\\{\frac{1}{\sqrt{\Lambda}}\ln\left(\sqrt{\frac{\Lambda
r^{2}}{M}}+\sqrt{\frac{\Lambda r^{2}}{M}-1}\right)\right\\}^{2}\,.$ (27)

Because $b=b_{1}$, the arguments around Eq. (22) tell that the radius where
$b(r)=b_{1}(r)=0$, that is, $r=\sqrt{\frac{M}{\Lambda}}$ when
$\frac{M}{\Lambda}>0$ corresponds to the horizon radius. Inside the horizon,
we find $\sqrt{\frac{\Lambda r^{2}}{M}-1}=\pm i\sqrt{1-\frac{\Lambda
r^{2}}{M}}$. Because $\frac{\Lambda r^{2}}{M}<1$ inside the horizon and
$\left(\sqrt{\frac{\Lambda r^{2}}{M}}\right)^{2}+\left(\sqrt{1-\frac{\Lambda
r^{2}}{M}}\right)^{2}=1$, we may write $\sqrt{\frac{\Lambda
r^{2}}{M}}=\cos\theta$ and we find

$\displaystyle,\phi(r)=-\frac{\theta^{2}}{2\Lambda}\,,$ (28)

which is surely real and negative. The expression in (27) tells $\phi$ is
positive outside the horizon. Therefore due to the factor $\omega(\phi)$, we
can connect the region inside the horizon with the region outside the horizon
smoothly and $\phi$ vanishes on the horizon. Equation (26) shows that when the
mimetic potential and the Lagrange multiplier are vanishing, we return to the
well-known BTZ BH of GR solution despite the non-vanishing of the mimetic
field. This result is consistent with the previous studies Nashed  _et al._
(2019). It is of interest to note that the mimetic field of Eq. (27) satisfies
the constraints (10). The relevant physics of the BTZ space-time is well
studied in the literature. From now on, we will put the cosmological constant
equal to zero, i.e., $\Lambda=0$ to show the effect of the Lagrange multiplier
and the mimetic potential on the geometry under consideration.

Case II: $V(r)=\mathrm{const}=c_{1}$.
By choosing $V(r)=\mathrm{const}=c_{1}$ and using $\omega(\phi)$ in (11), we
find

$\displaystyle b_{1}=$
$\displaystyle\,c_{2}-\frac{c_{3}r^{2}}{2}\,,\quad\lambda(r)=\frac{c_{2}}{r\left(c_{2}c_{3}\sqrt{c_{1}r^{2}-2c_{2}}-r\right)}\,,$
$\displaystyle b(r)=$ $\displaystyle
c_{4}\mathrm{e}^{\int\frac{2\left(2c_{2}-c_{1}r^{2}+c_{1}c_{2}c_{3}r\sqrt{c_{1}r^{2}-2c_{2}}\right)}{\left(c_{2}c_{3}\sqrt{c_{1}r^{2}-2c_{2}}-r\right)\left(c_{1}r^{2}-2c_{2}\right)}dr}$
$\displaystyle=$
$\displaystyle-\left\\{c_{4}\left(r{c_{3}}^{2}{c_{2}}^{2}c_{1}-r-\sqrt{2}c_{3}{c_{2}}^{\frac{3}{2}}C_{0}\right)\left(2{c_{3}}^{2}{c_{2}}^{3}c_{1}+c_{1}r\sqrt{2}c_{3}{c_{2}}^{\frac{3}{2}}C_{0}-2c_{2}-\sqrt{2}C_{0}\sqrt{c_{2}}\sqrt{c_{1}r^{2}-2c_{2}}\right)\right.$
$\displaystyle\left.\times\left({c_{3}}^{2}{c_{2}}^{2}c_{1}r^{2}-r^{2}-2{c_{3}}^{2}{c_{2}}^{3}\right)\right\\}\left[\left(c_{1}r\sqrt{2}c_{3}{c_{2}}^{\frac{3}{2}}C_{0}-2{c_{3}}^{2}{c_{2}}^{3}c_{1}+2c_{2}+\sqrt{2}C_{0}\sqrt{c_{2}}\sqrt{c_{1}r^{2}-2c_{2}}\right)\right.$
$\displaystyle\left.\times\left(r{c_{3}}^{2}{c_{2}}^{2}c_{1}-r+\sqrt{2}c_{3}{c_{2}}^{\frac{3}{2}}C_{0}\right)\right]^{-1}\,,$
$\displaystyle\phi(r)=$
$\displaystyle\frac{1}{\omega(r)}=\frac{1}{c_{1}}\left\\{\tan^{-1}\left(\frac{\sqrt{c_{1}}r}{\sqrt{2c_{2}-c_{1}r^{2}}}\right)\right\\}^{2}\,,$
(29)

where we have put $C_{0}=\sqrt{{c_{3}}^{2}{c_{2}}^{2}c_{1}-1}$. As Eq. (III)
shows that the constant $c_{1}$ cannot have a zero value which means that it
cannot reduce to the first case. In order that the radius where $b_{1}=0$,
that is, $r=\sqrt{\frac{2c_{2}}{c_{3}}}$ when $\frac{2c_{2}}{c_{3}}>0$, is the
horizon radius, by the arguments around Eq. (22), $b$ should also vanish when
$b$ vanishes, which gives a constraint on the parameters, as follows,

$\displaystyle
0=\left(c_{3}{c_{2}}^{2}c_{1}-1-{c_{3}}^{\frac{1}{2}}c_{2}C_{0}\right)\left({c_{3}}^{2}{c_{2}}^{2}c_{1}+c_{1}{c_{3}}^{\frac{1}{2}}c_{2}C_{0}-1-C_{0}\sqrt{\frac{c_{1}}{c_{3}}-1}\right)\left({c_{3}}^{2}{c_{2}}^{2}c_{1}-1-{c_{3}}^{3}{c_{2}}^{2}\right)\,.$
(30)

Here we have assumed $c_{2}$ and $c_{4}$ do not vanish. Other radii where $b$
vanishes corresponding to the curvature singularity and they are not the
horizon radius. We may choose the parameters so that the singular surfaces do
not appear or the singular surfaces are inside of the black hole horizon. In
order that the black hole exists, the radius $r=\sqrt{\frac{2c_{2}}{c_{3}}}$
to be the black hole horizon, which require $c_{2}$ and $c_{3}$ are negative.
Anyway, this solution can have only one horizon as in the BTZ black hole with
the angular momentum $J=0$.

The expression of $\phi(r)$ is valid when $r<\sqrt{\frac{2c_{2}}{c_{3}}}$.
When $r>\sqrt{\frac{2c_{2}}{c_{3}}}$, because
$\tan^{-1}\left(\frac{\sqrt{c_{1}}r}{\sqrt{2c_{2}-c_{1}r^{2}}}\right)=\pm
i\tanh^{-1}\left(\frac{\sqrt{c_{1}}r}{\sqrt{c_{1}r^{2}-2c_{2}}}\right)$, we
find $\phi(r)$ becomes negative and real

$\displaystyle\phi(r)=\frac{1}{\omega(r)}=-\frac{1}{c_{1}}\left\\{\tanh^{-1}\left(\frac{\sqrt{c_{1}}r}{\sqrt{c_{1}r^{2}-2c_{2}}}\right)\right\\}^{2}\,.$
(31)

Therefore the formalism using $\omega(\phi)$ in (8) is valid even in Case II.

Case III: $\lambda(r)\neq 0$ and $V(r)\neq 0$.

The differential equations (23), (24), and (25) possess many analytic
solutions, however, in this study, we are going to consider only the physical
one which yields metric potential $b$ asymptotes as AdS/dS and also the
mimetic potential behave asymptotically for $r\to\infty$ as constant as it
should for any physical model. Now we are going to use the constraint (10)
then we have three differential equations, (23), (24), and (25), in four
unknowns $b(r)$, $b_{1}(r)$, $\lambda(r)$ and $V(r)$. Therefore, we have to
fix anyone of the previous unknowns by hand. Here we will use the form of
$V(r)$ as,

$\displaystyle V(r)=$
$\displaystyle\,c_{5}+\frac{c_{6}}{r}+\frac{c_{7}}{r^{3}}\,,$ (32)

where $c_{5},c_{6}$, and $c_{7}$ are constants. The merit of the form of Eq.
(32) is that when $r\to\infty$ gives constant potential as it should be from
the view point of physics. Using Eq. (32) in Eqs. (23), (24), and (25) we get
the solution of the remanning unknowns in the form,

$\displaystyle b(r)=$
$\displaystyle\,\frac{r^{2}}{2}\left\\{\int\frac{c_{9}dr}{\sqrt{r^{3}\left(2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}\right)}}+c_{10}\right\\}^{2}\,,$
$\displaystyle b_{1}(r)=$
$\displaystyle\frac{2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}}{2rc_{8}}\,,$
$\displaystyle\lambda(r)=$
$\displaystyle\,\frac{1}{2}\left\\{-c_{9}\mathrm{e}^{2\int\frac{\left(2+5c_{6}c_{8}r+3c_{5}c_{8}r^{2}\right)dr}{2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}}}\frac{\left(c_{6}c_{8}r^{2}-3c_{7}c_{8}+r\right)}{\sqrt{r^{3}\left(2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}\right)}}-c_{10}\left(-3c_{7}c_{8}+r+c_{6}c_{8}r^{2}\right)\right.$
$\displaystyle\left.-\frac{rc_{9}\left(-2c_{7}c_{8}+r+2c_{6}c_{8}r^{2}+c_{5}c_{8}r^{3}\right)}{\sqrt{r^{3}\left(2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}\right)}}\right\\}\left(\frac{1}{\sqrt{r^{3}\left(2c_{7}c_{8}-r-2c_{6}c_{8}r^{2}-c_{5}c_{8}r^{3}\right)}}+c_{10}\right)r^{-3}{c_{8}}^{-1}\,.$
(33)

As discussed in the paragraph including Equations (21) and (22), $b_{1}(r)$
must vanish when $b(r)$ vanishes. First, let us consider the case that
$b_{1}(r)$ in (III) has only one zero $b_{1}(r)=0$ at $r=r_{\mathrm{h}}$. We
can always adjust the constant of the integration $c_{10}$ so that $b(r)$ also
vanishes at $r=r_{\mathrm{h}}$, $b\left(r_{\mathrm{h}}\right)=0$. Then, in
this case, the radius $r_{\mathrm{h}}$ can be regarded as the horizon radius.
We may also consider the case that $b_{1}(r)$ has two zeros at
$r=r_{\mathrm{h}}$ and $r=r_{\mathrm{h}1}$ and we assume
$r_{\mathrm{h}}>r_{\mathrm{h}1}$. In this case, $b_{1}(r)$ can be rewritten as
follows,

$\displaystyle b_{1}(r)=$
$\displaystyle\frac{-c_{5}\left(r-r_{\mathrm{h}}\right)\left(r-r_{\mathrm{h}1}\right)\left(r+r_{0}\right)}{2r}\,.$
(34)

Here $r_{0}$ is a positive constant. We should note the constant $c_{5}$ must
be negative. As in the above case that $b_{1}(r)$ has only one zero, by
choosing $c_{10}$ so that $b(r)$ also vanishes at $r=r_{\mathrm{h}}$, we find

$\displaystyle
b(r)=\frac{r^{2}}{2}\left\\{\int_{r_{\mathrm{h}}}^{r}\frac{c_{9}dr}{\sqrt{-c_{5}r^{3}\left(r-r_{\mathrm{h}}\right)\left(r-r_{\mathrm{h}1}\right)\left(r+r_{0}\right)}}\right\\}^{2}=\frac{r^{2}}{4}\left\\{\int_{r_{\mathrm{h}}}^{r}\frac{c_{9}dr}{\sqrt{rb_{1}(r)}}\right\\}^{2}\,.$
(35)

When $r_{\mathrm{h}1}<r<r_{\mathrm{h}}$, $b_{1}(r)$ is negative and therefore
$b(r)$ becomes negative,

$\displaystyle
b(r)=-\frac{r^{2}}{4}\left\\{\int_{r_{\mathrm{h}}}^{r}\frac{c_{9}dr}{\sqrt{-rb_{1}(r)}}\right\\}^{2}\,.$
(36)

In the region $r<r_{\mathrm{h}1}$, however, $b(r)$ becomes a complex number,

$\displaystyle b(r)=\frac{r^{2}}{4}\left\\{\pm
i\int_{r_{\mathrm{h}}}^{r_{\mathrm{h}1}}\frac{c_{9}dr}{\sqrt{-rb_{1}(r)}}+\int_{r_{\mathrm{h}1}}^{r}\frac{c_{9}dr}{\sqrt{rb_{1}(r)}}\right\\}^{2}\,.$
(37)

Therefore the model becomes physically inconsistent and irrelevant or there is
the end of space at $r=r_{\mathrm{h}1}$. The above arguments show that the
model can have only one horizon, even in Case III as in Cases I and II.

Using Eq. (III), we get the mimetic field $\phi$ and the function $\omega$ in
the form,

$\displaystyle\phi(r)=\frac{1}{\omega(r)}=\left(\int\\!{\frac{dr}{\sqrt{4\,{\frac{c_{7}}{r}}-4\,rc_{6}-2\,r^{2}c_{5}-2\,{c_{8}}^{-1}}}}\right)^{2}+2\,\int\\!{\frac{dr}{\sqrt{4\,{\frac{c_{7}}{r}}-4\,rc_{6}-2\,r^{2}c_{5}-2\,{c_{8}}^{-1}}}}c_{9}+{c_{9}}^{2}\,.$
(38)

We choose the constant of the integration $c_{9}$ so that $\phi$ vanishes when
$b(r)$ and $b_{1}(r)$ vanishes, that is

$\displaystyle\phi(r)=\left(\int_{r_{\mathrm{h}}}^{r}\frac{dr}{\sqrt{\frac{4c_{7}}{r}-4rc_{6}-2r^{2}c_{5}-2{c_{8}}^{-1}}}\right)^{2}=\frac{1}{2}\left(\int_{r_{\mathrm{h}}}^{r}\frac{dr}{\sqrt{b_{1}(r)}}\right)^{2}\,.$
(39)

The expression (39) is valid for $r$ where $b_{1}(r)>0$. For $r$ where
$b_{1}(r)<0$, $\phi$ becomes negative and real again,

$\displaystyle\phi(r)=-\frac{1}{2}\left(\int_{r_{\mathrm{h}}}^{r}\frac{dr}{\sqrt{-b_{1}(r)}}\right)^{2}\,.$
(40)

Therefore the formalism using $\omega(\phi)$ in (8) is valid, again.

##  IV Summary and Discussion

The topic of DM became a hot topic based on the data coming from observations
of astrophysics as well as cosmology. To explain the issue of dark matter,
scientists supposed two directions: The first one is to amend the field
equation of GR and the second direction is to amend the standard model by
building a new particle kind. Many studies have explained that the two
directions are equal Sebastiani  _et al._ (2017). Additionally, it is familiar
in the frame of amended theories of gravity that they have new degrees of
freedom in addition to the usual massless graviton of GR. As we discussed in
Section I that the target of the mimetic theory of gravity was mimetic DM, and
was a good theory to explain the presence of cold DM. We have shown, however,
that the mimetic theory with the constraint
$g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1$ cannot realize the
black hole geometry with the horizon(s). To bypass such problem, we changed
the mimetic constraint as
$\omega(\phi)g^{\rho\sigma}\partial_{\rho}\phi\partial_{\sigma}\phi=-1,$ where
$\omega(\phi)$ is a function of the scalar field $\phi$. Because we can choose
$\omega(\phi)$ changes its signature at the horizon, the constraint becomes
consistent with the existence of the horizon(s). As applications of the
formalism of the constraint, we considered the mimetic gravity in
$(2+1)$-dimensional space-time and found that the formalism works well.

In future work, we may consider the time-dependent black hole solution or the
black hole solution in the expanding universe, which is difficult to be found
even in other gravity theories. In the case of the mimetic theory, the vector
$(\partial_{\alpha}\phi)$ is time-like when we consider the cosmology and the
vector is space-like when we consider the static object like a black hole. If
we use the original mimetic constraint given by Eq. (2), it is difficult to
include both the time-like vector and the space-like vector for
$(\partial_{\alpha}\phi)$ but by introducing $\omega(\phi)$ as in Eq. (3), it
becomes possible to include both the time-like and the space-like vectors in a
solution, which may help to construct the time-dependent black hole solution
or the black hole solution in the expanding universe.

###### Acknowledgements.

This work is supported by the JSPS Grant-in-Aid for Scientific Research (C)
No. 18K03615 (S.N.).

## References

  * Hawking  _et al._ (2016) S. W. Hawking, M. J. Perry, and A. Strominger, Phys. Rev. Lett. 116, 231301 (2016), arXiv:1601.00921 [hep-th] .
  * Akiyama  _et al._ (2019) K. Akiyama _et al._ (Event Horizon Telescope), Astrophys. J. Lett. 875, L1 (2019), arXiv:1906.11238 [astro-ph.GA] .
  * Afshar  _et al._ (2016) H. Afshar, S. Detournay, D. Grumiller, W. Merbis, A. Perez, D. Tempo, and R. Troncoso, Phys. Rev. D 93, 101503 (2016), arXiv:1603.04824 [hep-th] .
  * Haco  _et al._ (2018) S. Haco, S. W. Hawking, M. J. Perry, and A. Strominger, JHEP 12, 098 (2018), arXiv:1810.01847 [hep-th] .
  * Haco  _et al._ (2019) S. Haco, M. J. Perry, and A. Strominger, (2019), arXiv:1902.02247 [hep-th] .
  * Grumiller  _et al._ (2020) D. Grumiller, A. Pérez, M. M. Sheikh-Jabbari, R. Troncoso, and C. Zwikel, Phys. Rev. Lett. 124, 041601 (2020), arXiv:1908.09833 [hep-th] .
  * Sheykhi (2021) A. Sheykhi, JHEP 01, 043 (2021), arXiv:2009.12826 [gr-qc] .
  * Chamseddine and Mukhanov (2013) A. H. Chamseddine and V. Mukhanov, JHEP 11, 135 (2013), arXiv:1308.5410 [astro-ph.CO] .
  * Chamseddine and Mukhanov (2017a) A. H. Chamseddine and V. Mukhanov, JCAP 1703, 009 (2017a), arXiv:1612.05860 [gr-qc] .
  * Chamseddine and Mukhanov (2017b) A. H. Chamseddine and V. Mukhanov, Eur. Phys. J. , 183 (2017b), arXiv:1612.05861 [gr-qc] .
  * Casalino  _et al._ (2018) A. Casalino, M. Rinaldi, L. Sebastiani, and S. Vagnozzi, Phys. Dark Univ. 22, 108 (2018), arXiv:1803.02620 [gr-qc] .
  * Casalino  _et al._ (2019) A. Casalino, M. Rinaldi, L. Sebastiani, and S. Vagnozzi, Class. Quant. Grav. 36, 017001 (2019), arXiv:1811.06830 [gr-qc] .
  * Vagnozzi (2017) S. Vagnozzi, Class. Quant. Grav. 34, 185006 (2017), arXiv:1708.00603 [gr-qc] .
  * Sheykhi and Grunau (2019) A. Sheykhi and S. Grunau, (2019), arXiv:1911.13072 [gr-qc] .
  * Chamseddine  _et al._ (2014) A. H. Chamseddine, V. Mukhanov, and A. Vikman, JCAP 1406, 017 (2014), arXiv:1403.3961 [astro-ph.CO] .
  * Dutta  _et al._ (2018) J. Dutta, W. Khyllep, E. N. Saridakis, N. Tamanini, and S. Vagnozzi, JCAP 1802, 041 (2018), arXiv:1711.07290 [gr-qc] .
  * Abbassi  _et al._ (2018) M. H. Abbassi, A. Jozani, and H. R. Sepangi, Phys. Rev. D 97, 123510 (2018), arXiv:1803.00209 [gr-qc] .
  * Zhong and Sáez-Chillón Gómez (2018) Y. Zhong and D. Sáez-Chillón Gómez, Symmetry 10, 170 (2018), arXiv:1805.03467 [gr-qc] .
  * Matsumoto (2016) J. Matsumoto, (2016), arXiv:1610.07847 [astro-ph.CO] .
  * Nojiri and Odintsov (2014) S. Nojiri and S. D. Odintsov,  (2014), 10.1142/S0217732314502113, [Erratum: Mod.Phys.Lett.A 29, 1450211 (2014)], arXiv:1408.3561 [hep-th] .
  * Odintsov and Oikonomou (2016a) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 93, 023517 (2016a), arXiv:1511.04559 [gr-qc] .
  * Nojiri  _et al._ (2016) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Rev. D 94, 104050 (2016), arXiv:1608.07806 [gr-qc] .
  * Sadeghnezhad and Nozari (2017) N. Sadeghnezhad and K. Nozari, Phys. Lett. , 134 (2017), arXiv:1703.06269 [gr-qc] .
  * Nashed (2021) G. G. L. Nashed, Astrophys. J. 919, 113 (2021), arXiv:2108.04060 [gr-qc] .
  * Gorji  _et al._ (2019) M. A. Gorji, S. Mukohyama, and H. Firouzjahi, JCAP 05, 019 (2019), arXiv:1903.04845 [gr-qc] .
  * Gorji  _et al._ (2018a) M. A. Gorji, S. Mukohyama, H. Firouzjahi, and S. A. Hosseini Mansoori, JCAP 08, 047 (2018a), arXiv:1807.06335 [hep-th] .
  * Bouhmadi-Lopez  _et al._ (2017) M. Bouhmadi-Lopez, C.-Y. Chen, and P. Chen, JCAP 1711, 053 (2017), arXiv:1709.09192 [gr-qc] .
  * Gorji  _et al._ (2018b) M. A. Gorji, S. A. Hosseini Mansoori, and H. Firouzjahi, JCAP 1801, 020 (2018b), arXiv:1709.09988 [astro-ph.CO] .
  * Firouzjahi  _et al._ (2018) H. Firouzjahi, M. A. Gorji, S. A. Hosseini Mansoori, A. Karami, and T. Rostami, JCAP 11, 046 (2018), arXiv:1806.11472 [gr-qc] .
  * Sebastiani  _et al._ (2017) L. Sebastiani, S. Vagnozzi, and R. Myrzakulov, Adv. High Energy Phys. 2017, 3156915 (2017), arXiv:1612.08661 [gr-qc] .
  * Chamseddine  _et al._ (2019) A. H. Chamseddine, V. Mukhanov, and T. B. Russ, Eur. Phys. J. C 79, 558 (2019), arXiv:1905.01343 [hep-th] .
  * Deruelle and Rua (2014) N. Deruelle and J. Rua, JCAP 1409, 002 (2014), arXiv:1407.0825 [gr-qc] .
  * Myrzakulov and Sebastiani (2015) R. Myrzakulov and L. Sebastiani, Gen. Rel. Grav. 47, 89 (2015), arXiv:1503.04293 [gr-qc] .
  * Myrzakulov  _et al._ (2016) R. Myrzakulov, L. Sebastiani, S. Vagnozzi, and S. Zerbini, Class. Quant. Grav. 33, 125005 (2016), arXiv:1510.02284 [gr-qc] .
  * Astashenok  _et al._ (2015) A. V. Astashenok, S. D. Odintsov, and V. K. Oikonomou, Class. Quant. Grav. 32, 185007 (2015), arXiv:1504.04861 [gr-qc] .
  * Odintsov and Oikonomou (2015) S. D. Odintsov and V. K. Oikonomou, Annals Phys. 363, 503 (2015), arXiv:1508.07488 [gr-qc] .
  * Nojiri  _et al._ (2017) S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Lett. B 775, 44 (2017), arXiv:1710.07838 [gr-qc] .
  * Odintsov and Oikonomou (2018) S. D. Odintsov and V. K. Oikonomou, Nucl. Phys. B 929, 79 (2018), arXiv:1801.10529 [gr-qc] .
  * Oikonomou (2016a) V. K. Oikonomou, Universe 2, 10 (2016a), arXiv:1511.09117 [gr-qc] .
  * Nashed (2014) G. G. L. Nashed, EPL 105, 10001 (2014), arXiv:1501.00974 [gr-qc] .
  * Gorji  _et al._ (2020) M. A. Gorji, A. Allahyari, M. Khodadi, and H. Firouzjahi, Phys. Rev. D 101, 124060 (2020), arXiv:1912.04636 [gr-qc] .
  * Nashed  _et al._ (2019) G. G. L. Nashed, W. El Hanafy, and K. Bamba, JCAP 01, 058 (2019), arXiv:1809.02289 [gr-qc] .
  * Chen  _et al._ (2018) C.-Y. Chen, M. Bouhmadi-Lopez, and P. Chen, Eur. Phys. J. , 59 (2018), arXiv:1710.10638 [gr-qc] .
  * Nashed and Nojiri (2021a) G. G. L. Nashed and S. Nojiri, Phys. Rev. D 104, 044043 (2021a), arXiv:2107.13550 [gr-qc] .
  * Nashed and Nojiri (2021b) G. Nashed and S. Nojiri,  2021, 007 (2021b).
  * Nashed and Nojiri (2021c) G. G. L. Nashed and S. Nojiri, (2021c), arXiv:2110.08560 [gr-qc] .
  * Nashed (2018) G. G. L. Nashed, International Journal of Geometric Methods in Modern Physics 15, 1850154-1056 (2018).
  * Nashed (2018) G. Nashed, Symmetry 10, 559 (2018).
  * Ben Achour  _et al._ (2018) J. Ben Achour, F. Lamy, H. Liu, and K. Noui, JCAP 05, 072 (2018), arXiv:1712.03876 [gr-qc] .
  * Zheng  _et al._ (2017) Y. Zheng, L. Shen, Y. Mou, and M. Li, JCAP 1708, 040 (2017), arXiv:1704.06834 [gr-qc] .
  * Elizalde  _et al._ (2020) E. Elizalde, G. G. L. Nashed, S. Nojiri, and S. D. Odintsov, Eur. Phys. J. C 80, 109 (2020), arXiv:2001.11357 [gr-qc] .
  * Shen  _et al._ (2019) L. Shen, Y. Zheng, and M. Li, JCAP 12, 026 (2019), arXiv:1909.01248 [gr-qc] .
  * Sheykhi (2020) A. Sheykhi, JHEP 07, 031 (2020), arXiv:2002.11718 [gr-qc] .
  * Cid and Labrana (2012) A. Cid and P. Labrana, Phys. Lett. B 717, 10 (2012), arXiv:1209.2339 [astro-ph.CO] .
  * Capozziello  _et al._ (2010) S. Capozziello, J. Matsumoto, S. Nojiri, and S. D. Odintsov, Phys. Lett. B 693, 198 (2010), arXiv:1004.3691 [hep-th] .
  * Oikonomou (2016b) V. K. Oikonomou, Mod. Phys. Lett. A 31, 1650191 (2016b), arXiv:1609.03156 [gr-qc] .
  * Nashed (2010) G. G. L. Nashed, Astrophys. Space Sci. 330, 173 (2010), arXiv:1503.01379 [gr-qc] .
  * Oikonomou (2016c) V. K. Oikonomou, Int. J. Mod. Phys. D 25, 1650078 (2016c), arXiv:1605.00583 [gr-qc] .
  * Nashed (2008) G. G. L. Nashed, Eur. Phys. J. C 54, 291 (2008), arXiv:0804.3285 [gr-qc] .
  * Myrzakulov and Sebastiani (2016) R. Myrzakulov and L. Sebastiani, Astrophys. Space Sci. 361, 188 (2016), arXiv:1601.04994 [gr-qc] .
  * Odintsov and Oikonomou (2016b) S. D. Odintsov and V. K. Oikonomou, Astrophys. Space Sci. 361, 236 (2016b), arXiv:1602.05645 [gr-qc] .
  * Odintsov and Oikonomou (2016c) S. D. Odintsov and V. K. Oikonomou, Phys. Rev. D 94, 044012 (2016c), arXiv:1608.00165 [gr-qc] .
  * Chen  _et al._ (2021) J. Chen, W.-D. Guo, and Y.-X. Liu, Eur. Phys. J. C 81, 709 (2021), arXiv:2011.03927 [gr-qc] .
  * Oikonomou (2015) V. K. Oikonomou, Phys. Rev. D 92, 124027 (2015), arXiv:1509.05827 [gr-qc] .
  * Nashed (2007) G. G. L. Nashed, Eur. Phys. J. C 49, 851 (2007), arXiv:0706.0260 [gr-qc] .
  * Zhong and Elizalde (2016) Y. Zhong and E. Elizalde, Mod. Phys. Lett. A 31, 1650221 (2016), arXiv:1612.04179 [gr-qc] .
  * Banados  _et al._ (1992) M. Banados, C. Teitelboim, and J. Zanelli, Phys. Rev. Lett. 69, 1849 (1992), arXiv:hep-th/9204099 .
  * Carlip (1995) S. Carlip, Class. Quant. Grav. 12, 2853 (1995), arXiv:gr-qc/9506079 .
  * Ashtekar  _et al._ (2003) A. Ashtekar, J. Wisniewski, and O. Dreyer, Adv. Theor. Math. Phys. 6, 507 (2003), arXiv:gr-qc/0206024 .
  * Sarkar  _et al._ (2006) T. Sarkar, G. Sengupta, and B. Nath Tiwari, JHEP 11, 015 (2006), arXiv:hep-th/0606084 .
  * Witten (1998) E. Witten, Adv. Theor. Math. Phys. 2, 505 (1998), arXiv:hep-th/9803131 .
  * Carlip (2005) S. Carlip, Class. Quant. Grav. 22, R85 (2005), arXiv:gr-qc/0503022 .
  * Birmingham  _et al._ (2002) D. Birmingham, I. Sachs, and S. N. Solodukhin, Phys. Rev. Lett. 88, 151301 (2002), arXiv:hep-th/0112055 .
  * Witten (2007) E. Witten, (2007), arXiv:0706.3359 [hep-th] .
  * Banados  _et al._ (1993) M. Banados, M. Henneaux, C. Teitelboim, and J. Zanelli, Phys. Rev. D 48, 1506 (1993), [Erratum: Phys.Rev.D 88, 069902 (2013)], arXiv:gr-qc/9302012 .
  * Martinez  _et al._ (2000) C. Martinez, C. Teitelboim, and J. Zanelli, Phys. Rev. D 61, 104013 (2000), arXiv:hep-th/9912259 .
  * Fernando and Mansouri (1997) S. Fernando and F. Mansouri, (1997), arXiv:gr-qc/9705016 .
  * Ren (2010) J. Ren, JHEP 11, 055 (2010), arXiv:1008.3904 [hep-th] .
  * Liu  _et al._ (2011) Y. Liu, Q. Pan, and B. Wang, Phys. Lett. B 702, 94 (2011), arXiv:1106.4353 [hep-th] .
  * Kord Zangeneh  _et al._ (2017) M. Kord Zangeneh, Y. C. Ong, and B. Wang, Phys. Lett. B 771, 235 (2017), arXiv:1704.00557 [hep-th] .
  * Mohammadi  _et al._ (2018) M. Mohammadi, A. Sheykhi, and M. Kord Zangeneh, Eur. Phys. J. C 78, 654 (2018), arXiv:1805.07377 [hep-th] .
  * Ghotbabadi  _et al._ (2018) B. B. Ghotbabadi, M. Kord Zangeneh, and A. Sheykhi, Eur. Phys. J. C 78, 381 (2018), arXiv:1804.05442 [hep-th] .
  * Panotopoulos and Rincón (2018) G. Panotopoulos and A. Rincón, Phys. Rev. D 97, 085014 (2018), arXiv:1804.04684 [hep-th] .
  * Rincón and Panotopoulos (2018) A. Rincón and G. Panotopoulos, Phys. Rev. D 97, 024027 (2018), arXiv:1801.03248 [hep-th] .
  * Clement (1996) G. Clement, Phys. Lett. B 367, 70 (1996), arXiv:gr-qc/9510025 .
  * Nojiri and Odintsov (1998) S. Nojiri and S. D. Odintsov, Mod. Phys. Lett. A 13, 2695 (1998), arXiv:gr-qc/9806034 .
  * Emparan  _et al._ (2000) R. Emparan, G. T. Horowitz, and R. C. Myers, JHEP 01, 021 (2000), arXiv:hep-th/9912135 .
  * Cadoni and Setare (2008) M. Cadoni and M. R. Setare, JHEP 07, 131 (2008), arXiv:0806.2754 [hep-th] .
  * Parsons and Ross (2009) J. Parsons and S. F. Ross, JHEP 04, 134 (2009), arXiv:0901.3044 [hep-th] .
  * Hendi (2012) S. H. Hendi, JHEP 03, 065 (2012), arXiv:1405.4941 [hep-th] .
  * Sheykhi  _et al._ (2014) A. Sheykhi, S. H. Hendi, and S. Salarpour, Phys. Scripta 10, 105003 (2014).
  * Xu (2014) W. Xu, Phys. Lett. B 738, 472 (2014), arXiv:1409.3368 [hep-th] .
  * Grumiller and Merbis (2020) D. Grumiller and W. Merbis, SciPost Phys. 8, 010 (2020), arXiv:1906.10694 [hep-th] .


# Chaos in Coupled Heteroclinic Cycles between Weak Chimeras

Artyom E. Emelin Lobachevsky University, Gagarin ave, 23, 603022 Nizhny
Novgorod, Russia Evgeny A. Grines Lobachevsky University, Gagarin ave, 23,
603022 Nizhny Novgorod, Russia Tatiana A. Levanova Lobachevsky University,
Gagarin ave, 23, 603022 Nizhny Novgorod, Russia

###### Abstract

Heteroclinic cycles are widely used in neuroscience in order to mathematically
describe different mechanisms of functioning of the brain and nervous system.
Heteroclinic cycles and interactions between them can be a source of different
types of nontrivial dynamics. For instance, as it was shown earlier, chaotic
dynamics can appear as a result of interaction via diffusive couplings between
two stable heteroclinic cycles between saddle equilibria. We go beyond these
findings by considering two rotating in opposite directions coupled stable
heteroclinic cycles between weak chimeras. Such an ensemble can be
mathematically described by a system of six phase equations. Using two
parameter bifurcation analysis we investigate the scenarios of appearance and
destruction of chaotic dynamics in the system under study.

## INTRODUCTION

One of the main goals of contemporary neurodynamics is to study the
fundamental principles of the brain and nervous system functioning [1]. Many
mechanisms in neurodynamics and biology are based on various interacting
periodic processes [2]. In order to describe such processes one can use
networks of different topologies with oscillators as network nodes [3]. In
this case, from the point of view of neurodynamics, a single element of such a
network can describe a certain group of neurons in the nervous system [4].

The functioning of described networks highly depends on the properties and
dynamics of its individual elements, as well as on the chosen topology and the
strength of couplings [5]. For the interacting populations of identical
oscillators a wide range of collective behaviour is possible, starting from
global synchronization [6, 7] to the patterns of localized synchrony [8, 9,
10]. In the ensembles of coupled oscillators with higher-order interactions
more sophisticated types of dynamics are possible, e.g., heteroclinic
switchings between patterns of localized frequency synchrony.

Studies of heteroclinic switchings between different metastable states in the
context of computational neuroscience have a long history [11, 12, 13, 14,
15]. Starting from pioneering works by V. Afraimovich [16, 17, 18], where
heteroclinic sequences consisting of saddles and heteroclinic orbits
connecting them were shown to be a mathematical images of sequential activity
in neural networks, to later research devoted to heteroclinic cycles based on
various saddle limit sets, such as saddle cycles [13, 19, 20, 21] and saddle
chaotic sets [22, 23], as well as to hierarchic heteroclinic structures of
different topology [24, 25].

However, the description of metastable states is not limited to the saddle
limit sets mentioned above. Chimera states, which consists of localized
patterns of frequency synchronization, are also one of the possible candidates
for describing such states. Dynamical transitions of the location of frequency
synchrony were observed in [26], where the system under study possessed
metastable chimera states coupled by heteroclinic transitions. In subsequent
papers [27] and [28], the stability of heteroclinic cycles between chimera
states was studied in detail. Also, previously obtained results were
generalized to cases of a larger number of elements and different topologies
of couplings between elements.

More complex structures based on coupled subsystems, each of which possesses a
stable heteroclinic cycle (here and further we will use the term "coupled
heteroclinic cycles"), have attracted particular interest recently. In [29] an
example of a such system with nearly symmetric coupling between its subsystems
was studied. Also, in the paper [30] several identical heteroclinic cycles
coupled by a diffusive term were studied. It was shown that in this case, with
different parameters of the coupling strength, the system is able to reproduce
various patterns of synchronous, quasi-periodic and chaotic dynamics.
Subsequently, in the work [31] a system was proposed, consisting of two
diffusively coupled subsystems, each of which is based on generalized a Lotka-
Volterra model. It was shown that the described model demonstrates the
predominant chaotic dynamics even for weak couplings between heteroclinic
cycles.

The goal of this paper is to generalize this result to systems consisting of
more complex heteroclinic cycles, in particular, to systems of coupled
heteroclinic cycles between chimeric states. In order to do this we consider a
system of six differential equations that describes the behavior of two
interconnected clusters with oppositely directed heteroclinic cycles. Each
cluster is described by a system of phase oscillators proposed in [27]. Based
on the fundamental principles of constructing networks of phase oscillators,
in our case, the most appropriate coupling function is selected, which is
necessary to obtain various complex patterns of behavior of a dynamic system.
Using various numerical methods, examples of the resulting collective dynamics
will be demonstrated.

The paper is organized as follows. We introduce the model in Section 1 paying
special attention to the following building blocks of the model. We discuss
the basic heteroclinic cycle between chimera states in Subsection 1.1. We
introduce couplings between two heteroclinic cycles and the resulting full
model in Subsection 1.2. The coupled system possesses several symmetries,
which we discuss in Subsection 1.3. Numerical exploration of described coupled
heteroclinic cycles is performed in Section 2. We discuss our findings and
draw conclusion in Section 3.

##  1 THE MODEL

###  1.1 Heteroclinic cycle between chimera states

To describe the dynamics of one heteroclinic cycle, we will use the system of
phase oscillators [27]. This model reproduces switching activity between
chimera states in the sense of Ashwin-Burylko [8]. The model consists of three
populations, each of which consist of two phase oscillators. The topology of
described ensemble is presented in Fig. 1

Figure 1: The topology of ensemble of phase oscillators that can reproduce
heteroclinic cycle between chimera states.

Initially, in order to describe the dynamics of the heteroclinic cycle of this
type one can use the following system of six ODEs:

$\left\\{\begin{array}[]{ll}\dot{\theta}_{\sigma,1}=\omega+g_{2}(\theta_{\sigma,2}-\theta_{\sigma,1})-KG_{4}(\theta_{\sigma-1};\theta_{\sigma,2}-\theta_{\sigma,1})+KG_{4}(\theta_{\sigma+1};\theta_{\sigma,2}-\theta_{\sigma,1})\\\
\dot{\theta}_{\sigma,2}=\omega+g_{2}(\theta_{\sigma,1}-\theta_{\sigma,2})-KG_{4}(\theta_{\sigma-1};\theta_{\sigma,1}-\theta_{\sigma,2})+KG_{4}(\theta_{\sigma+1};\theta_{\sigma,1}-\theta_{\sigma,2})\\\
i\in\overline{1,6}\end{array}\right.$ (1)

Here, interaction between oscillators within populations is pairwise and
determined by the function

$g_{2}(\nu)=\sin(\nu+\alpha_{2})-r\sin(2(\nu+\alpha_{2})),$

with parameters $\alpha_{2}$ and $r$. The interaction between populations is
determined by the non-pairwise interaction function

$G_{4}(\theta_{\tau};\phi)=\frac{1}{4}(g_{4}(\theta_{\tau,1}-\theta_{\tau,2}+\phi)+g_{4}(\theta_{\tau,2}-\theta_{\tau,1}+\phi)),$

where

$g_{4}(\nu)=\sin(\nu+\alpha_{4}).$

Here, the parameter $K$ is responsible for the strength of the connection
between the populations. Parameter $\alpha_{4}$ is responsible for the
interaction between populations.

Using reduction approach proposed in [27] we can reduce the dimension of the
system to three equations. In order to do this we introduce a phase difference
between phases of oscillators within each population, which allows us to
rewrite the system (1) in the following form:

$\dot{\psi_{i}}=\widehat{g_{2}}(\psi_{i})-\frac{K}{2}(\widehat{g_{4}}(\psi_{i-1}+\psi_{i})+\widehat{g_{4}}(\psi_{i}-\psi_{i-1}))+\frac{K}{2}(\widehat{g_{4}}(\psi_{i+1}+\psi_{i})+\widehat{g_{4}}(\psi_{i}-\psi_{i+1})),$
(2)

where

$\widehat{g}_{l}(\nu)=\frac{1}{2}(g_{l}(-\nu)-g_{l}(\nu)),\leavevmode\nobreak\
\leavevmode\nobreak\ l\in\\{2,4\\}$ (3)

It is easy to see that points of the phase space with any coordinate equal to
0 or $\pi$ are equilibrium states of the reduced system (2). Such equilibrium
states have the following interpretation: the elements of the $i$-th
population are in a synchronized state at the phase value $\phi_{i}=0$ (we
will denote it by the letter $S$), and a desynchronized state at the phase
value $\phi_{i}=\pi$ (we will denote it by the letter $D$). All such
combinations (except SSS and DDD) correspond to the saddle point equilibrium
of the system (2). In [27] it was shown that in the system (2) there is a
rough stable heteroclinic cycle that includes these metastable states (Fig.
2). It is noteworthy that a change in the sign of the parameter $K$ leads to a
change in the order of traversal of the heteroclinic cycle, but does not
affect its stability. Clearly, if $(x(t),\,y(t),\,z(t))$ is a trajectory of
system (2) for a certain $K$, then it can be shown by direct calculations that
$(x(t),\,z(t),\,y(t))$ is a solution of this system for $-K$ while the other
parameters’ values are fixed.

(a) (b)

Figure 2: Sequence of metastable states in a system (2)

###  1.2 Two coupled heteroclinic cycles

Let us start with two copies of system (1) that differ from each other only in
the sign of the parameter $K$. We will write it as follows:

$\left\\{\begin{array}[]{ll}\dot{\theta}_{\sigma,1}=\omega+g_{2}(\theta_{\sigma,2}-\theta_{\sigma,1})-\,KG_{4}(\theta_{\sigma-1};\theta_{\sigma,2}-\theta_{\sigma,1})+\,KG_{4}(\theta_{\sigma+1};\theta_{\sigma,2}-\theta_{\sigma,1})+\,\varepsilon\mathcal{C}(\theta_{\sigma,2}-\theta_{\sigma,1},\rho_{\sigma,2}-\rho_{\sigma,1}),\\\
\dot{\theta}_{\sigma,2}=\omega+g_{2}(\theta_{\sigma,1}-\theta_{\sigma,2})-\,KG_{4}(\theta_{\sigma-1};\theta_{\sigma,1}-\theta_{\sigma,2})+\,KG_{4}(\theta_{\sigma+1};\theta_{\sigma,1}-\theta_{\sigma,2})+\,\varepsilon\mathcal{C}(\theta_{\sigma,1}-\theta_{\sigma,2},\rho_{\sigma,2}-\rho_{\sigma,1}),\\\
\dot{\rho}_{\sigma,1}=\omega+g_{2}(\rho_{\sigma,2}-\rho_{\sigma,1})+\,KG_{4}(\rho_{\sigma-1};\rho_{\sigma,2}-\rho_{\sigma,1})-\,KG_{4}(\rho_{\sigma+1};\rho_{\sigma,2}-\rho_{\sigma,1})+\,\varepsilon\mathcal{C}(\rho_{\sigma,2}-\rho_{\sigma,1},\theta_{\sigma,2}-\theta_{\sigma,1}),\\\
\dot{\rho}_{\sigma,2}=\omega+g_{2}(\rho_{\sigma,1}-\rho_{\sigma,2})+\,KG_{4}(\rho_{\sigma-1};\rho_{\sigma,1}-\rho_{\sigma,2})-\,KG_{4}(\rho_{\sigma+1};\rho_{\sigma,1}-\rho_{\sigma,2})+\,\varepsilon\mathcal{C}(\rho_{\sigma,1}-\rho_{\sigma,2},\theta_{\sigma,2}-\theta_{\sigma,1}),\\\
\sigma\in\overline{1,3}.\end{array}\right.$ (4)

We will call each of the groups of oscillators whose dynamics is described by
the same letter a cluster. The coupling function $\mathcal{C}(x,y)$ governs
the interaction between populations of different clusters where $\varepsilon$
specifies the strength of this interaction. The topology of the resulting
ensemble is presented in Fig. 3.

Figure 3: The topology of the system (4).

Since equations (4) still depend only on phase differences, a reduction
similar to the one in the previous Subsection can be applied here. Namely,
introducing variables $\varphi_{\sigma}=\theta_{\sigma,2}-\theta_{\sigma,1}$
and $\psi_{\sigma}=\rho_{\sigma,2}-\rho_{\sigma,1}$ transforms equations (4)
into system

$\left\\{\begin{array}[]{l}\dot{\varphi}_{\sigma}=2\widehat{g_{2}}(\varphi_{\sigma})-\frac{K}{2}(\widehat{g_{4}}(\varphi_{\sigma-1}+\varphi_{\sigma})+\widehat{g_{4}}(\varphi_{\sigma}-\varphi_{\sigma-1}))+...\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
...+\frac{K}{2}(\widehat{g_{4}}(\varphi_{\sigma+1}+\varphi_{\sigma})+\widehat{g_{4}}(\varphi_{\sigma}-\varphi_{\sigma+1}))+\varepsilon
I(\varphi_{\sigma},\psi_{\sigma}).\par\\\
\dot{\psi}_{\sigma}=2\widehat{g_{2}}(\psi_{\sigma})+\frac{K}{2}(\widehat{g_{4}}(\psi_{\sigma-1}+\psi_{\sigma})+\widehat{g_{4}}(\psi_{\sigma}-\psi_{\sigma-1}))-...\\\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
...-\frac{K}{2}(\widehat{g_{4}}(\psi_{\sigma+1}+\psi_{\sigma})+\widehat{g_{4}}(\psi_{\sigma}-\psi_{\sigma+1}))+\varepsilon
I(\psi_{\sigma},\varphi_{\sigma}),\\\ \end{array}\right.$ (5)

where $I(x,y)=\mathcal{C}(-x,y)-\mathcal{C}(x,y)$ is a coupling function
between variables describing the behaviour of populations in reduced system.
While from the mathematical point of view the original inter-cluster coupling
function $\mathcal{C}(x,y)$ can be any periodic function of its arguments, the
inter-cluster coupling function between phase differences $I(x,y)$ clearly
must be odd with respect to its first argument.

Plugging the expressions for other coupling functions, we can rewrite Eqs. (5)
as

$\left\\{\begin{array}[]{ll}\dot{\varphi}_{1}=\left(-A\cos{\varphi_{2}}+A\cos{\varphi_{3}}+B\cos{\varphi_{1}}+C\right)\sin{\varphi_{1}}+\varepsilon
I(\varphi_{1},\psi_{1}),\\\
\dot{\varphi}_{2}=\left(-A\cos{\varphi_{3}}+A\cos{\varphi_{1}}+B\cos{\varphi_{2}}+C\right)\sin{\varphi_{2}}+\varepsilon
I(\varphi_{2},\psi_{2}),\\\
\dot{\varphi}_{3}=\left(-A\cos{\varphi_{1}}+A\cos{\varphi_{2}}+B\cos{\varphi_{3}}+C\right)\sin{\varphi_{3}}+\varepsilon
I(\varphi_{3},\psi_{3}),\\\
\dot{\psi}_{1}=\left(A\cos{\psi_{2}}-A\cos{\psi_{3}}+B\cos{\psi_{1}}+C\right)\sin{\psi_{1}}+\varepsilon
I(\psi_{1},\varphi_{1}),\\\
\dot{\psi}_{2}=\left(A\cos{\psi_{3}}-A\cos{\psi_{1}}+B\cos{\psi_{2}}+C\right)\sin{\psi_{2}}+\varepsilon
I(\psi_{2},\varphi_{2}),\\\
\dot{\psi}_{3}=\left(A\cos{\psi_{1}}-A\cos{\psi_{2}}+B\cos{\psi_{3}}+C\right)\sin{\psi_{3}}+\varepsilon
I(\psi_{3},\varphi_{3}),\\\ \end{array}\right.$ (6)

where $A=K\cos{\alpha_{4}}$, $B=4r\cos{2\alpha_{2}}$ and
$C=-2\cos{\alpha_{2}}$.

The system (6) possesses an important property: $\varphi_{i}=\pi k$ and
$\psi_{i}=\pi k$ are invariant hyperplanes for $k\in\mathbb{Z}$. In that case
the phase space is divided into invariant cubes. For example, the cube
$\\{\varphi_{i}\in[0,\,\pi],\;\psi_{i}\in[0,\,\pi]\\}$ is invariant: its faces
are invariant, and trajectories that start in the interior of this cube does
not leave it due to the invariance of cubes’ faces.

(a) (b)

Figure 4: Coupling functions. (a) $I(x)=\sin{x}$ and (b)
$I(z,w)=\sinh{z}\operatorname{sech}{w}-\tanh{w}$

While this is true for ’’natural’’ coupling $I(x,y)$ that comes from full
system (4) with coupling $\mathcal{C}(x,y)$, some of these properties are
preserved when we consider a "phase-like" system similar to Eqs. (6) where we
do not require $I(x,y)$ to be odd in its first argument. Namely, from now on
let us consider

$I(x,y)=\sin{(y-x)}.$

Such choice breaks the invariance of cube’s faces, but its interior is still
invariant: the vector field at the boundary now points inside or it is tangent
to it. Indeed, without loss of generality let us check $\dot{\varphi}_{1}$
when $\varphi_{1}=0$. Plugging that into Eqs. (6) gives
$\dot{\varphi}_{1}|_{\varphi_{1}=0}=\varepsilon\sin{\psi_{1}}$, which is non-
negative when $\psi_{1}\in[0,\;\pi]$.

When $\varphi_{1}=\pi$, we obtain
$\dot{\varphi}_{1}|_{\varphi_{1}=\pi}=\varepsilon\sin{(\psi_{1}-\pi)}=-\varepsilon\sin{\psi_{1}}$
which is non-positive when $\psi_{1}\in[0,\;\pi]$. Both these inequalities
entail that vector field points inside of the cube at hyperplanes
$\varphi_{1}=0$ and $\varphi_{1}=\pi$. The same reasoning could be applied to
all faces of the cube, hence we can conclude that trajectories starting in the
cube’s interior do not leave it.

In order to study dynamics of system (6) in this invariant cube we can utilize
a coordinate transformation

$z_{k}=\mathcal{F}(\varphi_{k}),\;w_{k}=\mathcal{F}(\psi_{k}),$ (7)

where

$\mathcal{F}(x)=\ln\sqrt{\frac{1-\cos{x}}{1+\cos{x}}}=\ln\tan{\frac{x}{2}}.$
(8)

Such transformation is a monotonous map from the interval $(0,\pi)$ to the
interval $(-\infty,+\infty)$: the values close to $0$ are mapped close to
$-\infty$, whereas values close to $\pi$ are mapped close to $+\infty$. This
transformation is inspired by the work of Pikovsky and Nepomnyaschy [31] and
serves the same purpose. In the absence of coupling (i.e., when
$\varepsilon=0$) the attractors of the system are likely to be situated on the
invariant faces of the cube. In that case a simulation of dynamics might
become very sensitive to unavoidable errors of approximation: numerical
trajectories might either stuck close to saddle equilibria or step outside of
invariant cube due to imprecise calculations. This transformation plays a role
of a "magnifying glass": when attractors are situated on the cube’s faces,
variables $z_{k}$ and $w_{k}$ tend to infinity, whereas they stay finite even
when an attractor is close to the previously invariant faces.

Let us write a transformed version of Eqs. (6) using the substitution (7).
Using tangent half-angle substitution

$\sin{x}=\frac{2u}{1+u^{2}},\;\cos{x}=\frac{1-u^{2}}{1+u^{2}},\;u=\tan{\frac{x}{2}},$

and utilizing reordered version of relation between angular and logarithmic
coordinates

$e^{{z_{k}}}=\tan{\frac{\varphi_{k}}{2}},\;e^{{w_{k}}}=\tan{\frac{\psi_{k}}{2}},$

we can express basic trigonometric functions of $\varphi_{k}$ and $\psi_{k}$
in terms of hyperbolic functions of $z_{k}$ and $w_{k}$:

$\cos{\varphi_{k}}=\frac{1-\tan^{2}{\frac{\varphi_{k}}{2}}}{1+\tan^{2}{\frac{\varphi_{k}}{2}}}=\frac{1-e^{2z_{k}}}{1+e^{2z_{k}}}=-\tanh{{z_{k}}},\;\sin{\varphi_{k}}=\frac{2\tan{\frac{\varphi_{k}}{2}}}{1+\tan^{2}{\frac{\varphi_{k}}{2}}}=\frac{2e^{{z_{k}}}}{1+e^{2z_{k}}}=\frac{1}{\cosh{{z_{k}}}};$
(9)

similar formulas can be written for $\psi_{k}$ and $w_{k}$.

Finally, combining formulas (9) and the fact that
$\mathcal{F}^{\prime}(x)=\frac{1}{\sin{x}}$, we arrive at the transformed
system of equations in logarithmic coordinates:

$\left\\{\begin{array}[]{ll}\dot{z}_{1}=\left(A\tanh{{z_{2}}}-A\tanh{{z_{3}}}-B\tanh{{z_{1}}}+C\right)+\varepsilon\bigl{(}\tanh{{w_{1}}}-\sinh{{z_{1}}}\;\operatorname{sech}{{w_{1}}}\bigr{)},\\\
\dot{z}_{2}=\left(A\tanh{{z_{3}}}-A\tanh{{z_{1}}}-B\tanh{{z_{2}}}+C\right)+\varepsilon\bigl{(}\tanh{{w_{2}}}-\sinh{{z_{2}}}\;\operatorname{sech}{{w_{2}}}\bigr{)},\\\
\dot{z}_{3}=\left(A\tanh{{z_{1}}}-A\tanh{{z_{2}}}-B\tanh{{z_{3}}}+C\right)+\varepsilon\bigl{(}\tanh{{w_{3}}}-\sinh{{z_{3}}}\;\operatorname{sech}{{w_{3}}}\bigr{)},\\\
\dot{w}_{1}=\left(-A\tanh{{w_{2}}}+A\tanh{{w_{3}}}-B\tanh{{w_{1}}}+C\right)+\varepsilon\bigl{(}\tanh{{z_{1}}}-\sinh{{w_{1}}}\;\operatorname{sech}{{z_{1}}}\bigr{)},\\\
\dot{w}_{2}=\left(-A\tanh{{w_{3}}}+A\tanh{{w_{1}}}-B\tanh{{w_{2}}}+C\right)+\varepsilon\bigl{(}\tanh{{z_{2}}}-\sinh{{w_{2}}}\;\operatorname{sech}{{z_{2}}}\bigr{)},\\\
\dot{w}_{3}=\left(-A\tanh{{w_{1}}}+A\tanh{{w_{2}}}-B\tanh{{w_{3}}}+C\right)+\varepsilon\bigl{(}\tanh{{z_{3}}}-\sinh{{w_{3}}}\;\operatorname{sech}{{z_{3}}}\bigr{)}.\end{array}\right.$
(10)

While all terms involving $\tanh$ and $\operatorname{sech}$ functions are
bounded, the only $\sinh$ term in each equation grows unbounded and
counteracts an unlimited growth of phase variable to either of infinities (at
least, theoretically).

###  1.3 Symmetries

Let us highlight the symmetry properties of the system (10). It is well known
that a presence of symmetries influences both possible geometries of
attractors (which can be asymmetric or posess some subgroup of self-
symmetries) and their possible bifurcations. This short section forms the
basis for an explanation in Section 2 of pitchfork bifurcation curves
appearing in parameter space.

Recall that a transformation of a phase space
$\mathcal{T}:\mathbb{R}^{n}\mapsto\mathbb{R}^{n}$ is called a symmetry of an
ODE $\dot{x}=F(x)$ if for any solution $\gamma(t)$ of this equation follows
that $\mathcal{T}(\gamma(t))$ is also a solution of this ODE. From this
definition follows a computable criterion: a map $\mathcal{T}$ is a symmetry
if $F(\mathcal{T}(x))\equiv\mathcal{T}^{\prime}(x)\cdot F(x)$ for any
$x\in\mathbb{R}^{n}$, where $\mathcal{T}^{\prime}(x)$ is a Jacobi matrix of
the map $\mathcal{T}$ computed at the point $x$. The symmetries of (10) are
the same as in [31] due to the same coupling principle between clusters and
the same symmetries of original equations (1) for a single cluster.

Using the notation from [31], we can say that system (10) is symmetric with
respect to following transformations:

  * •

renamings of variables inside each subsystem:

$\displaystyle\quad
R_{1}:\left(z_{1},z_{2},z_{3},w_{1},w_{2},w_{3}\right)\rightarrow\left(z_{2},z_{3},z_{1},w_{2},w_{3},w_{1}\right)$
$\displaystyle\text{ and
}R_{2}:\left(z_{1},z_{2},z_{3},w_{1},w_{2},w_{3}\right)\rightarrow\left(z_{3},z_{1},z_{2},w_{3},w_{1},w_{2}\right)\text{.
}$

  * •

exchange transformations between subsystems:

$\begin{gathered}T_{1}:\left(z_{1},z_{2},z_{3},w_{1},w_{2},w_{3}\right)\rightarrow\left(w_{1},w_{3},w_{2},z_{1},z_{3},z_{2}\right),\\\
T_{2}:\left(z_{1},z_{2},z_{3},w_{1},w_{2},w_{3}\right)\rightarrow\left(w_{3},w_{2},w_{1},z_{3},z_{2},z_{1}\right),\\\
\text{ and
}T_{3}:\left(z_{1},z_{2},z_{3},w_{1},w_{2},w_{3}\right)\rightarrow\left(w_{2},w_{1},w_{3},z_{2},z_{1},z_{3}\right).\end{gathered}$

There are following relations between these symmetries:
$R_{2}=R_{1}^{-1}=R_{1}^{2}$, $R_{1}=R_{2}^{-1}=R_{2}^{2}$,
$T_{1}^{2}=T_{2}^{2}=T_{3}^{2}=I$,$T_{2}=R_{2}T_{1}$, $T_{3}=R_{1}T_{1}$,
where $I$ is an identity map.

##  2 THE RESULTS

In order to investigate the phenomenon of appearance of chaotic dynamics in
the system (10) we conducted a numerical study. The results discussed below
were obtained by us using Julia software [32] and a package for numerical
bifurcation analysis MATCONT [33].

Our starting point for studying and describing dynamics of the system (10) is
a stable limit cycle found at parameter values

$\begin{array}[]{c}\alpha_{2}=\frac{\pi}{2},\;r=0.016850701746869687,\;\alpha_{4}=\pi,\\\
K=0.05977203237457928,\varepsilon=0.05206157664979711.\end{array}$

This is a comparatively large-period limit cycle (the period is roughly $200$
time units according to MATCONT), see Fig. 7. This limit cycle $\mathcal{L}$
was used as a starting point for an attractor continuation in parameters
$\varepsilon$ and $K$, leaving all other parameter values fixed.

Figure 5: The map of the Lyapunov spectrum of the system (10) and its enlarged
fragments. See detailed legend in subfigures

A map of the Lyapunov spectrum based on this attractor continuation is
presented on Fig. 5. Here the red marker corresponds to the regions with
hyperchaotic dynamics (i.e. regions with two positive Lyapunov exponents,
$\lambda_{1}>\lambda_{2}>0$), the green marker – to the regions with chaotic
dynamics (one positive Lyapunov exponents, $\lambda_{1}>0$), yellow marker –
to the regions with torus attractors ($\lambda_{1,2}\approx
0>\lambda_{3,4,5,6}$), blue marker – to the regions with limit cycle
attractors ($\lambda_{1}\approx 0>\lambda_{2,3,4,5,6}$), and white marker – to
the regions where trajectories of the system (10) diverge to infinity. A
detailed description of the calculation of maps of the largest Lyapunov
exponent can be found, e.g., in the following studies [34, 35].

While theoretically equations (10) preclude unbounded growth of solutions, in
certain cases such effect occurs even during attractor continuation. Here we
consider a trajectory going to infinity if any of the phase variable exceeds
20 in absolute value. It corresponds to a trajectory of the system (6) being
at distance of $4.12\times 10^{-9}$ (an approximate value of
$2\cdot\arctan{e^{-20}}$, see formula (8)) to $0$ or ${\pi}$ in some of its
phase variable. Let us describe some features of dynamics in the region of
bounded dynamics.

Continuing the limit cycle $\mathcal{L}$ with fixed $K$ and decreasing
$\varepsilon$ leads to its loss of stability in a Neimark-Sacker bifurcation.
This observation agrees with MATCONT computations, which locate the position
of $R_{3}$ and $R_{4}$ bifurcation points at the tips of periodicity windows
found at the map of Lyapunov exponents. Further change of parameters to the
left of Neimark-Sacker curve leads to a chaotization of attracting invariant
tori and formation of chaotic and hyperchaotic regimes.

Figure 6: A setup for a pitchfork bifurcation, $K=0.05977203237457928$. Two
unstable limit cycles $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ collide with a
stable limit cycle at $\varepsilon=0.05274407266311118$, which leads to its
loss of stability. (a) Projections of cycles $\mathcal{C}_{1}$ (green) and
$\mathcal{C}_{2}$ (orange) show that they are visibly different at
$\varepsilon\approx 0.0527226$ before the bifurcation; (b) however, they are
symmetry related: $\mathcal{C}_{2}$ (orange) and image of $\mathcal{C}_{1}$
(blue) under symmetry $T_{3}$ are the same set-wise; (c) limit cycle
$\mathcal{C}$ (gray) and its $T_{2}$ image (pink) are visibly different; (d)
meanwhile, $\mathcal{C}$ (gray) and its $T_{3}$ image (red) are the same set-
wise.

Continuing the limit cycle $\mathcal{L}$ to the right by increasing
$\varepsilon$ leads to the loss of its stability in a subcritical pitchfork
bifurcation. While in generic systems pitchfork bifurcation is usually a
codimension-2 phenomenon, here it is a codimension-1 bifurcation due to the
presence of symmetries. Fig. 6 illustrates the details of this bifurcation.
Here the initial stable limit cycle $\mathcal{L}$ already has a non-trivial
isotropy subgroup, i.e. a subgroup of symmetry group that leaves limit cycle
$\mathcal{L}$ invariant as a set. Namely, this limit cycle is mapped to itself
by a symmetry $T_{3}$ as illustrated on Fig. 6d, while other symmetries map it
to a completely different set (as an example, action of $T_{2}$ on it is
illustrated on Fig. 6c).

However, the unstable limit cycles $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$
that exist before bifurcation happens do not posess this self-symmetry, but
are mapped to each other by this $T_{3}$ symmetry. As parameter changes, they
both get close to a limit cycle $\mathcal{L}$ and collide with it, leading to
a stability loss; thus, a subcritical pitchfork bifurcation due to the system
being equivariant. After the loss of stability the trajectories go to a new
attractor that is not located near the old one. This repeats in other parts of
parameter space, which might explain a certain artifacts of the map of
Lyapunov exponents: in that case attractor and its characteristics might jump
abruptly.

A presence of multiple coexisting attractors can be further illustrated by the
right insert on Fig. 5. Here an attracting limit cycle found in that region
was used for numerical continuation in MATCONT. Fixing $K$ as in previous case
and varying $\varepsilon$ we find two bifurcation curves that bound stability
region of this regime: a curve of pitchfork bifurcation and a curve of period-
doubling bifurcation. The right insert also shows other bifurcations that
happen during this continuation, revealing a complex organization of parameter
and phase space. While we were unsuccessful in further numerical continuation
of large-period limit cycles (period is $\approx 400$) and discovering other
period doubling curves, we expect that in this region chaotic dynamics occurs
through a cascade of period doubling bifurcations.

Let us describe numerically typical dynamical regimes that can be observed in
the system (10) for different combinations of governing parameters
$\varepsilon$ and $K$ in the Fig. 5. Examples of regular and chaotic dynamics
are presented in Fig. 7 and Fig. 8, correspondingly.

Figure 7: Regular dynamics in a system (10). Left panel: time series. Right
panel: projections of 6-dimensional phase space of the system (10) on the
different phase planes. Parameter values: $\varepsilon=0.05206157664979711$,
$K=0.05977203237457928$. Figure 8: Chaotic dynamics in a system (10). Left
panel: time series. Right panel: projections of 6-dimensional phase space of
the system (10) on the different phase planes. Parameter values:
$\varepsilon=0.0527$, $K=0.06096$.

##  3 CONCLUSION

In this paper, we have studied coupled heteroclinic cycles between chimera
states rotating in opposite directions on the basis of the system of phase
equations. We generalize earlier result by Pikovsky and Nepomnyashchy [31],
where the appearance of chaotic dynamics was shown as a result of interaction
between oppositely rotating heteroclinic cycles between saddle equilibria for
the case of weak diffusive coupling between heteroclinics. Using two-parameter
bifurcation analysis we have identified regions on the plane of governing
parameters where different types of nontrivial dynamics exist, including
hyperchaotic, chaotic and regular spatio-temporal patterns. Note that in our
numerical simulations of the system under study, we have observed that chaotic
and hyperchaotic dynamics arise for small values of inter-cluster and intra-
cluster coupling strengths. The increase in these coupling strength leads to a
series of bifurcation transitions related to the change of symmetry of the
chaotic attractors.

We assume that the origin of chaotic and hyperchaotic dynamics in similar
systems is a general effect that can be observed in a range of governing
parameters for heteroclinic cycles between different types of attracting sets
(e.g. saddle cycles, saddle chaotic sets etc.) and various types of couplings
between them (e.g. mutual synaptic couplings). Nevertheless, the verification
of this hypothesis is a subject of future studies.

This paper is devoted to Prof. V. N. Belykh, who contributed greatly to the
various field of nonlinear dynamics, such as studies of synchronization and
clustering in ensembles of coupled dynamic systems [36, 37, 38, 39, 40, 7],
existence of homoclinic [41, 42, 43] and heteroclinic orbits and neurodynamics
[44, 45, 46, 47].

## ACKNOWLEDGMENTS

The authors thank Prof. A. S. Pikovsky and Prof. G.V. Osipov for useful
discussions.

## FUNDING

This work was supported by the Ministry of Science and Education of Russian
Federation Contract no. FSRW-2020-0036 (A.E.E. and E.A.G.) and RSF grant
22-12-00348 (T.A.L.).

## CONFLICT OF INTEREST

The authors declare that they have no conflicts of interest.

## References

  * [1] Christian Bick, Marc Goodfellow, Carlo R Laing, and Erik A Martens.  Understanding the dynamics of biological and neural oscillator networks through exact mean-field reductions: a review.  The Journal of Mathematical Neuroscience, 10:1–43, 2020.
  * [2] Arkady Pikovsky, Michael Rosenblum, and Jürgen Kurths.  Synchronization: a universal concept in nonlinear science, 2002.
  * [3] Steven H Strogatz.  Exploring complex networks.  Nature, 410(6825):268–276, 2001.
  * [4] Michael Breakspear.  Dynamic models of large-scale brain activity.  Nature neuroscience, 20(3):340–352, 2017.
  * [5] Arthur T Winfree.  The geometry of biological time, volume 2.  Springer, 1980.
  * [6] Vladimir N Belykh, Valentin S Petrov, and Grigory V Osipov.  Dynamics of the finite-dimensional kuramoto model: Global and cluster synchronization.  Regular and Chaotic Dynamics, 20(1):37–48, 2015.
  * [7] Nikita V Barabash, Vladimir N Belykh, Grigory V Osipov, and Igor V Belykh.  Partial synchronization in the second-order kuramoto model: An auxiliary system method.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 31(11), 2021.
  * [8] Peter Ashwin and Oleksandr Burylko.  Weak chimeras in minimal networks of coupled phase oscillators.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 25(1), 2015.
  * [9] Christian Bick and Peter Ashwin.  Chaotic weak chimeras and their persistence in coupled populations of phase oscillators.  Nonlinearity, 29(5):1468, 2016.
  * [10] Oleh E Omel’chenko.  The mathematics behind chimera states.  Nonlinearity, 31(5):R121, 2018.
  * [11] Valentin Afraimovich, Peter Ashwin, and Vivien Kirk.  Robust heteroclinic and switching dynamics.  Dynamical Systems, 25(3):285–286, 2010.
  * [12] Peter Ashwin, Özkan Karabacak, and Thomas Nowotny.  Criteria for robustness of heteroclinic cycles in neural microcircuits.  The Journal of Mathematical Neuroscience, 1(1):1–18, 2011.
  * [13] MA Komarov, GV Osipov, and JAK Suykens.  Sequentially activated groups in neural networks.  EPL (Europhysics Letters), 86(6):60006, 2009.
  * [14] VI Nekorkin, AS Dmitrichev, DV Kasatkin, and VS Afraimovich.  Relating the sequential dynamics of excitatory neural networks to synaptic cellular automata.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 21(4):043124, 2011.
  * [15] Vladimir Isaakovicih Nekorkin, AS Dmitrichev, DV Kasatkin, and Valentin Sendorovich Afraimovich.  Reducing the sequential dynamics of excitatory neural networks to cellular automata.  JETP letters, 95(9):492–496, 2012.
  * [16] Valentin S Afraimovich, Sze-Bi Hsu, and Huey-Er Lin.  Chaotic behavior of three competing species of may–leonard model under small periodic perturbations.  International Journal of Bifurcation and Chaos, 11(02):435–447, 2001.
  * [17] VS Afraimovich, VP Zhigulin, and MI Rabinovich.  On the origin of reproducible sequential activity in neural circuits.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 14(4):1123–1129, 2004.
  * [18] Valentin S Afraimovich, Mikhail I Rabinovich, and Pablo Varona.  Heteroclinic contours in neural ensembles and the winnerless competition principle.  International Journal of Bifurcation and Chaos, 14(04):1195–1208, 2004.
  * [19] MA Komarov, GV Osipov, and CS Zhou.  Heteroclinic contours in oscillatory ensembles.  Physical Review E, 87(2):022909, 2013.
  * [20] Tatiana A Levanova, Maxim A Komarov, and Grigory V Osipov.  Sequential activity and multistability in an ensemble of coupled van der pol oscillators.  The European Physical Journal Special Topics, 222(10):2417–2428, 2013.
  * [21] AO Mikhaylov, MA Komarov, TA Levanova, and GV Osipov.  Sequential switching activity in ensembles of inhibitory coupled oscillators.  EPL (Europhysics Letters), 101(2):20009, 2013.
  * [22] Michael Dellnitz, Michael Field, Martin Golubitsky, Andreas Hohmann, and Jun Ma.  Cycling chaos.  IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 42(10):821–823, 1995.
  * [23] Tatiana A Levanova, Grigory V Osipov, and Arkady Pikovsky.  Coherence properties of cycling chaos.  Communications in Nonlinear Science and Numerical Simulation, 19(8):2734–2739, 2014.
  * [24] Valentin S Afraimovich, Todd R Young, and Mikhail I Rabinovich.  Hierarchical heteroclinics in dynamical model of cognitive processes: chunking.  International Journal of Bifurcation and Chaos, 24(10):1450132, 2014.
  * [25] Valentin S Afraimovich, Michael A Zaks, and MI3799967 Rabinovich.  Mind-to-mind heteroclinic coordination: model of sequential episodic memory initiation.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 28(5), 2018.
  * [26] Christian Bick.  Heteroclinic switching between chimeras.  Physical Review E, 97(5):050201, 2018.
  * [27] Christian Bick.  Heteroclinic dynamics of localized frequency synchrony: heteroclinic cycles for small populations.  Journal of Nonlinear Science, 29(6):2547–2570, 2019.
  * [28] Christian Bick and Alexander Lohse.  Heteroclinic dynamics of localized frequency synchrony: Stability of heteroclinic cycles and networks.  Journal of Nonlinear Science, 29(6):2571–2600, 2019.
  * [29] Dong Li, MC Cross, Changsong Zhou, and Zhigang Zheng.  Quasiperiodic, periodic, and slowing-down states of coupled heteroclinic cycles.  Physical Review E, 85(1):016215, 2012.
  * [30] Maximilian Voit, Sara Veneziale, and Hildegard Meyer-Ortmanns.  Coupled heteroclinic networks in disguise.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8):083113, 2020.
  * [31] Arkady Pikovsky and Alexander Nepomnyashchy.  Chaos in coupled heteroclinic cycles and its piecewise-constant representation.  Physica D: Nonlinear Phenomena, page 133772, 2023.
  * [32] G Datseris.  Dynamicalsystems jl: A julia software library for chaos and nonlinear dynamics. j. open source softw. 3 (23): 598,(2018).
  * [33] Annick Dhooge, Willy Govaerts, Yu A Kuznetsov, Hil Gaétan Ellart Meijer, and Bart Sautois.  New features of the software matcont for bifurcation analysis of dynamical systems.  Mathematical and Computer Modelling of Dynamical Systems, 14(2):147–175, 2008.
  * [34] Alexey V Borisov, Alexey Yu Jalnine, Sergey P Kuznetsov, Igor R Sataev, and Julia V Sedova.  Dynamical phenomena occurring due to phase volume compression in nonholonomic model of the rattleback.  Regular and Chaotic Dynamics, 17:512–532, 2012.
  * [35] Alexey V Borisov, Alexey O Kazakov, and Igor R Sataev.  The reversal and chaotic attractor in the nonholonomic model of chaplygin’s top.  Regular and Chaotic Dynamics, 19:718–733, 2014.
  * [36] Vladimir N Belykh, Igor V Belykh, and Martin Hasler.  Connection graph stability method for synchronized coupled chaotic systems.  Physica D: nonlinear phenomena, 195(1-2):159–187, 2004.
  * [37] Igor V Belykh, Vladimir N Belykh, and Martin Hasler.  Blinking model and synchronization in small-world networks with a time-varying coupling.  Physica D: Nonlinear Phenomena, 195(1-2):188–206, 2004.
  * [38] NV Barabash and VN Belykh.  Synchronization thresholds in an ensemble of kuramoto phase oscillators with randomly blinking couplings.  Radiophysics and Quantum Electronics, 60(9):761–768, 2018.
  * [39] Igor V Belykh, Barrett N Brister, and Vladimir N Belykh.  Bistability of patterns of synchrony in kuramoto oscillators with inertia.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 26(9), 2016.
  * [40] Barrett N Brister, Vladimir N Belykh, and Igor V Belykh.  When three is a crowd: Chaos from clusters of kuramoto oscillators with inertia.  Physical Review E, 101(6):062206, 2020.
  * [41] Vladimir Nikolaevich Belykh.  Homoclinic and heteroclinic trajectories of a family of multidimensional dynamical systems.  Trudy Matematicheskogo Instituta imeni VA Steklova, 216:20–31, 1997.
  * [42] Vladimir N Belykh and Evgeniya Valerevna Pankratova.  Chaotic dynamics of two van der pol-duffing oscillators with huygens coupling.  Regular and Chaotic dynamics, 15:274–284, 2010.
  * [43] Vladimir Nikolaevich Belykh.  Bifurcation of separatrices of a saddle of the lorenz system.  Differentsial’nye Uravneniya, 20(10):1666–1674, 1984.
  * [44] VN Belykh, IV Belykh, M Colding-Jørgensen, and Erik Mosekilde.  Homoclinic bifurcations leading to the emergence of bursting oscillations in cell models.  The european physical journal E, 3:205–219, 2000.
  * [45] Vladimir Belykh, Igor Belykh, and Erik Mosekilde.  Hyperbolic plykin attractor can exist in neuron models.  International Journal of Bifurcation and Chaos, 15(11):3567–3578, 2005.
  * [46] Vladimir N Belykh, Evgeniya V Pankratova, and Erik Mosekilde.  Dynamics and synchronization of noise perturbed ensembles of periodically activated neuron cells.  International Journal of Bifurcation and Chaos, 18(09):2807–2815, 2008.
  * [47] Nikita V Barabash, Tatiana A Levanova, and Vladimir N Belykh.  Ghost attractors in blinking lorenz and hindmarsh–rose systems.  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(8), 2020.