a(t)\propto\exp\left(\left(\frac{2}{\gamma}-1\right)\sqrt{\frac{R_{s}}{12}}\left(t_{s}-t\right)^{-\gamma/2+1}\right)\,.$
(212)

When $1<\gamma<2$, $a(t)$ evolves as

$\displaystyle
a(t)\propto\exp\left(\frac{R_{s}}{6\gamma\left(\gamma-1\right)}\left(t_{s}-t\right)^{-\gamma}\right)\,.$
(213)

And if $\gamma<1$, we get

$\displaystyle
a(t)\propto\exp\left(H_{s}t+\frac{R_{s}}{6\gamma\left(\gamma-1\right)}\left(t_{s}-t\right)^{-\gamma}\right)\,.$
(214)

However, we see that a sudden future singularity appears at $t=t_{s}$
Barrow:2004xh ; Nojiri:2004pf ; Barrow:2004he ; Fernandez-Jambrina:2004yjt
when $\gamma<2$. Now, dealing with the case $\gamma<1$, as the second term in
Eq. (210) is smaller than the first term, one may solve the differential Eq.
(194) asymptotically in the following way

$\displaystyle P\sim
P_{s}\left(1+\frac{R_{s}}{3\beta(1-\beta)}\left(t_{s}-\phi\right)^{1-\beta}\right)\,,$
(215)

where $P_{s}$ is a constant, which finally leads to

$\displaystyle F(R)\sim F_{0}R+F_{1}R^{2\beta/\left(\beta+1\right)}\,.$ (216)

Now, since for $F(R)$ gravity, one can introduce the effective energy density
and effective pressure, see Eqs. (175), (176), (177), thus, for different
values of $\beta$ of the Hubble rate (114), the nature of the singularities
can be studied as follows. When $\beta>1$, as $t\to t_{s}$, we see that

$\displaystyle
a\sim\exp(h_{s}\left(t_{s}-t\right)^{1-\beta}/\left(\beta-1\right))\to\infty,\;\;\mbox{and
as a
consequence},\;\rho_{\mathrm{eff}}\to\infty,\,|p_{\mathrm{eff}}|\to\infty,$
(217)

which means that a Big Rip (Type I) singularity occurs. If $0<\beta<1$, $a$
goes to a constant but
$\rho_{\mathrm{eff}}\to\infty,\,|p_{\mathrm{eff}}|\to\infty$, that means a
Type III singularity occurs. If $-1<\beta<0$, then $a$ and
$\rho_{\mathrm{eff}}$ vanish but $|p_{\mathrm{eff}}|\to\infty$ which means
that a Type II singularity occurs. When $\beta<0$, instead of Eq. (114), as in
Eq. (210), one may assume, as in (115) with $H_{s}(t)$ is a constant,

$\displaystyle H\sim H_{s}+h_{s}\left(t_{s}-t\right)^{-\beta}\,.$ (218)

Now, if $-1<\beta<0$, for $t\to t_{s}$, $\rho_{\mathrm{eff}}$ approaches to
finite value $3H_{s}^{2}/\kappa^{2}$ but $|p_{\mathrm{eff}}|$ diverges, so we
have a Sudden singularity. If $\beta<-1$ but $\beta$ is not an integer, then
$a$ remains finite and $\rho_{\mathrm{eff}}$, $p_{\mathrm{eff}}$ vanish if
$H_{s}=0$ or $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are finite if
$H_{s}\neq 0$, however, higher derivatives of $H$ diverge. That means in this
case a Type IV singularity occurs. Thus, we see that $F(R)$ gravity may allow
various type of finite-time future singularities. This is not unnatural
because in the context of modified gravity, one can find an effective
phantom/quintessence phase Nojiri:2006ri  and it is well-known that a
phantom/quintessence-dominated Universe may end up with finite-time future
singularities of various types. Hence, with the reconstruction of the $F(R)$
gravity from a given cosmology one can find the possible functional forms of
$F(R)$ that may lead to finite-time future singularities. For example, from
the present discussions, one can see that $F(R)=R+\widetilde{\alpha}R^{n}$
with $n>2$ leads to a Type I singularity and $F(R)=R-\widetilde{\beta}R^{-n}$
with $n>0$ leads to a Type III singularity where $\widetilde{\alpha}$ and
$\widetilde{\beta}$ are any real numbers.

####  IV.5.1 Occurrence of Singularities in Different Frames: $F(R)$ Gravity

The choice of the physical frame in the context of $F(R)$ gravity is an
important point because cosmology of an $F(R)$ gravity in one frame could be
different in other frame. The accelerating phase of the Universe in one frame
may correspond to its decelerating phase Bahamonde:2017kbs . Also, the type of
a singularity changes from one frame to other frame Briscese:2006xu ;
Bahamonde:2016wmz . In this section, we shall discuss the second possibility
in detail, that means, how the choice of the frames affects the type of
finite-time future singularities appearing in $F(R)$ gravity theory. We start
with the vacuum $F(R)$ gravity in the Jordan frame whose action is given by
(174) when we neglect the contribution from matter by putting
${L}_{\mathrm{matter}}=0$ Bahamonde:2016wmz . Varying the action (174) with
${L}_{\mathrm{matter}}=0$ with respect to the metric $g_{\mu\nu}$ representing
a spatially flat FLRW geometry (3), one can obtain the following gravitational
field equations, which is equivalent to (175) with (176) and (177) when we
neglect the matter,

$\displaystyle 0=$
$\displaystyle\,-\frac{F(R)}{2}+3\left(H^{2}+\dot{H}\right)F^{\prime}(R)-18\left(4H^{2}\dot{H}+H\ddot{H}\right)F^{\prime\prime}(R)\,,$
(219) $\displaystyle 0=$
$\displaystyle\,\frac{F(R)}{2}-\left(\dot{H}+3H^{2}\right)F^{\prime}(R)+6\left(8H^{2}\dot{H}+4{\dot{H}}^{2}+6H\ddot{H}+\dddot{H}\right)F^{\prime\prime}(R)+36\left(4H\dot{H}+\ddot{H}\right)^{2}F^{\prime\prime\prime}(R)\,,$
(220)

where an overhead dot represents the differentiation with respect to the
cosmic time $t$ in (3) and the prime corresponds to the differentiation with
respect to the Ricci scalar $R$.

Now introducing the auxiliary fields $A$ and $B$, the action of Eq. (174) with
${L}_{\mathrm{matter}}=0$ written in the Jordan frame can be expressed into an
equivalent action as follows Maeda:1988ab  (see also Chiba:2003ir )

$\displaystyle S=\frac{1}{2\kappa^{2}}\int
d^{4}x\,\sqrt{-g}\Bigg{[}B(R-A)+F(A)\Bigg{]}\,.$ (221)

The action (221) can be varied with respect to the auxiliary scalar $B$
leading to the condition $A=R$, and hence, one can recover the action (174)
with ${L}_{\mathrm{matter}}=0$. Moreover, varying the action (221) with
respect to $A$, one can eliminate the auxiliary field $B$ from it ending up
$B=F^{\prime}(A)$. Hence, the action (221) takes the equivalent form

$\displaystyle S=\frac{1}{2\kappa^{2}}\int
d^{4}x\,\sqrt{-g}\Bigg{[}F^{\prime}(A)(R-A)+F(A)\Bigg{]}\,.$ (222)

Taking a conformal transformation of the metric tensor, one may obtain a
minimally coupled scalar-tensor theory, called the Einstein frame scalar-
tensor theory. We use a particular conformal factor taking the following
expression Bahamonde:2016wmz

$\displaystyle\hat{g}_{\mu\nu}=\frac{1}{F^{\prime}(A)}g_{\mu\nu}\,,$ (223)

which modifies the Ricci scalar as $R\rightarrow\hat{R}$. With the conformal
transformation (223), and then by defining a new scalar field $\sigma$ in
terms of the auxiliary scalar field $A$,

$\displaystyle\sigma=-\ln F^{\prime}(A)\,,$ (224)

the action (222) takes the form

$\displaystyle S=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-\hat{g}}\left\\{\hat{R}-\frac{3}{2}\hat{g}^{\mu\nu}\partial_{\mu}\sigma\partial_{\nu}\sigma-V(\sigma)\right\\}\,,$
(225)

where the potential $V(\sigma)$ takes the form

$\displaystyle
V(\sigma)=\frac{A}{F^{\prime}(A)}-\frac{F(A)}{F^{\prime}(A)^{2}}\,.$ (226)

Notice that with the use of Eq. (224), the potential of Eq. (226) can be
expressed in terms of the scalar field $\sigma$. Thus, corresponding to the
Jordan frame $F(R)$ gravity characterized by the action of Eq. (174) with
${L}_{\mathrm{matter}}=0$, the Einstein frame scalar-tensor theory can be
found as given by the action of Eq. (225). Finally, we note that with the use
of the following transformation

$\displaystyle\varphi=\sqrt{\frac{3}{2\kappa^{2}}}\sigma\,,$ (227)

the action (225) can be expressed into the canonical form Bahamonde:2016wmz

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left\\{\frac{R}{2\kappa^{2}}-\frac{1}{2}\partial_{\mu}\varphi\partial^{\mu}\varphi-V(\varphi)\right\\}\,,$
(228)

Alternatively, starting with the scalar-tensor canonical scalar field action
given by Bahamonde:2016wmz

$\displaystyle S=\int
d^{4}x\sqrt{-\hat{g}}\left\\{\frac{\hat{R}}{2\kappa^{2}}-\frac{1}{2}\hat{\partial}_{\mu}\varphi\hat{\partial}^{\mu}\varphi-V(\varphi)\right\\}\,,$
(229)

one can find its equivalent Jordan frame $F(R)$ gravity action. Assuming a
spatially flat FLRW line element (3), the FLRW equations corresponding to the
action (229) can be written as

$\displaystyle 3\widetilde{H}^{2}=\frac{1}{2}\dot{\varphi}^{2}+V\,,\quad
3\widetilde{H}^{2}+2\dot{\widetilde{H}}=-\frac{1}{2}\dot{\varphi}^{2}+V\,.$
(230)

where we have used ‘tilde’ in the physical quantities in this frame (i.e., the
Einstein frame) in order to differ the corresponding physical quantities in
the Jordan frame.

Now one can map the above action (229) to a modified $F(R)$ gravity theory. In
order to do this, we need to perform the conformal transformation given by
$g_{\mu\nu}\rightarrow\mathrm{e}^{\pm\sqrt{\frac{2}{3}}\kappa\varphi}\hat{g}_{\mu\nu}$,
and as a result, the FLRW metric becomes,

$\displaystyle
ds_{F(R)}^{2}=\mathrm{e}^{\pm\sqrt{\frac{2}{3}}\kappa\varphi}\left(-d\widetilde{t}^{2}+\widetilde{a}^{2}\left(~{}\widetilde{t}~{}\right)\sum_{i=1,2,3}\left(dx^{i}\right)^{2}\right),$
(231)

where we have introduced a new time coordinate $\widetilde{t}$, given by
$dt=\mathrm{e}^{\pm\frac{1}{2}\sqrt{\frac{2}{3}}\kappa\varphi}d\widetilde{t}$,
the solution of which, denoted by $t=f\left(~{}\widetilde{t}~{}\right)$, is an
increasing function. Let us note that if the function
$f\left(~{}\widetilde{t}~{}\right)$ allows singularities, then we may have
some problems. As we are considering two different frames, i.e., the Jordan
frame and the Einstein frame, therefore, following the relation between $t$
and $\widetilde{t}$, the range of the values of $\widetilde{t}$ could be
mapped to a different region in the $t$ coordinate. Let us consider an
interval $[\widetilde{t}_{1},\widetilde{t}_{2}]$ in the Einstein frame, with
the scale factor at $\widetilde{t}=\widetilde{t}_{1}$ being equal to
$\widetilde{a}\left(~{}\widetilde{t}_{1}~{}\right)=0$ and at
$\widetilde{t}=\widetilde{t}_{2}$, the scale factor is equal to,
$\widetilde{a}\left(~{}\widetilde{t}_{2}~{}\right)=0$ or
$\widetilde{a}\left(~{}\widetilde{t}_{2}~{}\right)=\infty$. As one can
identify that $\widetilde{t}=\widetilde{t}_{2}$, corresponds to a Big Crunch
or to a Big Rip singularity, respectively, in which case with potentially
$\widetilde{t}_{1}=-\infty$ and/or $\widetilde{t}_{2}=\infty$, so that these
singularities effectively do not occur, a fact which crucially depends on the
particular form of the scale factor at hand. The new range of $t$ coordinate
will be $[f\left(\widetilde{t}_{1}\right)\,,f\left(\widetilde{t}_{2}\right)]$,
assuming $\phi\left(~{}\widetilde{t}~{}\right)$ is regular everywhere in the
interval $[\widetilde{t}_{1},\widetilde{t}_{2}]$.

In the following we shall discuss the correspondence of the finite-time
singularities in the Jordan and Einstein frames using some simple examples. We
note that while discussing the singularities below, we have followed the same
unit system as in Ref. Bahamonde:2016wmz , that means, we have considered the
same unit system so that the gravitational coupling constant becomes unity,
i.e., $\kappa=1$.

(i) Power law cosmology

We consider the power law cosmology which is described by the following scale
factor,

$\displaystyle\widetilde{a}\left(\widetilde{t}\right)=\widetilde{a}_{c}\left(\frac{\widetilde{t}}{\widetilde{t}_{c}}\right)^{p}\,,$
(232)

where $\widetilde{t}_{c}$ is a fiducial value of the cosmic time, $p$ is a
positive real number and from the above relation we identify that
$\widetilde{a}(~{}\widetilde{t}_{c}~{})=\widetilde{a}_{c}$. Note that such a
power law cosmology described in Eq. (232) is a solution of the Friedmann
equation in the Einstein frame scalar-tensor theory when the potential has the
exponential form. In this case the scalar field evolves as

$\displaystyle\phi=\pm\sqrt{2p}\;\ln\left(\frac{\widetilde{t}}{\widetilde{t}_{c}}\right)\,.$
(233)

In this case $\widetilde{t}_{1}=0$ and $\widetilde{t}_{2}=\infty$ and in this
model the Hubble rate $\widetilde{H}$ diverges at $\widetilde{t}=0$, so we
have a Big Bang singularity in the Einstein frame.

Now in order to understand the behavior of the singularities in the Jordan
frame, we can follow the procedure as described in Section IV.5.1. The new
time variable $t$ in the Jordan frame can be found from the differential
equation

$\displaystyle\frac{dt}{d\widetilde{t}}=\left(\frac{\widetilde{t}}{\widetilde{t}_{c}}\right)^{\pm\sqrt{\frac{p}{3}}}\,,$
(234)

having the following solution,

$\displaystyle
t=\frac{3}{3\pm\sqrt{3p}}\widetilde{t}\left(\frac{\widetilde{t}}{\widetilde{t}_{c}}\right)^{\pm
2\sqrt{\frac{p}{3}}}\,.$ (235)

The corresponding scale factor as a function of the cosmic time $t$ takes the
form

$\displaystyle a(t)\sim t^{w}\quad\mbox{where}\quad w=\frac{\sqrt{3p}\pm
3p}{\sqrt{3p}\pm 3}\,.$ (236)

Now we have the following observations. When the minus sign is chosen in the
conformal factor, the cosmological evolution (236) has a Big Bang singularity
at $t=0$ if the power law parameter $p$ lies in the interval $1/3\leq p<3/4$.
If $p=1/3$, the Jordan frame metric becomes static and we do not have any
finite-time singularity. When $0<p<1/3$, the Big Bang singularity at
$\widetilde{t}=0$ in the Einstein frame becomes the beginning of a contracting
Universe in the Jordan frame. In the case $3/4<p<3$, from (235) we can see the
direction of time is reversed, meaning that this case must be disregarded.
Finally, for $p>3$ the time $\widetilde{t}=0$ correspond to $t=-\infty$, so
the Big Bang singularity disappears in the Jordan frame.

(ii) Cosmology generated by $R^{-n}$ gravity in the Jordan Frame

We now consider the cosmology driven by $R^{-n}$ gravity in the Jordan frame
(see Briscese:2006xu  for details). So, for $F(R)\sim R^{-n}$, from Eqs. (219)
and (220) one can see that the corresponding scale factor takes the power law
form

$\displaystyle a\sim\left(t_{s}-t\right)^{\frac{(n+1)(2n+1)}{n+2}}\,.$ (237)

Therefore, if either $n<-2$ or $-1<n<-1/2$, a Big Rip Type I singularity
appears at $t=t_{s}$ in the Jordan frame, and in the remaining cases, a Type
III Big Crunch singularity is present at this point. In this case, the
corresponding Einstein frame canonical scalar field takes the form

$\displaystyle\sigma\sim(n+1)\ln R\sim-2(n+1)\ln(t_{s}-t)\,,$ (238)

and the Ricci scalar takes the following expression

$\displaystyle R\sim\frac{6(n+1)(2n+1)(4n+5)n}{(n+2)^{2}(t_{s}-t)^{2}}\,.$
(239)

Now the time coordinate $\widetilde{t}$ in the corresponding Einstein frame
scalar-tensor theory is given by

$\displaystyle
d\widetilde{t}=\pm\mathrm{e}^{\frac{1}{2}\sigma}dt\;\;\sim\pm(t_{s}-t)^{-(n+1)}dt\,,$
(240)

which gives $\widetilde{t}=\pm(t_{s}-t)^{-n}$. Therefore, for $n>0$, $t$
approaches $t\rightarrow t_{s}$ in the Jordan frame corresponds to
$\widetilde{t}\rightarrow\pm\infty$ in the Einstein frame. As a consequence,
we observe that the singularity changes its structure abruptly, that means it
does not appear in finite-time in the Einstein frame scalar-tensor theory.
Nevertheless, a new additional singularity may be present, as when $t$
approaches infinity in the Einstein frame, it corresponds to the new time
coordinate $\widetilde{t}\to 0$, and thus any singularities at infinity can be
brought back to a finite-time. On the other hand, when $n<0$, the limit $t\to
t_{s}$ in the Jordan frame corresponds to $\widetilde{t}\to 0$ in the Einstein
frame. Finally, we note that the metric in the Einstein frame scalar-tensor
theory behaves as

$\displaystyle
ds^{2}=\mathrm{e}^{\sigma}\left(-dt^{2}+a^{2}(t)\sum_{i=1,2,3}(dx^{i})^{2}\right)\sim-d\widetilde{t}^{2}+\widetilde{a}^{2}\left(~{}\widetilde{t}~{}\right)\sum_{i=1,2,3}(dx^{i})^{2}\,,\quad\widetilde{a}^{2}\left(~{}\widetilde{t}~{}\right)\sim
a_{s}^{2}\;\widetilde{t}^{~{}\frac{2n(n^{2}-1)}{n+2}}\,,$ (241)

where the constant $a_{s}$ is a real number. In this case, the power of the
scale factor is negative only when $-2<n<-1$ or $0<n<1$, and thus a Big Rip
singularity is present. Thus, for the Big Rip singularity in the Jordan frame,
the scale factor in the Einstein frame behaves as
$\widetilde{a}^{2}\left(~{}\widetilde{t}~{}\right)\to 0$ when
$\widetilde{t}\to 0$, thus, it becomes a Type III Big Crunch singularity in
the Einstein frame.

(iii) A singular cosmological evolution

In this section, we examine how the simplest singular cosmology behaves in
different frames. We consider the following Hubble rate which describes the
simplest singular cosmology as in (114),

$\displaystyle
H\left(~{}\widetilde{t}~{}\right)=h_{s}\left(~{}\widetilde{t}-\widetilde{t_{s}}~{}\right)^{-\beta}\,,$
(242)

where $h_{s}$ is any positive real number and $\beta$ is a real number. From
the values of $\beta$ one can determine the singularity type. In particular,
we realize the following type of singularities based on the values of
$\alpha$:

  * •

For $\beta>1$, we realize a Type I singularity.

  * •

For $0<\beta<1$, a Type III singularity is found.

  * •

When $-1<\beta<0$, a Type II singularity occurs.

  * •

When $\beta<-1$, a Type IV singularity is realized.

Let us assume that the Hubble rate of Eq. (242) is given in the Einstein frame
and we aim to investigate the singular behavior captured in the Hubble rate of
Eq. (242) in the context of Jordan frame $F(R)$ theory. We use the conformal
factor given by $\mathrm{e}^{\sqrt{\frac{2}{3}}\varphi}$ which results in the
transformation of the metric as

$\displaystyle
ds_{F(R)}^{2}=\mathrm{e}^{\sqrt{\frac{2}{3}}\varphi}\left(-d\widetilde{t}^{2}+\widetilde{a}^{2}(~{}\widetilde{t}~{})\sum_{i=1,2,3}\left(dx^{i}\right)^{2}\right)\,,$
(243)

and the scale factor now becomes,

$\displaystyle
a(t)=\mathrm{e}^{\frac{1}{2}\sqrt{\frac{2}{3}}\varphi}\;\widetilde{a}\left(~{}\widetilde{t}~{}\right)\,,$
(244)

where the Jordan frame time parameter $t$ is defined as

$\displaystyle
dt=\mathrm{e}^{\frac{1}{2}\sqrt{\frac{2}{3}}\varphi}d\widetilde{t},$ (245)

the solution of which is an incomplete gamma function. However, for the Hubble
rate as given in Eq. (242), using the corresponding equations of motion, the
Einstein frame scalar field can be found to be

$\displaystyle\varphi=\frac{2\sqrt{2h_{s}\beta}\left(~{}\widetilde{t}-\widetilde{t}_{s}~{}\right)^{\frac{1-\beta}{2}}}{1-\beta}\,,$
(246)

where $h_{s}\beta>0$ is an essential criterion in order for the scalar field
to be canonical. We further notice that the transformation blows up at
$\widetilde{t}_{s}$ if $\beta>1$, hence, extra caution is required in this
case. in effect, the new Hubble rate in terms of our initial time coordinate
$\widetilde{t}$, i.e., $H(\tilde{t})\equiv\frac{1}{a}\frac{da}{d\tilde{t}}$,
is equal to,

$\displaystyle
H(\widetilde{t})=\frac{\sqrt{h_{s}\beta}\;(~{}\widetilde{t}-\widetilde{t}_{s}~{})^{-\frac{\beta+1}{2}}}{\sqrt{3}}+h_{s}(~{}\widetilde{t}-\widetilde{t}_{s}~{})^{-\beta}\,,$
(247)

and the question is what is the effect on the time coordinate $t$. We suppose
$t=f(~{}\widetilde{t}~{})$. Now, when, $f(~{}\widetilde{t}~{})\neq 0$, we
have,

$\displaystyle
f^{\prime}(~{}\widetilde{t}~{})=\mathrm{e}^{\frac{1}{2}\sqrt{\frac{2}{3}}\varphi}\,,$
(248)

and consequently the following relation

$\displaystyle\frac{dH}{dt}=\frac{d\widetilde{t}}{dt}\frac{d\widetilde{H}}{d\widetilde{t}}=\mathrm{e}^{-\frac{1}{2}\sqrt{\frac{2}{3}}\varphi}\frac{d\widetilde{H}}{d\widetilde{t}}\,.$
(249)

holds, that means, the expression $dH/dt$ diverges if and only if
$d\widetilde{H}/d\widetilde{t}$ diverges for $\beta<1$, as the conformal
factor is finite at $\widetilde{t}_{s}$. Now the Big Rip (Type I) singularity
occurs when the Hubble rate $H(t)$ diverges at a finite-time thus, we need to
examine whether $f(~{}\widetilde{t}_{s}~{})$ is finite or not. One can easily
derive $t_{s}$ given by Bahamonde:2016wmz

$\displaystyle
t_{s}=f(~{}\widetilde{t}_{s}~{})=c_{1}-c_{2}\Gamma\left(\frac{2}{\alpha+1}\right)\,,$
(250)

which is finite, provided $2/(1-\beta)$ is not a negative integer. This
readily implies that the singularity appears in the Hubble rate $H$ at a
finite-time, as long as $\beta\neq-(2/n)+1$ where $n\geq 2$ is an integer.
Accordingly, a Type II singularity is found if $dH/dt$ diverges, but $H$ does
not diverge. Combining together, these imply that $-3<\beta<-1$. On the other
hand, for $\beta<-3$, a Type IV singularity occurs. Now we investigate the
evolution of the scale factor when it is conformally transformed. The
evolution of the scale factor in the Einstein frame reads, as in (63),

$\displaystyle\widetilde{a}\left(~{}\widetilde{t}~{}\right)=\widetilde{a}_{s}\exp\left(\frac{h_{s}\left(~{}\widetilde{t}-\widetilde{t}_{s}~{}\right)^{1-\beta}}{1-\beta}\right)\,,$
(251)

which in the Jordan frame reads,

$\displaystyle
a\left(~{}\widetilde{t}~{}\right)=a_{s}\exp\left(\frac{3h_{s}\left(~{}\widetilde{t}-\widetilde{t}_{s}~{}\right)^{1-\beta}\pm
2\sqrt{3\beta
h_{s}}\left(~{}\widetilde{t}-\widetilde{t}_{s}~{}\right)^{\frac{1-\beta}{2}}}{3(1-\beta)}\right)\,.$
(252)

Now, depending on the values of $\beta$, from the scale factor (252), one can
witness the following type of singularities:

  1. 1.

For $\beta>1$, a Type I or no singularity occurs.

  2. 2.

For $-1<\beta<1$, a Type III singularity occurs.

  3. 3.

For $-3<\beta<-1$, a Type II singularity occurs.

  4. 4.

For $\beta<-3$, a Type IV singularity occurs.

In TABLE 1, we show the correspondence of the finite-time singularities
between the Jordan frame and the Einstein frame. TABLE 1 clearly shows that
the singularity in one frame may not be same in the other frame. For example,
as displayed in TABLE 1, the Type I (Big Rip) singularity in the Einstein
frame may correspond to a non-singular evolution in the Jordan frame. Further,
the Type II singularity in the Einstein frame could be modified to a more
severe Type III singularity in the Jordan frame, and the Type IV singularity
in the Einstein frame may correspond to a Type II singularity in the Jordan
frame.

Table 1: The table shows the correspondence of finite-time singularities in the Einstein and Jordan frames for the cosmological evolution given in terms of the Hubble rate (242) in the Einstein frame.  Singularity in the Einstein Frame | Singularity in the Jordan Frame
---|---
Type I | Type I or No Singularity
Type II | Type III
Type III | Type III
Type IV | Type IV or Type II

We close this section with a special case of the singular evolution (242)
describing the singular bounce cosmology, a special case of the symmetric
bounce Nojiri:2015sfd ; Nojiri:2016ppu ; Nojiri:2016ygo . One can rewrite the
scale factor of the cosmological evolution (242) as follows,

$\displaystyle a(t)=a_{s}\;\exp\left(h_{s}(t-t_{s})^{2(1+\epsilon)}\right)\,,$
(253)

for which the Hubble rate becomes

$\displaystyle H(t)=2(1+\epsilon)h_{s}(t-t_{s})^{2\epsilon+1},$ (254)

where $\epsilon>0$ and it has been chosen in such a way so that all the
quantities remain real. Now, in order to realize a bouncing scenario, for
$t<t_{s}$ the Hubble rate must become negative (i.e., $H<0$) and additionally,
in order for the bounce (253) to be a deformation of the symmetric bounce
described by $a(t)\sim\exp(\beta t^{2})$, the parameter $\epsilon$ must lie in
the interval $0<\epsilon<1$ having the following form

$\displaystyle\epsilon=\frac{2n}{2m+1}\,,$ (255)

where $m$ and $n$ are integers and they are chosen in such a way so that
$\epsilon<1$ is satisfied. Now, for this choice of $\epsilon$, the cosmology
described by the scale factor (253) and the Hubble rate (254) clearly depicts
a Type IV singular cosmology, in which case, the Hubble rate and its first
derivative with respect to the cosmic time, i.e., $\dot{H}$ remain finite, but
its second derivative with respect to the cosmic time, i.e., $\ddot{H}$
diverges. As demonstrated in Refs. Odintsov:2015jca ; Odintsov:2015gba  with
the use of reconstruction techniques, the pure $F(R)$ gravity realizing the
cosmological evolution (254) can approximately be given by

$\displaystyle F(R)=R+\frac{R^{2}}{4C}+\Lambda\,,$ (256)

where $C$ is a positive real number, near the bouncing point, which is
$t\simeq t_{s}$. Let us define a new parameter $x=t-t_{s}$ for simplicity.
Thus, the limit near the bouncing point, i.e., $t\simeq t_{s}$, corresponds to
the limit $x\rightarrow 0$. Now, in order to transform the theory in the
Einstein frame, we consider the following conformal transformation

$\displaystyle g_{\mu\nu}=\mathrm{e}^{-\sigma}\hat{g}_{\mu\nu}\,,$ (257)

where the scalar field $\sigma$ is equal to $\sigma=-\ln F^{\prime}(A)$. In
terms of the new parameter $x$, the Ricci scalar reads

$\displaystyle
R=12h_{s}(\epsilon+1)x^{2\epsilon}\Bigg{[}4h_{s}(\epsilon+1)x^{2\epsilon+2}+2\epsilon+1\Bigg{]},$
(258)

and if we are close to the singularity where $x$ is small, then the Ricci
scalar becomes,

$\displaystyle R\approx 12h_{s}(\epsilon+1)(2\epsilon+1)x^{2\epsilon}\,.$
(259)

Consequently, by combining Eqs. (256) and (259), one obtains

$\displaystyle F^{\prime}(R)\approx
1+\frac{6h_{s}(\epsilon+1)(2\epsilon+1)x^{2\epsilon}}{C}\equiv
1+Dx^{2\epsilon}\,,$ (260)

where $D$ is a positive real number. This means that the new time coordinate
of the Einstein frame FLRW metric will be given in terms of $x$ as follows

$\displaystyle
d\widetilde{t}=(1+Dx^{2\epsilon})^{\frac{1}{2}\sqrt{\frac{3}{2}}}dx\,,$ (261)

the solution of which is a hypergeometric function. The new scale factor in
terms of $x$ reads

$\displaystyle
a(~{}\widetilde{t}~{})=(1+Dx^{2\epsilon})^{\frac{1}{2}\sqrt{\frac{3}{2}}}\;\widetilde{a}(x)\,,$
(262)

and subsequently, the derivative of the scale factor is given by

$\displaystyle\frac{da}{d\widetilde{t}}=\frac{dx}{d\widetilde{t}}\frac{da}{dx}=\frac{dx}{d\widetilde{t}}\Bigg{[}(1+Dx^{2\epsilon})^{\frac{1}{2}\sqrt{\frac{3}{2}}}\frac{d\widetilde{a}}{dx}+2\epsilon
x^{2\epsilon}(1+Dx^{2\epsilon})^{\frac{1}{2}\sqrt{\frac{3}{2}}-1}\widetilde{a}(x)\Bigg{]}\,.$
(263)

As $dt=dx$ at $x\simeq 0$, thus, looking at the second derivative of the scale
factor $a(t)$, one may conclude that the scale factor diverges at $x=0$
provided $\epsilon<1/2$. Hence, the Type IV singularity in the Jordan frame
becomes a Type II singularity in the Einstein frame, also reflected from TABLE
1.

####  IV.5.2 Occurrence of Singularities in Different Frames: Unimodular
$F(R)$ Gravity

In this section, we discuss the correspondence between the frames in the
context of unimodular $F(R)$ gravity Nojiri:2015sfd ; Nojiri:2016ppu ;
Nojiri:2016ygo ; Odintsov:2016imq . The Jordan frame unimodular $F(R)$ gravity
is described by the action Nojiri:2015sfd

$\displaystyle S=\int
d^{4}x\bigg{[}\sqrt{-g}\;\bigg{(}F(R)-\lambda\bigg{)}+\lambda\bigg{]}+S_{\mathrm{matter}}\,,$
(264)

where we assume that $F(R)$ is a smooth function of the Ricci scalar $R$,
$\lambda$ denotes the Lagrange multiplier function and $S_{\mathrm{matter}}$
stands for the action of the matter fluids present in the Universe sector.
Notice that the variation of the action (264) with respect to $\lambda$ leads
to the unimodular constraint

$\displaystyle\sqrt{-g}$ $\displaystyle=1\,.$ (265)

The unimodular constraint given in Eq. (265) is the key point of unimodular
$F(R)$ gravity. Now if we vary the unimodular $F(R)$ gravity action (264) with
respect to the metric tensor $g_{\mu\nu}$, we obtain the unimodular $F(R)$
gravity field equations as follows

$\displaystyle
0=\frac{1}{2}g_{\mu\nu}\left(F(R)-\lambda\right)-R_{\mu\nu}F^{\prime}(R)+\nabla_{\mu}\nabla_{\nu}F^{\prime}(R)-g_{\mu\nu}\nabla^{2}F^{\prime}(R)+\frac{1}{2}T_{\mu\nu}\,.$
(266)

Now, in order to proceed with the cosmological evolution, one needs to be very
careful because the flat standard FLRW metric (3) does not satisfy the
unimodular constraint given in Eq. (265). Nevertheless, taking the following
coordinate transformation

$\displaystyle d\tau=a^{3}(t)dt\,,$ (267)

one can verify that the resulting metric satisfies the unimodular condition
(265). Now, with the use of the transformation (267), the FLRW metric in (3)
can be written as

$\displaystyle
ds^{2}=-a^{-6}\left(t\left(\tau\right)\right)d\tau^{2}+a^{2}\left(t\left(\tau\right)\right)\left(dx^{2}+dy^{2}+dz^{2}\right),$
(268)

which for the sake of brevity we call the unimodular FLRW metric. Using this
unimodular FLRW metric, the vacuum field equations become Nojiri:2015sfd ;
Nojiri:2016ppu ; Nojiri:2016ygo ; Odintsov:2016imq ,

$\displaystyle 0=$
$\displaystyle\,-\frac{a^{-6}}{2}\left(F(R)-\lambda\right)+\left(3\dot{K}+12K^{2}\right)F^{\prime}(R)-3K\frac{dF^{\prime}(R)}{d\tau}\,,$
(269) $\displaystyle 0=$
$\displaystyle\,\frac{a^{-6}}{2}\left(F(R)-\lambda\right)-\left(\dot{K}+6K^{2}\right)F^{\prime}(R)+5K\frac{dF^{\prime}(R)}{d\tau}+\frac{d^{2}F^{\prime}(R)}{d\tau^{2}}\,,$
(270)

where the function $K(\tau)$ is defined as the corresponding Hubble rate in
the “$\tau$” coordinate, that is,

$\displaystyle K=\frac{1}{a(\tau)}\frac{da(\tau)}{d\tau}\,.$ (271)

By using the unimodular FLRW metric of Eq. (268), the corresponding Ricci
scalar now becomes

$\displaystyle R$ $\displaystyle=a^{6}\left(6\dot{K}+30K^{2}\right)\,.$ (272)

Now, in the following sections we describe the correspondence of the Jordan
frame unimodular $F(R)$ gravity in the Einstein frame.

We begin with the unimodular $F(R)$ gravity action in the Jordan frame, i.e.,
Eq. (264) without the matter sector, and following the same approach as in the
ordinary $F(R)$ gravity, that means, introducing the auxiliary field $A$, we
rewrite the action (264) as follows

$\displaystyle S=\int
d^{4}x\Bigg{[}\sqrt{-g}\Bigg{(}F^{\prime}(A)(R-A)+F(A)-\lambda\Bigg{)}+\lambda\Bigg{]}\,.$
(273)

We note that the last term of the action (273) will remain unaffected by the
conformal transformations. Now to obtain a minimally coupled scalar-tensor
theory, we perform the conformal transformation as in the standard $F(R)$
gravity case, i.e., $\hat{g}_{\mu\nu}=\mathrm{e}^{\sigma}g_{\mu\nu}$, where
$\hat{g}$ denotes the metric in the Einstein frame, and $g$ is the metric in
the Jordan frame. Further, the scalar field $\sigma$ in terms of the auxiliary
field $A$ is given by, $\sigma=-\ln F^{\prime}(A)$, hence, the action (273)
becomes,

$\displaystyle S=\int
d^{4}x\left\\{\sqrt{-\hat{g}}\left(\hat{R}-\frac{3}{2}\hat{g}^{\mu\nu}\partial_{\mu}\sigma\partial_{\nu}\sigma-V(\sigma)-\lambda\mathrm{e}^{-2\sigma}\right)+\lambda\right\\}\,,$
(274)

which describes a canonical scalar field action, in the absence of any matter
sector. We note that the unimodular constraint is not unaffected by the
conformal transformation, so in the case it becomes

$\displaystyle\sqrt{-\hat{g}}=\mathrm{e}^{2\sigma}\,.$ (275)

In effect, the FLRW metric does not satisfy this constraint identically, so in
order to overcome this issue, similar to the earlier approach, we introduce a
new time coordinate $\tilde{\tau}$, which is related to the cosmic time
$\tilde{t}$ as follows

$\displaystyle
d\widetilde{\tau}=\widetilde{a}^{3}\left(~{}\widetilde{t}~{}\right)\mathrm{e}^{2\sigma\left(~{}\widetilde{t}~{}\right)}d\widetilde{t}\,,$
(276)

so that the conformally transformed unimodular constraint of Eq. (275) is
satisfied and the corresponding Einstein frame unimodular FLRW metric takes
the form

$\displaystyle
ds^{2}=-\frac{\mathrm{e}^{4\sigma\left(\widetilde{\tau}\right)}}{\widetilde{a}^{6}\left(~{}\widetilde{\tau}~{}\right)}d\widetilde{\tau}^{2}+\widetilde{a}^{2}\left(\widetilde{\tau}\right)\sum_{i=1}^{3}dx_{i}^{2}\,.$
(277)

In summary, in order to transform from one frame to other at the level of
metric, we have the following transformations. The scale factor transforms as

$\displaystyle\widetilde{a}\left(~{}\widetilde{\tau}~{}\right)$
$\displaystyle=\mathrm{e}^{\sigma/2}a(\tau)\,,$ (278)

where the parameter $\widetilde{\tau}$ is related to $\tau$ as follows,

$\displaystyle\left(\frac{\mathrm{e}^{4\sigma\left(\widetilde{\tau}\right)}}{\widetilde{a}^{6}\left(~{}\widetilde{\tau}~{}\right)}\right)d\widetilde{\tau}^{2}=\left(\frac{\mathrm{e}^{\sigma}}{a^{6}(\tau)}\right)d\tau^{2}\,.$
(279)

Combining Eqs. (278) and (279), one can easily observe that the new time
coordinate $\widetilde{\tau}$ is the same as the coordinate $\tau$, i.e.,
$\widetilde{\tau}=\tau$.

So far we discuss the unimodular $F(R)$ gravity scenario, however, we also
need to discuss the scalar-tensor unimodular gravity. To proceed, we begin
with the following minimally coupled scalar-tensor action

$\displaystyle S=\int
d^{4}x\left\\{\sqrt{-\hat{g}}\left(\frac{\hat{R}}{2\kappa^{2}}-\frac{1}{2}\hat{g}^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi-V(\varphi)-\lambda
h(\varphi)\right)+\lambda\right\\}\,,$ (280)

where $\lambda$ is a constant and in addition, for the time being we have
assumed that the determinant of the metric is given by an arbitrary function
of the scalar field but this will be determined later by the requirement that
the action has a Jordan frame. Now, considering the flat unimodular FLRW
metric, for the above action (280), one can write down the gravitational
equations which are nothing but the standard scalar field cosmological
equations with a modified scalar potential having the form $V(\varphi)+\lambda
h(\varphi)$. In particular, the gravitational equations read

$\displaystyle 3\widetilde{K}^{2}=$
$\displaystyle\,\frac{1}{2}\dot{\varphi}^{2}+\Big{(}V(\varphi)+\lambda
h(\varphi)\Big{)}\widetilde{a}^{-6}(~{}\widetilde{\tau}~{})\,,$ (281)
$\displaystyle 9\widetilde{K}^{2}+2\dot{\widetilde{K}}=$
$\displaystyle\,-\frac{1}{2}\dot{\varphi}^{2}+\Big{(}V(\varphi)+\lambda
h(\varphi)\Big{)}\widetilde{a}^{-6}(\widetilde{\tau})\,,$ (282)

where $\widetilde{K}(\widetilde{\tau})$ denotes the unimodular Hubble
parameter in the Einstein frame, explicitly given by

$\displaystyle\widetilde{K}\left(\widetilde{\tau}\right)$
$\displaystyle=\widetilde{K}(\tau)=\frac{1}{\widetilde{a}(\tau)}\frac{d\widetilde{a}(\tau)}{d\tau}\,.$
(283)

Now in order to conformally transform this action (280) to a Jordan frame, we
apply the following conformal transformation

$\displaystyle
g_{\mu\nu}=\mathrm{e}^{\pm\kappa\sqrt{\frac{2}{3}}\varphi}\hat{g}_{\mu\nu}\,.$
(284)

This rescaling eliminates the kinetic term from the action (280) and recast
into the following form

$\displaystyle S=\int
d^{4}x\Bigg{[}\sqrt{-g}\bigg{\\{}\mathrm{e}^{\pm\kappa\sqrt{\frac{2}{3}}\varphi}{2\kappa^{2}}R-\mathrm{e}^{\pm
2\kappa\sqrt{\frac{2}{3}}\varphi}\bigg{(}V(\varphi)+\lambda
h(\varphi)\bigg{)}\bigg{\\}}+\lambda\Bigg{]}\,.$ (285)

Now, varying the action (285) with respect to the scalar field $\phi$, which
is now just an auxiliary field, we get

$\displaystyle
R=\mathrm{e}^{\pm\kappa\sqrt{\frac{2}{3}}\varphi}\Bigg{[}4\kappa^{2}\bigg{(}V(\varphi)+\lambda
h(\varphi)\bigg{)}\pm
2\kappa\sqrt{\frac{2}{3}}\bigg{(}V^{\prime}(\varphi)+\lambda
h^{\prime}(\varphi)\bigg{)}\Bigg{]}\,.$ (286)

In order to find the Jordan frame unimodular $F(R)$ gravity, one needs to
invert Eq. (286) to get the scalar field $\phi(R)$, as a function of $R$ only,
and so independent of $\lambda$. This directs us to choose

$\displaystyle h(\varphi)=\mathrm{e}^{-2\kappa\sqrt{\frac{3}{2}}\varphi}\,.$
(287)

Hence, in this case we can invert to find $\varphi=\varphi(R)$, and
consequently, the resulting unimodular $F(R)$ gravity theory takes the form

$\displaystyle S=\int
d^{4}x\Bigg{[}\sqrt{-g}\bigg{(}F(R)-\lambda\bigg{)}+\lambda\Bigg{]},$ (288)

where the expression of $F(R)$ is given by

$\displaystyle
F(R)=\mathrm{e}^{\pm\kappa\sqrt{\frac{2}{3}}\varphi(R)}{2\kappa^{2}}R-\mathrm{e}^{\pm
2\kappa\sqrt{\frac{2}{3}}\varphi(R)}V(\varphi(R))\,.$ (289)

Now, having the correspondence between the Jordan and Einstein frames, in the
following we discuss the correspondence of the finite-time singularities in
the two frames. We again note that for all the physical quantities in the
Einstein frame we use tilde to differ it from the Jordan frame, and
additionally, we set the gravitational coupling constant to be unity
throughout the discussions, i.e., $\kappa=1$.

(i) Power law evolution

We consider the power law cosmology characterized by the following scale
factor expressed in terms of $\widetilde{\tau}$ (the time coordinate of the
unimodular Einstein frame FLRW metric (277)) as follows

$\displaystyle\widetilde{a}\left(\widetilde{\tau}\right)=\widetilde{a}_{c}\left(\frac{\widetilde{\tau}}{\widetilde{\tau}_{c}}\right)^{p}\,,$
(290)

where $\widetilde{\tau}_{c}$ is some fiducial time and $p$ is a positive real
number. For this scale factor, the unimodular Hubble parameter becomes
$\widetilde{K}(\tau)=p\widetilde{\tau}^{-1}$. Thus, one can see that at
$\widetilde{\tau}=0$, the unimodular Hubble rate diverges, and hence, for the
power law cosmology in the Einstein frame, a Big Bang singularity is realized.
The above scale factor is a solution of the Einstein frame Friedmann and
Raychaudhuri equations (281)-(282) when the potential has an exponential form.
In this case the scalar field takes the form

$\displaystyle\varphi\left(\widetilde{\tau}\right)=\pm\sqrt{2p(1-3p)}\,\log\left(\widetilde{\tau}/\widetilde{\tau}_{c}\right)\,.$
(291)

Now, using this solution (291), it is possible to change all the variables
from the Einstein frame to the Jordan frame with the use of the following
conformal transformation

$\displaystyle
g_{\mu\nu}=\mathrm{e}^{\pm\sqrt{\frac{2}{3}}\varphi}\hat{g}_{\mu\nu}\,.$ (292)

As discussed above, the time coordinate $\tau$ in the Jordan frame is
equivalent to the original time coordinate in the Einstein frame. Thus, using
the transformation (292), one can express scale factor in the Jordan frame as
follows

$\displaystyle
a(\tau)=\mathrm{e}^{\pm\frac{1}{2}\sqrt{\frac{2}{3}}\varphi}\widetilde{a}(\widetilde{\tau})\sim\tau^{p\pm\sqrt{\frac{p(1-3p)}{3}}}\,.$
(293)

Now, when the minus sign is taken in the conformal transformation, we see that
for $0<p<1/6$ the quantity $p-\sqrt{\frac{p(1-3p)}{3}}$ is negative, so at
$\tau=0$ the scale factor diverges, this means that in the Jordan frame we
have a contracting Universe. For $p=1/6$ the Universe becomes static. When
$1/6<p\leq 1/3$, in the Jordan frame, the Universe conserves the Big Bang
singularity at $\tau=0$, and finally, for $p>1/3$ the exponent in (293)
becomes complex, meaning that this case does not make physical sense.

(ii) The case $F(R)=R^{-n}$ in the Jordan Frame

Here we investigate a model starting in the Jordan frame and conformally
transform it to the Einstein frame. Let us consider the vacuum unimodular
$F(R)$ gravity having a typical form $F(R)\sim R^{-n}$. In this case, the
scale factor evolves as

$\displaystyle
a(\tau)\sim(\tau_{s}-\tau)^{\frac{1+3n+2n^{2}}{5+10n+6n^{2}}}\,,$ (294)

where we have used the $\tau$ coordinate of the unimodular FLRW metric. One
can see that a Big Rip singularity occurs in the $\tau$ coordinate, if the
parameter $n$ lies in the range $-1<n<-1/2$. Now, in terms of the original
FLRW cosmic time $t$, the scale factor takes the form

$\displaystyle a(t)\sim t^{\frac{(2n+1)(n+1)}{2+n}}\;,$ (295)

which shows that a Big Rip singularity in the $\tau$ coordinate, appears in
the $t$ coordinate if $n<-2$, or, $-1<n<-1/2$.

We shall now investigate what happens in the corresponding Einstein frame
scalar-tensor theory. The scalar field $\sigma$ of the conformal factor is
given by

$\displaystyle\sigma\sim(n+1)\ln
R\sim-\left(\frac{2(n+1)(2+n)}{5+10n+6n^{2}}\right)\ln(\tau_{s}-\tau)\,,$
(296)

where the Ricci scalar in the unimodular time parameter is given by

$\displaystyle
R\sim(\tau_{s}-\tau)^{-2\left(\frac{2+n}{5+10n+6n^{2}}\right)}\,.$ (297)

This means that under a conformal transformation, the scale factor transforms
as

$\displaystyle\widetilde{a}\left(\widetilde{\tau}\right)$
$\displaystyle=\mathrm{e}^{\sigma/2}a(\tau),$ (298)

where the parameter $\widetilde{\tau}$ is the same as the original time
coordinate $\tau$, that means, $\widetilde{\tau}=\tau$. Therefore, we get

$\displaystyle\widetilde{a}\sim(\tau_{s}-\tau)^{\frac{n^{2}-1}{6n^{2}+10n+5}}\,.$
(299)

These conditions show that we could have more values of $n$ for which the
power of the scale factor is negative. As one can see that for $-1<n<1$, the
power of the scale factor becomes negative. Thus, we conclude that a Type I
singularity will appear for more values of $n$.

(iii) A singular cosmological model

Here we consider a toy model where the unimodular Hubble parameter in the
Einstein frame is given by

$\displaystyle\widetilde{K}(~{}\widetilde{\tau}~{})$
$\displaystyle=f_{s}(\widetilde{\tau}-\widetilde{\tau}_{s})^{\alpha}\,,$ (300)

where $k_{s}$, $\alpha$ are real numbers and $\widetilde{\tau}_{s}$ refers to
a specific time instant of $\widetilde{\tau}$. For this Hubble parameter, the
scale factor can be given by

$\displaystyle\widetilde{a}\left(\widetilde{\tau}\right)$
$\displaystyle=a_{s}\exp\left[\frac{f_{s}\left(\widetilde{\tau}-\widetilde{\tau}_{s}\right)^{\alpha+1}}{\alpha+1}\right]\,.$
(301)

where $a_{s}$ is a constant.

Now, to convert this Einstein frame solution to the Jordan frame, one needs to
find the scalar field $\varphi$ giving rise to such a solution (300). To do
this it is essential to solve the differential equation obtained by
subtracting the unimodular Friedmann equations (281) and (282) as follows

$\displaystyle
0=\dot{\varphi}^{2}+6f_{s}^{2}\left(\widetilde{\tau}-\widetilde{\tau}_{s}\right)^{2\alpha}+2f_{s}\alpha\left(\widetilde{\tau}-\widetilde{\tau}_{s}\right)^{\alpha-1}\,.$
(302)

However, this equation cannot be solved in general for an arbitrary power of
$\alpha$, hence, one needs to find an approximate solution to the above
differential equation. Thus, in order to proceed, one needs to approximate the
solution around the singularity $\widetilde{\tau}=\widetilde{\tau}_{s}$. This
can be done for two cases as follows: when $\alpha<-1$ and the Hubble rate is
that of a Type I singularity, and when $\alpha>-1$ and the other types of
singularities are present. For $\alpha<-1$, and around the singularity, the
second term of (302) dominates over the third term and the differential
equation (302) is approximated to be

$\displaystyle\dot{\varphi}^{2}+6f_{s}^{2}\left(\widetilde{\tau}-\widetilde{\tau}_{s}\right)^{2\alpha}$
$\displaystyle\sim 0\,.$ (303)

In this case, the solution for $\varphi$ becomes imaginary, and hence, such a
Hubble rate could only be described by a phantom scalar field. For such a
phantom field, the corresponding Jordan frame $F(R)$ becomes complex and this
is unphysical.

However, for $\alpha>-1$, and around the singularity, the third term of (302)
dominates over the second term and therefore the scalar field behaves as

$\displaystyle\varphi\left(\widetilde{\tau}\right)$
$\displaystyle\sim\pm\frac{2\sqrt{-2f_{s}\alpha}}{\alpha+1}\left(\widetilde{\tau}-\widetilde{\tau}_{s}\right)^{\frac{\alpha+1}{2}}\,,$
(304)

which is real if $-2f_{s}\alpha\geq 0$ and nontrivial if $-2f_{s}\alpha\neq
0$. Thus, in this case one can continue further and conformally transform to
the Jordan frame. Applying the conformal transformation (284) and using the
fact that the time coordinate is unchanged, i.e., $\tau=\widetilde{\tau}$, the
scale factor in the Jordan frame reads

$\displaystyle a(\tau)\sim
a_{s}\exp\left(\frac{3f_{s}(\tau-\tau_{s})^{\alpha+1}\pm 2\sqrt{-3\alpha
f_{s}}\;(\tau-\tau_{s})^{\frac{\alpha+1}{2}}}{3(\alpha+1)}\right)\,.$ (305)

Now from this scale factor (305), one can describe a variety of finite-time
singularities as follows:

  * •

For $-1<\alpha<1$, a Type III singularity occurs.

  * •

For $1<\alpha<3$, a Type II singularity occurs.

  * •

For $3<\alpha$, a Type IV singularity occurs.

In TABLE 2 we show how the singularities change their types from one frame to
another444We note that in TABLE 2, we do not keep Type I singularity since it
only appears when $\alpha<-1$, and hence, we have a phantom scalar field and
consequently for such a phantom scalar field, the corresponding Jordan frame
$F(R)$ becomes complex and this is unphysical.. One can observe that the
unimodular $F(R)$ case behaves similarly to the standard $F(R)$ case. From
TABLE 2, we observe that a Type II singularity in the Einstein frame is
modified to the more severe Type III singularity in the Jordan frame. The Type
IV singularity in the Einstein frame may become a more severe Type II
singularity in the Jordan frame if the parameter $\alpha$ lies in the range
$1<\alpha<3$.

Table 2: We show the correspondence of finite-time singularities in the Einstein and Jordan frames, for the cosmological evolution (300) in the Einstein frame for the case $\alpha>-1$.  Singularity in the Einstein Frame | Singularity in the Jordan Frame
---|---
Type II | Type III
Type III | Type III
Type IV | Type IV or Type II

2.1. Qualitative Analysis of the Phase structure of Unimodular $F(R)$ gravity
near Finite-Time Singularities

Here we discuss the qualitative behavior of the dynamical system corresponding
to the vacuum unimodular $F(R)$ gravity near the finite-time singularities. To
begin with, we introduce the following variables Bahamonde:2016wmz

$\displaystyle
x_{1}=-\frac{1}{KF^{\prime}(R)}\frac{dF^{\prime}(R)}{d\tau}\,,\quad
x_{2}=-\frac{F(R)}{6K^{2}F^{\prime}(R)a^{6}}\,,\quad
x_{3}=\frac{R}{6K^{2}a^{6}}\,,\quad
x_{4}=\frac{\lambda}{6a^{6}K^{2}F^{\prime}(R)}\,,$ (306)

which allows one to write the Friedmann equation as

$\displaystyle x_{1}+x_{2}+x_{3}+x_{4}=1\,.$ (307)

This constraint tells us that in order to describe the dynamics one needs only
three variables as the fourth variable can be obtained using the remaining
three variables. Using the above variables, the Raychaudhuri equation becomes

$\displaystyle\frac{1}{F^{\prime}(R)K^{2}}\frac{d^{2}F^{\prime}(R)}{d\tau^{2}}=1+5x_{1}+3x_{2}+x_{3}+3x_{4}\,.$
(308)

Now, differentiating the variables as defined in Eq. (306) with respect to a
cosmological time $dN=K(\tau)d\tau$, one gets the following dynamical system
Bahamonde:2016wmz

$\displaystyle x_{1}^{\prime}=$
$\displaystyle\,-1+x_{1}^{2}-x_{1}x_{3}-3x_{2}-x_{3}-3x_{4}\,,$ $\displaystyle
x_{2}^{\prime}=$
$\displaystyle\,-m+4x_{2}+x_{1}x_{2}-2x_{2}x_{3}+50-16x_{3}\,,$ $\displaystyle
x_{3}^{\prime}=$ $\displaystyle\,m+20x_{3}-2x_{3}^{2}-50\,,$ $\displaystyle
x_{4}^{\prime}=$ $\displaystyle\,x_{4}(x_{1}-2x_{3}+4)\,,$ (309)

where $m=\ddot{K}/K^{3}$. Let us note that from the constraint Eq. (307), in
reality, one only has three dynamical equations in (IV.5.2). In fact, one can
disregard the fourth equation which is a combination of the other three. The
dynamical system, as one can notice, is non-autonomous due to the existence of
the term $m$. The critical points of the system (IV.5.2) have been shown in
TABLE 3. For an autonomous system $m$ should vanish. Now if we consider that
$K(\tau)=f_{s}(\tau-\tau_{s})^{\alpha}$ with $\alpha<-1$ which corresponds to
the Big Rip singularity as $\tau\to\tau_{s}$, then the parameter $m$ is equal
to, $m=\frac{(-1+\alpha)\alpha(\tau-\tau_{s})^{-2-2\alpha}}{f_{s}^{2}}\sim 0$,
and consequently, the dynamical system of Eq. (IV.5.2) becomes autonomous near
the Big Rip singularity. Thus, one may expect that finding the fixed points of
the dynamical system near the Big Rip singularity could offer some insights on
the non-autonomous system, however, from TABLE 3, one can clearly see that if
one sets $m=0$, then $x_{1}$ and $x_{2}$ become complex, that means the
critical points are not real. In fact, if we assume that $m$ is a constant,
then in order to have real critical points one needs $m\gtrsim 17$. This
indicates that the dynamics near the Big Rip singularity behaves very
strangely.

Table 3: The table shows the critical points of the dynamical system (IV.5.2).  Critical Point | $x_{1}$ | $x_{2}$ | $x_{3}$ | $x_{4}$
---|---|---|---|---
$P_{1}$ | $\frac{1}{2\sqrt{2}}\left(-\sqrt{m}-\sqrt{m+4\sqrt{2}\sqrt{m}-40}+2\sqrt{2}\right)$ | $\frac{1}{4}\left(3\sqrt{2}\sqrt{m}+\sqrt{2}\sqrt{m+4\sqrt{2}\sqrt{m}-40}-20\right)$ | $5-\frac{\sqrt{m}}{\sqrt{2}}$ | $0$
$P_{2}$ | $\frac{1}{2\sqrt{2}}\left(\sqrt{m}-\sqrt{m+4\sqrt{2}\sqrt{m}-40}+2\sqrt{2}\right)$ | $\frac{1}{4}\left(3\sqrt{2}\sqrt{m}-\sqrt{2}\sqrt{m+4\sqrt{2}\sqrt{m}-40}-20\right)$ | $5-\frac{\sqrt{m}}{\sqrt{2}}$ | $0$
$P_{3}$ | $\frac{1}{2\sqrt{2}}\left(\sqrt{m}+\sqrt{m-4\sqrt{2}\sqrt{m}-40}+2\sqrt{2}\right)$ | $\frac{1}{4}\left(-3\sqrt{2}\sqrt{m}-\sqrt{2}\sqrt{m-4\sqrt{2}\sqrt{m}-40}-20\right)$ | $5+\frac{\sqrt{m}}{\sqrt{2}}$ | $0$
$P_{4}$ | $\frac{1}{2\sqrt{2}}\left(\sqrt{m}-\sqrt{m-4\sqrt{2}\sqrt{m}-40}+2\sqrt{2}\right)$ | $\frac{1}{4}\left(-3\sqrt{2}\sqrt{m}+\sqrt{2}\sqrt{m-4\sqrt{2}\sqrt{m}-40}-20\right)$ | $5+\frac{\sqrt{m}}{\sqrt{2}}$ | $0$

###  IV.6 $F(G)$ Gravity

Apart from the popular modified gravity models, namely, $F(R)$ and $F(T)$,
another very interesting class of modified gravity models which can explain
the late-time cosmic acceleration is the string-inspired modified Gauss–Bonnet
gravity $-$ the so called $F(G)$-gravity Nojiri:2005jg ; Nojiri:2005am ;
Cognola:2006eg ; Cognola:2006sp ; Brevik:2006nh ; Li:2007jm ; Boehmer:2009fey
; Sadeghi:2009hi ; Elizalde:2010jx ; MohseniSadjadi:2010pc ; Garcia:2010xz ;
Banijamali:2012zz ; Rodrigues:2012qu ; Bamba:2014mya ; Momeni:2014bua ;
Huang:2015kca ; Abdolmaleki:2015iaa ; Astashenok:2015haa ; Kusakabe:2015ida ;
Oikonomou:2015qha ; Carloni:2017ucm ; Oikonomou:2017ppp ; Terrucha:2019jpm ;
Inagaki:2019rhm ; Bahamonde:2019swy  where $F(G)$ is any arbitrary function of
the Gauss-Bonnet invariant defined in (161). Even though this class of
modified gravity models can explain the late-time accelerating expansion of
the Universe, however, it has been found that they can lead to finite-time
future singularities Bamba:2008ut ; Bamba:2010wfw . In this section, we shall
describe the appearance of finite-time singularities in this class of modified
gravity models. Unfortunately this model includes ghosts DeFelice:2009ak
although the ghost-free modifications have been proposed Nojiri:2018ouv ;
Nojiri:2019dwl . In this section, however, we show the structure of the
singularities based on the original models for the purpose of the
illustration.

We start with the action of $F(G)$-gravity which is given by Nojiri:2005jg

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa^{2}}\bigg{(}R+F(G)\bigg{)}+{L}_{\mathrm{matter}}\right]\,.$
(310)

In the background of a spatially flat FLRW Universe in (3), one can write down
the equation of motion for this gravity as follows

$\displaystyle 24H^{3}\dot{F}^{\prime}(G)+6H^{2}+F(G)-GF^{\prime}(G)=$
$\displaystyle\,2\kappa^{2}\rho\,,$ (311) $\displaystyle
8H^{2}\ddot{F}^{\prime}(G)+16H\dot{F}^{\prime}(G)\left(\dot{H}+H^{2}\right)+\left(4\dot{H}+6H^{2}\right)+F(G)-GF^{\prime}(G)=$
$\displaystyle\,-2\kappa^{2}p\,,$ (312)

where the overhead dot stands for the derivative with respect to the cosmic
time and the prime denotes the differentiation with respect to $G$. The above
two equations (311) and (312) can also be alternatively represented as

$\displaystyle\rho_{\mathrm{eff}}=\frac{3}{\kappa^{2}}H^{2}\,,\quad
p_{\mathrm{eff}}=-\frac{1}{\kappa^{2}}\left(2\dot{H}+3H^{2}\right)\,,$ (313)

where $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are the effective energy
density and pressure of the Universe, respectively. We note that the energy
density and pressure of the matter sector described by ${L}_{\mathrm{matter}}$
in Eq. (310) are contained in $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$.
Now, using the expressions of $R$ and $G$ in the FLRW Universe, given by
$R=6\left(2H^{2}+\dot{H}\right)$ and $G=24H^{2}\left(H^{2}+\dot{H}\right)$,
one can visualize $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ of Eq. (313) as
follows:

$\displaystyle\rho_{\mathrm{eff}}=$
$\displaystyle\,\frac{1}{2\kappa^{2}}\left[-F(G)+24H^{2}\left(H^{2}+\dot{H}\right)F^{\prime}(G)-24^{2}H^{4}\left(2\dot{H}^{2}+H\ddot{H}+4H^{2}\dot{H}\right)F^{\prime\prime}(G)\right]+\rho\,,$
(314) $\displaystyle p_{\mathrm{eff}}=$
$\displaystyle\,\frac{1}{2\kappa^{2}}\Bigg{[}F(G)-24H^{2}\left(H^{2}+\dot{H}\right)F^{\prime}(G)+(8\times
24)\;H^{2}\Bigl{\\{}6\dot{H}^{3}+8H\dot{H}\ddot{H}+24\dot{H}^{2}H^{2}+6H^{3}\ddot{H}$
$\displaystyle\,+8H^{4}\dot{H}+H^{2}\dddot{H}\Bigr{\\}}F^{\prime\prime}(G)+(8\times
24^{2})\;H^{4}\left(2\dot{H}^{2}+H\ddot{H}+4H^{2}\dot{H}\right)^{2}F^{\prime\prime\prime}(G)\Bigg{]}+p\,,$
(315)

Assuming the matter sector with a constant EoS parameter $w\equiv p/\rho$, and
then combining the Friedmann and Raychaudhuri equations in Eq. (313), one
leads to the following equation:

$\displaystyle{G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)=-\frac{1}{\kappa^{2}}\left[2\dot{H}+3(1+w)H^{2}\right],$
(316)

where

$\displaystyle{G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)=p_{\mathrm{eff}}-w\rho_{\mathrm{eff}}.$
(317)

The Eq. (316) brings in a very crucial physical insight in terms of the new
function ${G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)$ which involves
$H$, $\dot{H}$ and the higher order derivatives of $H$. When a cosmological
model is prescribed in terms of the Hubble rate $H=H(t)$, the right-hand side
of Eq. (316) then reduces to a function, $f(t)$ of the cosmic time $t$. Now,
if the function ${G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)$ in Eq.
(317) is chosen in such a way so that
${G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)$ reproduces the above
function $f(t)$, then the aforementioned cosmology given by $H=H(t)$ can be
realized. Therefore, the function
${G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)$ plays a very crucial role
to judge the viability of a given cosmological scenario characterized by the
Hubble rate Bamba:2008ut .

The mathematical form of ${G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)$
can be determined from the prescribed gravitational theory. In the context of
$F(G)$-gravity, inserting Eqs. (314) and (315) into Eq. (317), one obtains,

$\displaystyle{G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)=\frac{1}{2\kappa^{2}}$
$\displaystyle\,\left[(1+w)F(G)-24(1+w)H^{2}\left(H^{2}+\dot{H}\right)F^{\prime}(G)+(8\times
24)H^{2}\left\\{6\dot{H}^{3}+8H\dot{H}\ddot{H}\right.\right.$
$\displaystyle\,\left.+6(4+w)\dot{H}^{2}H^{2}+3(2+w)H^{3}\ddot{H}+4(2+3w)H^{4}\dot{H}+H^{2}\dddot{H}\right\\}F^{\prime\prime}(G)$
$\displaystyle\,\left.+(8\times
24^{2})H^{4}\left(2\dot{H}^{2}+H\ddot{H}+4H^{2}\dot{H}\right)^{2}F^{\prime\prime\prime}(G)\right]\,.$
(318)

####  IV.6.1 Finite-time Singularities

In this section, we investigate the possibility of finite-time future
singularities in $F(G)$ gravity models. To start with we consider the
following expression of the Hubble parameter in (218). From Eq. (218) one can
clearly notice that if $\beta>0$, the $H(t)$ becomes singular in the limit
$t\rightarrow t_{s}$. Therefore, $t_{s}$ determines the time when a
singularity appears in the above model. While on the other hand, for
$\beta<0$, $H(t)$ does not exhibit any singular behavior, but for non-integer
negative values of $\beta$, some derivative of $H(t)$ may be singular and
consequently the curvature could exhibit singular behavior Bamba:2008ut .
Therefore, we see that for both positive and negative values of $\beta$,
singularities in $H(t)$ or in the curvature may appear. Therefore, we assume
$\beta\neq 0$ because the case $\beta=0$ corresponds to the de Sitter space,
which has no singularity.

Thus, for the above cosmology characterized by the Hubble rate in Eq. (218),
one can try to find the equivalent $F(G)$ models following the reconstruction
technique described in Refs. Nojiri:2006su ; Nojiri:2006gh ; Bamba:2008ut ;
Bamba:2010wfw . Now, with the choice of suitable functions $P(t)$ and $Q(t)$
of a scalar field $t$ which in this case is identified with the cosmic time,
the action integral in Eq. (310) can be expressed as

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa^{2}}\bigg{(}R+P(t)G+Q(t)\bigg{)}+{{L}}_{\mathrm{matter}}\right]\,.$
(319)

The variation of the action in Eq. (319) with respect to $t$ gives,

$\displaystyle\frac{dP(t)}{dt}G+\frac{dQ(t)}{dt}=0\,,$ (320)

from which one can in principle express $t=t(G)$. Now, substituting $t=t(G)$
into Eq. (319), one can express the action in terms of $F(G)$ where

$\displaystyle F(G)=P\left(t(G)\right)G+Q\left(t(G)\right)\,.$ (321)

Now, subtracting the equations of motion, namely, Eqs. (311) and (312), we
obtain a second order differential equation

$\displaystyle
8\frac{d}{dt}\left(H^{2}\frac{dP}{dt}\right)-8H^{3}\frac{dP}{dt}+4\dot{H}+2\kappa^{2}\rho_{0}(1+w)a^{-3(1+w)}=0\,,$
(322)

where $\rho_{0}$ refers to the present energy density of $\rho$. We note that
instead of considering only one fluid in Eqs. (314) and (315), if one would
consider several matter components with $w_{i}=p_{i}/\rho_{i}$ being the EoS
parameter of the $i$-th fluid, then in this case, the last term of the
equation (322) will be replaced by
$2\kappa^{2}\sum_{i}(1+w_{i})\rho_{i0}a^{-3(1+w_{i})}$. If one can solve for
$P(t)$ from Eq. (322), then from Eq. (311), one can find $Q(t)$ as follows

$\displaystyle
Q(t)=-24H^{3}\frac{dP}{dt}-6H^{2}+2\kappa^{2}\rho_{0}a^{-3(1+w)}\,.$ (323)

Thus, one can see that any cosmology which is prescribed in terms of any
Hubble rate, then this cosmology can be realized by some specific $F(G)$
gravity model. Just for simplicity, we neglect the matter term from the above
equations and focus on the specific $F(G)$ models described by the Hubble rate
given in Eq. (218). Therefore, under the absence of matter term, Eqs. (322)
and (323) become,

$\displaystyle
8\frac{d}{dt}\left(H\frac{dP}{dt}\right)-8H^{3}\frac{dP}{dt}+4\dot{H}=0\,,$
(324)

and

$\displaystyle Q(t)=-24H^{3}\frac{dP}{dt}-6H^{2}\,.$ (325)

In the following we show a correspondence between the finite-time
singularities and the associated $F(G)$ model.

(i) Big Rip singularity:

Let us focus on the Big Rip singularity which is realized by some of the
$F(G)$ gravity models. To begin with we consider $\beta=1$ and $H_{s}=0$ in
Eq. (218), for which the Hubble rate $H(t)$ takes the form of (63)
Bamba:2010wfw  and consequently one can derive $G$ as follows,

$\displaystyle G=\frac{24h_{s}^{3}}{(t_{s}-t)^{4}}(1+h_{s})\,.$ (326)

Thus, with the above choice of $H$, the most general solution of the second
order differential equation (324) is Bamba:2010wfw

$\displaystyle
P(t)=\frac{1}{4h_{s}(h_{s}-1)}(2t_{s}-t)t+c_{1}\frac{(t_{s}-t)^{3-h_{s}}}{3-h_{s}}+c_{2}\,,$
(327)

where $c_{1}$ and $c_{2}$ are arbitrary constants. Having the solution for
$P(t)$ as in Eq. (327), from Eq. (325), one can derive $Q(t)$ as Bamba:2010wfw

$\displaystyle
Q(t)=-\frac{6h_{s}^{2}}{(t_{s}-t)^{2}}-\frac{24h_{s}^{3}}{(t_{s}-t)^{3}}\left[\frac{(t_{s}-t)}{2h_{s}(h_{s}-1)}-c_{1}(t_{s}-t)^{2-h_{s}}\right]\,.$
(328)

Further, from Eq. (320) we obtain

$\displaystyle
t-t_{s}=\left[\dfrac{24(h_{s}^{3}+h_{s}^{4})}{G}\right]^{1/4}\,,$ (329)

which, as one can see, is consistent with Eq. (326). Thus, having all the
above information, using Eq. (321), one can write down the most general form
of $F(G)$ realizing the Big Rip singularity as follows Bamba:2010wfw :

$\displaystyle
F(G)=\left(\frac{\sqrt{6h_{s}^{3}(1+h_{s})}}{h_{s}(1-h_{s})}\right)\sqrt{G}+c_{1}G^{\frac{1+h_{s}}{4}}+c_{2}G\,.$
(330)

For $h_{s}=1$, that means when the Hubble rate becomes $H(t)=(t_{s}-t)^{-1}$,
one can also find another exact solution for $P(t)$ as follows,

$\displaystyle P(t)=\alpha(t_{s}-t)^{q}\ln\left[\gamma(t_{s}-t)^{z}\right]\,,$
(331)

where $\gamma$ is a positive real number, From Eq. (325), one can derive,

$\displaystyle
Q(t)=-\frac{12}{(t_{s}-t)^{2}}\ln\left[\gamma(t_{s}-t)\right]\,,$ (332)

Thus, using the expressions for $P(t)$ and $Q(t)$, the form of $F(G)$ is given
by Bamba:2010wfw

$\displaystyle F(G)=\frac{\sqrt{3}}{2}\sqrt{G}\ln(\gamma G)\,.$ (333)

This form of $F(G)$ corresponds to the Hubble rate (218) with $h_{s}=1$ and
$H_{s}=0$ and it realizes the Big Rip singularity. In general, for large
values of $G$, $F(G)\sim\alpha\sqrt{G}\ln(\gamma G)$ with $\alpha>0$,
$\gamma>0$, and the Big Rip singularity may appear. Moreover, the Big Rip
singularity may also appear for $F(G)\sim\alpha\sqrt{G}\ln(\gamma
G^{z}+G_{s})$ with $\alpha>0$, $\gamma>0$, $z>0$ where $G_{s}$ is a constant
such that $\gamma G^{z}+G_{s}>0$.

(ii) Other types of singularities:

Let us now discuss the appearance of other types of singularities in
$F(G)$-gravity models. We consider the Hubble rate of Eq. (218) with $H_{s}=0$
but $\beta\neq 1$. In this case, the evolution of the scale factor $a(t)$
turns out to be

$\displaystyle
a(t)=a_{0}\exp\left[\frac{h_{s}}{\beta-1}\left((t_{s}-t)^{1-\beta}-(t_{s}-t_{0})^{1-\beta}\right)\right]\,,$
(334)

and we observe the emergence of following singularities:

  * •

When $\beta>1$, then $H$ and $G$ are given by

$\displaystyle H=\frac{h_{s}}{(t_{s}-t)^{\beta}}\,,\quad
G\sim\frac{24h_{s}^{4}}{(t_{s}-t)^{4\beta}}\,.$ (335)

A solution of the Eq. (324) in the limit $t\rightarrow t_{s}$ can be found as

$\displaystyle P(t)\simeq\frac{\alpha}{(t_{s}-t)^{z}}\,,$ (336)

where $z=-2\beta$ and $\alpha=-1/4h_{s}^{2}$. The expression of $F(G)$ now
follows

$\displaystyle F(G)=-12\sqrt{\frac{G}{24}}\,.$ (337)

Hence, if for large values of $G$, $F(G)\sim-\alpha\sqrt{G}$ with $\alpha>0$,
then a Type I singularity could appear.

  * •

When $0<\beta<1$ for which the forms of $H$ and $G$ are given by

$\displaystyle H=\frac{h_{s}}{(t_{s}-t)^{\beta}}\,,\quad
G\sim\frac{24h_{s}^{3}\beta}{(t_{s}-t)^{3\beta+1}}\,.$ (338)

In this case, an asymptotic solution of Eq. (324) in the limit $t\rightarrow
t_{s}$ can be found as

$\displaystyle P(t)\simeq\frac{\alpha}{(t_{s}-t)^{z}}\;,$ (339)

where $z=-(1+\beta)$ and $\alpha=1/2h_{s}(1+\beta)$. The form of $F(G)$ now
becomes

$\displaystyle
F(G)=\frac{6h_{s}^{2}}{(\beta+1)}(3\beta+1)\left(\frac{|G|}{24h_{s}^{3}|\beta|}\right)^{2\beta/(3\beta+1)}\,.$
(340)

Hence, if for large values of $G$, $F(G)$ has the following form

$\displaystyle
F(G)\sim\alpha|G|^{\gamma}\,,\quad\gamma=\frac{2\beta}{3\beta+1}\,,$ (341)

where $\alpha>0$ and $0<\gamma<1/2$, then because we are assuming $0<\beta<1$
and in that case, a Type III singularity could emerge.

Finally, for a general value of $\beta$ we have the following:

    * –

If for $G\rightarrow-\infty$, $F(G)$ has the form as in Eq. (341) with
$\alpha>0$ and $-\infty<\gamma<0$, then we find $-1/3<\beta<0$ and a Type II
(sudden) singularity could appear.

    * –

If for $G\rightarrow 0^{-}$, $F(G)$ takes the form of Eq. (341) with
$\alpha<0$ and $1<\gamma<\infty$, then we obtain $-1<\beta<-1/3$ and a Type II
singularity could occur.

    * –

If for $G\rightarrow 0^{-}$, $F(G)$ assumes the form as in Eq. (341) with
$\alpha>0$ and $2/3<\gamma<1$, then we obtain $-\infty<\beta<-1$ and a Type IV
singularity could appear. We also require that $\gamma\neq 2n/(3n-1)$, i.e.,
$\beta\not=n$, where $n$ is a natural number.

Let us note that we can generate all the possible Type II singularities as
shown above except for $\beta=-1/3$, i.e., for $H=h_{s}(t_{s}-t)^{1/3}$
because $\gamma=0$. In this case, $G$ takes the form:

$\displaystyle G=24h_{s}^{3}\beta+24h_{s}^{4}(t_{s}-t)^{4/3}<0\,.$ (342)

Thus, in order to express $t$ in terms of $G$, it is essential to consider the
whole expression of $G$, and by considering the leading term involving
$(t_{s}-t)$ in (322) and (323) with $\rho_{0}=0$, we obtain

$\displaystyle
F(G)\simeq\frac{1}{4\sqrt{6}h_{s}^{3}}G(G+8h_{s}^{3})^{1/2}+\frac{2}{\sqrt{6}}(G+8h_{s}^{3})^{1/2}\,,$
(343)

which satisfies Eq. (316) in the limit $t\rightarrow t_{s}$. Consequently, the
specific model $F(G)=\sigma_{1}G(G+c_{3})^{1/2}+\sigma_{2}(G+c_{3})^{1/2}$,
where $\sigma_{1}$, $\sigma_{2}$ and $c_{3}$ are positive constants, can
generate a Type II singularity.

###  IV.7 $F(R,G)$ Gravity

The $F(R,G)$-gravity is a very generalized gravitational theory where $F$ is a
generalized function of the Ricci scalar $R$ and the Gauss-Bonnet invariant is
defined in Eq. (161). One can clearly see that the $F(R,G)$ gravity theory can
recover $F(R)$-gravity and $F(G)$-gravity as special cases. Moreover, with the
suitable choice for $F$ leading to $F(R,G)=R$, one recovers Einstein’s GR as
well. Thus, one can see that this modified gravity theory being the
generalized version of both $F(R)$ and $F(G)$ gravity theories can offer some
appealing consequences and very soon of its introduction, $F(R,G)$ gravity
theory received significant attention in the community Bamba:2010wfw ;
Elizalde:2010jx ; DeFelice:2010hg ; DeFelice:2011ka ; delaCruz-Dombriz:2011oii
; Makarenko:2012gm ; Atazadeh:2013cz ; SantosDaCosta:2018bbw ; Wu:2018jve ;
Odintsov:2018nch ; Camci:2018apx ; Barros:2019pvc ; KumarSanyal:2019rpa ;
Shekh:2019rib ; Inagaki:2019rhm ; Elizalde:2020zcb ; Navo:2020eqt ;
Odintsov:2021nim ; Lohakare:2021yuo ; Nojiri:2021mxf . In the present section
we shall investigate the finite-time future singularities appearing in this
generalized modified gravitational theory. This model includes ghosts
DeFelice:2009ak  and the ghost-free modifications have been proposed in Ref.
Nojiri:2021mxf . In this section, however, again for the illustration, we show
the structure of the singularities based on the original models.

The action of the ${F}(R,G)$-gravity is given by Bamba:2010wfw ; delaCruz-
Dombriz:2011oii ; Atazadeh:2013cz :

$\displaystyle S=\int
d^{4}x\sqrt{-g}\left[\frac{{F}(R,G)}{2\kappa^{2}}+{{L}}_{\mathrm{matter}}\right]\,,$
(344)

For the action (344), one can derive the gravitational equations as follows

$\displaystyle{{F}}^{\prime}_{R}\left(R_{\mu\nu}-\frac{1}{2}Rg_{\mu\nu}\right)=$
$\displaystyle\,\kappa^{2}T^{(\mathrm{matter})}_{\mu\nu}+\frac{1}{2}g_{\mu\nu}\left({F}-{{F}}^{\prime}_{R}R\right)+{\nabla}_{\mu}{\nabla}_{\nu}{{F}}^{\prime}_{R}-g_{\mu\nu}\Box{{F}}^{\prime}_{R}$
$\displaystyle\,+\left(-2RR_{\mu\nu}+4R_{\mu\rho}R_{\nu}{}^{\rho}-2R_{\mu}{}^{\rho\sigma\tau}R_{\nu\rho\sigma\tau}+4g^{\alpha\rho}g^{\beta\sigma}R_{\mu\alpha\nu\beta}R_{\rho\sigma}\right){{F}}^{\prime}_{G}$
$\displaystyle\,+2\left({\nabla}_{\mu}{\nabla}_{\nu}{{F}}^{\prime}_{G}\right)R-2g_{\mu\nu}\left(\Box{{F}}^{\prime}_{G}\right)R+4\left(\Box{{F}}^{\prime}_{G}\right)R_{\mu\nu}-4\left({\nabla}_{\rho}{\nabla}_{\mu}{{F}}^{\prime}_{G}\right)R_{\nu}{}^{\rho}$
$\displaystyle\,-4\left({\nabla}_{\rho}{\nabla}_{\nu}{{F}}^{\prime}_{G}\right)R_{\mu}{}^{\rho}+4g_{\mu\nu}\left({\nabla}_{\rho}{\nabla}_{\sigma}{{F}}^{\prime}_{G}\right)R^{\rho\sigma}-4\left({\nabla}_{\rho}{\nabla}_{\sigma}{{F}}^{\prime}_{G}\right)g^{\alpha\rho}g^{\beta\sigma}R_{\mu\alpha\nu\beta}\,,$
(345)

where ${\nabla}_{\mu}$ is the covariant derivative operator associated with
the metric tensor $g_{\mu\nu}$; $\Box\equiv
g^{\mu\nu}{\nabla}_{\mu}{\nabla}_{\nu}$ denotes the covariant d’Alembertian
for a scalar field;
$T^{(\mathrm{matter})}_{\mu\nu}=\mathrm{diag}\left(\rho,p,p,p\right)$ is the
energy-momentum tensor of the matter sector that includes all the ordinary
matter fluids with $\rho~{}(=\sum_{i}\rho_{i})$ and $p~{}(=\sum_{i}p_{i})$ are
the total energy density and pressure of all the fluids, respectively
($\rho_{i}$, $p_{i}$ being the representatives of the $i$-th fluid), and

$\displaystyle{{F}}^{\prime}_{R}=\frac{\partial{F}(R,G)}{\partial
R}\,,\quad{{F}}^{\prime}_{G}=\frac{\partial{F}(R,G)}{\partial G}\,.$ (346)

In the background of a spatially flat FLRW Universe (3), the gravitational
field equations in (345) can be expressed as

$\displaystyle\rho_{\mathrm{eff}}=\frac{3}{\kappa^{2}}H^{2}\,,\quad
p_{\mathrm{eff}}=-\frac{1}{\kappa^{2}}\left(2\dot{H}+3H^{2}\right)\,,$ (347)

where $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$, termed as the effective
energy density and pressure of the Universe, respectively, are given by

$\displaystyle\rho_{\mathrm{eff}}=$
$\displaystyle\,\frac{1}{{{F}}^{\prime}_{R}}\left[\rho+\frac{1}{2\kappa^{2}}\left\\{\left({{F}}^{\prime}_{R}R-{F}\right)-6H{\dot{{F}}}^{\prime}_{R}+G{{F}}^{\prime}_{G}-24H^{3}{\dot{{F}}}^{\prime}_{G}\right\\}\right]\,,$
(348) $\displaystyle p_{\mathrm{eff}}=$
$\displaystyle\,\frac{1}{{{F}}^{\prime}_{R}}\left[p+\frac{1}{2\kappa^{2}}\left\\{-\left({{F}}^{\prime}_{R}R-{F}\right)+4H{\dot{{F}}}^{\prime}_{R}+2{\ddot{{F}}}^{\prime}_{R}-G{{F}}^{\prime}_{G}+16H\left(\dot{H}+H^{2}\right){\dot{{F}}}^{\prime}_{G}+8H^{2}{\ddot{{F}}}^{\prime}_{G}\right\\}\right]\,.$
(349)

Let us now proceed towards the investigations of the finite-time singularities
in the ${F}(R,G)$-gravity. Similar to the earlier section IV.6 on $F(G)$, we
reconstruct the ${F}(R,G)$-gravity models that can produce finite-time
singularities. In order to do so, we consider the pure gravitation action of
${F}(R,G)$-gravity, that means action (344) without matter sector exactly what
we have considered in the earlier section IV.6. In this case, from Eqs. (348)
and (349), one can write down the gravitational equations as

$\displaystyle 0=$
$\displaystyle\,24H^{3}{\dot{{F}}}^{\prime}_{G}+6H^{2}{{F}}^{\prime}_{R}+6H{\dot{{F}}}^{\prime}_{R}+({F}-R{{F}}^{\prime}_{R}-G{{F}}^{\prime}_{G})\,,$
(350) $\displaystyle 0=$
$\displaystyle\,8H^{2}{\ddot{{F}}}^{\prime}_{G}+2{\ddot{{F}}}^{\prime}_{R}+4H{\dot{{F}}}^{\prime}_{R}+16H{\dot{{F}}}^{\prime}_{G}(\dot{H}+H^{2})+{{F}}^{\prime}_{R}(4\dot{H}+6H^{2})+{F}-R{{F}}^{\prime}_{R}-G{{F}}^{\prime}_{G}\,.$
(351)

We note that in the case of pure gravity (i.e., when no matter sector is
present in the gravitational action), the above equations (350) and (351) are
linearly dependent.

Now, following the similar fashion as adopted in the earlier section (IV.6),
with the choice of the proper functions $P(t)$, $Z(t)$ and $Q(t)$ of a scalar
field (where we identify the scalar field with time $t$), one can rewrite Eq.
(344) without ${{L}}_{\mathrm{matter}}$ as

$\displaystyle S=\frac{1}{2\kappa^{2}}\int
d^{4}x\sqrt{-g}\Bigg{(}P(t)R+Z(t)G+Q(t)\Bigg{)}\,.$ (352)

By varying the action in Eq. (319) with respect to $t$, we obtain

$\displaystyle P^{\prime}(t)R+Z^{\prime}(t)G+Q^{\prime}(t)=0\,,$ (353)

where the prime denotes differentiation with respect to $t$. From Eq. (353),
one can in principle express $t$ as a function of $R,G$, i.e., $t=t(R,G)$.
Now, with the substitution of $t=t(R,G)$ in Eq. (319), we can express
${F}(R,G)$ as

$\displaystyle{F}(R,G)=P\left(t(R,G)\right)R+Z\left(t(R,G)\right)G+Q\left(t(R,G)\right)\,.$
(354)

####  IV.7.1 Finite-time Singularities

We consider, once again, the following expression of the Hubble parameter in
(218) Bamba:2010wfw . From Eq. (218), one can express the scale factor as

$\displaystyle a(t)=\bar{a}\exp\left(g(t)\right)\,,$ (355)

where $\bar{a}$ is a constant and $g(t)$ is some differentiable function of
$t$ satisfying $\dot{g}(t)=H(t)$. Using the conservation law and Eq. (350), we
arrive at the following second order differential equation

$\displaystyle
P^{\prime\prime}(t)+4\dot{g}^{2}(t)Z^{\prime\prime}(t)-\dot{g}(t)P^{\prime}(t)+(8\dot{g}\ddot{g}-4\dot{g}^{3}(t))Z^{\prime}(t)+2\ddot{g}(t)P(t)=0\,,$
(356)

where we have implemented the expression of the scale factor as given in Eq.
(355). Now, from Eq. (350), one can derive $Q(t)$ as

$\displaystyle
Q(t)=-24\dot{g}^{3}(t)Z^{\prime}(t)-6\dot{g}^{2}(t)P(t)-6\dot{g}(t)P^{\prime}(t)\,.$
(357)

If $P(t)\neq 0$, then ${F}(R,G)$ can be written as

$\displaystyle{F}(R,G)=R\;\psi(R,G)+f(R,G)\,,$ (358)

where $\psi(R,G)\neq 0$ and $f(R,G)$ is any generic function of $R$ and $G$.
Now, from Eqs. (348) and (349), we obtain,

$\displaystyle\rho_{\mathrm{eff}}=-\frac{1}{2\kappa^{2}g(R,G)}\biggl{[}24H^{3}{\dot{{F}}}^{\prime}_{G}+6H^{2}\left(R\frac{d\psi(R,G)}{dR}+\frac{df(R,G)}{dR}\right)+6H{\dot{{F}}}^{\prime}_{R}+({F}-R{{F}}^{\prime}_{R}-G{{F}}^{\prime}_{G})\biggr{]},$
(359)

and

$\displaystyle p_{\mathrm{eff}}=$
$\displaystyle\,\frac{1}{2\kappa^{2}g(R,G)}\biggl{[}8H^{2}{\ddot{{F}}}^{\prime}_{G}+2{\ddot{{F}}}^{\prime}_{R}+4H{\dot{{F}}}^{\prime}_{R}+16H{\dot{{F}}}^{\prime}_{G}(\dot{H}+H^{2})$
$\displaystyle\,+\left(R\frac{d\psi(R,G)}{dR}+\frac{df(R,G)}{dR}\right)(4\dot{H}+6H^{2})+{F}-R{{F}}^{\prime}_{R}-G{{F}}^{\prime}_{G}\biggr{]}\,,$
(360)

and similar to the earlier Section IV.6, the function
${G}\left(H,\dot{H},\ddot{H},...\right)$ can be written as

$\displaystyle{G}\left(H,\dot{H},\ddot{H},\dddot{H},...\right)=-\frac{1}{\kappa^{2}}\left[2\dot{H}+3(1+w)H^{2}\right],$
(361)

where more explicitly,

$\displaystyle{G}(H,\dot{H},\ddot{H},...)=$
$\displaystyle\,p_{\mathrm{eff}}-w\rho_{\mathrm{eff}}$ $\displaystyle=$
$\displaystyle\,\frac{1}{2\kappa^{2}\psi(R,G)}\Bigg{[}(1+w)({F}-R{{F}}^{\prime}_{R}-G{{F}}^{\prime}_{G})+\left(R\frac{d\psi(R,G)}{dR}+\frac{df(R,G)}{dR}\right)\left(6H^{2}(1+w)+4\dot{H}\right)$
$\displaystyle\,+H{\dot{{F}}}^{\prime}_{R}(4+6w)+8H{\dot{{F}}}^{\prime}_{G}\left(2\dot{H}+H^{2}(2+3w)\right)+2{\ddot{{F}}}^{\prime}_{R}+8H^{2}{\ddot{{F}}}^{\prime}_{G}\Bigg{]}\,,$
(362)

where $w$ is the EoS parameter of the matter. It is important to note that the
use of the above equation (362) demands $\psi(R,G)\neq 0$ on the solution. Let
us now explicitly describe the finite-time singularities allowed in various
${F}(R,G)$-gravity models.

(i) Big Rip singularity

Let us first discuss the appearance of the Big Rip singularity in this
modified gravity model. We consider the Hubble rate of Eq. (218) with
$\beta=1$ and $H_{s}=0$ which gives the Hubble rate in (63) and consequently,
one can derive $R$ and $G$ for Eq. (63) as

$\displaystyle R=\frac{6h_{s}}{(t_{s}-t)^{2}}(2h_{s}+1)\,,\quad
G=\frac{24h_{s}^{3}}{(t_{s}-t)^{4}}(1+h_{s})\,.$ (363)

A very trivial solution of Eq. (356) can be given by

$\displaystyle P(t)=\alpha(t_{s}-t)^{z}\,,\quad Z(t)=\delta(t_{s}-t)^{x}\,,$
(364)

where $\alpha$ and $\delta$ be any real numbers; $x=3-h_{s}$, and $z$ can be
found from

$\displaystyle z_{\pm}=\frac{1-h_{s}\pm\sqrt{h_{s}^{2}-10h_{s}+1}}{2}\,.$
(365)

Therefore, the most general solution of $P(t)$ can be given as

$\displaystyle
P(t)=\alpha_{1}(t_{s}-t)^{z_{+}}+\alpha_{2}(t_{s}-t)^{z_{-}}\,,$ (366)

where $\alpha_{1}$ and $\alpha_{2}$ are any real numbers. Now, from Eq. (357),
one can derive

$\displaystyle
Q(t)=\dfrac{24h_{s}^{3}\delta(3-h_{s})}{(t_{s}-t)^{h_{s}+1}}+\dfrac{6h_{s}\alpha_{1}(z_{+}-h_{s})}{(t_{s}-t)^{2-z_{+}}}+\dfrac{6h_{s}\alpha_{2}(z_{-}-h_{s})}{(t_{s}-t)^{2-z_{-}}}\,.$
(367)

Now for the condition $0<h_{s}<5-2\sqrt{6}$ or $h_{s}>2+\sqrt{6}$, the
solution of ${F}(R,G)$ is given by

$\displaystyle{F}(R,G)=\bar{\alpha}_{1}R^{1-\frac{z_{+}}{2}}+\bar{\alpha}_{2}R^{1-\frac{z_{-}}{2}}+\bar{\delta}G^{\frac{h_{s}+1}{4}}\,,$
(368)

where some factors have been absorbed into the constants. Note also that
$z_{\pm}\neq 0$, otherwise $h_{s}$ vanishes which is a contradiction.

Another solution of Eq. (356) can be given by

$\displaystyle P(t)=\frac{\alpha}{(t_{s}-t)^{z}}\,,\quad
Z(t)=\frac{\delta}{(t_{s}-t)^{x}}\,,$ (369)

where $\delta$, $x$ are real numbers; $z=x+2$ and $\alpha$ is given by

$\displaystyle\alpha=\dfrac{4h_{s}^{2}\delta
x(h_{s}-x-3)}{x^{2}+(5-h_{s})x+6}\,.$ (370)

From Eq. (357), one finds

$\displaystyle Q(t)=-\frac{6h_{s}}{(t_{s}-t)^{x+4}}\left[4h_{s}^{2}\delta
x+\alpha(x+2+h_{s})\right]\,.$ (371)

The solution of the Eq. (353) is given by

$\displaystyle
t=t(R,G)=t_{s}-\left[\;\;\dfrac{-\alpha(x+2)R\pm\sqrt{\alpha^{2}(x+2)^{2}R^{2}+24h_{s}\bigl{(}4h_{s}^{2}\delta
x+\alpha(x+2+h_{s})\bigr{)}(x+4)\delta xG}}{2\delta xG}\;\;\right]^{1/2}\,,$
(372)

where $x\neq 0$ and $\delta\neq 0$. In order to get real solutions, one has to
ensure that first of all, we should have
$\Delta=\alpha^{2}(x+2)^{2}R^{2}+24h_{s}\bigl{(}4h_{s}^{2}\delta
x+\alpha(x+2+h_{s})\bigr{)}(x+4)\geq 0$ in Eq. (372) and secondly, the entire
term inside the third brace of Eq. (372) has to be positive. For $h_{s}>0$,
the principal cases with the real solutions of (372) are the following:

  * •

Case I: For $x>0$, $\delta>0$, $1+x\leq h_{s}<x+5+\frac{6}{x}$, we need to use
the sign $+$ in the r.h.s. of Eq. (372).

  * •

Case II: For $-\frac{3}{2}\leq x<0$, $\delta<0$, $h_{s}\geq x+1$, we need to
use the sign $+$ in the r.h.s. of Eq. (372).

  * •

Case III: For $-4<x<-\frac{3}{2}$, $\delta<0$, $h_{s}>x+5+\frac{6}{x}$, we
must use the sign $+$ in the r.h.s. of Eq. (372).

  * •

Case IV: For $x>0$, $\delta<0$, $x+5+\frac{6}{x}>h_{s}\geq 1+x$, we need to
use the sign $-$ in the r.h.s. of Eq. (372).

  * •

Case V: For $-\frac{3}{2}\leq x<0$, $\delta>0$, $h_{s}\geq x+1$, we must use
the sign $-$ in the r.h.s. of Eq. (372).

  * •

Case VI: For $-4<x<-\frac{3}{2}$, $\delta>0$, $h_{s}>x+5+\frac{6}{x}$, we need
to use the sign $-$ in the right hand side of Eq. (372).

  * •

Case VII: For $x=-4$, $\delta>0$, we must use the sign $-$ in the r.h.s. of
Eq. (372).

  * •

Case VIII: For $x=-4$, $\delta<0$, we must use the sign $+$ in the r.h.s. of
Eq. (372).

The solution of ${F}(R,G)$ is then given by

$\displaystyle{F}(R,G)=\left(\frac{\alpha}{(t_{s}-t(R,G))^{x+2}}\right)R+\left(\frac{\delta}{(t_{s}-t(R,G))^{x}}\right)G-\left(\frac{6h_{s}}{(t_{s}-t(R,G))^{x+4}}\right)\Bigg{[}4h^{2}\delta
x+\alpha(x+2+h_{s})\Bigg{]}\,,$ (373)

where $t(R,G)$ is given by Eq. (372). The expression for ${F}(R,G)$ in Eq.
(373) represents an exact solution of the equation of motion in Eqs. (348) and
(349) allowing the Big Rip singularity. Let us show some specific examples of
${F}(R,G)$-gravity models allowing the Big Rip singularity obtained for some
specific values of $\alpha$ and $x$.

  * •

For $\alpha=1$ and $x=-2$, we find

$\displaystyle{F}(R,G)=R+\left(\dfrac{\sqrt{6}\sqrt{h_{s}(1+h_{s})}}{(1-h_{s})}\right)\sqrt{G}\,,\quad
h\neq 1\,.$ (374)

  * •

If $\alpha=0$ and $x=h_{s}-3$ (this case corresponds to Case (I) $-$ Case
(VI), presented above), then one finds

$\displaystyle{F}(R,G)=\delta\;G^{\frac{h_{s}+1}{4}}\,,\quad\delta\neq 0\,,$
(375)

which is actually equivalent to Eq. (368) for $\alpha_{1}=\alpha_{2}=0$.

  * •

If $x=-4$, then the model reduces to

$\displaystyle{F}(R,G)=\dfrac{16h_{s}^{4}\delta}{(1+2h_{s}^{2})^{2}}\left[(9+21h_{s}+6h_{s}^{2})-(1+h_{s})^{2}\frac{R^{2}}{G}\right]\,,\quad\delta\neq
0\,.$ (376)

Thus, we see that for large values of $R$ and $G$,
${F}(R,G)\sim\pm\alpha\mp\delta(R^{2}/G)$ with $\alpha>0$ and $\delta>0$, then
the Big Rip singularity could appear.

  * •

For $x=h_{s}-1$, the model becomes (by absorbing some constants)

$\displaystyle{F}(R,G)=\delta\left(\dfrac{R}{G}\right)^{\frac{1-h_{s}}{2}}G\,,\quad\delta\neq
0\,,\quad h_{s}\neq 1\,.$ (377)

Thus, we see that for large values of $R$ and $G$, ${F}(R,G)\sim\delta
G^{\gamma}/R^{\gamma-1}$ with $\delta\neq 0$ and $\frac{1}{2}<\gamma<1$ or
$1<\gamma<+\infty$, and the Big Rip singularity could appear.

We close the discussion with a general model that allows Big Rip singularity.
One can verify that the model

$\displaystyle{F}(R,G)=\gamma\frac{G^{m}}{R^{n}}\,,$ (378)

with $\gamma$ being any real number, is a solution of Eqs. (348) and (349) in
the case of the Big Rip singularity for some value of $h_{s}$. In general, it
is possible to obtain solutions for $h_{s}>0$ if $m>0$, $n>0$ and $m>n$. For
example, the case $n=2$ and $m=3$ realizes the Big Rip singularity in
$h_{s}=5$; the case $n=1$ and $m=3$ realizes the big rip singularity in
$h_{s}=4+\sqrt{19}$ and so on. This is a generalization of Eq. (377). It is
important here to mention that for $m=-1$ and $n=-2$, we do not recover a
physical scenario because in this case we obtain $h_{s}=-3$. However, this
case (i.e., $m=-1$ and $n=-2$), has a similarity with Eq. (376) allowing the
Big Rip singularity. Finally, for the case with $m=0$ or $n=0$, we recover the
model in Eq. (368).

(ii) Other types of singularities

Here we discuss other types of singularities that may appear in various
${F}(R,G)$-gravity models. To begin with we consider the following expansion
of the Hubble rate in (114). For the above choice of the Hubble rate, an exact
solution of the differential Eq. (356) can be found as

$\displaystyle P(t)=-\lambda(4h_{s}^{2})(t_{s}-t)\,,\quad
Z(t)=\lambda(t_{s}-t)^{2\beta+1}\,,$ (379)

where $\lambda$ is a constant. The form of $Q(t)$ from Eq. (357) can now be
derived as

$\displaystyle
Q(t)=\dfrac{24h_{s}^{4}\lambda}{(t_{s}-t)^{2\beta-1}}+\dfrac{48h_{s}^{3}\beta}{(t_{s}-t)^{\beta}}\,.$
(380)

Now we examine the cases when $\beta>1$ or $\beta<1$.

  * •

For $\beta>1$, one may obtain the asymptotic real solution of Eq. (353):

$\displaystyle
t=t(R,G)=t_{s}-2^{1/2\beta}\left[\dfrac{h_{s}^{2}R+\sqrt{h_{s}^{4}R^{2}+6h_{s}^{4}(4\beta^{2}-1)G}}{(1+2\beta)G}\right]^{1/2\beta}\,.$
(381)

Thus, the mathematical form of ${F}(R,G)$ can now be expressed as

$\displaystyle{F}(R,G)=-4h_{s}^{2}\lambda(t_{s}-t(R,G))R+\lambda(t_{s}-t(R,G))^{1+2\beta}G+24h_{s}^{4}\lambda(t_{s}-t(R,G))^{1-2\beta}\,,\quad\beta>1\,.$
(382)

This is an asymptotic solution of Eq. (361) for $\beta>1$, when

$\displaystyle-\frac{1}{\kappa^{2}}\left[2\dot{H}+3(1+w)H^{2}\right]\sim-\dfrac{3(1+w)h_{s}^{2}}{\kappa^{2}}(t_{s}-t)^{-2\beta}\,.$
(383)

Interestingly, for $\beta\gg 1$, the mathematical form of ${F}(R,G)$ is given
by

$\displaystyle{F}(R,G)\simeq\lambda\left(\dfrac{\alpha G}{R+\sqrt{R^{2}+\gamma
G}}-R\right)\,,\quad\alpha>0\,,\quad\gamma>0\,,\quad\lambda\neq 0\,,$ (384)

and this is the asymptotic behavior of a ${F}(R,G)$ model where a “strong”
Type I singularity ($\beta\gg 1$) may appear.

We consider some other explicit cases where Type I singularity appears as
follows. Taking

$\displaystyle{F}(R,G)=\gamma\frac{G^{m}}{R^{n}}\,,$ (385)

and using Eqs. (348) and (349), one can verify that the function
${G}\left(H,\dot{H},\ddot{H},...\right)$ in Eq. (362) becomes,

$\displaystyle{G}(H,\dot{H},..)\simeq-\frac{3h_{s}^{2}(2m-n-1)(1+w)}{\kappa^{2}(t_{s}-t)^{2\beta}}\,,$
(386)

which under the restriction $2m-n-1>0$, represents an asymptotic solution of
Eq. (361) for $\beta>1$. Thus, we see that for the model ${F}(R,G)\simeq\gamma
G^{m}/R^{n}$ with $m>(n+1)/2$, the Type I singularity may appear. This
observation has some important consequences in the context of other modified
gravity models because one can see that the modified gravity models where
either $F(R)=R^{n}$ with $n>1$ or $F(G)=G^{m}$ with $m>1/2$ can allow
singularities.

One can also find some other models allowing the Type I singularities and we
can construct these models following section IV.6 as follows. The Type I
singularities actually correspond to the asymptotic limits for $R$ and $G$

$\displaystyle R\sim 12H^{2}\,,\quad G\sim 24H^{4}\,.$ (387)

Since $R$ and $H$ in the asymptotic limits are nothing but the functions of
the Hubble parameter only, therefore, one can write that,

$\displaystyle\lim_{t\rightarrow
t_{s}}24\left(\frac{R}{12}\right)^{2}=\lim_{t\rightarrow t_{s}}G\,.$ (388)

Now, if we substitute $G$ for $R$ in Eq. (337) considering Eq. (388), then we
obtain a zero function555The reason for obtaining the zero function is that,
in this case Eq. (337) becomes zero on the singularity solution. . Thus, if we
instead substitute $G$ for $G/R$, we therefore obtain the following model

$\displaystyle{F}(R,G)=R-\dfrac{6G}{R}\,,$ (389)

which is an asymptotic solution of Eq. (361) such as Eq. (383). Thus, we see
that Type I singularity appears in the model ${F}(R,G)\sim R-\alpha(G/R)$ with
$\alpha>0$.

  * •

For the Hubble rate in Eq. (114) with $\beta<1$, it is not possible to express
$G$ and $R$ as the functions of the same variable (i.e., $H$ or the same
combination of $H$ and $\dot{H}$). However, if we examine the asymptotic
behavior of $G$ and $R$, we have

$\displaystyle R\simeq\frac{6h_{s}\beta}{(t_{s}-t)^{\beta+1}}\,,\quad
G\simeq\frac{24h_{s}^{3}\beta}{(t_{s}-t)^{3\beta+1}}\,,$ (390)

and

$\displaystyle\frac{G}{R}\sim(t_{s}-t)^{-2\beta}\sim
G^{\frac{2\beta}{3\beta+1}}$ (391)

If we replace $G$ by $G/R$ (as given in Eq. (391)) in Eq. (340), then we find
that the asymptotic time dependence in Eq. (361) for $\beta<1$ is exactly same
as follows:

$\displaystyle-\frac{1}{\kappa^{2}}\left[2\dot{H}+3(1+w)H^{2}\right]\sim\frac{\alpha}{(t_{s}-t)^{\beta+1}}+\frac{\gamma}{(t_{s}-t)^{2\beta}}\,.$
(392)

With this consideration, one can derive a specific ${F}(R,G)$-gravity model
(by setting some parameters) from Eq. (340) as

$\displaystyle{F}(R,G)=R+\frac{3(3\beta+1)}{2(\beta+1)}\frac{G}{R}\,,$ (393)

which may encounter with other types of singularities. So, in the model
${F}(R,G)\sim R+\bar{\alpha}(G/R)$ with $\bar{\alpha}>0$, Type II, Type III
and Type IV singularities may appear. Then, by substituting $G$ for $R$ we get

$\displaystyle{F}(R,G)\simeq R+\bar{\delta}|R|^{\frac{2\beta}{1+\beta}}\,,$
(394)

which is a well-know result in cosmology. Now we have the following
observations:

    * –

In the model ${F}(R\rightarrow\infty)\sim
R+\bar{\delta}R^{\frac{2\beta}{1+\beta}}$, for $0<\beta<1$ and
$\bar{\delta}>0$, a Type III singularity may appear.

    * –

In the model ${F}(R\rightarrow-\infty)\sim
R+\bar{\delta}|R|^{\frac{2\beta}{1+\beta}}$, for $-1<\beta<0$ and
$\bar{\delta}>0$, a Type II singularity may appear.

    * –

In the model ${F}(R\rightarrow 0^{-})\sim
R+\bar{\delta}|R|^{\frac{2\beta}{1+\beta}}$, for $\beta<-1$ and $\beta\not=-n$
(being $n$ a natural number) and $\bar{\delta}<0$, a Type IV singularity may
appear.

###  IV.8 $F(T)$ Gravity

One of the alternative gravitational theories beyond Einstein’s GR is the
Teleparallel Equivalent of GR (TEGR) where instead of the curvature scalar $R$
defined by the Levi-Civita connection one uses the concept of torsion scalar
$T$ defined by the Weitzenböck connection Hehl:1976kj . The idea of this
gravitational theory was originally introduced by Einstein in 1928 under the
name “Fern-Parallelismus” or “distant parallelism” or “teleparallelism”
Einstein ; 2005physics…3046U . It was found that the modified teleparallel
gravity theory $F(T)$ in which $F$ is an arbitrary function of the torsion
scalar $T$, can explain the late-time accelerating expansion of the Universe
for some suitable choices of $F(T)$ Bengochea:2008gz ; Linder:2010py .
Moreover, it was also recognized that the $F(T)$ models can also take an
active role to describe the early inflationary phase Rezazadeh:2017edd ;
Awad:2017ign ; Bamba:2016wjm ; Rezazadeh:2015dza ; Nashed:2014lva . With such
appealing consequences, the modified teleparallel gravitation models got
massive attention in the astrophysical and cosmological community
Paliathanasis:2021uvd ; Paliathanasis:2017htk ; Paliathanasis:2016vsw ;
Bamba:2016gbu ; Haro:2014xtk ; Paliathanasis:2014iva ; Amoros:2013nxa ;
Basilakos:2013rua ; Bamba:2012vg . We refer to two recent reviews on $F(T)$
teleparallel gravity models Cai:2015emx ; Bahamonde:2021gfp . In Ref.
Izumi:2012qj ; Ong:2013qja , however, it has been proved that superluminal
modes appear in the $F(T)$ gravity and therefore the $F(T)$ gravity model is
physically inconsistent. In this section, however, just for the illustration,
we show the structure of the singularities in the framework of the $F(T)$
gravity.

In order to discuss the teleparallelism, we use the orthonormal tetrad
components $e_{A}(x^{\mu})$. We assume that the index $A$ runs over $0$, $1$,
$2$, $3$ for the tangent space at each point $x^{\mu}$ of the manifold and
$e_{A}^{\mu}$ forms the tangent vector of the manifold. The relation between
the metric $g^{\mu\nu}$ and the tetrad components is given by
$g_{\mu\nu}=\eta_{AB}e^{A}_{\mu}\mathrm{e}^{B}_{\nu}$ where
$\eta_{AB}=\mathrm{diag}(-1,1,1,1)$. Compared to the General theory of
Relativity where the torsionless Levi-Civita connection is used, in
Teleparallelism the curvatureless Weitzenböck connection is used Weitzenbock ,
hence, the gravitational field is described by the non-null torsion tensor

$\displaystyle T^{\rho}_{\
\mu\nu}\equiv\mathrm{e}^{\rho}_{A}\left(\partial_{\mu}e^{A}_{\nu}-\partial_{\nu}e^{A}_{\mu}\right)\,,$
(395)

The torsion scalar $T$, one of the main ingredients of the teleparallel
gravity, is defined by Hayashi:1979qx ; Maluf:2013gaa :

$\displaystyle T\equiv S_{\rho}^{\ \mu\nu}T^{\rho}_{\ \mu\nu}\,,$ (396)

where

$\displaystyle S_{\rho}^{\ \mu\nu}\equiv\frac{1}{2}\left(K^{\mu\nu}_{\ \
\rho}+\delta^{\mu}_{\rho}T^{\alpha\nu}_{\ \
\alpha}-\delta^{\nu}_{\rho}T^{\alpha\mu}_{\ \ \alpha}\right)\,,$ (397)

and $K^{\mu\nu}_{\verb||\rho}$, the contorsion tensor is given by

$\displaystyle K^{\mu\nu}_{\ \ \rho}\equiv-\frac{1}{2}\left(T^{\mu\nu}_{\ \
\rho}-T^{\nu\mu}_{\ \ \rho}-T_{\rho}^{\ \ \mu\nu}\right)\,.$ (398)

The action of the modified teleparallel gravity as Linder:2010py

$\displaystyle S=\int
d^{4}x|e|\left[\frac{F(T)}{2{\kappa}^{2}}+{{L}}_{\mathrm{matter}}\right]\,,$
(399)

$|e|=\det\left(e^{A}_{\mu}\right)=\sqrt{-g}$. Now, varying the action in Eq.
(399) with respect to the vierbein vector field $e_{A}^{\mu}$ one can obtain
the gravitational field equations Bengochea:2008gz

$\displaystyle\frac{1}{e}\partial_{\mu}\left(eS_{A}^{\
\mu\nu}\right)F^{\prime}-e_{A}^{\lambda}T^{\rho}_{\ \mu\lambda}S_{\rho}^{\
\nu\mu}F^{\prime}+S_{A}^{\
\mu\nu}\partial_{\mu}TF^{\prime\prime}+\frac{1}{4}e_{A}^{\nu}F=\frac{{\kappa}^{2}}{2}e_{A}^{\rho}{T^{(\mathrm{M})}}_{\rho}^{\
\nu}\,,$ (400)

where prime denotes the differentiation with respect to the torsion scalar $T$
and ${T^{(\mathrm{M})}}_{\rho}^{\ \nu}$ is the energy-momentum tensor
describing the entire matter sector.

Now, let us assume that the background manifold is described well by the
spatially flat FLRW Universe (3). Then, vierbein takes the
form$e^{A}_{\mu}=\left(1,a(t),a(t),a(t)\right)$ where $a(t)$ is the scale
factor of the FLRW Universe. One can quickly derive the line element of the
FLRW Universe (3). For this tetrad, using Eqs. (395), (397), and (398), one
can derive the torsion scalar $T=-6H^{2}$. In this flat FLRW background, one
can write down the gravitational equations

$\displaystyle H^{2}=$
$\displaystyle\,\frac{\kappa^{2}}{3}(\rho+\rho_{\mathrm{GDE}})=\frac{\kappa^{2}}{3}\rho_{\mathrm{eff}},$
(401) $\displaystyle\dot{H}=$
$\displaystyle\,-\frac{\kappa^{2}}{2}(\rho+p+\rho_{\mathrm{GDE}}+p_{\mathrm{GDE}})=-\frac{\kappa^{2}}{2}(\rho_{\mathrm{eff}}+p_{\mathrm{eff}})\,,$
(402)

where $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are the total energy
density and total pressure of the Universe, respectively, $\rho$ and $p$
denote the energy density and pressure of all the perfect fluids comprising
the entire matter sector where they satisfy the conservation law
$\dot{\rho}+3H(p+\rho)=0$. Further, $\rho_{\mathrm{GDE}}$ and
$p_{\mathrm{GDE}}$ are given by

$\displaystyle\rho_{\mathrm{GDE}}=\frac{1}{2{\kappa}^{2}}\left(-T-F+2TF^{\prime}\right)\,,\quad
p_{\mathrm{GDE}}=-\frac{1}{2{\kappa}^{2}}\left[\left(4-4F^{\prime}-8TF^{\prime\prime}\right)\dot{H}-T-F+2TF^{\prime}\right]\,.$
(403)

Notice that the effective fluid characterized by $\rho_{\mathrm{GDE}}$ and
$p_{\mathrm{GDE}}$ satisfies the continuity equation
$\dot{\rho}_{\mathrm{GDE}}+3H(p_{\mathrm{GDE}}+\rho_{\mathrm{GDE}})=0$. As one
can clearly notice that the choice of $F(T)$ highly influences the
gravitational equations. In the following we shall focus on the finite-time
future singularities that may appear in the context of modified teleparallel
gravity.

We consider, as in the previous sections, the following Hubble rate in (63)
and (218) Bamba:2012vg

$\displaystyle H\sim$
$\displaystyle\,\frac{h_{\mathrm{s}}}{\left(t_{\mathrm{s}}-t\right)^{\beta}}\,,\quad\mathrm{for}\,\,\,\beta>0\,,$
(404) $\displaystyle H\sim$
$\displaystyle\,H_{\mathrm{s}}+\frac{h_{\mathrm{s}}}{\left(t_{\mathrm{s}}-t\right)^{\beta}}\,,\quad\mathrm{for}\
\beta<-1\,,\ \mbox{and}\ -1<\beta<0\,.$ (405)

One can clearly notice that as $t\rightarrow t_{s}$, finite-time singularities
may appear depending on the nature of $\beta$. This can be illustrated as
follows: when $t\to t_{\mathrm{s}}$, for $\beta>0$, both $H\sim
h_{\mathrm{s}}\left(t_{\mathrm{s}}-t\right)^{-\beta}$ and $\dot{H}\sim\beta
h_{\mathrm{s}}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta+1\right)}$ diverge
to infinity; for $-1<\beta<0$, $H$ is finite, but $\dot{H}$ becomes infinity;
for $\beta<-1$, but $\beta$ is not any integer, both $H$ and $\dot{H}$ are
finite, but the higher derivatives of $H$ can diverge to infinity.

From Eq. (404) one can derive the evolution of the scale factor as

$\displaystyle a\sim$
$\displaystyle\,a_{\mathrm{s}}\exp\left[\frac{h_{\mathrm{s}}}{\beta-1}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}\right]\,,\quad\mathrm{for}\,\,\,0<\beta<1,\;\mbox{and}\,\,\,\beta>1\,,$
(406) $\displaystyle a\sim$
$\displaystyle\,a_{\mathrm{c}}\left(\frac{t_{\mathrm{s}}-t_{\mathrm{c}}}{t_{\mathrm{s}}-t}\right)^{h_{\mathrm{s}}},\quad\mathrm{for}\,\,\,\beta=1\,,$
(407)

where $a_{\mathrm{s}}$ and $a_{\mathrm{c}}$ are some positive constants. Now,
from Eq. (406) one can see that, when $t\to t_{\mathrm{s}}$, for $\beta\geq
1$, $a\to\infty$, whilst for $\beta<0$ and $0<\beta<1$, $a\to a_{\mathrm{s}}$.
Additionally, one can further notice that from
$\rho_{\mathrm{eff}}=3H^{2}/\kappa^{2}$ (Eq. (401)) and (404) that, for
$\beta>0$, $H\to\infty$ and as a result,
$\rho_{\mathrm{eff}}=3H^{2}/\kappa^{2}\to\infty$, while on the other hand, for
$\beta<0$, $H$ asymptotically becomes finite and $\rho_{\mathrm{eff}}$
asymptotically approaches to a finite constant value $\rho_{\mathrm{s}}$.
Furthermore, from $\dot{H}\sim\beta
h_{\mathrm{s}}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta+1\right)}$ and Eq.
(405) one can see that for $\beta>-1$, $\dot{H}\to\infty$ and as a
consequence,
$P_{\mathrm{eff}}=-\left(2\dot{H}+3H^{2}\right)/\kappa^{2}\to\infty$. While
for $\beta<-1$, but $\beta$ is not any integer, $a$, $\rho_{\mathrm{eff}}$ and
$P_{\mathrm{eff}}$ are finite since both $H$ and $\dot{H}$ are finite,
however, the higher derivatives of $H$ diverges.

Thus, we see that for the Hubble parameter described in Eqs. (404) and (405)
the $F(T)$ gravity models can encounter with the finite-time singularities for
the following values of $\beta$:

  * •

For $\beta\geq 1$, Type I (“Big Rip”) singularity appears.

  * •

For $-1<\beta<0$, Type II (“sudden”) singularity appears.

  * •

For $0<\beta<1$, Type III singularity appears.

  * •

For $\beta<-1$, but $\beta$ is not any integer, Type IV singularity appears.

Table 4: The table summarizes the emergence of the finite-time future singularities for different values of $\beta$ in the expressions of the Hubble parameter given in Eqs. (404) and (405) along with the behavior of $H$ and $\dot{H}$ in the limit of $t\to t_{\mathrm{s}}$. Table courtesy Ref. Bamba:2012vg .  Region of $\beta~{}(\neq 0,\,-1)$ | Type of the Singularity |  $H$ ($t\to t_{\mathrm{s}}$) |  $\dot{H}$ ($t\to t_{\mathrm{s}}$)
---|---|---|---
$\beta\geq 1$ |  Type I | $H\to\infty$ | $\dot{H}\to\infty$
$-1<\beta<0$ |  Type II | $H\to H_{\mathrm{s}}$ | $\dot{H}\to\infty$
$0<\beta<1$ |  Type III | $H\to\infty$ | $\dot{H}\to\infty$
$\beta<-1$, but $\beta$ is not any integer |  Type IV | $H\to H_{\mathrm{s}}$ | $\dot{H}\to 0$
|  |  | (and higher order
|  |  | derivatives of $H$ diverge)

In TABLE 4, we summarize the appearance of finite-time singularities out of
the expressions of the Hubble parameter as given in Eqs. (404) and (405).
Thus, one can see that the above expressions of the Hubble parameter can be
used to reconstruct the $F(T)$ gravity models which may allow the finite-time
singularities. We focus on the finite-time singularities when the geometrical
sector dominates over the matter sector, that means when
$\rho_{\mathrm{eff}}\simeq\rho_{\mathrm{GDE}}$ and $p_{\mathrm{eff}}\simeq
p_{\mathrm{GDE}}$, that means
$w_{\mathrm{eff}}=p_{\mathrm{eff}}/\rho_{\mathrm{eff}}\simeq p_{\rm
GDE}/\rho_{\mathrm{GDE}}$. Following this, one can have

$\displaystyle w_{\mathrm{eff}}\simeq
w_{\mathrm{GDE}}=-\frac{\left(4-4F^{\prime}-8TF^{\prime\prime}\right)\dot{H}-T-F+2TF^{\prime}}{-T-F+2TF^{\prime}}\,.$
(408)

On the other hand, from Eq. (403) one can write

$\displaystyle p_{\mathrm{GDE}}=-\rho_{\mathrm{GDE}}+I(H,\dot{H})\,,$ (409)

where
$I(H,\dot{H})=-\frac{1}{\kappa^{2}}\left[2-2F^{\prime}-4TF^{\prime\prime}\right]\dot{H}$.
From Eq.(402) one can quickly write
$p_{\mathrm{GDE}}=-\rho_{\mathrm{GDE}}-2\dot{H}/\kappa^{2}$, where we have
using that the geometrical sector dominates over the matter one. Now,
comparing this with (409), we obtain the following differential equation

$\displaystyle\dot{H}+\frac{\kappa^{2}}{2}I(H,\dot{H})=0\quad\Longrightarrow\quad\dot{H}\left[F^{\prime}+2TF^{\prime\prime}\right]=0\,.$
(410)

In the following we shall illustrate the the appearance of the finite-time
singularities in a specific $F(T)$ gravity model having the power law form:

$\displaystyle F(T)=AT^{\alpha}\,,$ (411)

where $A~{}(\neq 0)$ and $\alpha~{}(\neq 0)$ are constants. With this choice
of $F(T)$, Eq. (410) becomes

$\displaystyle\left(2\alpha-1\right)H^{2\left(\alpha-1\right)}=0\,.$ (412)

Notice that for both Eqs. (404) and (405), $\dot{H}\neq 0$. In the limit of
$t\rightarrow t_{s}$, Eq. (412) needs to be satisfied. Now, from Eqs. (404)
and (405), one can see that for $\beta>0$ [here two different cases may arise
in this way: (i) $\beta\geq 1$ (i.e., Type I singularity) and $0<\beta<1$
(Type III singularity)], $\alpha<1$, so that Eq. (412) is asymptotically
satisfied. While on the other hand, for $\beta<0$ [here one can similarly
explore two different cases: $-1<\beta<0$ (Type II singularity and $\beta<-1$
(Type IV singularity)], $\alpha=1/2$, in which Eq. (412) is always satisfied.
At this point it is very crucial to mention that the conditions: (i) for
$\beta>0$ (Type I singularity corresponds to $\beta\geq 1$ and Type III
singularity corresponds to $0<\beta<1$), $\alpha<0$, while (ii) for $\beta<0$
(Type II singularity corresponds to $-1<\beta<0$ and Type IV singularity
corresponds to $\beta<-1$), $\alpha=1/2$ are “necessary conditions” for the
appearance of the finite-time future singularities but they are not the
sufficient conditions. Because if $\alpha<0$, then the Type I singularity with
$\beta\geq 1$ appears rather than the Type III singularity with $0<\beta<1$,
since in the limit of $t\to t_{\mathrm{s}}$, both $H$ and $\dot{H}$ with
$\beta\geq 1$ diverge more rapidly than those with $0<\beta<1$. This is
because of the absolute value of the power $\beta$, namely, the absolute value
of $\beta$ for the Type I singularity ($\beta\geq 1$) is larger than the
absolute value of $\beta$ for the Type III singularity ($0<\beta<1$). As a
consequence, the Type I singularity is realized faster than the Type III
singularity, and this results in the appearance of the Type I singularity. In
a similar fashion, if $\alpha=1/2$, then Type IV singularity with $\beta<-1$
appears compared to the Type II singularity with $-1<\beta<0$ because again in
the limit of $t\to t_{\mathrm{s}}$, $H\to H_{\mathrm{s}}$ and $\dot{H}\to 0$
with $\beta<-1$ are realized more quickly than $H\to H_{\mathrm{s}}$ and
$\dot{H}\to\infty$ with $-1<\beta<0$. This is again because of the absolute
value of the power $\beta$, namely, the absolute value of $\beta$ for the Type
IV singularity ($\beta<-1$) is larger than the absolute value of $\beta$ for
the Type II singularity ($-1<\beta<0$). Hence, the Type IV singularity appears
faster than the Type II singularity, and Type IV singularity appears as a
result.

Finally, we note that Type V (“$w$”) singularity can also appear in this
model. In the Type V singularity, a specific choice of the scale factor can be
considered Dabrowski:2009kg

$\displaystyle a(t)=$
$\displaystyle\,a_{\mathrm{s}}\left(1-\frac{3\sigma}{2}\left\\{\frac{n-1}{n-\left[2/\left(3\sigma\right)\right]}\right\\}^{n-1}\right)^{-1}+\frac{1-2/\left(3\sigma\right)}{n-2/\left(3\sigma\right)}\times
na_{\mathrm{s}}\left(1-\frac{2}{3\sigma}\left\\{\frac{n-\left[2/\left(3\sigma\right)\right]}{n-1}\right\\}^{n-1}\right)^{-1}\left(\frac{t}{t_{\mathrm{s}}}\right)^{2/\left(3\sigma\right)}$
$\displaystyle\,+a_{\mathrm{s}}\left(\frac{3\sigma}{2}\left\\{\frac{n-1}{n-\left[2/\left(3\sigma\right)\right]}\right\\}^{n-1}-1\right)^{-1}\left[1-\frac{1-2/\left(3\sigma\right)}{n-2/\left(3\sigma\right)}\frac{t}{t_{\mathrm{s}}}\right]^{n}\,,$
(413)

where $n$ $\sigma$ are any arbitrary real numbers. In the limit of $t\to
t_{\mathrm{s}}$, we have $H(t\to t_{\mathrm{s}})\to 0$ and $\dot{H}(t\to
t_{\mathrm{s}})\to 0$. The effective EoS
$w_{\mathrm{eff}}=\left(1/3\right)\left(2q-1\right)\to\infty$ where
$q\equiv-\ddot{a}a/\dot{a}^{2}$ is the deceleration parameter of the Universe.
Thus, for the power law model (411) with $\alpha>1$, Eq. (412) can be
satisfied asymptotically. Hence, the Type V (“$w$”) singularity can appear in
this model.

###  IV.9 Non-local Gravity

Along with the known modified gravitational theories presented in the earlier
sections, an addition to the class of modified gravity theories is the non-
local gravitational theory, i.e., non-local additions to GR Deser:2007jk .
These non-local corrections naturally arise due to quantum effects and with
such new corrections, it is also possible to explain the accelerating
expansion of the Universe Deser:2007jk . In this section we shall describe how
the finite-time future singularities may appear in this gravitational theory.

The action of the non-local gravity is given by Bamba:2012ky

$\displaystyle S=\int
d^{4}x\sqrt{-g}\Bigg{[}\frac{1}{2\kappa^{2}}\left\\{R\left(1+f(\Box^{-1}R)\right)-2\Lambda\right\\}+{L}_{\mathrm{matter}}\left(Q;g\right)\Bigg{]}\,,$
(414)

where $R$ is the Ricci scalar, $g$ is the determinant of the metric tensor
$g_{\mu\nu}$, $\Box\equiv g^{\mu\nu}{\nabla}_{\mu}{\nabla}_{\nu}$ with
${\nabla}_{\mu}$ being the covariant derivative is the covariant d’Alembertian
for a scalar field, and $f$ is any arbitrary function of $\Box^{-1}R$,
$\Lambda$ is a cosmological constant, and
${L}_{\mathrm{matter}}\left(Q;g\right)$ denotes the matter Lagrangian in which
$Q$ stands for the matter fields.

Introducing two scalar fields $\eta$ and $\xi$, we can rewrite the above
action (414) as

$\displaystyle S=$ $\displaystyle\,\int
d^{4}x\sqrt{-g}\Bigg{[}\frac{1}{2\kappa^{2}}\Bigg{\\{}R\left(1+f(\eta)\right)+\xi\left(\Box\eta-R\right)-2\Lambda\Bigg{\\}}+{L}_{\mathrm{matter}}\Bigg{]}$
$\displaystyle=$ $\displaystyle\,\int
d^{4}x\sqrt{-g}\Bigg{[}\frac{1}{2\kappa^{2}}\Bigg{\\{}R\left(1+f(\eta)\right)-\partial_{\mu}\xi\partial^{\mu}\eta-\xi
R-2\Lambda\Bigg{\\}}+{L}_{\mathrm{matter}}\Bigg{]}\,.$ (415)

In the context of a spatially flat FLRW metric (3), one can write down the
gravitational field equations as

$\displaystyle 0=$
$\displaystyle\,-3H^{2}\left(1+f(\eta)-\xi\right)+\frac{1}{2}\dot{\xi}\dot{\eta}-3H\left(f^{\prime}(\eta)\dot{\eta}-\dot{\xi}\right)+\Lambda+\kappa^{2}\rho\,,$
(416) $\displaystyle 0=$
$\displaystyle\,\left(2\dot{H}+3H^{2}\right)\left(1+f(\eta)-\xi\right)+\frac{1}{2}\dot{\xi}\dot{\eta}+\left(\frac{d^{2}}{dt^{2}}+2H\frac{d}{dt}\right)\left(f(\eta)-\xi\right)-\Lambda+\kappa^{2}p\,,$
(417)

where we have considered that the scalar fields $\eta$ and $\xi$ depend only
on time Bamba:2012ky . Note that here $\rho$ and $p$ are respectively the
energy density and pressure of the matter sector. One can further write down
the equations of motion for the scalar fields $\eta$ and $\xi$ as

$\displaystyle 0=$
$\displaystyle\,\ddot{\eta}+3H\dot{\eta}+6\dot{H}+12H^{2}\,,$ (418)
$\displaystyle 0=$
$\displaystyle\,\ddot{\xi}+3H\dot{\xi}-\left(6\dot{H}+12H^{2}\right)f^{\prime}(\eta)\,,$
(419)

where we have used $R=6\dot{H}+12H^{2}$. Let us now focus on the emergence of
the finite-time singularities in this gravitational theory.

####  IV.9.1 Finite-time Singularities

We start with the following expression of the Hubble parameter in (63) or
(114) (see Bamba:2012ky ). The scale factor is also given in (63). Now with
the use of $\ddot{\eta}+3H\dot{\eta}=a^{-3}d\left(a^{3}\dot{\eta}\right)/dt$
and Eq. (418), one can express $\eta$ as

$\displaystyle\eta=-\int^{t}\frac{1}{a^{3}}\left(\int^{\bar{t}}Ra^{3}d\bar{t}\right)dt\,.$
(420)

We note that in the limit $t\to t_{s}$: for $\beta>1$, $\dot{H}\ll H^{2}$, and
hence, $R\sim 12H^{2}$, whilst for $-1<\beta<0$ and $0<\beta<1$, $\dot{H}\gg
H^{2}$, and hence, $R\sim 6\dot{H}$. By applying these relations to Eq. (420)
and then taking the leading term in terms of $\left(t_{\mathrm{s}}-t\right)$,
one can express $\eta$ for different regions of $\beta$ as follows
Bamba:2012ky :

$\displaystyle\eta\sim$
$\displaystyle\,-\frac{4h_{\mathrm{s}}}{\beta-1}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}+{\eta}_{\mathrm{c}}\quad(\beta>1)\,,$
(421) $\displaystyle\eta\sim$
$\displaystyle\,-\frac{6h_{\mathrm{s}}}{\beta-1}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}+{\eta}_{\mathrm{c}}\quad(-1<\beta<0\,,\,0<\beta<1)\,,$
(422)

where ${\eta}_{\mathrm{c}}$ is an integration constant. As discussed in
Bamba:2012ky , for ${\eta}_{\mathrm{c}}=0$, the finite-time future
singularities depicted by the Hubble parameter in Eq. (63) or (114) do not
appear, however, for ${\eta}_{\mathrm{c}}\neq 0$, finite time future
singularities can occur.

For ${\eta}_{\mathrm{c}}\neq 0$, note that, if the power in terms of
$\left(t_{\mathrm{s}}-t\right)$ is negative, i.e. $-\left(\beta-1\right)<0$,
then the first term proportional to
$\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}$ is the leading one. On
the other hand, if the power in terms of $\left(t_{\mathrm{s}}-t\right)$ is
positive, i.e. $-\left(\beta-1\right)>0$, then the second constant term
becomes the leading one. Thus, for $\beta>1$, the first term is the leading
one, i.e., $\eta\propto\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}$,
while for $-1<\beta<0$ and $0<\beta<1$, the second term is the leading one,
i.e., $\eta\sim{\eta}_{\mathrm{c}}$. For $\beta=1$ in Eq. (63) or (114), it
follows from Eq. (418) that,

$\displaystyle\eta\sim
6h_{\mathrm{s}}\left[\left(1+2h_{\mathrm{s}}\right)/\left(1+3h_{\mathrm{s}}\right)\right]\ln\left(t_{\mathrm{s}}-t\right)+{\eta}_{\mathrm{c}}.$
(423)

We continue by taking a particular form of $f(\eta)$ as Bamba:2012ky

$f(\eta)=f_{\mathrm{s}}\eta^{\sigma}\,,$ (424)

where $f_{\mathrm{s}}$ and $\sigma$ are non-zero constants. Now, using
$\ddot{\xi}+3H\dot{\xi}=a^{-3}d\left(a^{3}\dot{\xi}\right)/dt$ and Eq. (419),
$\xi$ can be expressed as

$\xi=\int^{t}\frac{1}{a^{3}}\left(\int^{\bar{t}}\frac{df(\eta)}{d\eta}Ra^{3}d\bar{t}\right)dt\,.$
(425)

Now we apply the conditions $R\sim 12H^{2}$ (for $\beta>1$) and $R\sim
6\dot{H}$ (for $\beta<1$) to Eq. (425) and taking into account the leading
term in terms of $\left(t_{\mathrm{s}}-t\right)$, we obtain Bamba:2012ky

$\displaystyle\xi\sim-
f_{\mathrm{s}}\left(-\frac{4h_{\mathrm{s}}}{\beta-1}\right)^{\sigma}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)\sigma}+{\xi}_{\mathrm{c}}~{},\quad(\beta>1)\,,$
(426)
$\displaystyle\xi\sim\frac{6f_{\mathrm{s}}h_{\mathrm{s}}\sigma{\eta}_{\mathrm{c}}^{\sigma-1}}{\beta-1}\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}+{\xi}_{\mathrm{c}}~{},\quad(-1<\beta<0\,~{},\,0<\beta<1)\,,$
(427)

where ${\xi}_{\mathrm{c}}$ is the constant of integration. Now, if the power
in terms of $\left(t_{\mathrm{s}}-t\right)$ is negative, i.e.
$-\left(\beta-1\right)\sigma<0$, the first term proportional to
$\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)\sigma}$ is the leading
one. On the other hand, if the power in terms of
$\left(t_{\mathrm{s}}-t\right)$ is positive, i.e.
$-\left(\beta-1\right)\sigma>0$, then the second constant term becomes the
leading one. Therefore, we see that for $\beta>1$, $\sigma>0$,
$\xi\propto\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)\sigma}$, while
for ($\beta>1$, $\sigma<0$) and ($-1<\beta<0\,,\,0<\beta<1$),
$\xi\sim{\xi}_{\mathrm{c}}$. Therefore one arrives at the following three
cases Bamba:2012ky :

  * •

For $\beta>1$, $\sigma>0$:
$\eta\propto\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}$ and
$\xi\propto\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)\sigma}$.

  * •

For $\beta>1$, $\sigma<0$:
$\eta\propto\left(t_{\mathrm{s}}-t\right)^{-\left(\beta-1\right)}$ and
$\xi\sim{\xi}_{\mathrm{c}}$.

  * •

For $-1<\beta<0\,,\,0<\beta<1$: $\eta\sim{\eta}_{\mathrm{c}}$ and
$\xi\sim{\xi}_{\mathrm{c}}$.

Now, having all of these, following Bamba:2012ky  one can further investigate
the possibility of the finite time future singularities in this gravitational
theory when the Hubble rate is given in Eq. (63) or (114) and the results are
as follows

  * •

For $\sigma<0$, with ${\eta}_{\mathrm{c}}\neq 0$ and ${\xi}_{\mathrm{c}}=1$,
Type I (Big Rip) singularity can occur for $\beta>1$.

  * •

If ${\eta}_{\mathrm{c}}\neq 0$, and
$f_{\mathrm{s}}{\eta}_{\mathrm{c}}^{\sigma-1}\left(6\sigma-{\eta}_{\mathrm{c}}\right)+{\xi}_{\mathrm{c}}-1=0$,
then Type II (sudden) singularity can occur for $-1<\beta<0$.

  * •

If ${\eta}_{\mathrm{c}}\neq 0$, and
$f_{\mathrm{s}}{\eta}_{\mathrm{c}}^{\sigma-1}\left(6\sigma-{\eta}_{\mathrm{c}}\right)+{\xi}_{\mathrm{c}}-1=0$,
there Type III singularity can occur for $0<\beta<1$.

###  IV.10 Non-minimal Maxwell-Einstein Gravity

In this section, we discuss the finite-time future singularities appearing in
the non-minimal Maxwell-Einstein gravity with general gravitational coupling.
To begin with, we consider the following action Bamba:2008ja ; Bamba:2008ut :

$\displaystyle S_{\mathrm{GR}}=\int
d^{4}x\sqrt{-g}\left[\frac{1}{2\kappa^{2}}R-\frac{1}{4}I(R)F_{\mu\nu}F^{\mu\nu}\right]\,,$
(428)

where $I(R)$ is any function of the Ricci scalar,
$F_{\mu\nu}={\partial}_{\mu}A_{\nu}-{\partial}_{\nu}A_{\mu}$ is the
electromagnetic field-strength tensor in which $A_{\mu}$ denotes the $U(1)$
gauge field. It is well-known that the coupling between the scalar curvature
and the Lagrangian of the electromagnetic field, as shown in action (428),
arises in curved spacetime due to one-loop vacuum-polarization effects in
Quantum Electrodynamics Drummond:1979pp .

Now. taking the variations of the action in Eq. (428) with respect to the
metric $g_{\mu\nu}$ and the $U(1)$ gauge field $A_{\mu}$, one can find the
gravitational field equations and the equation of motion of $A_{\mu}$,
respectively, as Bamba:2008ja ; Bamba:2008ut

$\displaystyle
R_{\mu\nu}-\frac{1}{2}g_{\mu\nu}R=\kappa^{2}T^{(\mathrm{EM})}_{\mu\nu}\,,$
(429)

and

$\displaystyle-\frac{1}{\sqrt{-g}}{\partial}_{\mu}\left(\sqrt{-g}I(R)F^{\mu\nu}\right)=0\,,$
(430)

where $T^{(\mathrm{EM})}_{\mu\nu}$ in Eq. (429) denotes the contribution to
the energy-momentum tensor from the electromagnetic field:

$\displaystyle T^{(\mathrm{EM})}_{\mu\nu}=$
$\displaystyle\,I(R)\left(g^{\alpha\beta}F_{\mu\beta}F_{\nu\alpha}-\frac{1}{4}g_{\mu\nu}F_{\alpha\beta}F^{\alpha\beta}\right)$
$\displaystyle\,+\frac{1}{2}\left[I^{\prime}(R)F_{\alpha\beta}F^{\alpha\beta}R_{\mu\nu}+g_{\mu\nu}\Box\left(I^{\prime}(R)F_{\alpha\beta}F^{\alpha\beta}\right)-{\nabla}_{\mu}{\nabla}_{\nu}\left(I^{\prime}(R)F_{\alpha\beta}F^{\alpha\beta}\right)\right]\,,$
(431)

in which the prime refers to the derivative with respect to $R$;
$\nabla_{\mu}$ denotes the covariant derivative operator associated with
$g_{\mu\nu}$, and $\equiv g_{\mu\nu}\nabla_{\mu}\nabla_{\nu}$ is the covariant
d’Alembertian for a scalar field. To understand the dynamics of the Universe
within this gravitational theory, we consider the flat FLRW metric and
following Ref. Bamba:2008ut  we consider the case where only the effects of
magnetic fields are present that means the effects of electric fields are
negligible. Additionally, we further consider that only one component of the
magnetic field $\vec{B}$ is non-zero, that means other two components are
zero. Thus, it follows from $\mathrm{div}\vec{B}=0$ that the off-diagonal
components of the last term of the r.h.s. of Eq. (431) for
$T^{(\mathrm{EM})}_{\mu\nu}$, i.e.,
${\nabla}_{\mu}{\nabla}_{\nu}\left[I^{\prime}(R)F_{\alpha\beta}F^{\alpha\beta}\right]$
are zero. Hence, all of the off-diagonal components of
$T^{(\mathrm{EM})}_{\mu\nu}$ are zero. Furthermore, since we have assumed that
we only have the magnetic fields as background quantities at the
$0^{\mathrm{th}}$ order, therefore, the magnetic fields remain independent of
the space components $\vec{x}$.

In the FLRW spacetime, the equation of motion for the $U(1)$ gauge field in
the Coulomb gauge, ${\partial}^{j}A_{j}(t,\vec{x})=0$, and the case of
$A_{0}(t,\vec{x})=0$, becomes

$\displaystyle{\ddot{A}}_{i}(t,\vec{x})+\left(H+\frac{\dot{I}}{I}\right){\dot{A}}(t,\vec{x})-\frac{1}{a^{2}}\mathop{\Delta}\limits^{(3)}\,A_{i}(t,\vec{x})=0\,,$
(432)

where $\mathop{\Delta}\limits^{(3)}={\partial}^{i}{\partial}_{i}$ is the flat
3-dimensional Laplacian. From Eq. (432) it follows that the Fourier mode
$A_{i}(k,t)$ satisfies the equation

$\displaystyle{\ddot{A}}_{i}(k,t)+\left(H+\frac{\dot{I}}{I}\right){\dot{A}}_{i}(k,t)+\frac{k^{2}}{a^{2}}A_{i}(k,t)=0\,,$
(433)

which can alternatively be expressed (by replacing the cosmic time $t$ in
terms of the conformal time $\eta$) as

$\displaystyle\frac{\partial^{2}A_{i}(k,\eta)}{\partial\eta^{2}}+\frac{1}{I(\eta)}\frac{dI(\eta)}{d\eta}\frac{\partial
A_{i}(k,\eta)}{\partial\eta}+k^{2}A_{i}(k,\eta)=0\,.$ (434)

We note that finding the exact solution of Eq. (434) for any arbitrary
coupling function $I(\eta)$ is not possible while with the use of the
Wentzel–Kramers–Brillouin (WKB) approximation on subhorizon scales and the
long-wavelength approximation on superhorizon scales, and finally matching
these solutions at the horizon crossing Bamba:2006ga ; Bamba:2007sx , an
approximate solution can be found as

$\displaystyle\left|A_{i}(k,\eta)\right|^{2}=|\bar{C}(k)|^{2}=\frac{1}{2kI(\eta_{k})}\left|1-\left[\frac{1}{2}\frac{1}{kI(\eta_{k})}\frac{dI(\eta_{k})}{d\eta}+i\right]k\int_{\eta_{k}}^{{\eta}_{\mathrm{f}}}\frac{I(\eta_{k})}{I\left(\widetilde{\widetilde{\eta}}\right)}d\widetilde{\widetilde{\eta}}\,\right|^{2}\,,$
(435)

where $\eta_{k}$ and ${\eta}_{\mathrm{f}}$ denotes the conformal time at the
horizon-crossing and the conformal time at the end of inflation, respectively.
From Eq. (435),one may obtain the amplitude of the proper magnetic fields in
the position space

$\displaystyle|{B_{i}}^{(\mathrm{proper})}(t)|^{2}=\frac{k|\bar{C}(k)|^{2}}{\pi^{2}}\frac{k^{4}}{a^{4}}\,,$
(436)

on a comoving scale $L=2\pi/k$. From Eq. (436) one can see that the proper
magnetic fields evolve as
$|{B_{i}}^{(\mathrm{proper})}(t)|^{2}=|\bar{B}|^{2}/a^{4}$, where $|\bar{B}|$
is a constant Bamba:2008ut . This means that the impact of the coupling
function $I$ on the value of the proper magnetic fields exists only during the
inflationary stage. On the other hand, the conductivity of the Universe,
${\sigma}_{\mathrm{c}}$, is extremely small during the period of inflation
since at that time only a few charged particles exist. After the reheating
stage, a number of charged particles are produced, so that the conductivity
immediately jumps to a large value ${\sigma}_{\mathrm{c}}\gg H$. Consequently,
for a large enough conductivity at the reheating stage, the proper magnetic
fields evolve in proportion to $a^{-2}(t)$ in the radiation-dominated stage
and in the subsequent matter-dominated stage Ratra:1991bn .

Now, from (431), one can express the effective energy density and effective
pressure of the Universe as follows

$\displaystyle\rho_{\mathrm{eff}}=$
$\displaystyle\,\left\\{\frac{I(R)}{2}+3\left[-\left(5H^{2}+\dot{H}\right)I^{\prime}(R)+6H\left(4H\dot{H}+\ddot{H}\right)I^{\prime\prime}(R)\right]\right\\}\frac{|\bar{B}|^{2}}{a^{4}}\,,$
(437) $\displaystyle p_{\mathrm{eff}}=$
$\displaystyle\,\biggl{[}-\frac{I(R)}{6}+\left(-H^{2}+5\dot{H}\right)I^{\prime}(R)-6\left(-20H^{2}\dot{H}+4\dot{H}^{2}-H\ddot{H}+\dddot{H}\right)I^{\prime\prime}(R)$
$\displaystyle\,-36\left(4H\dot{H}+\ddot{H}\right)^{2}I^{\prime\prime\prime}(R)\biggr{]}\frac{|\bar{B}|^{2}}{a^{4}}\,,$
(438)

where the following relations under the assumption of negligible electric
fields have been used:
$g^{\alpha\beta}F_{0\beta}F_{0\alpha}-\left(1/4\right)g_{00}F_{\alpha\beta}F^{\alpha\beta}=|{B_{i}}^{(\mathrm{proper})}(t)|^{2}/2,$
and $F_{\alpha\beta}F^{\alpha\beta}=2|{B_{i}}^{(\mathrm{proper})}(t)|^{2}$.

We recall that in terms of the effective energy density and pressure, the
Friedmann and Raychaudhuri equations for this gravitational theory can be
written as

$\displaystyle\frac{3}{\kappa^{2}}H^{2}\ =$
$\displaystyle\,\rho_{\mathrm{eff}}\,,$ (439)
$\displaystyle-\frac{1}{\kappa^{2}}\left(2\dot{H}+3H^{2}\right)=$
$\displaystyle\,p_{\mathrm{eff}}\,,$ (440)

where $\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are given in Eqs. (437) and
(438), respectively. Since the Friedmann and Raychaudhuri equations (439) and
(440) includes $I(R)$, therefore, depending on the functional forms of $I(R)$,
the Maxwell-Einstein gravity may allow finite-time future singularities. In
the following we discuss the possible finite-time future singularities.

(i) Big Rip singularity

We consider the Hubble rate in (185). For this Hubble rate, one can derive the
scale factor and the Ricci scalar as in (64) and (363). We assume that for the
large curvature, $I(R)$ behaves as

$\displaystyle I(R)\sim I_{s}R^{\alpha},$ (441)

where $I_{s}$ and $\alpha$ are constants. Thus, with the choice of $I(R)$ in
Eq. (441), $\rho_{\mathrm{eff}}$ in Eq. (437) behaves as
$\left(t_{s}-t\right)^{-2\alpha+4h_{s}}$ while from the l.h.s. of the
Friedmann equation, i.e., Eq. (440) $\rho_{\mathrm{eff}}$ behaves as
$\left(t_{s}-t\right)^{-2}$. From the consistency in $\rho_{\mathrm{eff}}$, we
have,

$\displaystyle-2=-2\alpha+4h_{s},\quad\mbox{i.e., }\quad\alpha=1+2h_{s}\,.$
(442)

Now, from Eq. (440) one can write,

$\displaystyle\frac{3h_{s}^{2}}{\kappa^{2}}=$
$\displaystyle\,I_{s}\left(12h_{s}^{2}+6h_{s}\right)^{\alpha-2}$
$\displaystyle\,\times\left\\{\frac{\left(12h_{s}^{2}+6h_{s}\right)^{2}}{2}+3\left[-\alpha\left(12h_{s}^{2}+6h_{s}\right)\left(h_{s}+5h_{s}^{2}\right)+6\alpha\left(\alpha-1\right)h_{s}\left(2h_{s}+4h_{s}^{2}\right)\right]\right\\}\frac{|\bar{B}|^{2}}{a_{s}^{4}}~{},$
(443)

which can be simplified with the use of $\alpha$ from (442) as

$\displaystyle\frac{3h_{s}^{2}}{\kappa^{2}}=-\frac{I_{s}h_{s}\left(12h_{s}^{2}+6h_{s}\right)^{\alpha}|\bar{B}|^{2}}{2a_{s}^{4}}\,,$
(444)

which clearly infers that $I_{s}$ must be negative. This tells that $I(R)$ is
also negative and therefore the gauge field becomes ghost and the model
becomes physically inconsistent.

Thus, it follows from Eqs. (441) and (442) that the Big rip singularity as
described through the Hubble rate in Eq. (185) may emerge when for the large
curvature, $I(R)$ behaves as $R^{1+2h_{s}}$. However, for other choices of
$I(R)$, we do not get the Big rip singularity. We in fact note that for
$I(R)=I_{s}R^{\alpha}$, $H=h_{s}/\left(t_{s}-t\right)$ is an exact solution.

(ii) Other types of singularities

We now focus on other types of singularities that may arise for a more general
choice of $I(R)$. Thus, we consider the general case where the Hubble
parameter takes the form of (185). For this choice of the Hubble rate, one can
find the scale factor and the scalar curvature as follows

$\displaystyle
a=a_{s}\exp\left[\frac{h_{s}}{\beta-1}\left(t_{s}-t\right)^{-\left(\beta-1\right)}\right]\,,\quad
R=6h_{s}\left[\beta+2h_{s}\left(t_{s}-t\right)^{-\left(\beta-1\right)}\right]\left(t_{s}-t\right)^{-\left(\beta+1\right)}\,,$
(445)

where $a_{s}$ is a positive constant. Similar to the previous case, we assume
that in the large curvature, $I(R)$ behaves as in Eq. (441). We note that the
case with $\beta>1$ is unphysical because, for $\beta>1$, $a\to\infty$ in the
limit $t\rightarrow t_{s}$, and as a consequence, $\rho_{\mathrm{eff}}\to 0$
and $p_{\mathrm{eff}}\to 0$ since $\rho_{\mathrm{eff}}\propto a^{-4}$ and
$p_{\mathrm{eff}}\propto a^{-4}$. On the other hand, $H\to\infty$ as
$t\rightarrow t_{s}$. Thus, we see that for $\beta>1$, Eqs. (439) and (440)
are not satisfied leading to a unphysical case. Therefore, we restrict
ourselves to $\beta<1$ and explore the possible singularities as follows.

  * •

If $0<\beta<1$ and $\alpha>0$, then $\rho_{\mathrm{eff}}$ in Eq. (439) evolves
as
$\rho_{\mathrm{eff}}\propto\left(t_{s}-t\right)^{-\alpha\left(\beta+1\right)}$,
however, looking at Eq. (440), we see that $\rho_{\mathrm{eff}}$ evolves as
$\rho_{\mathrm{eff}}\propto\left(t_{s}-t\right)^{-2\beta}$. Therefore, from
the consistency, we must have,

$\displaystyle-2\beta=-\alpha\left(\beta+1\right)\ \Longrightarrow\
\beta=\frac{\alpha}{2-\alpha},\quad\mbox{or}\quad\alpha=\frac{2\beta}{\beta+1}\,.$
(446)

Now, from Eq. (440), one can find,

$\displaystyle\frac{3h_{s}^{2}}{\kappa^{2}}=-\frac{I_{s}\left(6h_{s}\beta\right)^{\alpha}\left(1-\beta\right)|\bar{B}|^{2}}{2a_{s}^{4}\left(\beta+1\right)}\,,$
(447)

where we have used Eq. (446). Notice that Eq. (447) demands $I_{s}$ to be
negative. Thus, in this case (i.e., $\alpha>0$ and $0<\beta<1$), we see that,
in the limit $t\to t_{s}$, $a\to a_{s}$, $R\to\infty$,
$\rho_{\mathrm{eff}}\to\infty$, and $|p_{\mathrm{eff}}|\to\infty$. Hence, the
Type III singularity appears.

  * •

If $-1<\beta<0$ and $\left(\beta-1\right)/\left(\beta+1\right)<\alpha<0$, then
in the limit $t\to t_{s}$, $a\to a_{s}$, $R\to\infty$, $\rho_{\mathrm{eff}}\to
0$, and $|p_{\mathrm{eff}}|\to\infty$. Although we can see that the final
value of $\rho_{\mathrm{eff}}$ is not a finite value but it vanishes
($\rho_{\mathrm{eff}}\rightarrow 0$ as $t\rightarrow t_{s}$), however, this
singularity can be considered to be the singularity of Type II. Because, when
$I$ and $H$ are given by $I=1+I_{s}R^{\alpha}$ and
$H=H_{s}+h_{s}\left(t_{s}-t\right)^{-\beta}$ (where $H_{s}$ is a constant),
respectively, then in the limit $t\rightarrow t_{s}$,
$\rho_{\mathrm{eff}}\to\rho_{s}$, a constant, where $\rho_{s}$ can be found
from Eqs. (437) and (439) as
$\rho_{s}=3H_{s}^{2}/\kappa^{2}=|\bar{B}|^{2}/\left(2_{s}^{4}\right)$.

  * •

If $\beta<-1$, then in the limit $t\to t_{s}$, $a\to a_{s}$ and $R\to 0$. Now,
for $\alpha\geq\left(\beta-1\right)/\left(\beta+1\right)$, in the limit $t\to
t_{s}$, $\rho_{\mathrm{eff}}\to 0$, $|p_{\mathrm{eff}}|\to 0$, and higher
derivatives of $H$ diverge. That means, Type IV singularity appears.

We further note here that for $-1<\beta<0$ and $\alpha>0$,
$\rho_{\mathrm{eff}}\to\infty$, but $H\to 0$ in the limit $t\rightarrow
t_{s}$, which shows that Eq. (439) is not satisfied. That means, the case with
$\alpha>0$ and $-1<\beta<0$ is not physical in the sense that the Friedmann
equation is inconsistent. We also find that If $\alpha<0$ and $0<\beta<1$,
then $\rho_{\mathrm{eff}}\to 0$, but $H\to\infty$. Hence, Eq. (439) is again
not satisfied. Further, if $-1<\beta<0$ and
$\alpha\leq\left(\beta-1\right)/\left(\beta+1\right)$, then we have unphysical
scenario because in this case, in the limit $t\to t_{s}$, $a\to a_{s}$,
$R\to\infty$, $\rho_{\mathrm{eff}}\to 0$, and $|p_{\mathrm{eff}}|\to 0$, but
$\dot{H}\to\infty$. Thus, we see that Eq. (440) is not satisfied.

We mention that if $I(R)$ is a constant666Let us note that the case $I(R)=1$
corresponds to the ordinary Maxwell theory., we do not find any finite-time
future singularity. We close this section with two examples of $I(R)$ as
follows. Suppose $I(R)$ takes the Hu–Sawicki form given by Hu:2007nk

$\displaystyle
I(R)=I_{\mathrm{HS}}(R)\equiv\frac{c_{1}\left(R/m^{2}\right)^{n}}{c_{2}\left(R/m^{2}\right)^{n}+1}\,,$
(448)

where $c_{1}$ and $c_{2}$ are dimensionless constants, $n$ is a positive
constant, and $m$ denotes a mass scale. We see that the model in Eq. (448)
satisfies the following conditions: (i)
$\lim_{R\to\infty}I_{\mathrm{HS}}(R)=c_{1}/c_{2}=\mbox{const}$ and (ii)
$\lim_{R\to 0}I_{\mathrm{HS}}(R)=0$. The second example of $I(R)$ having the
same features as with the Hu–Sawicki form of $I(R)$ in Eq. (448) takes the
form Nojiri:2007as

$\displaystyle
I(R)=I_{\mathrm{NO}}(R)\equiv\frac{\left[\left(R/M^{2}\right)-\left(R_{\mathrm{c}}/M^{2}\right)\right]^{2q+1}+{\left(R_{\mathrm{c}}/M^{2}\right)}^{2q+1}}{c_{3}+c_{4}\left\\{\left[\left(R/M^{2}\right)-\left(R_{\mathrm{c}}/M^{2}\right)\right]^{2q+1}+{\left(R_{\mathrm{c}}/M^{2}\right)}^{2q+1}\right\\}}\,,$
(449)

where $c_{3}$ and $c_{4}$ are dimensionless constants, $q$ is a positive
integer, $M$ denotes a mass scale, and $R_{\mathrm{c}}$ is current curvature.
One can check that $I(R)$ given in Eq. (449) satisfies the conditions: (i)
$\lim_{R\to\infty}I_{\mathrm{NO}}(R)=1/c_{4}=\mbox{const}$ and (ii)
$\lim_{R\to 0}I_{\mathrm{NO}}(R)=0.$ If $\beta<-1$ and $I(R)$ is given either
by $I_{\mathrm{HS}}(R)$ in Eq. (448) or $I_{\mathrm{NO}}(R)$ in Eq. (449), in
the limit $t\to t_{s}$, $a\to a_{s}$, $R\to 0$, $\rho_{\mathrm{eff}}\to 0$,
and $|p_{\mathrm{eff}}|\to 0$. In addition, higher derivatives of $H$ diverge.
Thus, the Type IV singularity emerges. Thus, we see that the Maxwell theory
which is non-minimally coupled to Einstein gravity may produce finite-time
future singularities depending on the form of non-minimal gravitational
coupling.

We conclude that the general conditions on $I(R)$ for which the finite-time
future singularities characterized by Eq. (185) cannot emerge are that in the
limit $t\to t_{s}$, $I(R)\to\bar{I}$ (where $\bar{I}(\neq 0)$ is a finite
constant), $I^{\prime}(R)\to 0$, $I^{\prime\prime}(R)\to 0$, and
$I^{\prime\prime\prime}(R)\to 0$.

J.1: Influence of non-minimal gravitational coupling on the finite-time future
singularities in modified gravity

In this section, we discuss the influence of non-minimal gravitational
coupling of the electromagnetic field on the $F(R)$ gravity theory
Bamba:2008ut . This is a generalization of the previous section IV.10. In this
case, the total energy density and pressure of the Universe are contributed
as, $\rho_{\mathrm{tot}}=\rho_{\mathrm{eff}}+\rho_{\mathrm{MG}}$ and
$p_{\mathrm{tot}}=p_{\mathrm{eff}}+p_{\mathrm{MG}}$, respectively, where
$\rho_{\mathrm{eff}}$ and $p_{\mathrm{eff}}$ are given by Eqs. (437) and
(438), respectively. And $\rho_{\mathrm{MG}}$ and $p_{\mathrm{MG}}$ can be
found from Eqs. (176) and (177) as follows

$\displaystyle\rho_{\mathrm{MG}}=$
$\displaystyle\,\frac{1}{\kappa^{2}}\left[-\frac{1}{2}f(R)+3\left(H^{2}+\dot{H}\right)f^{\prime}(R)-18\left(4H^{2}\dot{H}+H\ddot{H}\right)f^{\prime\prime}(R)\right]\,,$
(450) $\displaystyle p_{\mathrm{MG}}=$
$\displaystyle\,\frac{1}{\kappa^{2}}\left[\frac{1}{2}f(R)-\left(3H^{2}+\dot{H}\right)f^{\prime}(R)+6\left(8H^{2}\dot{H}+4{\dot{H}}^{2}+6H\ddot{H}+\dddot{H}\right)f^{\prime\prime}(R)+36\left(4H\dot{H}+\ddot{H}\right)^{2}f^{\prime\prime\prime}(R)\right]\,.$
(451)

In this case, it follows from Eqs. (111), (437), (438), (450) and (451) that
the flat Friedmann and Raychaudhuri equations can be given by Bamba:2008ut

$\displaystyle\frac{3}{\kappa^{2}}H^{2}=\rho_{\mathrm{tot}}=$
$\displaystyle\,\left\\{\frac{I(R)}{2}+3\left[-\left(5H^{2}+\dot{H}\right)I^{\prime}(R)+6H\left(4H\dot{H}+\ddot{H}\right)I^{\prime\prime}(R)\right]\right\\}\frac{|\bar{B}|^{2}}{a^{4}}$
$\displaystyle\,+\frac{1}{\kappa^{2}}\left[-\frac{1}{2}\left(F(R)-R\right)+3\left(H^{2}+\dot{H}\right)\left(F^{\prime}(R)-1\right)-18\left(4H^{2}\dot{H}+H\ddot{H}\right)F^{\prime\prime}(R)\right]\,,$
(452)
$\displaystyle-\frac{1}{\kappa^{2}}\left(2\dot{H}+3H^{2}\right)=p_{\mathrm{tot}}=$
$\displaystyle\,\biggl{[}-\frac{I(R)}{6}+\left(-H^{2}+5\dot{H}\right)I^{\prime}(R)-6\left(-20H^{2}\dot{H}+4\dot{H}^{2}-H\ddot{H}+\dddot{H}\right)I^{\prime\prime}(R)$
$\displaystyle\,-36\left(4H\dot{H}+\ddot{H}\right)^{2}I^{\prime\prime\prime}(R)\biggr{]}\frac{|\bar{B}|^{2}}{a^{4}}+\frac{1}{\kappa^{2}}\biggl{[}\frac{1}{2}\left(F(R)-R\right)-\left(3H^{2}+\dot{H}\right)\left(F^{\prime}(R)-1\right)$
$\displaystyle\,+6\left(8H^{2}\dot{H}+4{\dot{H}}^{2}+6H\ddot{H}+\dddot{H}\right)F^{\prime\prime}(R)+36\left(4H\dot{H}+\ddot{H}\right)^{2}F^{\prime\prime\prime}(R)\biggr{]}\,.$
(453)

Now, using Eqs. (452) and (453), we find

$\displaystyle\Bigg{[}\frac{I(R)}{3}+2\left(-8H^{2}+\dot{H}\right)I^{\prime}(R)+6\left(32H^{2}\dot{H}-4\dot{H}^{2}+4H\ddot{H}-\dddot{H}\right)I^{\prime\prime}(R)-36\left(4H\dot{H}+\ddot{H}\right)^{2}I^{\prime\prime\prime}(R)\Bigg{]}\frac{|\bar{B}|^{2}}{a^{4}}$
$\displaystyle+\frac{1}{\kappa^{2}}\left[2\dot{H}F^{\prime}(R)+6\left(-4H^{2}\dot{H}+4\dot{H}^{2}+3H\ddot{H}+\dddot{H}\right)F^{\prime\prime}(R)+36\left(4H\dot{H}+\ddot{H}\right)^{2}F^{\prime\prime\prime}(R)\right]=0\,.$
(454)

Recalling that in modified gravity with the ordinary Maxwell theory, if $F(R)$
behaves as in Eq. (216), i.e.,

$\displaystyle F(R)\sim F_{0}R+F_{1}R^{q}\,,\quad\mbox{where $q$ is a
constant}\,,$ (455)