# Tachyonic cosmology with cubic contractions of the Riemann tensor

Mihai Marciu  Faculty of Physics, University of Bucharest
<EMAIL_ADDRESS>
###### Abstract

A new cosmological theory is proposed in the theoretical framework of modified
gravity theories which is based on a tachyonic field non–minimally coupled
with a specific topological invariant constructed with third order
contractions of the Riemann tensor. After proposing the action of the
cosmological model, the modified Friedmann relations and the Klein–Gordon
equations are obtained, describing the corresponding geometrical corrections
to the Einstein–Hilbert action. The physical features of the cosmological
system are investigated by adopting the dynamical system analysis in the case
of an exponential function for the geometrical coupling term. The
investigation revealed that the cosmological system can explain the current
accelerated expansion of the Universe and the matter dominated epoch, showing
a high compatibility to the recent history of our Universe for various values
of the coupling coefficients.

##  I Introduction

The general relativity theory represents an important success in describing
the behavior of the Universe at the large scale structure dynamics. This
theory represents the fundamental structure considered for describing the
large scale dynamics of the Universe, later embedded into the $\Lambda$CDM
model [1]. Although the $\Lambda$CDM simple model can explain various
phenomena at the large scale [1], it suffers from various pathological
inconsistencies [2, 3, 4, 5]. The modified gravity approach [6, 7, 8, 9, 10]
represents an important direction which aims to further correct the basic
theory of general relativity, by taking into account various possible
couplings with different specific invariant components which are embedded into
the corresponding action [10]. To this regard, many studies have proposed
different modifications of the basic Einstein–Hilbert action, by including or
extending the geometrical sector [10, 11]. The main aim of these theories is
related to the description of the accelerated expansion at a consistent level,
offering some possible solution to various fundamental problems in the modern
cosmology [12, 7, 13]. From an experimental point of view, the accelerated
expansion represents a cryptic phenomena at the large scale, affecting the
large scale evolution of the Universe [12, 14]. The consistent description and
analysis of this phenomena is expected to offer new insights and revolutions
in the modern cosmology, having specific effects in the development of various
physical theories and technologies. Since the discovery, the accelerated
expansion have been probed through various astrophysical studies [15, 16, 17,
18, 19].

In the modified gravity landscape the Einsteinian cubic gravity [20]
represents an interesting and possible viable theory studied in the recent
years [21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31]. This theory is based on a
specific invariant, constructed using third order contractions of the Riemann
tensor [20]. The non–linear extension of the Einsteinian cubic gravity was
proposed in Ref. [32], a theory capable of explaining the late time
acceleration of the Universe. The dynamical analysis of the later theory was
performed in Refs. [33, 34] for specific parameterizations in the action.
Recently, various black hole solutions have been investigated in the
Einsteinian cubic gravity theory [29, 35, 36, 37, 38, 39]. Furthermore, the
properties of the inflationary era have been analyzed in some specific
theories [40, 41, 42, 43, 44, 45, 46]. The coupling of a scalar field with a
cubic term has been investigated by considering a dynamical system analysis
[47]. Recently, it has been shown that the cubic gravity theory can be
associated to the developments of various types of pathological instabilities
[48] which have to be addressed in order to construct a viable theory [49].

In the scalar tensor theories, a special class of dark energy models is
represented by the tachyonic cosmologies, a novel approach originated from
string theory [50, 51]. The development of tachyonic models has been
considered in the past years, leading to various theoretical constructions
[52, 53, 54, 55, 56, 57, 58] which can explain the recent accelerated
expansion at the large scale structure. In the scalar tensor theories based on
general relativity the study of tachyonic fields have been considered in
various cosmological applications [59, 60, 61, 62, 63]. Moreover, the study of
cosmological models containing tachyonic fields have been applied in
teleparallel gravity [64, 65, 66, 67, 68, 69, 70], a viable alternative theory
[71]. In these theories, the choice of the potential energy term play a
fundamental role, dictating the future and past dynamics of the aforementioned
models [72, 61]. The study of tachyonic cosmological scenarios for various
classes of potential energies have been considered [73, 74, 65, 72], leading
to viable models which can explain various cosmological features.

In the present paper we shall further extend the tachyonic cosmology by
considering a non–minimal coupling with a novel topological invariant, based
on cubic contractions of the Riemann tensor. After we deduce the corresponding
field equation, we shall study the physical implications by adopting the
dynamical system analysis. The study takes into consideration two specific
cases associated to the behavior of the potential energy term. In the first
case we consider an exponential representation, while in the second one an
inverse hyperbolic sine function is studied. For all of the previous mentioned
cases we have analyzed the structure and properties of the phase space,
discussing possible physical effects.

The present manuscript is organized as follows: in Sec. II we discuss the
action and the corresponding field equations for the dark energy model. Then,
in Sec. III we analyze the physical features for an exponential coupling and
potential by considering the dynamical system analysis. In Sec. IV we discuss
the phase space structure where the potential energy term is beyond the usual
exponential case, considering an inverse hyperbolic function. Finally, in Sec.
V we have a short summary of our analysis, discussing the main conclusions
which are applicable to the present study.

##  II The description of the field equations

In the present study we shall consider a tachyonic cosmological model
non–minimally coupled with a topological invariant constructed from the cubic
contractions of the Riemann tensor. The action corresponding to the present
study is the following:

$S=S_{m}+\int
d^{4}x\sqrt{-g}\Big{[}\frac{R}{2}-V(\phi)\sqrt{1+\epsilon\frac{g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi}{V(\phi)}}+f(\phi)P\Big{]},$
(1)

where the topological invariant is based on specific contractions of the
Riemann tensor in the third order [32],

$P=\beta_{1}R_{\mu\quad\nu}^{\quad\rho\quad\sigma}R_{\rho\quad\sigma}^{\quad\gamma\quad\delta}R_{\gamma\quad\delta}^{\quad\mu\quad\nu}+\beta_{2}R_{\mu\nu}^{\rho\sigma}R_{\rho\sigma}^{\gamma\delta}R_{\gamma\delta}^{\mu\nu}\\\
+\beta_{3}R^{\sigma\gamma}R_{\mu\nu\rho\sigma}R_{\quad\quad\gamma}^{\mu\nu\rho}+\beta_{4}RR_{\mu\nu\rho\sigma}R^{\mu\nu\rho\sigma}+\beta_{5}R_{\mu\nu\rho\sigma}R^{\mu\rho}R^{\nu\sigma}\\\
+\beta_{6}R_{\mu}^{\nu}R_{\nu}^{\rho}R_{\rho}^{\mu}+\beta_{7}R_{\mu\nu}R^{\mu\nu}R+\beta_{8}R^{3},$
(2)

with $\beta_{j},(j=\\{1,...,8\\})$ constant parameters.

In this case the potential energy is $V(\phi)$, and $\epsilon$ is a constant
parameter which describes the canonical representation of the tachyonic field.
For a canonical tachyonic field $\epsilon=+1$, while for the non–canonical
case we have $\epsilon=-1$. The matter part in the action is denoted by
$S_{m}$, describing a barotropic fluid which characterizes the dark matter
sector, having the corresponding density $\rho_{m}$ and pressure $p_{m}$,
satisfying the equation of state $p_{m}=\rho_{m}w_{m}$, with $w_{m}$ a
constant coefficient describing a non–relativistic behavior. In this case, the
dark matter fluid satisfies the standard continuity equation.

Next, the large scale structure dynamics in the Universe is described by the
following Robertson–Walker metric,

$ds^{2}=-dt^{2}+a^{2}(t)\delta_{uv}dx^{u}dx^{v},$ (3)

with $a(t)$ the scale factor, $H(t)$ the Hubble parameter, and $t$ the cosmic
time. If we adopt the following relations between the constant parameters
$\beta_{j},(j=\\{1,...,8\\})$ [32],

$\beta_{7}=\frac{1}{12}\big{[}3\beta_{1}-24\beta_{2}-16\beta_{3}-48\beta_{4}-5\beta_{5}-9\beta_{6}\big{]},$
(4)
$\beta_{8}=\frac{1}{72}\big{[}-6\beta_{1}+36\beta_{2}+22\beta_{3}+64\beta_{4}+3\beta_{5}+9\beta_{6}\big{]},$
(5) $\beta_{6}=4\beta_{2}+2\beta_{3}+8\beta_{4}+\beta_{5},$ (6)
$\bar{\beta}=(-\beta_{1}+4\beta_{2}+2\beta_{3}+8\beta_{4}),$ (7)

then the third order tensor polynomial constructed from the cubic contractions
of the Riemann tensor is equal to [32]:

$P=6\bar{\beta}H^{4}(2H^{2}+3\dot{H}).$ (8)

Here, the dot represents the differentiation with respect to the cosmic time
$t$, while the prime denotes the differentiation with respect to the argument
of the specific function. The variation of the action described in the Eq. 1
with respect to the tachyonic field $\phi(t)$ gives the corresponding
Klein–Gordon equation,

$\Big{[}1-\epsilon\frac{\dot{\phi}^{2}}{V(\phi)}\Big{]}^{-\frac{3}{2}}\Big{[}3\epsilon\frac{V^{\prime}(\phi)\dot{\phi}^{2}}{V(\phi)}-2(V^{\prime}(\phi)+\epsilon\ddot{\phi})\Big{]}\\\
+12\beta H^{4}f^{\prime}(\phi)(2H^{2}+3\dot{H})-\frac{6\epsilon
H\dot{\phi}}{\sqrt{1-\epsilon\frac{\dot{\phi}^{2}}{V(\phi)}}}=0.$ (9)

Furthermore, the variation of the action (1) with respect to the inverse
metric leads to the modified Friedmann relations which have the following form
[32]:

$3H^{2}=\rho_{m}+\rho_{\phi},$ (10) $3H^{2}+2\dot{H}=-p_{m}-p_{\phi},$ (11)

where the energy density of the tachyonic field is

$\rho_{\phi}=\frac{V(\phi)}{\sqrt{1-\epsilon\frac{\dot{\phi}^{2}}{V(\phi)}}}+6\beta
f(\phi)H^{6}-18\beta H^{5}\frac{df(\phi)}{d\phi}\dot{\phi},$ (12)

with the pressure

$p_{\phi}=-V(\phi)\sqrt{1-\epsilon\frac{\dot{\phi}^{2}}{V(\phi)}}-6\beta
f(\phi)H^{6}-12\beta f(\phi)H^{4}\dot{H}\\\ +12\beta
H^{5}\frac{df(\phi)}{d\phi}\dot{\phi}+24\beta
H^{3}\frac{df(\phi)}{d\phi}\dot{H}\dot{\phi}\\\ +6\beta
H^{4}\dot{\phi}^{2}\frac{d^{2}f(\phi)}{d\phi^{2}}+6\beta
H^{4}\frac{df(\phi)}{d\phi}\ddot{\phi}.$ (13)

Finally, we can define the barotropic parameter associated to the dark energy
field,

$w_{\bf{\phi}}=\frac{p_{\phi}}{\rho_{\phi}},$ (14)

and the effective (total) equation of state for our cosmological model,

$w_{\bf{eff}}=\frac{p_{m}+p_{\phi}}{\rho_{m}+\rho_{\phi}}=-1-\frac{2}{3}\frac{\dot{H}}{H^{2}}.$
(15)

If we introduce the matter density parameter,

$\Omega_{m}=\frac{\rho_{m}}{3H^{2}},$ (16)

and the density parameter corresponding to the tachyonic field,

$\Omega_{\phi}=\frac{\rho_{\phi}}{3H^{2}},$ (17)

we have the following constraint,

$\Omega_{m}+\Omega_{\phi}=1.$ (18)

##  III Dynamical effects in the case of an exponential potential

In this section we shall discuss the main physical features of the present
cosmological model by applying the dynamical system analysis, an important
tool in the study of various modified gravity theories. Analyzing the
Friedmann constraint equation (10) we introduce the following dimension–less
variables:

$x=\frac{\dot{\phi}}{\sqrt{V(\phi)}},$ (19)

$y=\frac{\sqrt{V(\phi)}}{H\sqrt{3}},$ (20)

$z=\beta f(\phi)H^{4},$ (21)

in the case of an exponential coupling,

$f(\phi)=f_{0}e^{\alpha\phi},$ (22)

and potential energy,

$V(\phi)=e^{-\lambda\phi}.$ (23)

Hence, in this case we have

$V^{\prime}(\phi)=-\lambda V(\phi),$ (24)

where $\lambda$ is a positive constant coefficient which characterizes the
steepness of the potential energy.

The Friedmann constraint equation (10) can be written as:

$1=\frac{y^{2}}{\sqrt{1-\epsilon x^{2}}}+2z-6\alpha zxy\sqrt{3}+\Omega_{m}.$
(25)

Next, if we introduce another variable $N=log(a)$ and change the dependence of
the dimension–less components to $N$ (the e-fold variable), we obtain the
following autonomous dynamical system, approximating the evolution of the
cosmological system:

$x^{\prime}=\frac{\ddot{\phi}}{H^{2}}\frac{1}{y\sqrt{3}}+\frac{\lambda}{2}\sqrt{3}x^{2}y,$
(26)

$y^{\prime}=-\lambda\frac{\sqrt{3}}{2}xy^{2}-y\frac{\dot{H}}{H^{2}},$ (27)

$z^{\prime}=\alpha zxy\sqrt{3}+4z\frac{\dot{H}}{H^{2}},$ (28)

where the prime ′ describe the derivative with respect to $N$.

In this case the Klein–Gordon equation (9) can be written in the following
way:

$-\frac{6\sqrt{3}H^{2}xy\epsilon}{\sqrt{1-x^{2}\epsilon}}+\frac{-27H^{4}\lambda
x^{2}y^{4}\epsilon-6H^{2}y^{2}\left(\epsilon\ddot{\phi}-3H^{2}\lambda
y^{2}\right)}{3H^{2}y^{2}\left(1-x^{2}\epsilon\right)^{3/2}}+12\alpha
z\left(2H^{2}+3\dot{H}\right)=0,$ (29)

while the acceleration equation (11) is equal to:

$-3H^{2}-2\dot{H}=3H^{2}w_{m}\Omega_{m}+18\alpha^{2}H^{2}x^{2}y^{2}z-3H^{2}y^{2}\sqrt{1-x^{2}\epsilon}+12\sqrt{3}\alpha
H^{2}xyz-6H^{2}z+24\sqrt{3}\alpha xyz\dot{H}-12z\dot{H}+6\alpha z\ddot{\phi}.$
(30)

In this way we obtain the final expression of the autonomous dynamical system:

$x^{\prime}=\frac{1}{y\left(54\alpha^{2}z^{2}\left(1-x^{2}\epsilon\right)^{3/2}+6z\epsilon\left(2\sqrt{3}\alpha
xy-1\right)+\epsilon\right)}\cdot\\\ \Big{(}-9\sqrt{3}\alpha
x^{2}y^{2}z\epsilon
w_{m}-162\alpha^{2}xyz^{2}w_{m}\sqrt{1-x^{2}\epsilon}-18\sqrt{3}\alpha
x^{2}z^{2}\epsilon w_{m}\sqrt{1-x^{2}\epsilon}+18\sqrt{3}\alpha
z^{2}w_{m}\sqrt{1-x^{2}\epsilon}+9\sqrt{3}\alpha x^{2}z\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}\\\ -9\sqrt{3}\alpha
zw_{m}\sqrt{1-x^{2}\epsilon}+162\alpha^{2}x^{3}yz^{2}\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}+9\sqrt{3}\alpha y^{2}zw_{m}+45\sqrt{3}\alpha
x^{4}y^{2}z\epsilon^{2}-36\alpha\lambda
x^{3}y^{3}z\epsilon-18x^{3}yz\epsilon^{2}+3x^{3}y\epsilon^{2}\\\
-54\sqrt{3}\alpha^{3}x^{2}y^{2}z^{2}\sqrt{1-x^{2}\epsilon}+27\sqrt{3}\alpha^{2}\lambda
x^{2}y^{2}z^{2}\sqrt{1-x^{2}\epsilon}-54\sqrt{3}\alpha
x^{2}y^{2}z\epsilon+6\sqrt{3}\lambda x^{2}y^{2}z\epsilon-\sqrt{3}\lambda
x^{2}y^{2}\epsilon+36\alpha^{2}xyz^{2}\sqrt{1-x^{2}\epsilon}\\\
+6\sqrt{3}\alpha x^{2}z^{2}\epsilon\sqrt{1-x^{2}\epsilon}-6\sqrt{3}\alpha
z^{2}\sqrt{1-x^{2}\epsilon}+5\sqrt{3}\alpha
x^{2}z\epsilon\sqrt{1-x^{2}\epsilon}-5\sqrt{3}\alpha
z\sqrt{1-x^{2}\epsilon}+54\sqrt{3}\alpha^{3}x^{4}y^{2}z^{2}\epsilon\sqrt{1-x^{2}\epsilon}\\\
-27\sqrt{3}\alpha^{2}\lambda
x^{4}y^{2}z^{2}\epsilon\sqrt{1-x^{2}\epsilon}-36\alpha^{2}x^{3}yz^{2}\epsilon\sqrt{1-x^{2}\epsilon}+36\alpha\lambda
xy^{3}z+18xyz\epsilon-3xy\epsilon+\sqrt{3}\lambda y^{2}+9\sqrt{3}\alpha
y^{2}z-6\sqrt{3}\lambda y^{2}z\Big{)},$ (31)

$y^{\prime}=\frac{1}{2\sqrt{1-x^{2}\epsilon}\left(54\alpha^{2}z^{2}\left(1-x^{2}\epsilon\right)^{3/2}+6z\epsilon\left(2\sqrt{3}\alpha
xy-1\right)+\epsilon\right)}\cdot\\\ \Big{(}18\sqrt{3}\alpha xy^{2}z\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}-6yz\epsilon w_{m}\sqrt{1-x^{2}\epsilon}+3y\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}-3y^{3}\epsilon w_{m}-54\sqrt{3}\alpha^{2}\lambda
x^{5}y^{2}z^{2}\epsilon^{2}+72\alpha^{2}x^{4}yz^{2}\epsilon^{2}\\\
+108\sqrt{3}\alpha^{2}\lambda
x^{3}y^{2}z^{2}\epsilon+18\alpha^{2}x^{2}y^{3}z\epsilon\sqrt{1-x^{2}\epsilon}-63\alpha\lambda
x^{2}y^{3}z\epsilon\sqrt{1-x^{2}\epsilon}+18\alpha\lambda
y^{3}z\sqrt{1-x^{2}\epsilon}+3x^{2}y^{3}\epsilon^{2}-6\sqrt{3}\alpha
xy^{2}z\epsilon\sqrt{1-x^{2}\epsilon}\\\ +6\sqrt{3}\lambda
xy^{2}z\epsilon\sqrt{1-x^{2}\epsilon}-\sqrt{3}\lambda
xy^{2}\epsilon\sqrt{1-x^{2}\epsilon}-144\alpha^{2}x^{2}yz^{2}\epsilon-6yz\epsilon\sqrt{1-x^{2}\epsilon}+3y\epsilon\sqrt{1-x^{2}\epsilon}+18\sqrt{3}\alpha
x^{3}y^{2}z\epsilon^{2}\sqrt{1-x^{2}\epsilon}\\\ -54\sqrt{3}\alpha^{2}\lambda
xy^{2}z^{2}-3y^{3}\epsilon+72\alpha^{2}yz^{2}\Big{)},$ (32)

$z^{\prime}=\frac{1}{\sqrt{1-x^{2}\epsilon}\left(54\alpha^{2}z^{2}\left(1-x^{2}\epsilon\right)^{3/2}+6z\epsilon\left(2\sqrt{3}\alpha
xy-1\right)+\epsilon\right)}\cdot\\\ \Big{(}-36\sqrt{3}\alpha xyz^{2}\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}+12z^{2}\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}-6z\epsilon
w_{m}\sqrt{1-x^{2}\epsilon}+6y^{2}z\epsilon
w_{m}+54\sqrt{3}\alpha^{3}x^{5}yz^{3}\epsilon^{2}-144\alpha^{2}x^{4}z^{3}\epsilon^{2}\\\
-108\sqrt{3}\alpha^{3}x^{3}yz^{3}\epsilon+54\alpha\lambda
x^{2}y^{2}z^{2}\epsilon\sqrt{1-x^{2}\epsilon}-36\alpha\lambda
y^{2}z^{2}\sqrt{1-x^{2}\epsilon}-6x^{2}y^{2}z\epsilon^{2}+6\sqrt{3}\alpha
xyz^{2}\epsilon\sqrt{1-x^{2}\epsilon}+\sqrt{3}\alpha
xyz\epsilon\sqrt{1-x^{2}\epsilon}\\\
+288\alpha^{2}x^{2}z^{3}\epsilon+12z^{2}\epsilon\sqrt{1-x^{2}\epsilon}-6z\epsilon\sqrt{1-x^{2}\epsilon}-36\sqrt{3}\alpha
x^{3}yz^{2}\epsilon^{2}\sqrt{1-x^{2}\epsilon}+54\sqrt{3}\alpha^{3}xyz^{3}+6y^{2}z\epsilon-144\alpha^{2}z^{3}\Big{)}.$
(33)

Figure 1: The figure describes a region in the phase space structure where the
A cosmological solution is saddle, in the case where $w_{m}=-0.0001$. Figure
2: The figure describes the value of the matter density parameter $\Omega_{m}$
as a function of dark matter equation of state parameter $w_{m}$ and the
strength of the potential energy, encoded into the value of the $\lambda$
coefficient. Figure 3: The figure describes the value of the matter density
parameter $\Omega_{m}$ as a function of the strength of the potential energy,
encoded into the value of the $\lambda$ coefficient. In this figure we have
considered $w_{m}=-0.0001$. Figure 4: The figure describes the evolution
towards the A critical point in the $\\{x,y\\}$ plane
($\alpha=1,w_{m}=-0.01,\lambda=8.9,\epsilon=1$). Figure 5: The evolution
towards the A cosmological solution in the $\\{x,z\\}$ plane. Figure 6: The
evolution towards A in the $\\{y,z\\}$ plane Figure 7: The figure describes
the evolution in the 3D space for various initial conditions, towards the A
cosmological solution.  Figure 8: The figure describes a specific region where
the de–Sitter critical point B is stable. The figure takes into account also
the existence conditions, implying that the $y$ component is real and
positive.

For our cosmological system we have obtained two classes of critical points.
The first critical point is located at the following coordinates:

$A=\Big{[}x=\sqrt{1+w_{m}},y=\frac{\sqrt{3}\sqrt{1+w_{m}}}{\lambda},z=0\Big{]}.$
(34)

At this solution we note that the dark matter equation of state dictates the
value of the kinetic term, while for the potential component we have an
influence due to the steepness of the potential, embedded into the value of
the $\lambda$ coefficient. The effective equation of state is equal to

$w_{\bf{eff}}=w_{m},$ (35)

with the matter density parameter

$\Omega_{m}=1-\frac{3(1+w_{m})}{\sqrt{-w_{m}}\lambda^{2}}.$ (36)

Due to the existence conditions the matter density parameter should be real in
the $[0,1]$ interval, implying that the $w_{m}$ parameter is slightly negative
and close to zero, a value not ruled out by astrophysical observations. This
cosmological solution represents a matter dominated era, a critical point
which appeared also in various dynamical studies [65]. This cosmological
solution is similar to the one found in the minimal coupling case [73]. For
this cosmological solution we have obtained the following eigenvalues:

$\Big{[}\frac{3(\alpha-2\lambda)\left(w_{m}+1\right)}{\lambda},\frac{3}{4}\left(w_{m}-1\pm\frac{\Pi}{\Psi}\right)\Big{]},$
(37)

where

$\Pi=\sqrt{-\lambda^{4}w_{m}^{3}\left(w_{m}+1\right)\left(\lambda^{2}+17\lambda^{2}w_{m}^{2}+14\lambda^{2}w_{m}+48\left(-w_{m}\right){}^{5/2}+48\sqrt{-w_{m}}+96\sqrt{-w_{m}}w_{m}\right)},$
(38) $\Psi=\lambda^{3}w_{m}\sqrt{-w_{m}\left(w_{m}+1\right)}.$ (39)

From the above expressions we can note that this solution has a high
sensitivity to the values of various coefficients, $w_{m},\lambda,\alpha$. For
this solution, we have displayed in Fig. 1 some possible intervals where the
corresponding behavior is associated to a saddle dynamics. In the
implementation of this figure we have used an interpolation method based on a
neural network approach, with an initial grid considered as the training set.
The architecture of the neural network is based on four internal linear
layers, separated through a specific activation function. After training, the
neural network can be considered for the interpolation, obtaining possible
regions in the phase space structure where the $A$ cosmological solution has a
saddle behavior. For example, in the case where $\Omega_{m}\approx 0.70$, if
we set $\lambda=32$, $w_{m}=-0.0001$, and $\alpha=100$, we obtain a saddle
cosmological solution which can mimic the matter dominated era.

In Fig. 2 we have the dependence of the matter density parameter $\Omega_{m}$
in the $\\{w_{m},\lambda\\}$ plane. It can be seen that the matter density
parameter can span easily by fine–tuning the $[0,1]$ interval, being
compatible to various astrophysical observations. Furthermore, in Fig. 3 we
have displayed the value of the matter density parameter $\Omega_{m}$ as a
function of $\lambda$, the strength of the potential energy. In this case we
have fixed $w_{m}=-0.0001$. Lastly, for this solution we have analyzed the
phase space structure, displaying various aspects of the numerical evolution
in Figs. 4–7.

The second cosmological solution is located at the coordinates:

$B=\Big{[}x=0,y=\frac{\sqrt{2\alpha}}{\sqrt{2\alpha-\lambda}},z=\frac{\lambda}{2(\lambda-2\alpha)}\Big{]},$
(40)

a de–Sitter epoch where the effective equation of state corresponds to a
cosmological constant,

$w_{\bf{eff}}=-1.$ (41)

The matter density parameter is $\Omega_{m}=0$, a solution where the scalar
field completely dominates in terms of density parameters. For this solution,
we can note that the field is at rest, without any kinetic energy. The
potential energy variable $y$ and the coupling component which is encoded into
the $z$ variable are both influenced by the strength of the coupling function,
displayed by the $\alpha$ coefficient, and the steepness of the exponential
potential. From a dynamical point of view we have obtained the following
eigenvalues,

$\Bigg{[}-3\left(w_{m}+1\right),\frac{1}{\Theta}(\Delta\pm\sqrt{3}\sqrt{\Xi})\Bigg{]},$
(42)

where

$\Theta=2(2\alpha-\lambda)^{3}\left(\alpha^{2}\left(27\lambda^{2}+8\right)+4\alpha\lambda-4\lambda^{2}\right)^{2},$
(43)

$\Delta=24\alpha^{7}\left(27\lambda^{2}+8\right)^{2}-12\alpha^{6}\lambda\left(2187\lambda^{4}+864\lambda^{2}+64\right)+18\alpha^{5}\lambda^{2}\left(729\lambda^{4}-288\lambda^{2}-128\right)-3\alpha^{4}\lambda^{3}\left(729\lambda^{4}-3456\lambda^{2}-640\right)\\\
+\alpha^{3}\left(480\lambda^{4}-4536\lambda^{6}\right)+72\alpha^{2}\lambda^{5}\left(9\lambda^{2}-14\right)+384\alpha\lambda^{6}-48\lambda^{7},$
(44)

$\Xi=-(\lambda-2\alpha)^{6}\left(\alpha^{2}\left(27\lambda^{2}+8\right)+4\alpha\lambda-4\lambda^{2}\right)^{3}\left(32\alpha^{3}\lambda-\alpha^{2}\left(161\lambda^{2}+24\right)+4\alpha\lambda\left(8\lambda^{2}-3\right)+12\lambda^{2}\right).$
(45)

For the B critical point we have analyzed the values of the corresponding
eigenvalues and the associated dynamical behavior. In Fig. 8 we have
considered an interpolation method based on a neural network approach,
determining a possible region where the $B$ cosmological solution has a stable
dynamical behavior. Note that the corresponding dynamical behavior is
influenced by the strength of the coupling with the cubic component (which
encodes geometrical effects due to the third order contractions of the Riemann
tensor), and the steepness of the potential energy.

##  IV  Beyond the exponential potential

In this section we shall investigate the structure and properties of the phase
space in the case where the potential energy term is beyond the usual
exponential case. To this regard, we shall consider that the potential is
represented by an inverse hyperbolic sine,

$V(\phi)=V_{0}sinh^{-\xi}(\chi\phi),$ (46)

where $V_{0},\xi,\chi$ are constant parameters. This potential has been
considered in various cosmological models [75, 65, 76, 77], representing a
viable function. In order to study such a potential energy, we need to
introduce another variable,

$\lambda=-\frac{1}{V(\phi)}\frac{dV(\phi)}{d\phi}.$ (47)

Then, the dynamical system is a four dimensional system, where we have to add
the following differential equation,

$\lambda^{\prime}=-\sqrt{3}\lambda^{2}xy(\Gamma-1),$ (48)

where $\Gamma$ is defined as follows:

$\Gamma=\frac{V(\phi)\frac{d^{2}V(\phi)}{d\phi^{2}}}{(\frac{dV(\phi)}{d\phi})^{2}}.$
(49)

Taking into account the inverse hyperbolic sine decomposition, we have:

$\Gamma=1+\frac{1}{\xi}-\frac{\xi\chi^{2}}{\lambda^{2}},$ (50)

and the last dynamical equation reduces to:

$\lambda^{\prime}=-xy\sqrt{3}(-\xi\chi^{2}+\frac{\lambda^{2}}{\xi}).$ (51)

In what follows we shall present the corresponding critical points and the
specific fundamental properties obtained in the case where the potential is
represented by a hyperbolic sine function. The first critical point found in
our analysis is located at the following coordinates:

$C=\Big{[}x=0,y,z=\frac{1}{2}-\frac{y^{2}}{2},\lambda=\frac{2\alpha(y^{2}-1)}{y^{2}}\Big{]}.$
(52)

We can note that for this solution the field is at rest, without any kinetic
energy, while the $y$ variable associated to the potential energy represents
an independent parameter. The specific value of the last two variables
$(z,\lambda)$, which are associated to the non–minimal coupling function and
the steepness of the potential, respectively, are influenced by the potential
component. From a physical point of view this cosmological solution describes
a de–Sitter epoch $(\Omega_{m}=0,w_{\bf{eff}}=-1)$. For this solution, the
general form of the eigenvalues is too complex to be written here. However, if
we set $y=1$, we obtain a simpler form of the corresponding eigenvalues:

$\Big{[}0,3\left(w_{m}+1\right),\frac{1}{2}\left(-\sqrt{12\xi\chi^{2}+9}-3\right),\frac{1}{2}\left(\sqrt{12\xi\chi^{2}+9}-3\right)\Big{]},$
(53)

describing a non hyperbolic solution which is saddle. In this case the
dynamical aspects are influenced by the values of the $\xi$ and $\chi$
coefficients.

The second cosmological solution is represented by the following critical
point,

$D^{\pm}=\Big{[}x=\pm\sqrt{1+w_{m}},y=\frac{\sqrt{3}\sqrt{w_{m}+1}}{\xi\chi},z=0,\lambda=\xi\chi\Big{]}.$
(54)

This solution corresponds to a matter dominated era
$(\Omega_{m}=1-\frac{3\left(w_{m}+1\right)}{\xi^{2}\chi^{2}\sqrt{-w_{m}}},w_{\bf{eff}}=w_{m})$,
having the following eigenvalues (for the $D^{+}$ solution):

$\Bigg{[}-\frac{6\left(w_{m}+1\right)}{\xi},\frac{3\left(w_{m}+1\right)(\alpha-2\xi\chi)}{\xi\chi},\\\
\frac{3}{4}\left(\pm\frac{\sqrt{\xi^{6}\chi^{4}w_{m}^{2}\left(w_{m}+1\right)\left(17\xi^{2}\chi^{2}w_{m}^{2}+14\xi^{2}\chi^{2}w_{m}+48\left(-w_{m}\right){}^{5/2}+48\sqrt{-w_{m}}+96\sqrt{-w_{m}}w_{m}+\xi^{2}\chi^{2}\right)}}{\xi^{4}\chi^{3}w_{m}\sqrt{w_{m}+1}}+w_{m}-1\right)\Bigg{]}.$
(55)

Hence, this physical point is similar to the $B$ solution discussed earlier in
the case where the potential is represented by an exponential function. From a
dynamical perspective if we consider that the dark matter pressure is negative
but very close to zero and $\xi<0$, then the first eigenvalue becomes
positive, implying that this solution cannot be stable. Hence, it is either
saddle or unstable, depending of the values of the specific coefficients
$\\{\alpha,\xi,\chi\\}$. By imposing that the second eigenvalue is negative
and $w_{m}=-0.0001$, we have obtained some constraints to the specific
coefficients,

$\xi<0\land\left(\chi<0\lor\alpha>0\lor\chi>\frac{0.5\alpha}{\xi}\right)\land\left(\chi<\frac{0.5\alpha}{\xi}\lor\alpha\leq
0\lor\chi>0\right),$ (56)

where the matter dominated era appears as a saddle cosmological solution in
the phase space structure. As can be noted, the structure of the phase space
in the case where the potential energy term is represented by an inverse
hyperbolic sine is similar to the exponential potential case.

##  V Summary and Conclusions

In the present manuscript we have proposed a novel cosmological model, by
adding to the Einstein–Hilbert Lagrangian a tachyonic field non–minimally
coupled with a topological invariant constructed with specific third order
contractions of the Riemann tensor. In this scenario the dark energy component
is represented by the tachyonic field which depends on the cosmic time. After
proposing the action corresponding to this cosmological system, we have
obtained the basic equations which describe the evolution of such a
theoretical model, obtained by applying the variational principle. To this
regard, the Klein–Gordon equation is obtained by varying the action with
respect to the tachyonic field, while the modified Friedmann relations are
deduced by varying the inverse metric. Since the present model does not take
into account an interaction between the tachyonic field and the matter
component, the standard continuity equation is also satisfied.

The physical features of the current cosmological system are investigated by
adopting the dynamical system analysis, an important tool used in the study of
various modified gravity theories. In the present paper we have analyzed the
phase space structure and properties in the case of an exponential coupling
function. Furthermore, for the potential energy we have considered two
specific cases, the exponential potential and a distinct potential, the
hyperbolic potential. In the case of an exponential potential the structure of
the phase space has three dimensions, having two types of cosmological
solutions. The first type is represented by the de–Sitter epoch, a
cosmological solution where the tachyonic field acts as a cosmological
constant, with a constant equation of state. At this critical point the value
of the coupling coefficient which encodes specific interactions with the
topological cubic invariant is affecting the dynamical consequences.

To this regard, we have obtained specific constraints for the coupling
coefficients where the attractor behavior is attained in the distant future.
The second type of cosmological solutions is represented by the matter
dominated era, an epoch where the total (effective) equation of state of the
cosmological system corresponds to the matter component. From a physical point
of view this solution is viable only if we take into account that the dark
matter fluid has a negative equation of state, slightly close to zero. This
implies that the dark matter component corresponds to an exotic fluid with a
negative pressure. From a theoretical perspective such a solution has been
also found in different tachyonic dark energy models [65]. In the case of the
second solution the dynamical features have been investigated, revealing some
values of the coupling coefficients where the matter epoch corresponds to a
saddle behavior, compatible with the recent evolution at the large scale
structure. Lastly, we have also considered that the potential energy term is
represented by an inverse hyperbolic sine function, discussing the phase space
structure and the corresponding dynamical effects. Finally, due to the
presented arguments we can note that the present cosmological setup represents
a viable alternative theory which can explain the evolution of the Universe –
the matter dominated epoch and the dark energy phenomenon, constituting a
feasible theoretical framework, at least at the level of background dynamics.

###### Acknowledgements.

This work was supported by a grant of the Romanian Ministry of Research,
Innovation and Digitalization, CNCS - UEFISCDI, project number PN-III-P4-ID-
PCE-2020-1142, within PNCDI III. For this work we have considered various
computations in Wolfram Mathematica [78] and xAct [79].

## References

  * Will [2014] C. M. Will, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
  * Riess  _et al._ [2021] A. G. Riess, S. Casertano, W. Yuan, J. B. Bowers, L. Macri, J. C. Zinn, and D. Scolnic, Astrophys. J. Lett. 908, L6 (2021), arXiv:2012.08534 [astro-ph.CO] .
  * Asgari  _et al._ [2021] M. Asgari _et al._ (KiDS), Astron. Astrophys. 645, A104 (2021), arXiv:2007.15633 [astro-ph.CO] .
  * Di Valentino  _et al._ [2021] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, Class. Quant. Grav. 38, 153001 (2021), arXiv:2103.01183 [astro-ph.CO] .
  * Vagnozzi [2020] S. Vagnozzi, Phys. Rev. D 102, 023518 (2020), arXiv:1907.07569 [astro-ph.CO] .
  * Nojiri and Odintsov [2006] S. Nojiri and S. D. Odintsov, eConf C0602061, 06 (2006), arXiv:hep-th/0601213 .
  * Bamba  _et al._ [2012] K. Bamba, S. Capozziello, S. Nojiri, and S. D. Odintsov, Astrophys. Space Sci. 342, 155 (2012), arXiv:1205.3421 [gr-qc] .
  * Nojiri and Odintsov [2003] S. Nojiri and S. D. Odintsov, Phys. Rev. D 68, 123512 (2003), arXiv:hep-th/0307288 .
  * Nojiri and Odintsov [2011] S. Nojiri and S. D. Odintsov, Phys. Rept. 505, 59 (2011), arXiv:1011.0544 [gr-qc] .
  * Nojiri  _et al._ [2017] S. Nojiri, S. D. Odintsov, and V. K. Oikonomou, Phys. Rept. 692, 1 (2017), arXiv:1705.11098 [gr-qc] .
  * Frusciante and Perenon [2020] N. Frusciante and L. Perenon, Phys. Rept. 857, 1 (2020), arXiv:1907.03150 [astro-ph.CO] .
  * Copeland  _et al._ [2006] E. J. Copeland, M. Sami, and S. Tsujikawa, Int. J. Mod. Phys. D 15, 1753 (2006), arXiv:hep-th/0603057 .
  * Capozziello and De Laurentis [2011] S. Capozziello and M. De Laurentis, Phys. Rept. 509, 167 (2011), arXiv:1108.6266 [gr-qc] .
  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * González-Morán  _et al._ [2021] A. L. González-Morán, R. Chávez, E. Terlevich, R. Terlevich, D. Fernández-Arenas, F. Bresolin, M. Plionis, J. Melnick, S. Basilakos, and E. Telles, Mon. Not. Roy. Astron. Soc. 505, 1441 (2021), arXiv:2105.04025 [astro-ph.CO] .
  * Bermejo-Climent  _et al._ [2021] J. R. Bermejo-Climent, M. Ballardini, F. Finelli, D. Paoletti, R. Maartens, J. A. Rubiño Martín, and L. Valenziano, Phys. Rev. D 103, 103502 (2021), arXiv:2106.05267 [astro-ph.CO] .
  * Cao  _et al._ [2021] S. Cao, J. Ryan, and B. Ratra, Mon. Not. Roy. Astron. Soc. 504, 300 (2021), arXiv:2101.08817 [astro-ph.CO] .
  * Yang  _et al._ [2021] W. Yang, E. Di Valentino, S. Pan, Y. Wu, and J. Lu, Mon. Not. Roy. Astron. Soc. 501, 5845 (2021), arXiv:2101.02168 [astro-ph.CO] .
  * Muccino  _et al._ [2021] M. Muccino, L. Izzo, O. Luongo, K. Boshkayev, L. Amati, M. Della Valle, G. B. Pisani, and E. Zaninoni, Astrophys. J. 908, 181 (2021), arXiv:2012.03392 [astro-ph.CO] .
  * Bueno and Cano [2016a] P. Bueno and P. A. Cano, Phys. Rev. D 94, 104005 (2016a), arXiv:1607.06463 [hep-th] .
  * Hennigar and Mann [2017] R. A. Hennigar and R. B. Mann, Phys. Rev. D 95, 064055 (2017), arXiv:1610.06675 [hep-th] .
  * Bueno and Cano [2016b] P. Bueno and P. A. Cano, Phys. Rev. D 94, 124051 (2016b), arXiv:1610.08019 [hep-th] .
  * Hennigar  _et al._ [2018] R. A. Hennigar, M. B. J. Poshteh, and R. B. Mann, Phys. Rev. D 97, 064041 (2018), arXiv:1801.03223 [gr-qc] .
  * Mehdizadeh and Ziaie [2019] M. R. Mehdizadeh and A. H. Ziaie, Mod. Phys. Lett. A 35, 2050017 (2019), arXiv:1903.10907 [gr-qc] .
  * Cisterna  _et al._ [2020] A. Cisterna, N. Grandi, and J. Oliva, Phys. Lett. B 805, 135435 (2020), arXiv:1811.06523 [hep-th] .
  * Poshteh and Mann [2019] M. B. J. Poshteh and R. B. Mann, Phys. Rev. D 99, 024035 (2019), arXiv:1810.10657 [gr-qc] .
  * Bueno  _et al._ [2018] P. Bueno, P. A. Cano, and A. Ruipérez, JHEP 03, 150, arXiv:1802.00018 [hep-th] .
  * Konoplya  _et al._ [2020] R. A. Konoplya, A. F. Zinhailo, and Z. Stuchlik, Phys. Rev. D 102, 044023 (2020), arXiv:2006.10462 [gr-qc] .
  * Adair  _et al._ [2020] C. Adair, P. Bueno, P. A. Cano, R. A. Hennigar, and R. B. Mann, Phys. Rev. D 102, 084001 (2020), arXiv:2004.09598 [gr-qc] .
  * Giri and Rudra [2021] K. Giri and P. Rudra, (2021), arXiv:2107.12417 [astro-ph.CO] .
  * Marciu [2021] M. Marciu, (2021), arXiv:2103.08420 [gr-qc] .
  * Erices  _et al._ [2019] C. Erices, E. Papantonopoulos, and E. N. Saridakis, Phys. Rev. D 99, 123527 (2019), arXiv:1903.11128 [gr-qc] .
  * Marciu [2020a] M. Marciu, Phys. Rev. D 101, 103534 (2020a), arXiv:2003.06403 [gr-qc] .
  * Quiros  _et al._ [2020] I. Quiros, R. García-Salcedo, T. Gonzalez, J. L. M. Martínez, and U. Nucamendi, Phys. Rev. D 102, 044018 (2020), arXiv:2003.10516 [gr-qc] .
  * Kord Zangeneh and Kazemi [2020] M. Kord Zangeneh and A. Kazemi, Eur. Phys. J. C 80, 794 (2020), arXiv:2003.04458 [hep-th] .
  * Frassino and Rocha [2020] A. M. Frassino and J. V. Rocha, Phys. Rev. D 102, 024035 (2020), arXiv:2002.04071 [hep-th] .
  * Burger  _et al._ [2020] D. J. Burger, W. T. Emond, and N. Moynihan, Phys. Rev. D 101, 084009 (2020), arXiv:1910.11618 [hep-th] .
  * Cano and Pereñiguez [2020] P. A. Cano and D. Pereñiguez, Phys. Rev. D 101, 044016 (2020), arXiv:1910.10721 [hep-th] .
  * Emond and Moynihan [2019] W. T. Emond and N. Moynihan, JHEP 12, 019, arXiv:1905.08213 [hep-th] .
  * Arciniega  _et al._ [2020a] G. Arciniega, J. D. Edelstein, and L. G. Jaime, Phys. Lett. B 802, 135272 (2020a), arXiv:1810.08166 [gr-qc] .
  * Arciniega  _et al._ [2020b] G. Arciniega, P. Bueno, P. A. Cano, J. D. Edelstein, R. A. Hennigar, and L. G. Jaime, Phys. Lett. B 802, 135242 (2020b), arXiv:1812.11187 [hep-th] .
  * Edelstein  _et al._ [2020] J. D. Edelstein, D. Vázquez Rodríguez, and A. Vilar López, JCAP 12, 040, arXiv:2006.10007 [hep-th] .
  * Arciniega  _et al._ [2019] G. Arciniega, P. Bueno, P. A. Cano, J. D. Edelstein, R. A. Hennigar, and L. G. Jaime, Int. J. Mod. Phys. D 28, 1944008 (2019).
  * Quiros  _et al._ [2021] I. Quiros, R. De Arcia, R. García-Salcedo, T. Gonzalez, F. X. Linares Cedeño, and U. Nucamendi, Phys. Rev. D 103, 064043 (2021), arXiv:2007.06111 [gr-qc] .
  * Cano  _et al._ [2021] P. A. Cano, K. Fransen, and T. Hertog, Phys. Rev. D 103, 103531 (2021), arXiv:2011.13933 [hep-th] .
  * Edelstein  _et al._ [2021] J. D. Edelstein, R. B. Mann, D. V. Rodríguez, and A. Vilar López, JHEP 01, 029, arXiv:2007.07651 [hep-th] .
  * Marciu [2020b] M. Marciu, Phys. Rev. D 102, 023517 (2020b), arXiv:2004.07120 [gr-qc] .
  * Pookkillath  _et al._ [2020] M. C. Pookkillath, A. De Felice, and A. A. Starobinsky, JCAP 07, 041, arXiv:2004.03912 [gr-qc] .
  * Jiménez and Jiménez-Cano [2021] J. B. Jiménez and A. Jiménez-Cano, JCAP 01, 069, arXiv:2009.08197 [gr-qc] .
  * Mazumdar  _et al._ [2001] A. Mazumdar, S. Panda, and A. Perez-Lorenzana, Nucl. Phys. B 614, 101 (2001), arXiv:hep-ph/0107058 .
  * Sen [2005] A. Sen, Phys. Scripta T 117, 70 (2005), arXiv:hep-th/0312153 .
  * Piao  _et al._ [2002] Y.-S. Piao, R.-G. Cai, X.-m. Zhang, and Y.-Z. Zhang, Phys. Rev. D 66, 121301 (2002), arXiv:hep-ph/0207143 .
  * Bagla  _et al._ [2003] J. S. Bagla, H. K. Jassal, and T. Padmanabhan, Phys. Rev. D 67, 063504 (2003), arXiv:astro-ph/0212198 .
  * Shchigolev and Rotova [2012a] V. K. Shchigolev and M. P. Rotova, Mod. Phys. Lett. A 27, 1250086 (2012a), arXiv:1203.5030 [gr-qc] .
  * Shchigolev and Rotova [2012b] V. K. Shchigolev and M. P. Rotova, Grav. Cosmol. 18, 88 (2012b), arXiv:1105.4536 [gr-qc] .
  * Avelino  _et al._ [2011] P. P. Avelino, L. Losano, and J. J. Rodrigues, Phys. Lett. B 699, 10 (2011), arXiv:1103.1384 [astro-ph.CO] .
  * Gibbons [2002] G. W. Gibbons, Phys. Lett. B 537, 1 (2002), arXiv:hep-th/0204008 .
  * Mukohyama [2002] S. Mukohyama, Phys. Rev. D 66, 024009 (2002), arXiv:hep-th/0204084 .
  * Teixeira  _et al._ [2019] E. M. Teixeira, A. Nunes, and N. J. Nunes, Phys. Rev. D 100, 043539 (2019), arXiv:1903.06028 [gr-qc] .
  * Gibbons [2003] G. W. Gibbons, Class. Quant. Grav. 20, S321 (2003), arXiv:hep-th/0301117 .
  * Copeland  _et al._ [2005] E. J. Copeland, M. R. Garousi, M. Sami, and S. Tsujikawa, Phys. Rev. D 71, 043003 (2005), arXiv:hep-th/0411192 .
  * Banijamali and Solbi [2017] A. Banijamali and M. Solbi, Gen. Rel. Grav. 49, 103 (2017).
  * Sen [2002] A. Sen, Mod. Phys. Lett. A 17, 1797 (2002), arXiv:hep-th/0204143 .
  * Banijamali and Fazlpour [2012] A. Banijamali and B. Fazlpour, Astrophys. Space Sci. 342, 229 (2012), arXiv:1206.3580 [physics.gen-ph] .
  * Bahamonde  _et al._ [2019] S. Bahamonde, M. Marciu, and J. L. Said, Eur. Phys. J. C 79, 324 (2019), arXiv:1901.04973 [gr-qc] .
  * Rezaei Akbarieh and Izadi [2019] A. Rezaei Akbarieh and Y. Izadi, Eur. Phys. J. C 79, 366 (2019), arXiv:1812.06649 [gr-qc] .
  * Motavalli  _et al._ [2016] H. Motavalli, A. R. Akbarieh, and M. Nasiry, J. Exp. Theor. Phys. 123, 33 (2016).
  * Fazlpour and Banijamali [2015] B. Fazlpour and A. Banijamali, Adv. High Energy Phys. 2015, 283273 (2015), arXiv:1408.0203 [gr-qc] .
  * Banijamali [2014] A. Banijamali, Adv. High Energy Phys. 2014, 631630 (2014).
  * Otalora [2013] G. Otalora, Phys. Rev. D 88, 063505 (2013), arXiv:1305.5896 [gr-qc] .
  * Bahamonde  _et al._ [2021] S. Bahamonde, K. F. Dialektopoulos, C. Escamilla-Rivera, G. Farrugia, V. Gakis, M. Hendry, M. Hohmann, J. L. Said, J. Mifsud, and E. Di Valentino, (2021), arXiv:2106.13793 [gr-qc] .
  * Abramo and Finelli [2003] L. R. W. Abramo and F. Finelli, Phys. Lett. B 575, 165 (2003), arXiv:astro-ph/0307208 .
  * Quiros  _et al._ [2010] I. Quiros, T. Gonzalez, D. Gonzalez, and Y. Napoles, Class. Quant. Grav. 27, 215021 (2010), arXiv:0906.2617 [gr-qc] .
  * Fang and Lu [2010] W. Fang and H.-Q. Lu, Eur. Phys. J. C 68, 567 (2010), arXiv:1007.2330 [hep-th] .
  * Roy and Bhadra [2018] N. Roy and N. Bhadra, JCAP 06, 002, arXiv:1710.05968 [gr-qc] .
  * Urena-Lopez and Matos [2000] L. A. Urena-Lopez and T. Matos, Phys. Rev. D 62, 081302 (2000), arXiv:astro-ph/0003364 .
  * Sahni and Starobinsky [2000] V. Sahni and A. A. Starobinsky, Int. J. Mod. Phys. D 9, 373 (2000), arXiv:astro-ph/9904398 .
  * [78] Wolfram Research, https://www.wolfram.com/mathematica, Mathematica.
  * [79] J. M. Martin-Garcia, xAct: Efficient tensor computer algebra for the Wolfram Language.