# Geometric Perfect Fluids and Dark Side of the Universe

Metin Gürses a, Yaghoub Heydarzade a, and Çetin Şentürk b

aDepartment of Mathematics, Faculty of Sciences, Bilkent University, 06800
Ankara, Turkey
bDepartment of Aeronautical Engineering, University of Turkish Aeronautical
Association, 06790 Ankara, Turkey email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
(August 28, 2024)

###### Abstract

Recently we showed that in FLRW cosmology, the contribution from higher
curvature terms in any generic metric gravity theory to the energy-momentum
tensor is of the perfect fluid form. Such a geometric perfect fluid can be
interpreted as a fluid remaining from the beginning of the universe where the
string theory is thought to be effective. Just a short time after the
beginning of the Universe, it is known that the Einstein-Hilbert action is
assumed to be modified by adding all possible curvature invariants. We propose
that the observed late-time accelerating expansion of the Universe can be
solely driven by this geometric fluid. To support our claim, we specifically
study the quadratic gravity field equations in $D$-dimensions. We show that
the field equations of this theory for the FLRW metric possess a geometric
perfect fluid source containing two critical parameters $\sigma_{1}$ and
$\sigma_{2}$. To analyze this theory concerning its parameter space
$(\sigma_{1},\sigma_{2})$, we obtain the general second-order nonlinear
differential equation governing the late-time dynamics of the deceleration
parameter $q$. Hence using some present-day cosmological data as our initial
conditions, our findings for the $\sigma_{2}=0$ case are as follows: $(i)$ In
order to have a positive energy density for the geometric fluid $\rho_{g}$,
the parameter $\sigma_{1}$ must be negative for all dimensions up to $D=11$,
$(ii)$ For a suitable choice of $\sigma_{1}$, the deceleration parameter
experiences signature changes in the past and future, and in the meantime it
lies within a negative range which means that the current observed accelerated
expansion phase of the Universe can be driven solely by the curvature of the
spacetime, $(iii)$ $q$ experiences a signature change and as the dimension $D$
of spacetime increases, this signature change happens at earlier and later
times, in the past and future, respectively.

###### Contents

  1. 1 Introduction
  2. 2 Generic Gravity Field Equations in Perfect Fluid Form
  3. 3 Cosmological Parameters in Generic Gravity Theories
  4. 4 Quadratic Gravity and Criticality in $D$-Dimensions
  5. 5 Some Particular Cosmological Solutions in Quadratic Gravity
    1. 5.1 Solutions with $\sigma_{2}=0$
      1. 5.1.1 Exponential Solutions
      2. 5.1.2 A More General Solution
      3. 5.1.3 Approximate Solutions Near Criticality with $\sigma_{2}=0$
    2. 5.2 Solutions with $\sigma_{1}=0$
  6. 6 Conclusion

##  1 Introduction

Although General Relativity (GR) has been immensely successful in explaining a
wide range of gravitational phenomena, there are certain observations that
have motivated researchers to consider modifications to the theory. Two of
these motivations are: (i) Theoretical consistency of GR: Modifications to
gravity theories can arise from attempts to reconcile GR with other
fundamental theories, such as quantum gravity or string theory. GR is assumed
to be a low-energy approximation of a more complete theory where the effective
action includes higher-curvature terms in addition to the usual Einstein-
Hilbert term. Hence, modified gravity can be seen as an exploration of how
gravity might behave at very low or high energy limits where the effects of
quantum physics become significant. See for instance [1, 2, 3, 4, 5], and (ii)
Dark matter and dark energy: The need for dark matter and dark energy to
explain the observed motion of galaxies and the accelerated expansion of the
universe, respectively, has led to questions about whether our understanding
of gravity is complete. Modified gravity theories also seek to address these
phenomena without introducing the need for dark matter or dark energy [6, 7,
8, 9, 10, 11, 12, 13, 14, 15, 16, 17].

Higher order curvature corrections to Einstein’s field equations have been
considered by many authors [18, 19, 20, 21, 22, 23, 24, 25, 26]. Recently we
showed that, in Friedman-Lemaitre-Robertson-Walker (FLRW) cosmology, the
contribution of higher curvature terms in any generic theory of gravity to the
energy momentum tensor is of the perfect fluid form [27]. This is the reason
that some authors have observed this fact [28, 29, 30] and verified it in some
particular modified gravity theories such as $f(R)$ gravity [31], Gauss-Bonnet
gravity [32], and other higher order gravities [33, 34]. In [27], the cases of
general Lovelock and $F(R,G)$ theories are given as examples.

The FLRW metric is the most known and most studied metric in general
relativity. This metric is mainly used to describe the Universe as a
homogeneous, isotropic fluid distribution [35]. It is known that FLRW
cosmology in Einstein’s theory is not sufficient to explain the accelerating
expansion of the universe. To explain this phenomena, staying in Einstein
theory it is claimed that, in addition to the known ordinary matter
distribution, the Universe should contain a dark component driving the
accelerating expansion, the so called Dark Energy. If we consider the low
energy limit of the string theory, then the Einstein equations are modified by
adding all possible curvature invariants. We call such a theory as generic
gravity where the action takes the form

$I=\int\,d^{D}\,x\,\sqrt{-g}\,\left[\frac{1}{\kappa}\left(R-2\Lambda\right)+\mathcal{F}(g,\,\mbox{Riem},\,\nabla\mbox{Riem},\,\nabla\nabla\mbox{Riem},\dots)+\mathcal{L}_{M}\right],$
(1)

here ${\cal F}$ is a function of all combinations of the metric tensor,
curvature tensor, and the covariant derivatives of the curvature tensor of any
order. Then the field equations take the form

$\frac{1}{\kappa}\left(G_{\mu\nu}+\Lambda g_{\mu\nu}\right)+{\cal
E}_{\mu\nu}=T^{M}_{\mu\nu},$ (2)

where $G_{\mu\nu}$ is the Einstein tensor, $\Lambda$ is the cosmological
constant, $T^{M}_{\mu\nu}$ is the matter energy momentum tensor of perfect
fluid distribution (with the energy density $\rho$ and pressure $p$), and
$\mathcal{E}_{\mu\nu}$ resulting from the higher order curvature terms is any
second rank tensor obtained from the metric tensor, Riemann tensor, Ricci
tensor, Ricci scalar and their covariant derivatives at any order. In [27], we
have shown that $\mathcal{E}_{\mu\nu}$ can be written as the combination of
the metric tensor $g_{\mu\nu}$ and $u_{\mu}u_{\nu}$; that is,

${\cal E}_{\mu\nu}=Ag_{\mu\nu}+Bu_{\mu}u_{\nu},$ (3)

where $A$ and $B$ are functions of scale factor $a(t)$ and its derivatives
with respect to time $t$. This implies that in the Einstein field equations,
in addition to the matter fluid energy momentum tensor $T^{M}_{\mu\nu}$, there
exists another fluid distribution which we call it as geometric fluid
distribution $T^{g}_{\mu\nu}$ with the energy density $\rho_{g}=A-B$ and
pressure $p_{g}=-A$. Here we adopt the idea that the source of the dark
matter/energy is the geometrical fluid distribution. Hence we conjecture that
higher curvature modifications of the Einstein theory are complete in the
sense that all cosmological observations can be explained by choosing
appropriate modified theories studied by several authors in [18, 19, 20, 21,
22, 23, 24, 25, 26]. To support our claim, we specifically study the quadratic
gravity field equations in $D$-dimensions. We show that the field equations of
this theory for the FLRW metric possess a geometric perfect fluid source
containing two critical parameters $\sigma_{1}$ and $\sigma_{2}$. To analyze
this theory concerning its parameter space $(\sigma_{1},\sigma_{2})$, we
obtain the general second-order nonlinear differential equation governing the
late-time dynamics of the deceleration parameter $q$. Hence using some
present-day cosmological data as our initial conditions, our findings for the
$\sigma_{2}=0$ case are as follows: $(i)$ In order to have a positive energy
density for the geometric fluid $\rho_{g}$, the parameter $\sigma_{1}$ must be
negative for all dimensions up to $D=11$, $(ii)$ For a suitable choice of
$\sigma_{1}$, the deceleration parameter experiences signature changes in the
past and future, and in the meantime it lies within a negative range which
means that the current observed accelerated expansion phase of the Universe
can be driven solely by the curvature of the spacetime, $(iii)$ $q$
experiences a signature change and as the dimension $D$ of spacetime
increases, this signature change happens at earlier and later times, in the
past and future respectively, $(iv)$ The geometric equation of state parameter
$w_{g}$ is negative for all integers $4\leq D\leq 10$, specifically, for
$D=4$, $w_{g}=-0.85$. For $\sigma_{1}=0$ (Critical Gravity), we find that
there are two cases both representing a possibility of having an accelerating
expansion. Furthermore, we present some particular cosmological solutions in
quadratic gravity depending upon the parameters $\sigma_{1}$ and $\sigma_{2}$.

In the next section we summarize our theorem for the $D-$dimensional FLRW
metrics. For illustration, we shall study the quadratic gravity theory in
detail.

##  2 Generic Gravity Field Equations in Perfect Fluid Form

Using the covariant decomposition, one can write the $D$-dimensional FLRW
metric as

$g_{\mu\nu}=-u_{\mu}\,u_{\nu}+a^{2}\,h_{\mu\nu},$ (4)

where $\mu,\,\nu=0...D-1$, $a=a(t),\leavevmode\nobreak\
u_{\mu}=\delta^{0}_{\mu}$, and $h_{\mu\nu}$ reads as

$h_{\mu\nu}=\begin{pmatrix}0&0&\ldots&0\\\ 0&&&\\\ \vdots&&h_{ij}&\\\ 0&&&\\\
\end{pmatrix},$ (5)

where $h_{ij}=h_{ij}(x^{a})$ is the metric of the spatial section of the
spacetime possessing the constant curvature $k$, and $i,j=1,...,D-1$. One
notes that

$\displaystyle u^{\mu}\,h_{\mu\nu}=u_{\mu}\,h^{\mu\nu}=0,$ $\displaystyle
h^{\mu}_{\alpha}=h^{\mu\alpha}\,h_{\alpha\nu}=\delta^{\mu}_{\nu}+u^{\mu}\,u_{\nu}.$
(6)

The corresponding Christoffel symbols can be obtained as

$\Gamma^{\mu}_{\alpha\beta}=\gamma^{\mu}_{\alpha\beta}-a\,\dot{a}\,u^{\mu}\,h_{\alpha\beta}+H\,\left(2u_{\alpha}\,u^{\mu}\,u_{\beta}+u_{\beta}\,\delta^{\mu}_{\alpha}+u_{\alpha}\,\delta^{\mu}_{\beta}\right),$
(7)

where the dot sign represents the derivative with respect to time $t$,
$H=\dot{a}/a$ and $\gamma^{\mu}_{\alpha\beta}$ is defined as

$\gamma^{\mu}_{\alpha\beta}=\frac{1}{2}a^{2}\,h^{\mu\gamma}\,\left(h_{\gamma\alpha,\beta}+h_{\gamma\beta,\alpha}-h_{\alpha\beta,\gamma}\right).$
(8)

The Riemann curvature tensor can be written in the following linear form in
terms of the metric $g_{\alpha\gamma}$ and the four vectors $u_{\alpha}$ and
$u^{\mu}$

$R^{\mu}_{\alpha\beta\gamma}=\left[\delta^{\mu}_{\beta}g_{\alpha\gamma}-\delta^{\mu}_{\gamma}g_{\alpha\beta}\right]\rho_{1}+\left[u^{\mu}\left(g_{\alpha\gamma}u_{\beta}-g_{\alpha\beta}u_{\gamma}\right)-u_{\alpha}\,\left(\delta^{\mu}_{\gamma}u_{\beta}-\delta^{\mu}_{\beta}u_{\gamma}\right)\right]\rho_{2},$
(9)

where $\rho_{1}$ and $\rho_{2}$ are defined as

$\displaystyle\rho_{1}=H^{2}+\frac{k}{a^{2}},$ (10)
$\displaystyle\rho_{2}=H^{2}+\frac{k}{a^{2}}-\frac{\ddot{a}}{a}=\dot{H}+\frac{k}{a^{2}}.$
(11)

The contractions of the Riemann tensor (9) gives the Ricci tensor and Ricci
scalar, respectively, as

$\displaystyle R_{\mu\nu}=Pg_{\mu\nu}+Qu_{\mu}u_{\nu},$ $\displaystyle
R=DP-Q,$ (12)

where

$\displaystyle P$ $\displaystyle=$ $\displaystyle(D-1)\rho_{1}-\rho_{2},$ (13)
$\displaystyle Q$ $\displaystyle=$ $\displaystyle(D-2)\rho_{2}.$ (14)

One can also verify that the Weyl tensor

$C^{\mu}_{\alpha\beta\gamma}=R^{\mu}_{\alpha\beta\gamma}+\frac{1}{D-2}\left(\delta^{\mu}_{\gamma}R_{\alpha\beta}-\delta^{\mu}_{\beta}R_{\alpha\gamma}+g_{\alpha\beta}R^{\mu}_{\gamma}-g_{\alpha\gamma}R^{\mu}_{\beta}\right)+\frac{1}{(D-1)(D-2)}\left(\delta^{\mu}_{\beta}g_{\alpha\gamma}-\delta^{\mu}_{\gamma}g_{\alpha\beta}\right)R$
(15)

vanishes for the metric (4). Since the conformal tensor is zero, the curvature
tensor is expressed in terms of the Ricci tensor. This means that, for the
FLRW spacetime, the basic geometrical tensors are the metric and Ricci
tensors. The covariant derivatives of the Ricci tensor are given by

$\displaystyle
R_{\mu\nu;\alpha}=\dot{P}\,u_{\alpha}\,g_{\mu\nu}-QH\left(u_{\nu}\,g_{\mu\alpha}+u_{\mu}\,g_{\nu\alpha}\right)+\left(\dot{Q}-2QH\right)u_{\mu}\,u_{\nu}\,u_{\alpha},$
(16) $\displaystyle{\square
R}_{\mu\nu}=-\left[\ddot{P}+(D-1)H\dot{P}-2QH^{2}\right]\,g_{\mu\nu}+\left[2DQH^{2}-\ddot{Q}-(D-1)H\dot{Q}\right]u_{\mu}u_{\nu},$
(17) $\displaystyle\square
R=-D\ddot{P}-D(D-1)H\dot{P}+\ddot{Q}+(D-1)H\dot{Q}.$ (18)

(NOTE: The above equations are from [27] (Eq.(19)) where some typos are
corrected.)

The field equations of any generic gravity theory in D dimensions with the
action, together with the matter fields,

$I=\int\,d^{D}\,x\,\sqrt{-g}\,\left[\frac{1}{\kappa}\left(R-2\Lambda\right)+\mathcal{F}(g,\,\mbox{Riem},\,\nabla\mbox{Riem},\,\nabla\nabla\mbox{Riem},\cdots)+\mathcal{L}_{M}\right],$
(19)

take the form

$\frac{1}{\kappa}(G_{\mu\nu}+\Lambda g_{\mu\nu})+{\cal
E}_{\mu\nu}=T^{M}_{\mu\nu},$ (20)

where $G_{\mu\nu}$ is the Einstein tensor, $\Lambda$ is the cosmological
constant, $T^{M}_{\mu\nu}$ is the energy momentum tensor coming from the
matter fields denoted by $\mathcal{L}_{M}$, and $\cal{E}_{\mu\nu}$ resulting
from the higher order curvature terms contained in the function $\mathcal{F}$
is any second rank tensor obtained from the metric tensor, Riemann tensor,
Ricci tensor, and their covariant derivatives at any order. Hence we arrive at
the following theorem [27]:

Theorem 1 : For the $D$-dimensional FLRW spacetimes given in (4), any second
rank symmetric tensor obtained from the metric tensor, Riemann tensor, Ricci
tensor, and their covariant derivatives at any order becomes a combination of
the metric tensor $g_{\mu\nu}$ and $u_{\mu}u_{\nu}$; that is,

${\cal{E}}_{\mu\nu}=Ag_{\mu\nu}+Bu_{\mu}u_{\nu},$ (21)

where $A$ and $B$ are functions of the scale factor $a(t)$ and its time
derivatives at any order and they depend on the underlying gravity theory.

And we have the following corollary of this theorem:

Corollary 2: The field equations of any generic gravity theory given in (20)
take the form

$\frac{1}{\kappa}(G_{\mu\nu}+\Lambda
g_{\mu\nu})=T^{M}_{\mu\nu}+T^{g}_{\mu\nu},$ (22)

where $T^{M}_{\mu\nu}$ is the energy momentum tensor of perfect fluid
distribution representing the baryonic matter fields,

$T^{M}_{\mu\nu}=(\rho+p)u_{\mu}u_{\nu}+pg_{\mu\nu},$ (23)

with $\rho$ and $p$ being respectively the energy density and pressure of the
fluid, and $T^{g}_{\mu\nu}$ is the tensor coming from the higher order
curvature terms in (19) which is also in the perfect fluid form,

$T^{g}_{\mu\nu}=(\rho_{g}+p_{g})u_{\mu}u_{\nu}+p_{g}g_{\mu\nu},$ (24)

with $\rho_{g}=A-B$ and $p_{g}=-A$, due to (21). Hence defining an effective
energy-momentum tensor as

$T^{eff}_{\mu\nu}\equiv
T^{M}_{\mu\nu}+T^{g}_{\mu\nu}=(\rho_{eff}+p_{eff})u_{\mu}u_{\nu}+p_{eff}g_{\mu\nu},$
(25)

the generic gravity field equations (22) for the FLRW spacetime with a perfect
fluid source become

$\displaystyle\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\rho_{1}-\Lambda\right]=\rho+A-B\equiv\rho_{eff},$
(26)
$\displaystyle-\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\rho_{1}-(D-2)\rho_{2}-\Lambda\right]=p-A\equiv
p_{eff}.$ (27)

Thus the interpretations of the functions $A$ and $B$ appearing in the above
formulation can be given as follows: the combination “$A-B$” is the energy
density and “$-A$” is the pressure of a perfect fluid of purely geometric
origin. The functions $A$ and $B$ differ in different modified theories. In
each modification it is possible to arrange parameters of the theories to meet
the observations. In particular, we shall analyze the quadratic gravity theory
in $D$ dimensions in Section 4 and show that these purely geometric terms
solely drive the late time accelerating expansion of the Universe consistently
with the current observations. In the next section, we shall give the
cosmological parameters for generic gravity theories described by the action
in (19) in $D$ dimensions.

##  3 Cosmological Parameters in Generic Gravity Theories

Using (10) and (26), one can obtain

$\displaystyle H^{2}$ $\displaystyle=$
$\displaystyle\frac{2\kappa}{(D-1)(D-2)}\left[\rho+A-B+\frac{\Lambda}{\kappa}\right]-\frac{k}{a^{2}}$
(28) $\displaystyle=$
$\displaystyle\frac{2\kappa}{(D-1)(D-2)}\sum_{i}\rho_{i}-\frac{k}{a^{2}},$

where the label ‘$i$’ denotes ‘$m$’, ‘$r$’, ‘$\Lambda$’, or ‘$g$’ representing
matter, radiation, cosmological constant (or dark energy), and “dark geometric
fluid”, respectively. Defining the corresponding dimensionless density
parameter for each of the mentioned components of the Universe as

$\Omega_{i}=\frac{\rho_{i}}{\rho_{crit}}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\text{with}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \rho_{crit}=\frac{(D-1)(D-2)H^{2}}{2\kappa},$ (29)

we can write (28) as

$1=\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}+\Omega_{g}+\Omega_{k},$ (30)

where

$\Omega_{m}=\frac{\rho_{m}}{\rho_{crit}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Omega_{r}=\frac{\rho_{r}}{\rho_{crit}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Omega_{\Lambda}=\frac{\rho_{\Lambda}}{\rho_{crit}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Omega_{g}=\frac{A-B}{\rho_{crit}},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \Omega_{k}=-\frac{k}{a^{2}H^{2}}.$
(31)

and note that $\rho_{\Lambda}=\Lambda/\kappa$. The first observation here is
the contribution of the dark geometric fluid in determining the spatial
curvature of the Universe. One observes that

$\displaystyle\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}+\Omega_{g}<1\leavevmode\nobreak\
\Leftrightarrow\ \textnormal{open universe ($k=-1$)}$ (32)
$\displaystyle\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}+\Omega_{g}=1\leavevmode\nobreak\
\Leftrightarrow\ \textnormal{flat universe ($k=0$)}$ (33)
$\displaystyle\Omega_{m}+\Omega_{r}+\Omega_{\Lambda}+\Omega_{g}>1\leavevmode\nobreak\
\Leftrightarrow\ \textnormal{closed universe ($k=1$)}.$ (34)

On the other hand, using (10), (11), and (28), the equation (27) can be
written equivalently as

$-\frac{\ddot{a}}{aH^{2}}=\frac{\kappa}{(D-1)(D-2)H^{2}}\left[(D-3)\rho+(D-1)p-2A-B-2\frac{\Lambda}{\kappa}\right].$
(35)

Now defining the left hand side as the deceleration parameter $q$ and
considering the barotropic equation of state $p_{j}=\omega_{j}\rho_{j}$, where
$j=(m,r)$, for the matter ($\omega_{m}=0$) and the radiation
($\omega_{r}=1/3$), the equation (35) can be rewritten, with the help of (29),
as

$\displaystyle q\equiv-\frac{\ddot{a}}{aH^{2}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\left[\sum_{j}[(D-3)+(D-1)w_{j}]\Omega_{j}-2\Omega_{\Lambda}-\Omega^{*}_{g}\right]$
(36) $\displaystyle=$
$\displaystyle\frac{1}{2}\left[(D-3)\Omega_{m}+\frac{2}{3}(2D-5)\Omega_{r}-2\Omega_{\Lambda}-\Omega^{*}_{g}\right],$

where

$\Omega^{*}_{g}=\frac{2A+B}{\rho_{crit}}.$ (37)

Here we see that the geometric fluid has negative contribution to the
deceleration (or, equivalently, positive contribution to the acceleration) of
the Universe, if $2A+B>0$. It can be observed that, setting $D=4$ and
neglecting higher curvature modifications, one recovers the deceleration
parameter in GR as

$q=\frac{1}{2}\left(\Omega_{m}+2\Omega_{r}-2\Omega_{\Lambda}\right).$ (38)

At this point, before proceeding further, it is appropriate to introduce the
cosmological parameters

$j(t)\equiv\frac{\dddot{a}}{aH^{3}},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
s(t)\equiv\frac{\ddddot{a}}{aH^{4}},$ (39)

which are called “jerk” and “snap”, respectively. These parameters, together
with $H$ and $q$, are defined by expanding the scale factor in a Taylor series
in the vicinity of the current time $t_{0}$ [36]:

$\frac{a(t)}{a_{0}}=1+H_{0}(t-t_{0})-\frac{1}{2}q_{0}H_{0}^{2}(t-t_{0})^{2}+\frac{1}{3!}j_{0}H_{0}^{3}(t-t_{0})^{3}+\frac{1}{4!}s_{0}H_{0}^{4}(t-t_{0})^{4}+O((t-t_{0})^{5}),$
(40)

where $H_{0}$, $q_{0}$, $j_{0}$, and $s_{0}$ are the present-day values of the
Hubble, deceleration, jerk, and snap parameters and they can be used to
determine the evolutionary behavior of the Universe.

##  4 Quadratic Gravity and Criticality in $D$-Dimensions

The action of the quadratic gravity theory [37, 38] is

$\displaystyle I$ $\displaystyle=$ $\displaystyle\int
d^{D}x\,\sqrt{-g}\left[\frac{1}{\kappa}\left(R-2\Lambda\right)+\alpha
R^{2}+\beta
R_{\mu\nu}^{{}^{2}}+\gamma\left(R_{\mu\nu\sigma\rho}^{2}-4R_{\mu\nu}^{2}+R^{2}\right)+\mathcal{L}_{M}\right],$
(41)

giving the field equations

$\frac{1}{\kappa}\left(R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R+\Lambda
g_{\mu\nu}\right)+{\cal{E}}_{\mu\nu}=(\rho+p)\,u_{\mu}\,u_{\nu}+p\,g_{\mu\nu},$
(42)

where $\rho$ and $p$ are the energy density and pressure of the matter perfect
fluid. Considering the FLRW metric (4), we find

$\displaystyle{\cal{E}}_{\mu\nu}$ $\displaystyle\equiv$ $\displaystyle 2\alpha
R\left(R_{\mu\nu}-\frac{1}{4}g_{\mu\nu}R\right)+\left(2\alpha+\beta\right)\left(g_{\mu\nu}\square-\nabla_{\mu}\nabla_{\nu}\right)R$
(43)
$\displaystyle+2\gamma\left[RR_{\mu\nu}-2R_{\mu\sigma\nu\rho}R^{\sigma\rho}+R_{\mu\sigma\rho\tau}R_{\nu}^{\phantom{\nu}\sigma\rho\tau}-2R_{\mu\sigma}R_{\nu}^{\phantom{\nu}\sigma}-\frac{1}{4}g_{\mu\nu}\left(R_{\tau\lambda\sigma\rho}^{2}-4R_{\sigma\rho}^{2}+R^{2}\right)\right]$
$\displaystyle+\beta\,\square\left(R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R\right)+2\beta\left(R_{\mu\sigma\nu\rho}-\frac{1}{4}g_{\mu\nu}R_{\sigma\rho}\right)R^{\sigma\rho}$
$\displaystyle=$ $\displaystyle Ag_{\mu\nu}+Bu_{\mu}u_{\nu},$

where $A$ and $B$ are given by

$\displaystyle A$ $\displaystyle=$
$\displaystyle\frac{\alpha}{2}\left[\left(D-1\right)\left(D\rho_{1}-2\rho_{2}\right)\right]\left[-(D-1)(D-4)\rho_{1}+2(D-3)\rho_{2}\right]$
(44)
$\displaystyle+(2\alpha+\beta)\left[(D-1)(D-2)H\left(-D\dot{\rho}_{1}+2\dot{\rho}_{2}\right)-D(D-1)\ddot{\rho}_{1}+2(D-1)\ddot{\rho}_{2}\right]$
$\displaystyle+\frac{\gamma}{2}\left[(D-2)(D-3)(D-4)\left(-(D-1)\rho_{1}^{2}+4\rho_{1}\rho_{2}\right)\right]$
$\displaystyle+\frac{\beta}{2}\left[(D-2)\left[(D-1)(\ddot{\rho}_{1}+(D-1)H\dot{\rho}_{1})\right]-2(\ddot{\rho}_{2}+(D-1)H\dot{\rho}_{2}-2H^{2}\rho_{2})\right]$
$\displaystyle+\frac{\beta}{2}\left[-(D-1)(D-4)\left((D-1)\rho_{1}^{2}+\rho_{2}^{2}\right)+4(D^{2}-5D+5)\rho_{1}\rho_{2}\right],$
$\displaystyle B$ $\displaystyle=$ $\displaystyle
2\alpha\left[\left(D-1\right)(D-2)\left(D\rho_{1}-2\rho_{2}\right)\rho_{2}\right]$
(45)
$\displaystyle+(2\alpha+\beta)\left[(D-1)\left[H\left(-D\dot{\rho}_{1}-2\dot{\rho}_{2}\right)-(D\ddot{\rho}_{1}-2\ddot{\rho}_{2})\right]\right]$
$\displaystyle+2\gamma\left[(D-2)(D-3)(D-4)\rho_{1}\rho_{2}\right]$
$\displaystyle-\beta\left[(D-2)\left[\ddot{\rho}_{2}+(D-1)H\dot{\rho}_{2}-2DH^{2}\rho_{2})\right]\right]$
$\displaystyle+2\beta\left[(D-2)^{2}\rho_{1}\rho_{2}\right].$

Hence, we can write (42) as

$\displaystyle\frac{1}{\kappa}\left[\frac{1}{2}(D-2)(D-1)\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)-\Lambda\right]=\rho+\rho_{g_{1}}+\rho_{g_{2}}\equiv\rho_{eff},$
(46)
$\displaystyle-\frac{1}{\kappa}\left[\frac{(D-2)(D-3)}{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)+(D-2)\frac{\ddot{a}}{a}-\Lambda\right]=p+p_{g_{1}}+p_{g_{2}}\equiv
p_{eff},$ (47)

where $\rho_{g_{1}},\rho_{g_{2}},p_{g_{1}}$, and $p_{g_{2}}$ have the
geometric origin given by

$\displaystyle\rho_{g_{1}}=-\frac{(D-1)\sigma_{1}}{2Da^{4}}\Biggl{\\{}k^{2}(D-2)^{2}+\dot{a}^{2}\left[2D(D-3)a\ddot{a}-4k(D-2)\right]$
$\displaystyle\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\
\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\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
-\left(D^{2}-4\right)\dot{a}^{4}-Da^{2}\ddot{a}^{2}+2Da^{2}\dot{a}\dddot{a}\Biggr{\\}},$
(48)
$\displaystyle\rho_{g_{2}}=-(D-1)\sigma_{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)^{2},$
(49) $\displaystyle
p_{g_{1}}=-\frac{\sigma_{1}}{2Da^{4}}\Biggl{\\{}(D-2)(D-5)(D+2)\dot{a}^{4}-k^{2}(D-2)^{2}(D-5)+a\ddot{a}\left(8k(D-2)-3D(D-3)a\ddot{a}\right)$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+2\dot{a}^{2}\left[2k(D-2)(D-5)-\left((D-9)D^{2}+12D+8\right)a\ddot{a}\right]-2D\ddddot{a}a^{3}-4(D-3)D\dddot{a}a^{2}\dot{a}\Biggr{\\}},$
(50) $\displaystyle
p_{g_{2}}=\sigma_{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)\left[(D-5)\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)+4\frac{\ddot{a}}{a}\right],$
(51)

with

$\displaystyle\sigma_{1}=4(D-1)\alpha+D\beta,$ (52)
$\displaystyle\sigma_{2}=\frac{(D-2)(D-4)}{2D}\left[(D-1)(D-2)\alpha+D(D-3)\gamma\right].$
(53)

Remark 3: When the parameters $\sigma_{1}$ and $\sigma_{2}$ vanish together,
the geometric contributions given in (48)-(51), resulting from the higher
curvature terms in the action (42), vanish identically and the equations (46)
and (47) reduce to the equations in pure Einstein’s gravity. These so-called
critical points; i.e. $\sigma_{1}=0$ and $\sigma_{2}=0$, were identified and
studied in higher curvature gravity theories first in four dimensions (where
$\sigma_{2}$ identically vanishes) in [39] and later in $D$ dimensions in
[40]. In these works, it is shown that the linearized excitations around these
critical points have vanishing energies and the mass and corresponding entropy
of the usual Schwarzschild-AdS black hole solution turn out to be zero at
criticality.

Remark 4: It can also be observed that the above expressions are invariant
under the scale transformations $a\to\eta a$ and $k\to\eta^{2}k$, where $\eta$
is a constant.

Remark 5: From the positiveness of $\rho_{g_{2}}$ in (49) it follows that
$\sigma_{2}<0$.

From the equations (46) and (47), we can also write

$\frac{\ddot{a}}{a}=-\frac{\kappa}{(D-1)(D-2)}\left[(D-3)\rho_{eff}+(D-1)p_{eff}-2\frac{\Lambda}{\kappa}\right].$
(54)

Thus, when $\Lambda=0$, to have a universe expanding in an accelerating
fashion ($\ddot{a}>0$), it must be that

$w_{eff}=\frac{p_{eff}}{\rho_{eff}}<-\frac{D-3}{D-1}.$ (55)

Specifically, for $D=4$, $w_{eff}<-1/3$. Using (46), along with the definition
of the Hubble parameter, $H=\dot{a}/a$, one can also write (54) as

$\dot{H}=-\frac{\kappa}{D-2}(\rho_{eff}+p_{eff})+\frac{k}{a^{2}},$ (56)

which is more convenient than (54) because it does not involve the
cosmological constant $\Lambda$ explicitly. This equation relates the
acceleration of the Universe to its energy and momentum content. Together with
the expressions (48)-(51), the equation (56) becomes highly nonlinear and
therefore it is not possible to give a compact analytical solution for $a(t)$.
However, we can use (56) to investigate the late-time accelerated expansion of
the Universe during which the higher curvature terms in the action (41) are
assumed to be dominant. Indeed, recently it was observed that the Universe had
entered into an accelerated expansion phase [41, 42] and a possible cause of
this late time acceleration is the curvature of the spacetime itself related
to the combination of the higher curvature terms that may appear in the action
of a generic gravity theory. This observation implies that the acceleration,
defined by $\ddot{a}(t)$, of the Universe has changed its sign from negative
to positive, or in other words, the deceleration parameter $q(t)$ has recently
experienced a sign change from positive to negative. Therefore, for the
purpose of investigating the late time acceleration or deceleration behavior
of the Universe, we shall work with the deceleration parameter $q(t)$, instead
of $a(t)$, and convert (56) as a differential equation for $q(t)$ which can be
solved numerically for given initial values of $q$ and its derivatives. Let us
first write $\rho_{eff}$ and $p_{eff}$ in terms of the Hubble ($H$),
deceleration ($q$), jerk ($j$), and snap ($s$) parameters. From (48)-(51),

$\displaystyle\rho_{eff}=\rho-\frac{\sigma_{1}(D-1)}{2D}\Biggl{[}\frac{k^{2}(D-2)^{2}}{a^{4}}-4\frac{k}{a^{2}}(D-2)H^{2}$
$\displaystyle\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\
\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\
-2D(D-3)qH^{4}-(D^{2}-4)H^{4}-Dq^{2}H^{4}+2DjH^{4}\Biggl{]}$
$\displaystyle\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\ -\sigma_{2}(D-1)\left(\frac{k}{a^{2}}+H^{2}\right)^{2},$
(57) $\displaystyle
p_{eff}=p-\frac{\sigma_{1}}{2D}\Biggl{[}-\frac{k^{2}}{a^{4}}(D-2)^{2}(D-5)+4\frac{k}{a^{2}}(D-2)H^{2}[D-5-2q]+(D-2)(D+2)(D-5)H^{4}$
$\displaystyle\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\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-3D(D-3)q^{2}H^{4}+2(D^{3}-9D^{2}+12D+8)qH^{4}-2DsH^{4}-4(D-3)DjH^{4}\Biggr{]},$
$\displaystyle\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\
+\sigma_{2}\left(\frac{k}{a^{2}}+H^{2}\right)\left[(D-5)\left(\frac{k}{a^{2}}+H^{2}\right)-4qH^{2}\right],$
(58)

where

$j=-\frac{\dot{q}}{H}+q(1+2q),\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
s=-\frac{\ddot{q}}{H^{2}}+2\frac{\dot{q}}{H}(1+3q)-q(1+2q)(1+3q).$ (59)

In getting these expressions, use has been made of

$q=-\frac{\ddot{a}}{aH^{2}}=-\frac{\dot{H}}{H^{2}}-1.$ (60)

Now we shall set the ordinary matter and curvature to zero in the formulation,
i.e. $\rho=0=p$ and $k=0$, to both comply with the observations and
investigate the late-time acceleration of the Universe resulting purely from
the geometric terms related to curvature of the spacetime. Then, (56) becomes
an equation involving $H$, $q$, and the first and second derivatives of $q$
with respect to time; that is,

$\displaystyle\left[1+\frac{4H^{2}\kappa\sigma_{2}}{D-2}\right](q+1)+\frac{\kappa\sigma_{1}}{(D-2)D}\Biggl{\\{}-2H^{2}\left(D^{2}-4\right)+H(D-7)D\dot{q}+D\ddot{q}$
$\displaystyle\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\
+Hq\left[H\left(-6D^{2}+15D+8\right)-6D\dot{q}\right]-4H^{2}(D-5)Dq^{2}+6H^{2}Dq^{3}\Biggr{\\}}=0.$
(61)

However, this form is not appropriate for solving the equation numerically. To
obtain an appropriate form, we will change the time derivatives to derivatives
with respect to $H$, since both $H$ and $q$ are functions of time only and
related to each other by (60). That is, using

$\displaystyle\dot{q}=-H^{2}(1+q)q^{\prime},$
$\displaystyle\ddot{q}=H^{4}(1+q)^{2}q^{\prime\prime}+2H^{3}(1+q)^{2}q^{\prime}+H^{4}(1+q)q^{\prime
2},$ (62)

where the prime denotes derivatives with respect to $H$, we can bring (4) into
the following second-order nonlinear differential equation for $q$ in $H$:

$\displaystyle\left[1+\frac{4H^{2}\kappa\sigma_{2}}{D-2}\right](q+1)+\frac{H^{2}\kappa\sigma_{1}}{(D-2)D}\Biggl{\\{}8-2D^{2}+6Dq^{3}-H(D-9)Dq^{\prime}+H^{2}Dq^{\prime
2}+H^{2}Dq^{\prime\prime}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+Dq^{2}(20-4D+8Hq^{\prime}+H^{2}q^{\prime\prime})+q[8+15D-6D^{2}-H(D-17)Dq^{\prime}+H^{2}Dq^{\prime
2}+2H^{2}Dq^{\prime\prime}]\Biggr{\\}}=0.$

Once the initial values $H_{0}$, $q_{0}$, and $q^{\prime}_{0}$ are given, this
equation can be numerically solved for the evolution of the deceleration
parameter $q$ in $H$. This kind of analysis was previously exploited in the
context of $f(R)$ gravity in [15, 23, 24]. In the following analysis we will
take the present-day values of the cosmological parameters $H$, $q$, $j$, and
$s$ as [43, 44]

$H_{0}=67.4\leavevmode\nobreak\
\text{km}\cdot\text{s}^{-1}\cdot\text{Mpc}^{-1},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
q_{0}=-0.71,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ j_{0}=1.26,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ s_{0}=0.04.$ (64)

Here, we shall consider two particular cases regarding our parameter space
$(\sigma_{1},\sigma_{2})$ corresponding to $\sigma_{2}=0$ (Case I) and
$\sigma_{1}=0$ (Case II).

Case I: When $\sigma_{2}$=0, the deceleration equation (4) becomes

$\displaystyle(q+1)+\frac{h^{2}\kappa\sigma_{1}}{(D-2)D}\Biggl{\\{}8-2D^{2}+6Dq^{3}-h(D-9)Dq^{\prime}+h^{2}Dq^{\prime
2}+h^{2}Dq^{\prime\prime}$ $\displaystyle\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+Dq^{2}(20-4D+8hq^{\prime}+h^{2}q^{\prime\prime})+q[8+15D-6D^{2}-h(D-17)Dq^{\prime}+h^{2}Dq^{\prime
2}+2h^{2}Dq^{\prime\prime}]\Biggr{\\}}=0.$

Note that in writing this equation we have rescaled the Hubble parameter $H$
in (4) as

$H\equiv(100\leavevmode\nobreak\
\text{km}\cdot\text{s}^{-1}\cdot\text{Mpc}^{-1})\leavevmode\nobreak\ h$ (66)

to construct a dimensionless parameter $h$ and assumed all the numerical
constants are absorbed into the constant $\sigma_{1}$.

Remark 6: Before presenting the solution of (4), it would be useful to look at
the energy density $\rho_{g}$ related to the geometry: Since it must be
positive, we can determine the sign of the constant $\sigma_{1}$ by evaluating
the energy density at the present time. From (48), we can infer that the
energy density at the present epoch is

$\rho_{g}=\frac{(D-1)h_{0}^{2}\sigma_{1}}{2D}\left[2D(D-3)q_{0}+(D^{2}-4)+Dq_{0}^{2}-2Dj_{0}\right].$
(67)

Remark 7: Inserting $h_{0}$, $q_{0}$, and $j_{0}$ from (64) into $\rho_{g}$
and graphing with respect to $\sigma_{1}$ and $D$, we obtain the plot shown in
Fig. 1. As is obvious from the figure, to have a positive $\rho_{g}$, the
parameter $\sigma_{1}$ must be negative for all dimensions up to $D=11$.

Figure 1: The plot of $\rho_{g}$ vs $\sigma_{1}$ and $D$ for $h_{0}=0.674$,
$q_{0}=-0.71$, and $j_{0}=1.26$.

Now we can solve the equation (4) numerically: First, to observe the effect of
the value of $\sigma_{1}$ on the solution, we plot the solution of $q(h)$ for
$D=4$ and for different values of $\sigma_{1}$, this is shown in Fig 2.

Figure 2: The plot of $q$ as a function of $h$ for different values of
$\sigma_{1}$ and for $D=4$, $\kappa=1$, $\sigma_{2}=0$, $h_{0}=0.674$,
$q_{0}=-0.71$, and $q_{0}^{\prime}=4.92$.

Remark 8: Since $h$ is related to the inverse of the cosmic time $t$, in
Fig.2, $h>h_{0}$ defines the past and $h<h_{0}$ defines the future in the
cosmic evolution of the Universe. As is obvious from the figure, the
deceleration parameter $q$ experiences two signature changes: one is at the
past and the other is at the future.

Remark 9: There is some time interval in which the deceleration parameter is
negative (i.e. $q<0$). In particular, close to the present value of the Hubble
parameter $h_{0}=0.674$, the deceleration parameter $q$ is negative. This
means that the observed accelerated expansion phase of the Universe can be
driven solely by the curvature of the spacetime.

Remark 10: It can also be observed from Fig.2 that, as the value of
$\sigma_{1}$ decreases in negative (or increases in magnitude), the curves are
opening out and approaching each other, and for very small values (or large
magnitudes) they are becoming almost identical. This stems from the fact that,
as the magnitude of the parameter $\sigma_{1}$ increases in (4), the first
term in the equation can be neglected and the equation is effectively reduced
to the one in which the curly bracketed expression equals zero.

Figure 3: Behavior of $q$ for $4\leq D\leq 10$, $\kappa=1$, $\sigma_{1}=-10$,
$\sigma_{2}=0$, $h_{0}=0.674$, $q_{0}=-0.71$, $q^{\prime}_{0}=4.92$.

Remark 11: Additionally, we can also investigate the behavior of $q(h)$ with
$D$. This is shown in Fig. 3. Here it is explicitly seen that, as $D$
increases, the signature change of $q$ occurs at later times.

On the other hand, one can also study the behavior of the equation of state
parameter of the dark fluid $w_{g}$ stemmed from the terms proportional to
$\sigma_{1}$ in (57) and (58); that is, since $w_{g}=p_{g}/\rho_{g}$, at the
present time $t_{0}$ we obtain

$w_{g}=\frac{(D-2)(D+2)(D-5)-3D(D-3)q_{0}^{2}+2(D^{3}-9D^{2}+12D+8)q_{0}-2Ds_{0}-4(D-3)Dj_{0}}{(D-1)[-2D(D-3)q_{0}-(D^{2}-4)-Dq_{0}^{2}+2Dj_{0}]},$
(68)

the plot of which is given in Fig. 4 with respect to $D$, where the dotted
lines represent the upper bounds in (55) decreasing in value as $D$ increases.
Remark 12: As is obvious, $w_{g}$ is negative for all integers $4\leq D\leq
10$, consistently with Fig. 3, and satisfying the condition (55) for all
dimensions. Specifically, for $D=4$, $w_{g}=-0.85$.

Figure 4: Behavior of $w_{g}$ for $\sigma_{2}=0$, $q_{0}=-0.71$, $j_{0}=1.26$,
$s_{0}=0.04$.

Case II: When $\sigma_{1}$=0 and $D\neq 4$, the deceleration equation (4)
becomes

$\left[1+\frac{4H^{2}\kappa\sigma_{2}}{D-2}\right](q+1)=0.$ (69)

Recalling that $\sigma_{2}<0$ follows from the positiveness of $\rho_{g2}$,
there are two possibilities:

$(i)$ $\left[1+\frac{4H^{2}\kappa\sigma_{2}}{D-2}\right]\neq 0$ and $q=-1$,
independently of the number of dimensions. Thus the deceleration parameter is
always negative, representing an accelerating Universe driven by the curvature
of the spacetime. Also, one can show that for the equation of state of the
geometric fluid $w_{g}$ stemmed from the terms proportional to $\sigma_{2}$ in
(57) and (58) is

$w_{g}=-\frac{D-5-4q}{D-1}=-1,$ (70)

for $q=-1$ independently of the number of dimensions. Hence, the geometric
fluid can derive the accelerating expansion of the Universe playing the role
of an effective cosmological constant.

$(ii)$ When $H^{2}=-\frac{(D-2)}{4\kappa\sigma_{2}}$, it follows that $q=-1$
automatically. This case represents an exact exponential solution for the
scale factor, i.e $a(t)=a_{0}e^{\pm\lambda t}$ where
$\lambda=\sqrt{-(D-2)/4\sigma_{2}}$.

##  5 Some Particular Cosmological Solutions in Quadratic Gravity

In this section, we shall investigate some particular solutions to the general
equations presented in (46) and (47).

###  5.1 Solutions with $\sigma_{2}=0$

When $\sigma_{2}$ vanishes, the field equations (46) and (47) become

$\displaystyle\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)-\Lambda\right]$
$\displaystyle\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\
=\rho-\frac{(D-1)\sigma_{1}}{2Da^{4}}\Biggl{[}k^{2}(D-2)^{2}+\dot{a}^{2}\left(2D(D-3)a\ddot{a}-4k(D-2)\right)$
$\displaystyle\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\
\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\
\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\
-\left(D^{2}-4\right)\dot{a}^{4}-Da^{2}\ddot{a}^{2}+2Da^{2}\dot{a}\dddot{a}\Biggr{]},$
(71)
$\displaystyle-\frac{1}{\kappa}\left[\frac{(D-2)(D-3)}{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)+(D-2)\frac{\ddot{a}}{a}-\Lambda\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
=p+\frac{\sigma_{1}}{2Da^{4}}\Biggl{\\{}(D-2)(D-5)(D+2)\dot{a}^{4}-k^{2}(D-2)^{2}(D-5)+a\ddot{a}\left[8k(D-2)-3D(D-3)a\ddot{a}\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
+2\dot{a}^{2}\left[2k(D-2)(D-5)-\left(D((D-9)D+12)+8\right)a\ddot{a}\right]-2D\ddddot{a}a^{3}-4(D-3)D\dddot{a}a^{2}\dot{a}\Biggr{\\}}.$
(72)

It must be noted that these equations are valid either when $D=4$ and
$(D-1)(D-2)\alpha+D(D-3)\gamma\neq 0$, or when $D\neq 4$ and
$(D-1)(D-2)\alpha+D(D-3)\gamma=0$, or when $D=4$ and
$(D-1)(D-2)\alpha+D(D-3)\gamma=0$. In four dimensions, these reduce to

$\displaystyle\frac{1}{\kappa a^{4}}\left[3a^{2}(\dot{a})^{2}+3ka^{2}-\Lambda
a^{4}\right]=\rho-\frac{3\sigma_{1}}{2a^{4}}\left[2a^{2}\dot{a}\dddot{a}-a^{2}(\ddot{a})^{2}+2a(\dot{a})^{2}\,\ddot{a}-3(\dot{a})^{4}-2k(\dot{a})^{2}+k^{2}\right],$
(73) $\displaystyle-\frac{1}{\kappa
a^{4}}\left[2a^{3}\ddot{a}+a^{2}(\dot{a})^{2}+ka^{2}-\Lambda a^{4}\right]$
$\displaystyle\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\
=p+\frac{\sigma_{1}}{2a^{4}}\left[2a^{3}\ddddot{a}+4a^{2}\dot{a}\dddot{a}+3a^{2}(\ddot{a})^{2}-12a(\dot{a})^{2}\ddot{a}+3(\dot{a})^{4}-4ka\ddot{a}+2k(\dot{a})^{2}-k^{2}\right].$
(74)

where $\sigma_{1}=4(3\alpha+\beta$). This case is known as the critical
gravity where (73) and (74) reduces to the ones in the usual Einstein’s
gravity when $\sigma_{1}=0$ [37, 38]. Here we have the following immediate
Remarks:

Remark 13: When the critical parameter $\sigma_{1}=0$ [39, 40] it is quite
interesting that for this special case the above energy density and pressure
expressions for the FLRW metric reduce to the corresponding expressions in
pure Einstein theory. This means that the highly nontrivial tensor field
${\mathcal{E}}_{\mu\nu}$ given in (21) reduces to

${\mathcal{E}}_{\mu\nu}=\alpha\left(2RR_{\mu\nu}-\frac{1}{2}R^{2}g_{\mu\nu}+\frac{1}{2}g_{\mu\nu}\square
R+\nabla_{\mu}\nabla_{\nu}R-3\,\square
R_{\mu\nu}-6R_{\mu\rho\nu\sigma}R^{\rho\sigma}+\frac{3}{2}g_{\mu\nu}R^{\rho\sigma}R_{\rho\sigma}\right),$
(75)

which vanishes identically for $\alpha\neq 0$.

Remark 14: The case when $k=0$ and vanishing of the coefficient of
$\sigma_{1}$ in (73) and (74) corresponds to the work of Barrow et al. [21].
They found a power law solution for $a$ and studied the stability of the
solution.

####  5.1.1 Exponential Solutions

Now we shall study a special case which may correspond to the late-time
accelerating era of the Universe. Let $a(t)=a_{0}\,e^{H_{0}t}$, where $a_{0}$
and $H_{0}$ are the scale factor and the Hubble constant at the time when the
accelerating era begins, respectively. It is interesting that, for an
exponentially expanding flat universe ($k=0$), the contributions of the
quadratic gravity terms related to $\sigma_{1}$ in (71) and (72) vanish
identically. This means that the presence of the bare cosmological constant
$\Lambda$ is crucial for having exponential solutions in arbitrary $D$
dimensions. When $k\neq 0$, the field equations (71) and (72) become

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{(D-1)(D-2)}{2D\kappa
a_{0}^{4}}\left[D\beta_{0}\,a_{0}^{4}+k\,a_{0}^{2}\beta_{2}\,e^{-2H_{0}t}+(D-2)k^{2}\kappa\sigma_{1}e^{-4H_{0}t}\right],$
(76) $\displaystyle p$ $\displaystyle=$ $\displaystyle-\frac{D-2}{2D\kappa
a_{0}^{4}}\left[D\beta_{0}\,a_{0}^{4}+(D-3)k\,a_{0}^{2}\,\beta_{2}\,e^{-2H_{0}t}+(D-2)(D-5)k^{2}\kappa\sigma_{1}e^{-4H_{0}t}\right],$
(77)

where $\beta_{0}=H_{0}^{2}-\frac{2\Lambda}{(D-1)(D-2)}$,
$\beta_{2}=D-4\kappa\sigma_{1}H_{0}^{2}$, and
$\sigma_{1}=4(D-1)\alpha+D\beta$. We have the following consequences.

(a) If $\beta_{0}>0$ and $\beta_{2}>0$, the energy density remains positive
for all $t$. At late times when $t\to\infty$, we have
$\rho\to\frac{(D-1)(D-2)\beta_{0}}{2\kappa}$ and
$p\to-\frac{(D-2)\beta_{0}}{2\kappa}$. Hence, the equation of state is of dark
energy type, i.e $p\to-\frac{\rho}{D-1}$.

(b) When $\beta_{0}=0$ or $H_{0}=\sqrt{\frac{2\Lambda}{(D-1)(D-2)}}$ and as
$t\to\infty$, then we obtain

$\displaystyle\rho$ $\displaystyle=$ $\displaystyle\frac{(D-1)(D-2)}{2D\kappa
a_{0}^{4}}\left[k\,a_{0}^{2}\beta_{2}\,e^{-2H_{0}t}+(D-2)k^{2}\kappa\sigma_{1}e^{-4H_{0}t}\right],$
(78) $\displaystyle p$ $\displaystyle=$ $\displaystyle-\frac{D-2}{2D\kappa
a_{0}^{4}}\left[(D-3)k\,a_{0}^{2}\,\beta_{2}\,e^{-2H_{0}t}+(D-2)(D-5)k^{2}\kappa\sigma_{1}e^{-4H_{0}t}\right].$
(79)

In these expressions, since the last terms decay faster than the first terms,
one can deduce an equation of state at late times

$p=-\frac{D-3}{D-1}\rho.$ (80)

It can be seen that, for $D\geq 4$ and $\rho>0$, the pressure is always
negative. In $D=4$, this gives $p=-\rho/3$ which corresponds to the equation
of state of cosmic strings [45].

(c) To have positive pressure, we let $\beta_{0}=0$ and $\beta_{2}=0$ together
so that $H_{0}=\sqrt{\frac{2\Lambda}{(D-1)(D-2)}}$ and
$4\kappa\sigma_{1}H_{0}^{2}=D$. We find

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{(D-1)(D-2)^{2}}{2Da_{0}^{4}}k^{2}\sigma_{1}e^{-4H_{0}t},$
(81) $\displaystyle p$ $\displaystyle=$
$\displaystyle-\frac{(D-2)^{2}(D-5)}{2Da_{0}^{4}}k^{2}\sigma_{1}e^{-4H_{0}t}.$
(82)

This special solution can only be obtained when both $k$ and $\sigma_{1}$ are
nonvanishing. Hence it can be obtained neither in Einstein theory nor in the
work of Barrow et al [21]. In other words, it is the effect of the quadratic
gravity which predicts a de Sitter era at late times with the acceleration of
the expansion being constant, i.e. the square of the Hubble constant
($H_{0}=\frac{\dot{a}}{a}$) at the beginning of the accelerating phase.
Furthermore, the above expressions for the energy density and pressure of the
fluid provides an equation of state

$p=-\frac{D-5}{D-1}\rho,$ (83)

which, for $D\geq 6$ and $\rho>0$, the pressure is always negative. It can
also be seen that in $D=5$ the pressure vanishes corresponding to the equation
of state of dust matter, and in $D=4$ it becomes $p=\frac{1}{3}\rho$ which
corresponds to the equation of state of radiation. This solution is valid for
both closed ($k=1$) and open ($k=-1$) universes.

At this point, it should be stressed that all the above equations of state can
mimic a variety of sources of “geometric” origin filling the Universe which
accelerates like the de Sitter universe at late times.

####  5.1.2 A More General Solution

Let us assume that $a(t)$ satisfies the following differential equation

$\dot{a}^{2}=f(a),$ (84)

where $f(a)$ is any arbitrary function of the scale factor $a(t)$. Hence, the
energy density and pressure expressions in (71) and (72) become

$\displaystyle\rho=\frac{(D-1)(D-2)}{2\kappa
a^{4}}\Biggl{\\{}(f+k)a^{2}-\frac{2\Lambda}{(D-1)(D-2)}a^{4}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{\sigma_{1}\kappa}{D}\left[(D-2)k^{2}-(D+2)f^{2}-\frac{D}{4(D-2)}a^{2}f^{\prime
2}-4kf+\frac{D(D-3)}{(D-2)}aff^{\prime}+\frac{D}{(D-2)}a^{2}ff^{\prime\prime}\right]\Biggr{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (85) $\displaystyle
p=-\frac{(D-1)(D-2)}{2\kappa
a^{4}}\Biggl{\\{}\frac{D-3}{D-1}(f+k)a^{2}-\frac{2\Lambda}{(D-1)(D-2)}a^{4}+\frac{1}{(D-1)}a^{3}f^{\prime}$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{\sigma_{1}\kappa}{D(D-1)}\Biggl{[}(D+2)(D-5)f^{2}-(D-2)(D-5)k^{2}+4kaf^{\prime}-\frac{3D(D-3)}{4(D-2)}a^{2}f^{\prime
2}+4(D-5)kf$ $\displaystyle\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{D^{3}-9D^{2}+12D+8}{(D-2)}aff^{\prime}-\frac{2D(D-3)}{(D-2)}a^{2}ff^{\prime\prime}-\frac{D}{2(D-2)}a^{3}(f^{\prime}f^{\prime\prime}+2ff^{\prime\prime\prime})\Biggr{]}\Biggr{\\}},$
(86)

where $f^{\prime}=\frac{df}{da}$. Now we shall give some examples:

Example 1: Let $f=a_{0}k+\frac{a_{1}}{a}+a_{3}a^{2}$. For this choice, the
acceleration of the Universe, $\ddot{a}=f^{\prime}/2$, becomes positive when
$a>(\frac{a_{1}}{2a_{3}})^{1/3}$. Taking $a_{0}=-1$ and
$a_{3}=\frac{2\Lambda}{(D-1)(D-2)}$, we get

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{2(D-4)\sigma_{1}\Lambda^{2}}{D(D-2)^{2}}-\frac{a_{1}(D-1)}{2D\kappa
a^{6}}\Biggl{\\{}D(D-2)a^{3}$ (87)
$\displaystyle+\frac{\kappa\sigma_{1}}{4}\biggl{[}(8D^{2}-19D-16)a_{1}+12D(D-3)ka-\frac{8(D+2)(D-4)\Lambda}{(D-1)(D-2)}a^{3}\biggr{]}\Biggr{\\}},$
$\displaystyle p$ $\displaystyle=$
$\displaystyle-\frac{2(D-4)\sigma_{1}\Lambda^{2}}{D(D-2)^{2}}-\frac{a_{1}}{2D\kappa
a^{6}}\Biggl{\\{}D(D-2)(D-4)a^{3}$ (88)
$\displaystyle+\frac{\kappa\sigma_{1}}{4}\biggl{[}(D-7)(8D^{2}-19D-16)a_{1}+12D(D-3)(D-6)ka-\frac{8(D+2)(D-4)\Lambda}{(D-1)(D-2)}a^{3}\biggr{]}\Biggr{\\}},$

where $a_{1}$ is an arbitrary constant. From these expressions, it can readily
be observed that, when $a_{1}=0$, the equation of state of the fluid reduces
to the form

$p=-\rho,$ (89)

which corresponds to a cosmological constant equation of state for $D\neq 4$.
For $D=4$ and $a_{1}=0$, the energy density and pressure vanish identically.
This means that, for $f=-k+\frac{\Lambda}{3}\,R^{2}$, the FLRW metric is the
vacuum solution of the quadratic gravity field equations. In particular, when
$k=0$, this is the usual de Sitter solution with the scale factor
$a=a_{0}e^{\sqrt{\Lambda/3}t}$.

Example 2: Let $f=a_{0}k+\frac{a_{2}}{a^{2}}+a_{3}a^{2}$. This time, the
acceleration of the Universe, $\ddot{a}=f^{\prime}/2$, becomes positive when
$a>(\frac{a_{2}}{a_{3}})^{1/4}$. Taking $a_{0}=-1$, we obtain

$\displaystyle\rho=\frac{D-1}{2D\kappa}\left[(D-1)(D-4)\kappa\sigma_{1}a_{3}^{2}+D(D-2)\left(a_{3}-\frac{2\Lambda}{(D-1)(D-2)}\right)\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{a_{2}(D-1)}{2D\kappa
a^{8}}\Biggl{\\{}D(D-2)a^{4}-\kappa\sigma_{1}\biggl{[}(D-4)(1+3D)a_{2}-4D(D-4)ka^{2}+2(D^{2}-5D-4)a_{3}a^{4}\biggr{]}\Biggr{\\}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (90) $\displaystyle
p=-\frac{D-1}{2D\kappa}\left[(D-1)(D-4)\kappa\sigma_{1}a_{3}^{2}+D(D-2)\left(a_{3}-\frac{2\Lambda}{(D-1)(D-2)}\right)\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\frac{a_{2}}{2D\kappa
a^{8}}\Biggl{\\{}D(D-2)(D-5)a^{4}-\kappa\sigma_{1}\biggl{[}(D-4)(D-9)(1+3D)a_{2}-4D(D-4)(D-7)ka^{2}$
$\displaystyle\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\
\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\
\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\
\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\ +2(D-5)(D^{2}-5D-4)a_{3}a^{4}\biggr{]}\Biggr{\\}},$ (91)

where $a_{2}$ and $a_{3}$ are arbitrary constants. Here we have the following
special case: When $D=4$,

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{3a_{3}-\Lambda}{\kappa}+\frac{3a_{2}}{\kappa
a^{4}}(1+2\kappa\sigma_{1}a_{3}),$ (92) $\displaystyle p$ $\displaystyle=$
$\displaystyle-\frac{3a_{3}-\Lambda}{\kappa}+\frac{a_{2}}{\kappa
a^{4}}(1+2\kappa\sigma_{1}a_{3}).$ (93)

Now choosing $a_{3}=\Lambda/3$, the equation of state becomes

$p=\frac{1}{3}\rho,$ (94)

which corresponds to an radiation equation of state for arbitrary $a_{2}$.
This means that, for $f=-k+\frac{a_{2}}{a^{2}}+a_{3}a^{2}$, the spatially flat
($k=0$) FLRW metric solves the quadratic gravity field equations with the
scale factor

$a(t)=\frac{1}{\sqrt{2a_{3}}}e^{\sqrt{a_{3}}t}\sqrt{1-a_{2}a_{3}e^{-4\sqrt{a_{3}}t}}.$
(95)

As it can be seen, $a(t)\to\frac{1}{\sqrt{2a_{3}}}e^{\sqrt{a_{3}}t}$ as
$t\to\infty$, which means that, although the equation of state represents a
radiation field, there is a de Sitter like expansions at late times.

Example 3: Let $f=a_{0}e^{ma}$ where $a_{0}$ and $m$ are two real arbitrary
constants. For this function, the scale factor is

$a=-\frac{2}{m}\log\left(1-\frac{m\sqrt{a_{0}}}{2}t\right)$ (96)

where $a>0$ requires $m<0$. Thus, this corresponds to a decelerating universe.
In this case, the energy density and pressure take the form

$\displaystyle\rho$ $\displaystyle=$
$\displaystyle\frac{a_{0}^{2}(D-1)\sigma_{1}\left[3m^{2}Da^{2}+4m(D-3)Da-4\left(D^{2}-4\right)\right]}{8Da^{4}}e^{2ma}+\frac{a_{0}(D-1)(D-2)}{2\kappa
a^{2}}e^{ma}-\frac{\Lambda}{\kappa},$ (97) $\displaystyle p$ $\displaystyle=$
$\displaystyle-\frac{a_{0}^{2}\sigma_{1}\left[6m^{3}Da^{3}+11m^{2}(D-3)Da^{2}+4m\left(D^{3}-9D^{2}+12D+8\right)a-4(D-5)(D-2)(D+2)\right]}{8Da^{4}}e^{2ma}$
(98) $\displaystyle-\frac{a_{0}(D-2)(ma+D-3)}{2\kappa
a^{2}}e^{ma}+\frac{\Lambda}{\kappa}.$

Example 4: Let $f=a_{0}a^{m}$ where $a_{0}$ and $m$ are two arbitrary real
constants. This function produces the scale factor as

$a=\left[\left(\frac{2-m}{2}\right)\sqrt{a_{0}}\right]^{\frac{2}{2-m}}t^{\frac{2}{2-m}}$
(99)

where $a_{0}>0$ and $m\neq 2$. From this scale factor, one obtains the
acceleration as $\ddot{a}=\frac{a_{0}m}{2}a^{m-1}$. To assure the positiveness
of $a$ and $\ddot{a}$, it must be that $0<m<2$. In this case, the energy
density and pressure are

$\displaystyle\rho=\frac{1}{\kappa}\left[\frac{(D-2)(D-1)\left(a_{0}a^{m}+k\right)}{2a^{2}}-\Lambda\right]$
$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
+\frac{(D-1)\sigma_{1}}{8Da^{4}}\Biggl{\\{}a_{0}^{2}\left[3m^{2}D+4m(D-4)D-4(D-2)(D+2)\right]a^{2m}$
$\displaystyle\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\
\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\
\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\
-16a_{0}k(D-2)a^{m}+4k^{2}(D-2)^{2}\Biggr{\\}},$ (100) $\displaystyle
p=-\frac{1}{\kappa}\left[\frac{1}{2}a_{0}(D-2)(m+D-3)a^{m-2}+\frac{k(D-1)(D-5)}{2a^{2}}-\Lambda\right]$
$\displaystyle-\frac{\sigma_{1}}{8Da^{4}}\Biggl{\\{}a_{0}^{2}\left[6m^{3}D+m^{2}D(11D-47)+4m\left(D^{3}-11D^{2}+20D+8\right)-4(D+2)(D-2)(D-5)\right]a^{2m}$
$\displaystyle\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\
\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\
-16a_{0}k(D-2)(m+D-5)a^{m}+4k^{2}(D-5)(D-2)^{2}\Biggr{\\}}.$ (101)

It can be observed that, when $t\to\infty$, the contributions of the higher
curvature terms decay faster than the Einstein terms, and hence

$\displaystyle\rho$ $\displaystyle\to$
$\displaystyle\frac{1}{\kappa}\left[\frac{(D-2)(D-1)\left(a_{0}a^{m}+k\right)}{2a^{2}}-\Lambda\right]$
$\displaystyle p$ $\displaystyle\to$
$\displaystyle-\frac{1}{\kappa}\left[\frac{1}{2}a_{0}(D-2)(m+D-3)a^{m-2}+\frac{k(D-1)(D-5)}{2a^{2}}-\Lambda\right].$

Now, for the flat universe and $\Lambda=0$, we can deduce that

$w\equiv\frac{p}{\rho}=-\frac{D-3+m}{D-1}$ (102)

which reduces to $w=-\frac{1}{3}(1+m)$ in $D=4$. Since $0<m<2$, it must be
that $-1<w<-\frac{D-3}{D-1}$. One can see that as $D$ increases, the range
squeezes and the upper bound approaches to -1.

On the other hand, at early times (as $t\to 0$), the only conditions on the
parameters of the scale factor are $a_{0}>0$ and $m<2$. There are the
following cases:

  * •

For $m<0$, we have

$\displaystyle\rho$ $\displaystyle\to$
$\displaystyle\frac{(D-1)\sigma_{1}a_{0}^{2}a^{2m}}{8Da^{4}}\left[3m^{2}D+4m(D-4)D-4(D-2)(D+2)\right],$
(103) $\displaystyle p$ $\displaystyle\to$
$\displaystyle-\frac{(D-5+2m)\sigma_{1}a_{0}^{2}a^{2m}}{8Da^{4}}\left[3m^{2}D+4m(D-4)D-4(D-2)(D+2)\right].$
(104)

Then one can deduces the following equation of state

$w\equiv\frac{p}{\rho}=-\frac{D-5+2m}{D-1},$ (105)

which becomes $w=\frac{1-2m}{3}>\frac{1}{3}$ in $D=4$. For $m=-1$, it gives
the stiff fluid equation of state, $w=1$.

  * •

For $m=0$, we have

$\displaystyle\rho$ $\displaystyle\to$
$\displaystyle\frac{(D-1)(D-2)\sigma_{1}}{2Da^{4}}\left[k^{2}(D-2)-4a_{0}k-(D+2)a_{0}^{2}\right],$
(106) $\displaystyle p$ $\displaystyle\to$
$\displaystyle-\frac{\sigma_{1}(D-2)(D-5)}{2Da^{4}}\left[k^{2}(D-2)-4a_{0}k-(D+2)a_{0}^{2}\right].$
(107)

Then

$w\equiv\frac{p}{\rho}=-\frac{D-5}{D-1}.$ (108)

This becomes $w=\frac{1}{3}$ in $D=4$ and $w=0$ in $D=5$. The former
represents the radiation and the latter represents the dust matter.

  * •

For $0<m<2$, we have

$\displaystyle\rho$ $\displaystyle\to$
$\displaystyle\frac{(D-1)(D-2)^{2}\sigma_{1}k^{2}}{2Da^{4}}$ (109)
$\displaystyle p$ $\displaystyle\to$
$\displaystyle-\frac{(D-2)^{2}(D-5)\sigma_{1}k^{2}}{2Da^{4}}.$ (110)

Then

$w\equiv\frac{p}{\rho}=-\frac{D-5}{D-1}.$ (111)

Again this represents the radiation ($w=\frac{1}{3}$) in $D=4$ and the dust
matter ($w=0$) in $D=5$.

####  5.1.3 Approximate Solutions Near Criticality with $\sigma_{2}=0$

Now, we assume that, in generic gravity theories in (19), the action contains
some number of coupling constants and hence $A$ and $B$ in (26) and (27) are
functions of these coupling constants and/or some combinations $\alpha_{i}$ of
them. Assuming that these coupling constants are relatively smaller than the
other parameters in these functions, we can expand the scale factor $a$ in
terms of these parameters

$a(t,\alpha_{i})=a_{0}(t)+\sum_{i}\alpha_{i}a_{i}(t)+\mathcal{O}(\alpha_{i}^{2}).$
(112)

Following this approach, we obtain the functions $A$ and $B$ as

$\displaystyle A=\sum_{i}\alpha_{i}A_{i}+\mathcal{O}(\alpha_{i}^{2}),$ (113)
$\displaystyle B=\sum_{i}\alpha_{i}B_{i}+\mathcal{O}(\alpha_{i}^{2}),$ (114)

where $A_{i}$ and $B_{i}$ are functions depending on the explicit gravity
theory. In what follows we shall keep only the terms linear in $\alpha_{i}$.
Hence Eqs.(26) and (27) reduce to

$\displaystyle\rho-\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\rho_{10}-\Lambda\right]+\sum_{i}\alpha_{i}\left[A_{i}-B_{i}-\frac{(D-1)(D-2)}{2\kappa}\rho_{1i}\right]=0,$
(115) $\displaystyle
p+\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\rho_{10}-(D-2)\rho_{20}-\Lambda\right]-\sum_{i}\alpha_{i}\left[A_{i}-\frac{(D-2)}{\kappa}\left(\frac{D-1}{2}\rho_{1i}-\rho_{2i}\right)\right]=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (116)

where

$\displaystyle\rho_{10}=H_{0}^{2}+\frac{k}{a_{0}^{2}},$ (117)
$\displaystyle\rho_{20}=H_{0}^{2}+\frac{k}{a_{0}^{2}}-\frac{\ddot{a}_{0}}{a_{0}}=\dot{H_{0}}+\frac{k}{a_{0}^{2}},$
(118)
$\displaystyle\rho_{1i}=\frac{6}{a_{0}}\left(\dot{a}_{i}H_{0}-a_{i}H_{0}^{2}-\frac{ka_{i}}{a_{0}^{2}},\right)$
(119)
$\displaystyle\rho_{2i}=\frac{2}{a_{0}}\left[\dot{a}_{i}H_{0}-a_{i}H_{0}^{2}-\frac{ka_{i}}{a_{0}^{2}}+\frac{1}{2}\left(\frac{\ddot{a}_{0}a_{i}}{a_{0}}+\ddot{a}_{i}\right)\right],$
(120)

where $H_{0}=\dot{a}_{0}/a_{0}$. Now one can find $a_{i}$ from (115) by
assuming that the extra terms in $\alpha_{i}$ vanish, i.e.,

$\sum_{i}\alpha_{i}\left[A_{i}-B_{i}-\frac{(D-1)(D-2)}{2\kappa}\rho_{1i}\right]=0.$
(121)

Remark 15: Then the expression for $\rho$ becomes exactly the energy density
without the modification of the generic gravity, but the expression for $p$
has contributions from the generic gravity theory.

Now we will specifically consider the critical quadratic gravity theory given
by the action (41) with $\sigma_{2}=0$ in (53). Taking $i=1$ and
$\alpha_{i}=\sigma_{1}$ in the above formulation, from (112) the scale factor
becomes

$a(t,\sigma_{1})=a_{0}(t)+\sigma_{1}a_{1}(t)+\mathcal{O}(\sigma_{1}^{2}).$
(122)

and the field equation (115) reads as

$\displaystyle\rho-\frac{1}{\kappa}\left[\frac{(D-1)(D-2)}{2}\left(H_{0}^{2}+\frac{k}{a_{0}^{2}}\right)-\Lambda\right]-\frac{(D-1)(D-2)\sigma_{1}}{\kappa
a_{0}}\Biggl{\\{}\dot{a}_{1}H_{0}-a_{1}H_{0}^{2}-\frac{ka_{1}}{a_{0}^{2}}$
$\displaystyle\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{\kappa}{2D(D-2)a_{0}^{3}}\Biggl{[}k^{2}(D-2)^{2}+\dot{a}_{0}^{2}\left[2D(D-3)a_{0}\ddot{a}_{0}-4k(D-2)\right]$
$\displaystyle\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\
\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\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\left(D^{2}-4\right)\dot{a}_{0}^{4}-Da_{0}^{2}\ddot{a}_{0}^{2}+2Da_{0}^{2}\dot{a}_{0}{\dddot{a}}_{0}\Biggr{]}\Biggr{\\}}=0.$
(123)

To determine $a_{1}$, we assume that the coefficient terms of $\sigma_{1}$
vanish, ie.,

$\displaystyle\dot{a}_{1}H_{0}-a_{1}\left(H_{0}^{2}+\frac{k}{a_{0}^{2}}\right)+\frac{\kappa}{2D(D-2)a_{0}^{3}}\Biggl{[}k^{2}(D-2)^{2}+\dot{a}_{0}^{2}\left[2D(D-3)a_{0}\ddot{a}_{0}-4k(D-2)\right]$
$\displaystyle\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\
\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\
\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\
-\left(D^{2}-4\right)\dot{a}_{0}^{4}-Da_{0}^{2}\ddot{a}_{0}^{2}+2Da_{0}^{2}\dot{a}_{0}{\dddot{a}}_{0}\Biggr{]}=0,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (124)

which can also be obtained from (121). This equation can be rewritten in the
following 1st order linear differential equation form

$\dot{a}_{1}+R(t)a_{1}=S(t),$ (125)

where

$\displaystyle R(t)=-\left(H_{0}+\frac{k}{a_{0}^{2}H_{0}}\right),$ (126)
$\displaystyle
S(t)=-\frac{\kappa}{2D(D-2)a_{0}^{2}\dot{a_{0}}}\Biggl{\\{}k^{2}(D-2)^{2}+\dot{a}_{0}^{2}\left[2D(D-3)a_{0}\ddot{a}_{0}-4k(D-2)\right]$
$\displaystyle\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\
\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\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
-\left(D^{2}-4\right)\dot{a}_{0}^{4}-Da_{0}^{2}\ddot{a}_{0}^{2}+2Da_{0}^{2}\dot{a}_{0}{\dddot{a}}_{0}\Biggr{\\}}.$
(127)

The equation (125) admits the general solution

$a_{1}(t)=\frac{C}{\lambda(t)}+\frac{\mu(t)}{\lambda(t)},$ (128)

where $C$ is an integration constant and

$\lambda(t)=e^{\int R(t)dt},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\mu(t)=\int\lambda(t)S(t)dt.$ (129)

Taking $D=4$, we consider the following cases:

  1. 1.

Letting $a_{0}(t)=b_{0}e^{c_{0}t}$, where $b_{0}$ and $c_{0}$ are constants,
we find $a_{1}(t)$ as

$a_{1}(t)=e^{c_{0}t}\left[Ce^{-d_{0}e^{-2c_{0}t}}+\frac{\kappa}{4b_{0}}\left(ke^{-2c_{0}t}-4c_{0}^{2}b_{0}^{2}\right)\right],\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ d_{0}=\frac{k}{2c_{0}^{2}b_{0}^{2}}$ (130)

Hence $a(t,\sigma_{1})$ in (122) reads as

$a(t,\sigma_{1})=b_{0}e^{c_{0}t}+\sigma_{1}e^{c_{0}t}\left[Ce^{-d_{0}e^{-2c_{0}t}}+\frac{\kappa}{4b_{0}}\left(ke^{-2c_{0}t}-4c_{0}^{2}b_{0}^{2}\right)\right].$
(131)

For a flat universe $(k=0)$, it reduces to

$a(t,\sigma_{1})=\chi_{0}e^{c_{0}t},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\chi_{0}=b_{0}+\sigma_{1}\left(C-\kappa c_{0}^{2}b_{0}\right).$ (132)

Then, the contribution of the higher curvature terms in the scale factor for a
flat universe has the same exponential form $e^{c_{0}t}$ as in GR. For
$\chi_{0}>0$, the total acceleration, i.e $\ddot{a}(t,\sigma_{1})$, is
positive.

  2. 2.

Letting $a_{0}(t)=m_{0}t^{n}$, where $m_{0}$ and $n$ are constants, we find
$a_{1}(t)$ in (128) with

$\displaystyle\lambda(t)=t^{-n}\cdot e^{bt^{2(1-n)}},$ (133)
$\displaystyle\mu(t)=\int e^{bt^{2(1-n)}}\left(-\frac{3}{4}m_{0}\kappa
n(1-2n)t^{-3}+\frac{\kappa nk}{2m_{0}}t^{-2n-1}-\frac{\kappa
k^{2}}{4nm_{0}^{3}}t^{1-4n}\right)dt,$ (134)

where $b=-\frac{k}{2nm_{0}^{2}(1-n)}$. Here we consider the following two
cases.

     * •

$k=0,\leavevmode\nobreak\ \leavevmode\nobreak\ n=\frac{1}{2}$ corresponding to
the radiation era in a flat universe in the context of GR. For this case, we
obtain

$a_{1}(t)=Ct^{\frac{1}{2}},$ (135)

and then

$a(t,\sigma_{1})=\chi_{0}t^{\frac{1}{2}},\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \chi_{0}=m_{0}+\sigma_{1}C$ (136)

Then, the contribution of the higher curvature terms in the total scale factor
for a flat universe has the same power law form $t^{\frac{1}{2}}$ as in GR. In
contrast to the previous case of the scale factor with an exponential form,
here we see that the total acceleration
$\ddot{a}(t,\sigma_{1})=-\frac{1}{4}\chi_{0}t^{-\frac{3}{2}}$ is positive if
$\chi_{0}<0$. For $\chi_{0}>0$, the higher curvature modifications cannot
support the acceleration of the universe.

     * •

$k=0,\leavevmode\nobreak\ \leavevmode\nobreak\ n=\frac{2}{3}$ corresponding to
the dust matter era in a flat universe in the context of GR. For this case, we
obtain

$a_{1}(t)=Ct^{\frac{2}{3}}-\frac{1}{12}\kappa m_{0}t^{-\frac{4}{3}},$ (137)

and then

$a(t,\sigma_{1})=\chi_{0}t^{\frac{2}{3}}-\frac{1}{12}\kappa
m_{0}t^{-\frac{4}{3}},\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \chi_{0}=m_{0}+\sigma_{1}C$ (138)

In contrast to two previous cases, here we see that the contribution of the
higher curvature terms in the scale factor are not of the same kind of GR, and
there is an extra $t^{-\frac{4}{3}}$ type of contribution. The total
acceleration has the form
$\ddot{a}(t,\sigma_{1})=-\frac{2}{9}\chi_{0}t^{-\frac{4}{3}}-\frac{28}{108}\sigma_{1}\kappa
m_{0}t^{-\frac{10}{3}}$ which is negative for $(m_{0},C,\sigma_{1},\kappa)>0$.
Then, the higher curvature modifications cannot support the acceleration of
the universe in a matter dominated era.

###  5.2 Solutions with $\sigma_{1}=0$

When $\sigma_{1}=0$ and $D\neq 4$, the field equations (46) and (47) become

$\displaystyle\frac{1}{\kappa}\left[\frac{1}{2}(D-2)(D-1)\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)-\Lambda\right]=\rho-(D-1)\sigma_{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)^{2},$
(139)
$\displaystyle-\frac{1}{\kappa}\left[\frac{(D-2)(D-3)}{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)+(D-2)\frac{\ddot{a}}{a}-\Lambda\right]$
$\displaystyle\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\
\leavevmode\nobreak\ \leavevmode\nobreak\
=p+\sigma_{2}\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)\left[(D-5)\left(\frac{k}{a^{2}}+\frac{\dot{a}^{2}}{a^{2}}\right)+4\frac{\ddot{a}}{a}\right].$
(140)

In this case, when $[(D-1)(D-2)\alpha+D(D-3)\gamma]=0$ (i.e. $\sigma_{2}=0$),
one again recovers the expressions in the Einstein gravity in $D\neq 4$
dimensions.

To see the effects of quadratic gravity terms on the expansion of the universe
at late times, neglecting $\rho$ and $p$ of baryonic matter in (139) and
(140). Then, from (139), one can immediately obtain the following solutions
for the scale factor:

$a(t)=\begin{cases}\frac{\sinh[\sqrt{h_{0}}(t-t_{0})]}{\sqrt{h_{0}}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for $k=-1$},\\\ \\\
e^{\sqrt{h_{0}}(t-t_{0})}\leavevmode\nobreak\ \leavevmode\nobreak\ \text{for
$k=0$},\\\ \\\
\frac{\cosh[\sqrt{h_{0}}(t-t_{0})]}{\sqrt{h_{0}}}\leavevmode\nobreak\
\leavevmode\nobreak\ \text{for $k=+1$},\end{cases}$ (141)

where $t_{0}$ is an arbitrary integration constant and

$h_{0}\equiv-\frac{D-2}{4\kappa\sigma_{2}}\left[1\mp\sqrt{1+\frac{16\sigma_{2}\Lambda}{(D-1)(D-2)^{2}}}\right]$
(142)

with $\sigma_{2}\Lambda\geq-\frac{(D-1)(D-2)^{2}}{16}$. One can verify that,
with this solution, (140) is identically satisfied. The energy density and
pressure of the geometric fluid becomes

$\displaystyle\rho_{g2}$ $\displaystyle=$
$\displaystyle-\sigma_{2}(D-1)h_{0}^{2},$ (143) $\displaystyle p_{g2}$
$\displaystyle=$ $\displaystyle\sigma_{2}(D-1)h_{0}^{2}.$ (144)

Remark 16: These give an equation of state $p_{g2}=-\rho_{g2}$ which
corresponds to the vacuum equation of state. It must be observed that the
positiveness of the energy density requires $\sigma_{2}<0$. Also it should be
stated that in the absence of cosmological constant ($\Lambda=0$), the higher
curvature terms in the theory behaves like en effective cosmological constant
driving the late time exponentially accelerating expansion. These results are
consistent with the ones discussed in Case II in section 4.

##  6 Conclusion

We have given all our findings both in abstract and in the introduction. Here
we give a short summary of this work. We consider FLRW cosmology in the
context of generic gravity theories in which the action includes all the
combinations of the metric tensor, curvature tensor, and covariant derivatives
of the curvature tensor of any order. Very recently we showed that in such
theories with FLRW geometry, contributions of all higher order terms reduce to
a perfect fluid form which we call it now as the geometric fluid. Hence all
generic theories of gravity in FLRW geometry are equivalent to Einstein’s
theory of general relativity where the source term contains both matter and
geometric fluids. We propose that the source of dark energy/matter is this
geometrical fluid arising from higher order gravity theories. Choosing any
higher order gravity, the parameters of the theory can be suitably arranged
that the corresponding geometric fluid contributes to the accelerated
expansion of the universe. We verified our assertion by taking the quadratic
gravity as an example. Furthermore we have given some particular exact
cosmological solutions of quadratic gravity theory with matter and geometrical
fluids.

## References

  * [1] T. Damour, and A. M. Polyakov, String theory and gravity, General Relativity and Gravitation 26, 1171 (1994).
  * [2] GT. Horowitz, M. Kolanowski, G.N. Remnen, and J.E. Santos, Extremal Kerr Black Holes as Amplifiers of New Physics, Phys. Rev. Lett. 131, 091402 (2023).
  * [3] D.G. Boulware, S. Deser, Classical general relativity derived from quantum gravity, Annals of Physics, 89, 193 (1975).
  * [4] C. Teitelboim, Quantum mechanics of the gravitational field, Phys. Rev. D 25, 3159 (1982).
  * [5] C. Rovelli, Quantum Gravity, Cambridge university press (2004).
  * [6] J.D. Bekenstein, Alternatives to dark matter: Modified gravity as an alternative to dark matter, arXiv:1001.3876 (2010).
  * [7] S. Capozziello, and M. De Laurentis, The dark matter problem from $f(R)$ gravity viewpoint, Annalen der Physik 524, 545 (2012).
  * [8] C.G.Böhmer, T. Harko, and F.S.N. Lobo, Dark matter as a geometric effect in $f(R)$ gravity, Astroparticle Physics 29, 386 (2008).
  * [9] R. Zaregonbadi, M. Farhoudi, and N. Riazi, Dark matter from $f(R,T)$ gravity, Phys. Rev. D 94, 084052 (2016).
  * [10] S. Nojiri and S.D. Odintsov, Introduction to modified gravity and gravitational alternative for dark energy, International Journal of Geometric Methods in Modern Physics 4, 115 (2007).
  * [11] V. Sahni, and Y. Shtanov, Braneworld models of dark energy, Journal of Cosmology and Astroparticle Physics 11, 014 (2003).
  * [12] R. Maartens, and K. Koyama, Brane-world gravity, Living Reviews in Relativity 13, 1 (2010).
  * [13] E. Bertschinger and P. Zukin, Distinguishing modified gravity from dark energy, Phys. Rev. D 78, 024015 (2008).
  * [14] K. Koyama, Cosmological tests of modified gravity, Rep. Prog. Phys. 79, 046902 (2016).
  * [15] S. M. Carroll, A. De Felice, V. Duvvuri, D. A. Easson, M. Trodden, and M. S. Turner, Cosmology of generalized modified gravity models, Phys. Rev. D 71, 063513 (2005).
  * [16] J. Harada, Gravity at cosmological distances: Explaining the accelerating expansion without dark energy, Phys. Rev. D 108, 044031 (2023).
  * [17] S. Nojiri, and S.D. Odintsov, Modified gravity and its reconstruction from the universe expansion history, Journal of Physics: Conference Series 66, 1, IOP Publishing (2007).
  * [18] D. Wands, Extended gravity theories and the Einstein-Hilbert action, Class. Quantum Grav. 11, 269 (1994).
  * [19] H.-J. Schmidt, Variational derivatives of arbitrarily high order and multi-inflation cosmological models. Class. Quantum Grav. 7, 1023 (1990).
  * [20] S. Capozziello, M. De Laurentis, Extended theories of gravity, Phys. Rep. 509, 167 (2011).
  * [21] J. D. Barrow and S. Hervik, Evolution of Universes in Quadratic Theories of Gravity, Phys. Rev. D 74, 124017 (2006).
  * [22] S. Nojiri, S.D. Odintsov, Modified Gauss-Bonnet theory as gravitational alternative for dark energy, Phys. Lett. B631, 1 (2005).
  * [23] S. Das, N. Banerjee, and N. Dadhich, Curvature-driven acceleration: a utopia or reality?, Class. Quantum Grav. 23, 4159 (2006).
  * [24] B. C. Paul, P. S. Debnath, and S. Ghose, Accelerating universe in modified theories of gravity, Phys. Rev. D 79, 083534 (2009).
  * [25] D. Lovelock, The Einstein tensor and its generalizations, Journal of Mathematical Physics, 12, 498 (1971).
  * [26] A. Casalino, A. Colleaux, M. Rinaldi, and S. Vicentini, Regularized Lovelock Gravity, Physics of the Dark Universe 31, 100770 (2021).
  * [27] M. Gürses and Y. Heydarzade, FLRW-Cosmology in generic gravity theories, Eur. Phys. J. C. 80, 1061 (2020).
  * [28] B. Y. Chen, A simple characterization of generalized Robertson Walker space-times, Gen. Relativ. Gravit. 46, 1833 (2014).
  * [29] C. A. Mantica and L.G. Molinari, On the Weyl and Ricci tensors of Generalized Robertson-Walker space-times, Journal of Mathematical Physics, 57 102502 (2016).
  * [30] C. A. Mantica, L.G. Molinari and U. Chand De, A condition for a perfect-fluid space-time to be a generalized Robertson-Walker space-time, Journal of Mathematical Physics, 57 022508 (2016).
  * [31] S. Capozziello, C.A. Mantica, L.G. Molinari, Cosmological perfect fluids in $f(R)$ gravity, Int. J. Geom. Methods Mod. Phys. 16, 1950008 (2019).
  * [32] S. Capozziello, C. A. Mantica and L.G. Molinari, Cosmological Perfect fluids in Gauss-Bonnet Gravity, International Journal of Geometrical Methods in Modern Physics, 16, 1950133 (2019).
  * [33] S. Capozziello, C.A. Mantica, and L.G. Molinari, Cosmological perfect fluids in higher-order gravity", General Relativity and Gravitation 52, 36 (2020).
  * [34] S. Capozziello, C.A. Mantica, L.G. Molinari, Geometric perfect fluids from extended gravity, Europhysics Letters, 137, 19001 (2022).
  * [35] S. Weinberg, Cosmology, Oxford University Press, Oxford (2008).
  * [36] M. Visser, Jerk, snap and the cosmological equation of state, Class. Quantum Grav. 21, 2603 (2004).
  * [37] I. Güllü and B. Tekin, Massive higher derivative gravity in D-dimensional anti-de Sitter spacetimes, Phys. Rev. D 80, 064033 (2009).
  * [38] S. Deser and B. Tekin, Gravitational energy in quadratic-curvature gravities, Phys. Rev. Lett. 89, 101101 (2002); Energy in generic higher curvature gravity theories, Phys. Rev. D 67, 084009 (2003).
  * [39] H. Lü and C. N. Pope, Critical Gravity in Four Dimensions, Phys. Rev. Lett. 106, 181302 (2011).
  * [40] S. Deser, H. Liu, H. Lü, C. N. Pope, T. Ç. Şişman, and B. Tekin, Critical points of $D$-dimensional extended gravities, Phys. Rev. D 83, 061502(R) (2011).
  * [41] A. G. Riess et al., Observational evidence from supernovae for an accelerating universe and a cosmological constant, Astron. J. 116, 1009 (1998).
  * [42] S. Perlmutter et al., Measurements of $\omega$ and $\Lambda$ from 42 high-redshift supernovae, Astrophys. J. 517, 565 (1999).
  * [43] N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641, A6 (2020).
  * [44] L. Liu, L. Qiao, B. Chang, and L. Xu, Revisiting cosmography via Gaussian process, Eur. Phys. J. C. 83, 374 (2023).
  * [45] A. Vilenkin, Cosmic strings, Phys. Rev. D 24, 2082 (1981).