# The little rip in classical and quantum $f(R)$ cosmology

Teodor Borislavov Vasilev<EMAIL_ADDRESS>Departamento de Física Teórica and
IPARCOS, Universidad Complutense de Madrid,
E-28040 Madrid, Spain.  Mariam Bouhmadi-López<EMAIL_ADDRESS>IKERBASQUE, Basque Foundation for Science, 48011, Bilbao, Spain  Department of
Physics, University of the Basque Country, UPV/EHU, P.O. Box 644, 48080
Bilbao, Spain.  Prado Martín-Moruno<EMAIL_ADDRESS>Departamento de Física
Teórica and IPARCOS, Universidad Complutense de Madrid,
E-28040 Madrid, Spain.

###### Abstract

The little rip is a cosmological abrupt event predicted by some phantom dark
energy models that could describe the future evolution of our Universe. This
event can be interpreted as a big rip singularity delayed indefinitely,
although in those models bounded structures will be destroyed in a finite
cosmic time in the future. In this work, we analyse the little rip cosmology
from a classical and quantum point of view within the scheme of alternative
metric $f(R)$ theories of gravity. The quantum analysis is performed in the
framework of $f(R)$ quantum geometrodynamics by means of the modified Wheeler-
DeWitt equation. In this context, we show that the DeWitt criterion can be
satisfied. Similar to what happens in general relativity, this result points
towards the avoidance of the little rip in $f(R)$ quantum cosmology.

##  I Introduction

Understanding the mechanisms involved in the accelerated expansion of our
Universe has become one of the greatest milestones in modern cosmology. Since
the moment this acceleration was first discovered [1, 2], a large number of
models have been proposed to explain the origin of this phase. In general
relativity (GR), the accelerated expansion is attributed to an exotic form of
energy with negative pressure that pushes the cosmos further into expansion.
This extra content is designated as dark energy (DE). However, the nature of
this energy is still unknown. Among the DE models that best fit the existing
observational data is the standard $\Lambda$CDM, where the cosmological
constant $\Lambda$ plays the role of DE. Nevertheless, the theoretically
expected and the observed values for $\Lambda$ differ in several orders of
magnitude (for a review on the topic, see Ref. [3]) and, therefore, one can
expect that this is just a useful effective model. Thus, various alternative
models for DE describing the late-time cosmology without a cosmological
constant have been proposed. Some examples are scalar fields in the form of
quintessence [4, 5], k-essence [6], DE of phantom nature [7], tachyonic matter
[8, 9], Chaplygin gas [10, 11] and holographic DE [12], among others. On the
other hand, it is also possible to describe the same cosmological evolution,
which in GR is attributed to DE, without the introduction of new material
content. In this fashion, alternative theories of gravity have acquired a
renewed interest, since they can provide a different framework for explaining
the accelerated expansion of our Universe. Accordingly, the late-time cosmic
expansion can be understood as a consequence of suitable modifications to the
underlining theory of gravity rather than due to some exotic matter content.
Some examples of this approach are $f(R)$ theories of gravity [13] (see also
Refs. [14, 15]), Gauss-Bonnet gravity [16, 17], $f(R,\mathcal{T})$ gravity
[18], where $\mathcal{T}$ stands for the trace of the energy momentum tensor,
$f(T)$ modified teleparallel gravity [19], being $T$ the torsion scalar, and
Horndeski theories [20] (see also, e.g., Ref. [21]). For a review on the state
of the art of DE cosmology see, for example, Refs. [22, 23] and references
therein.

From a practical point of view, whatever the origin of DE may be, it can be
described effectively by an equation of state parameter $w$. This parameter
represents the ratio between the pressure and energy density of DE, or the
effective energy density and pressure in the case of alternative theories of
gravity. Latest cosmological data suggests that $w$ lies on a very narrow band
around -1 [24, 25]. Therefore, the possibility of DE being of phantom nature
($w<-1$) is not observationally excluded. On the contrary, this is even
suggested by some data [26] and could be a solution to alleviate the $H_{0}$
tension (see, for example, Refs. [27, 28]). However, considering phantom DE
rises back ancient questions regarding the final fate of our Universe, since
phantom DE is well-known for possibly leading the cosmic expansion towards a
future singularity. All bounded structures in the universe and space-time
itself might be ripped apart at a final big rip (BR) singularity [29], where
the Hubble rate, its cosmic time derivative and the scale factor diverge at a
finite cosmic time. The cosmos could also reach an infinite expansion rate at
a finite size of the observable universe, freezing its evolution at a big
freeze singularity [30, 31]. Nonetheless, if $w<-1$ but $w\to-1$ sufficiently
rapidly, then, the occurrence of future singularities may be infinitely
delayed in time [32], in which case the singularity is called an abrupt cosmic
event. Indeed, the little rip (LR) abrupt event is just a big rip that would
take place at the infinite asymptotic future, although bounded structures will
be destroyed in a finite time from the present, see Refs. [32, 33]. Another
abrupt cosmic event appearing in phantom DE is the little sibling of the big
rip (LSBR). This event is characterized by the divergence of the Hubble rate
and the scale factor at an infinite future cosmic time, while the cosmic time
derivative of the Hubble rate remains finite [34]. For observational
constraints on these rip-like cosmic catastrophes see Refs. [35, 36]. (See
also Refs. [37, 38, 39] for other examples of cosmic singularities and Refs.
[40, 41, 42] for a detailed classification of DE singularities.) Nevertheless,
it is commonly believed that quantum gravity effects may smooth or avoid these
(classical) cosmic catastrophes, see Refs. [43, 44] (see also Refs. [40, 42,
45, 46, 47, 48]).

Previous works on quantum cosmology have shown that the aforementioned phantom
rip-like doomsdays, namely the BR, the LR and the LSBR, can be avoided due to
quantum effects rising up as the universe approaches the classical singularity
[43, 49, 50]. Since the background late-time cosmology can be equivalently
described in the context of GR or by alternative theories of gravity, it is
natural to wonder whether these singularities are still avoided in the quantum
realm for a different underlying theory of gravity. Following this line of
thought, it has already been established that the BR and the LSBR can be
avoided due to quantum effects in $f(R)$ cosmology [51, 52]. However, up-to
our knowledge, no results for the LR abrupt event in quantum $f(R)$ cosmology
have been published so far. Thus, in this work we address the missing quantum
fate of the LR in metric $f(R)$ theories of gravity. To do so, we consider a
group of $f(R)$ theories of gravity that predict the occurrence of the LR
abrupt event at the classical level. (For previous works on the LR cosmology
in alternative theories of gravity see Refs. [53, 54, 55].) The quantum
analysis is performed in the framework of $f(R)$ quantum geometrodynamics,
with the Wheeler-DeWitt equation [56] being adapted to the $f(R)$ gravity case
[57]. Consequently, we explore the possibility of avoiding the LR cosmic
doomsday in $f(R)$ quantum cosmology.

This paper is organized as follows. In Sec. II, we review some basic results
and observational constraints on the LR abrupt event in the framework of GR.
In Sec. III, we consider that the LR could take place when the description of
gravity is that provided by alternative $f(R)$ theories of gravity. For that
aim, we briefly introduce the reconstruction method for metric $f(R)$ gravity
in Sec. III.1. Thereafter, in Sec. III.2, we apply the reconstruction method
to the LR model reviewed in Sec. II. Thus, we obtain the group of metric
$f(R)$ theories of gravity predicting the LR abrupt cosmic event, recovering
the $f(R)$ function first presented in Ref. [54]. In Sec. IV, we study the
quantum fate of the LR predicted in one of the theories obtained in the past
section. We summarize previous results on the formulation of quantum cosmology
in the framework of general $f(R)$ quantum geometrodynamics in Sec. IV.1.
Next, in Secs. IV.2 and IV.3, we take different approaches to solve the
modified Wheeler-DeWitt equation and we analyse the avoidance of the
singularity by means of the DeWitt criterion. We summarize our results in Sec.
V. Finally, we discuss the validity of the approximations carried out to solve
the modified Wheeler-DeWitt equation in Appendix A.

##  II The little rip

Let us begin with a brief summary on the phenomenology of the LR abrupt event
in GR. Throughout this work we limit ourselves to homogeneous and isotropic
cosmological scenarios, which are described by the Friedmann-Lemaître-
Robertson-Walker (FLRW) metric given by the line element

$ds^{2}=-dt^{2}+a(t)^{2}ds^{2}_{3},$ (1)

where $a(t)$ stands for the scale factor, $ds_{3}^{2}$ represents the three-
dimensional metric and we have used the geometric unit system $8\pi G=c=1$.
Then, for the content of the universe represented by a perfect fluid with
energy density $\rho$ and pressure $p$, the Einstein’s field equations reduce
to the Friedman and Raychaudhuri equations

$\displaystyle\frac{\dot{a}^{2}}{a^{2}}$
$\displaystyle=H^{2}=\frac{1}{3}\rho-\frac{k}{a^{2}},$ (2)
$\displaystyle\frac{\ddot{a}}{a}$
$\displaystyle=\dot{H}+H^{2}=-\frac{1}{2}\left(p+\frac{\rho}{3}\right),$ (3)

respectively, where the dot represents the derivative with respect to the
cosmic time, $H$ denotes the Hubble rate and $k$ indicates the spatial
curvature of the Universe (not fixed at this point). In the common
interpretation of cosmological data within GR’s framework, this perfect fluid
is constituted with three different species: radiation, (dark and barionic)
matter (M) and DE. According to the latest observations [24, 25] the current
density parameters for M and DE are $\Omega_{M,0}\sim 0.315$ and
$\Omega_{DE,0}\sim 0.685$, respectively, whereas radiation is negligible at
the present time. Thus, DE is the dominant cosmic ingredient today.
Furthermore, it will be even more dominant in the future since matter tends to
dilute (faster). Thence, from a practical perspective, the contribution of
matter and radiation can be neglected when analysing the asymptotic future
evolution of these cosmological models. Consequently, we consider that $p$ and
$\rho$ in the above equations are those corresponding to DE.

The LR abrupt event can be understood as a BR singularity that has been
delayed indefinitely. This abrupt event was first discovered in Ref. [58] (see
also Refs. [40, 59]) and, thereafter, named as “little rip” in Ref. [32]. (See
also Ref. [60] where the LR was found in the context of brane cosmology.) For
the occurrence of this cosmic event, DE can be modelled by the following
equation of state (EoS)

$p=-\rho-A\sqrt{\rho},$ (4)

being $A$ a positive parameter. From the conservation of the energy momentum
tensor it follows that the energy density $\rho$ evolves with the scale factor
as

$\rho=\rho_{0}\left[1+\frac{3A}{2\sqrt{\rho_{0}}}\ln\left(\frac{a}{a_{0}}\right)\right]^{2},$
(5)

being $\rho_{0}$ and $a_{0}$ the current values of the DE density and the
scale factor, respectively. Then, combining the Friedman equation (2) and Eq.
(5), the time dependence of the scale factor reads

$\displaystyle a(t)=$ $\displaystyle
a_{0}\exp\Bigg{\\{}\frac{2\sqrt{\rho_{0}}}{3A}\Bigg{[}\left(1+\frac{3A}{2\sqrt{\rho_{0}}}\ln\frac{a_{\star}}{a_{0}}\right)$
$\displaystyle\times\exp\left(\frac{\sqrt{3}A}{2}\left(t-t_{\star}\right)\right)-1\Bigg{]}\Bigg{\\}},$
(6)

where we have denoted by $t_{\star}$ some arbitrary (future) moment in the
expansion history of the Universe from which we can safely assume DE is the
only content of the cosmos and $a_{\star}$ represents the corresponding scale
factor. The EoS parameter $w$ for the DE follows from Eqs. (4) and (5). That
is

$w=-1-\frac{A}{\sqrt{\rho_{0}}+\frac{3A}{2}\ln\left(\frac{a}{a_{0}}\right)}.$
(7)

Note that the preceding EoS parameter is clearly less than $-1$ (phantom DE)
and asymptotically approaches $-1$ as the universe expands. Nevertheless, the
behaviour of the universe is not that of a de Sitter model since the DE
density (5) and pressure are not constant. Furthermore, they even tend to
explode with the super-accelerated expansion of the universe, as found in Ref.
[58] (see also Ref. [32] and references therein). On the other hand, the
corresponding Hubble rate and its cosmic time derivative are

$\displaystyle H(t)$
$\displaystyle=\sqrt{\frac{\rho_{0}}{3}}\left(1+\frac{3A}{2\sqrt{\rho_{0}}}\ln\frac{a_{\star}}{a_{0}}\right)\exp\left[\frac{\sqrt{3}}{2}A\left(t-t_{\star}\right)\right],$
(8) $\displaystyle\dot{H}(t)$
$\displaystyle=\frac{A}{2}\sqrt{\rho_{0}}\left(1+\frac{3A}{2\sqrt{\rho_{0}}}\ln\frac{a_{\star}}{a_{0}}\right)\exp\left[\frac{\sqrt{3}}{2}A\left(t-t_{\star}\right)\right].$
(9)

Therefore, in this model the scale factor $a$, the Hubble parameter $H$ and
its cosmic derivative $\dot{H}$ diverge in the infinite distant future.
Nonetheless, bounded structures are shown to be disintegrated in a finite time
from present [32]. Moreover, disintegration can take place even before than
for the BR scenario [32].

From an observational point of view, the DE model described by the EoS given
in Eq. (4) has been shown to be compatible with the current expansion history
of our Universe, see Ref. [36]. In fact, authors in Ref. [36] found a
preferred value for the parameter $A$ when binding the model with current
observational data. Here we have denoted this value as $\bar{A}$. This is (in
geometric units $8\pi G=c=1$)

$\displaystyle\bar{A}=2.75\times 10^{-28}\ m^{-1},$ (10)

where $m$ stands for metres. Note that the small value obtained there suggests
that tiny deviations form the $\Lambda$CDM scenario are, indeed, the
observational preferred situation [36].

Even though in the present work we will focus on the LR abrupt event modelled
by the EoS given in Eq. (4), it should be stressed that this is not the only
scenario where a fate à la LR occurs. For example a more general EoS, still of
the form of Eq. (4), was considered in Ref. [40, 45]. This is

$p=-\rho-A\rho^{\alpha}.$ (11)

This EoS was thoroughly discussed in terms of singularity occurrence in Ref.
[58]. As there concluded, depending on the value of $\alpha$ several types of
DE-driven singularities (at finite time from present epoch) may occur.
Nevertheless, if $\alpha\leq 1/2$ the singularity is infinitely delayed in
time, thus, an abrupt event takes place. In this case, a straightforward
examination of the Hubble rate and its time derivative reveals that $\dot{H}$
remains finite and $H$ diverge when $\alpha\leq 0$. This corresponds to a
final fate à la LSBR. Whereas, both quantities diverge for $0<\alpha\leq 1/2$.
Hence, for the latter range of values for $\alpha$, the DE described by EoS
(11) leads the universe towards a LR abrupt event. For other examples of LR
cosmology see Refs. [32, 33] (see also Refs. [53, 54, 55] for LR cosmology in
$f(R)$ gravity).

##  III The LR in $f(R)$ theories of gravity

Given some cosmological background evolution it is possible to find an
alternative theory of gravity that leads to the same expansion history. The
group of techniques used to perform such a “reconstruction” task are commonly
known as “reconstruction methods” (for a review on the topic see, for example,
Ref. [61] and references therein). In this section we focus on reconstruction
methods within the framework of metric $f(R)$ alternative theories of gravity.
Thence, we shall look for an $f(R)$ theory of gravity able to reproduce the
super-accelerated expansion of the relativistic model filled with phantom DE
described by the EoS (4). Subsequently, as the LR abrupt event is inevitable
in the latter case, then, the reconstructed $f(R)$ theory will suffer the same
classical fate. For previous works on reconstruction techniques in $f(R)$
gravity see, for instance, Refs. [61, 62, 63, 64, 65, 66, 67]. See also Refs.
[51, 52, 68] for successful reconstruction of phantom DE-driven rip-like
events in metric $f(R)$ theories of gravity.

###  III.1 The reconstruction method

Following a line of reasoning similar to that presented in Ref. [68] (see also
Ref. [54]) we consider two cosmological evolutions to be equivalent at the
background level if the corresponding geometrical variables $H$, $\dot{H}$,
$R$ and $\dot{R}$ are identical.

Within GR, the expansion of the isotropic and homogeneous relativistic
universe is ruled by the Friedmann and Raychaudhuri equations (2) and (3),
respectively. Accordingly, the scalar curvature reads

$R=6\left(\dot{H}+2H^{2}+\frac{k}{a^{2}}\right)=\rho-3p.$ (12)

From the continuity equation for the perfect fluid, i.e.

$\dot{\rho}+3H(p+\rho)=0,$ (13)

and the Friedmann equation (2), it follows

$\displaystyle\dot{\rho}$
$\displaystyle=-3(p+\rho)\left(\frac{1}{3}\rho-\frac{k}{a^{2}}\right)^{\frac{1}{2}},$
(14) $\displaystyle\dot{p}$
$\displaystyle=-3(p+\rho)\left(\frac{1}{3}\rho-\frac{k}{a^{2}}\right)^{\frac{1}{2}}\frac{dp}{d\rho},$
(15)

where we have assumed $p=p(\rho)$. Therefore, the cosmic time derivative of
the scalar curvature $R$ is given by

$\dot{R}=-3\,(p+\rho)\left(\frac{1}{3}\rho-\frac{k}{a^{2}}\right)^{\frac{1}{2}}\left(1-3\frac{dp}{d\rho}\right).$
(16)

On the other hand, in the framework of metric $f(R)$ theories of gravity, the
evolution of the universe is described by the action

$\displaystyle S=\frac{1}{2}\int d^{4}x\sqrt{-g}f(R)+S_{m},$ (17)

where $S_{m}$ stands for the minimally coupled matter fields. Under these
circumstances, the field equations are no longer Eqs. (2) and (3). In fact the
first field equation, the so-called modified Friedmann equation, reads

$3H^{2}\frac{df}{dR}=\frac{1}{2}\left(R\frac{df}{dR}-f\right)-3H\dot{R}\frac{d^{2}f}{dR^{2}}-3\frac{k}{a^{2}}+\rho_{m},$
(18)

being $\rho_{m}$ the energy density of the minimally coupled matter fields. In
the reconstruction method, the above expression is considered as a
differential equation for some, a priori unknown, function $f(R)$, where the
coefficients are already fixed. That is, when the geometrical quantities
involved in Eq. (18) are set to be equal to those of the GR model that we want
to reproduce, then the background cosmological expansion of the resulting
metric $f(R)$ theory of gravity will be equivalent to that provided by the
general relativistic model. Therefore, in the next section we solve Eq. (18)
for $f$ when considering that the behaviour of $H$, $R$ and $\dot{R}$ is given
by the relativistic formulas (2) to (16). As we are interested on the
asymptotic behaviour of the universe, we can neglect the matter part, which
will be quickly red-shifted in the future. Furthermore, we shall focus our
attention to flat FLRW for the same reason, as the contribution of the spatial
curvature $k$ dilutes ($\propto a^{-2}$) with the cosmic expansion and,
ultimately, will become unimportant when compared with the phantom DE density.

###  III.2 $f(R)$ theories predicting the LR

We are mainly interested in $f(R)$ functions that mimics the behaviour given
in Eq. (II); i.e. on $f(R)$ functions that describe the same asymptotic future
behaviour of a FLRW universe to that of GR whose matter content corresponds to
a DE fluid whose EoS reads as in Eq. (4). Under those ansätze, the matter
content corresponding to barionic and dark matter can be neglected.

Note that the Hubble parameter in Eq. (8) is an exponential function of the
cosmic time. Thus, its time derivative is proportional to itself. For the sake
of simplicity we denoted by $\beta$ that proportionality constant; i.e.

$\dot{H}=\beta H,$ (19)

where we have defined

$\displaystyle\beta\coloneqq\frac{\sqrt{3}}{2}A.$ (20)

Consequently, it is straightforward to express the curvature scalar and its
time derivative in terms of the Hubble parameter. These are

$\displaystyle R$ $\displaystyle=6H\left(\beta+2H\right),$ (21)
$\displaystyle\dot{R}$ $\displaystyle=6\beta H\left(\beta+4H\right).$ (22)

As a result of these expressions, Eq. (18) simplifies when rewritten in terms
of the Hubble rate. Thus, the modified Friedmann equation, with $\rho_{m}=0$
and $k=0$, transforms into

$\displaystyle\beta H^{2}\left(\beta+4H\right)f_{HH}-\left[4\beta
H^{2}+H\left(\beta+H\right)\left(\beta+4H\right)\right]f_{H}$
$\displaystyle+\left(\beta+4H\right)^{2}f=0.$ (23)

The general solution to the preceding second order differential equation was
already found in Ref. [54], that is

$\displaystyle f(H)=$ $\displaystyle C_{1}\left(H^{4}-5\beta
H^{3}+2\beta^{2}H^{2}+2\beta^{3}H\right)$ $\displaystyle+C_{2}\bigg{[}\beta
H\left(\beta^{2}+4\beta H-H^{2}\right)e^{\frac{H}{\beta}}$
$\displaystyle+\left(H^{4}-5\beta
H^{3}+2\beta^{2}H^{2}+2\beta^{3}H\right)\textup{Ei}\left(\frac{H}{\beta}\right)\bigg{]},$
(24)

being $C_{1}$ and $C_{2}$ integration constants and Ei the exponential
integral function (see definition 5.1.2 in Ref. [69]).

On the other hand, to obtain the expression for $f(R)$ we need to calculate
the inverse of Eq. (21). This is

$H=\frac{1}{12}\left(-3\beta\pm\sqrt{9\beta^{2}+12R}\right).$ (25)

However, since the scalar curvature increases along with the DE density, see
(12) and (5), and the universe goes through an eternal expansion ($H>0$), only
the positive branch of the solution is compatible with our model. Thus, we
choose the positive sign in expression (25). Therefore, we get [54]

$\displaystyle f(R)$
$\displaystyle=c_{1}\left[27\beta^{4}+150\beta^{2}R-\beta\left(9\beta^{2}+12R\right)^{\frac{3}{2}}+2R^{2}\right]$
$\displaystyle+c_{2}\Bigg{\\{}\beta\left(-3\beta^{2}-2R+9\beta\sqrt{9\beta^{2}+12R}\right)$
$\displaystyle\times\left(-3\beta+\sqrt{9\beta+12R}\right)\exp\left(-\frac{1}{4}+\frac{1}{4}\sqrt{1+\frac{4R}{3\beta}}\right)$
$\displaystyle+\left[27\beta^{4}+150\beta^{2}R-\beta\left(9\beta^{2}+12R\right)^{\frac{3}{2}}+2R^{2}\right]$
$\displaystyle\times\textup{Ei}\left(-\frac{1}{4}+\frac{1}{4}\sqrt{1+\frac{4R}{3\beta}}\right)\Bigg{\\}},$
(26)

where $c_{1}$ and $c_{2}$ are arbitrary constants. We emphasize that the
metric $f(R)$ theory of gravity presented in Eq. (III.2) leads to the same
cosmological expansion as the DE fluid described by the EoS (4) in the context
of GR. However, in the framework of the modified theory of gravity, the
evolution of the DE content is mimicked by the modifications appearing in the
Friedmann equation as a result of $f(R)\neq R$. Consequently, since the
relativistic model is doomed to evolve towards the LR abrupt event, then the
reconstructed theory in Eq. (III.2) will suffer the same classical fate. We
further discuss on the viability of the theory given by Eq. (III.2) in Sec.
IV.2.

##  IV The LR in quantum $f(R)$ cosmology

The quantum fate of classical singularities can be addressed in the framework
of quantum cosmology: the application of quantum theory to the Universe as a
whole (see Refs. [70, 71] for a review on the topic). Although there are
different approaches to quantum cosmology, we focus on one of the first
attempts to quantize cosmological backgrounds [56, 72]. This quantum cosmology
is based on a canonical quantization with the Wheeler-DeWitt equation playing
a central role [56, 72, 73]. In addition, following the ideas presented in
Ref. [57], we adapt the Wheeler-DeWitt equation to $f(R)$ theories of gravity.
The resulting scheme is known as $f(R)$ quantum geometrodynamics. Within this
framework, we analyse the singularity avoidance by means of the DeWitt (DW)
criterion [56]. This is, the classical singularity might be avoided due to
quantum effects if the wave function of the Universe vanishes in the
configuration space close to the event. Therefore, this criterion is based on
a probabilistic interpretation of the wave function. Note that this criterion
has been successfully applied in several cosmological scenarios, see, e.g.,
Refs. [42, 48, 43, 49, 50, 51, 52, 74, 75, 76] among others.

###  IV.1 Modified Wheeler–DeWitt equation

Cosmological models in $f(R)$ theories of gravity can be described, in the so-
called Jordan frame, by the action

$S=\frac{1}{2}\int d^{4}x\sqrt{-g}f(R).$ (27)

When considering a FLRW background geometry, the above action can be cast into

$S=\frac{1}{2}\int dt\ \mathcal{L}(a,\dot{a},\ddot{a}),$ (28)

where the Lagrangian reads

$\mathcal{L}(a,\dot{a},\ddot{a})=\mathcal{V}_{(3)}\,a^{3}f(R),$ (29)

denoting by $\mathcal{V}_{(3)}$ the spatial $3$-dimensional volume. It is
well-known that metric $f(R)$ theories of gravity carry an extra degree of
freedom in comparison with GR (Einstein and Jordan formulation of $f(R)$
gravity can be found in [77, 78, 79, 80] and references therein). Therefore,
for the canonical quantization of alternative $f(R)$ theories of gravity a new
variable can be introduced to emphasize the existence of an extra degree of
freedom. Furthermore, for a suitable choice of this new variable the second
derivatives of the scale factor appearing in Eq. (28) can be removed.
Following the line of reasoning presented in Ref. [57], we select $R$ to be
the new variable. Thus, the action in Eq. (28) becomes

$\displaystyle S=\frac{1}{2}\int dt\ \mathcal{L}(a,\dot{a},R,\dot{R}).$ (30)

Nevertheless, since $R$ and $a$ are not independent (at the classical level),
their relation needs to be properly introduced in the theory via a Lagrange
multiplier, $\mu$, for the constraint $R=R(a,\dot{a},\ddot{a})$. Thence,

$\mathcal{L}=\mathcal{V}_{(3)}a^{3}\left\\{f(R)-\mu\left[R-6\left(\frac{\ddot{a}}{a}+\frac{\dot{a}^{2}}{a^{2}}+\frac{k}{a^{2}}\right)\right]\right\\}.$
(31)

The Lagrange multiplier can be determined varying the action with respect $R$.
This is

$\displaystyle\mu=f_{R}(R),$ (32)

with the notation $f_{R}\coloneqq df/dR$. Then, the Lagrangian can be
reformulated as

$\displaystyle\mathcal{L}(a,\dot{a},R,\dot{R})=$
$\displaystyle\mathcal{V}_{(3)}\Big{\\{}a^{3}\Big{[}f(R)-Rf_{R}(R)\Big{]}$
$\displaystyle-6a^{2}f_{RR}(R)\dot{a}\dot{R}+6af_{R}(R)(k-\dot{a}^{2})\Big{\\}},$
(33)

where $f_{RR}\coloneqq d^{2}f/dR^{2}$. Thereafter, in order to ease the
application of the quantization procedure, we proceed to diagonalize the
derivative part by the introduction of a new set of variables alike in Ref.
[57]. These are

$\displaystyle q$ $\displaystyle\coloneqq
a\sqrt{R_{\star}}\left(\frac{f_{R}}{f_{R_{\star}}}\right)^{\frac{1}{2}},$
(34a) $\displaystyle x$
$\displaystyle\coloneqq\frac{1}{2}\ln\left(\frac{f_{R}}{f_{R_{\star}}}\right),$
(34b)

being $R_{\star}$ a constant (not fixed at this point) needed for the above
change of variables to be well-defined. In the new variables the Lagrangian
(IV.1) transforms into

$\displaystyle\mathcal{L}(x,\dot{x},q,\dot{q})=$
$\displaystyle\mathcal{V}_{(3)}\left(\frac{R_{\star}f_{R}}{f_{R_{\star}}}\right)^{-\frac{3}{2}}q^{3}\left[f-6f_{R}\frac{\dot{q}^{2}}{q^{2}}\right.$
$\displaystyle\left.-Rf_{R}+6f_{R}\dot{x}^{2}+6k\frac{R_{\star}}{f_{R_{\star}}}\frac{f^{2}_{R}}{q^{2}}\right],$
(35)

where $f$ and $f_{R}$ are now understood as functions of $x$. This form of the
Lagrangian is already suitable for the quantization procedure.

Since the kinetic part has been diagonalized, the derivation of the
Hamiltonian is straightforward. The conjugate momenta read

$\displaystyle P_{q}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\dot{q}}=-12\mathcal{V}_{(3)}R_{\star}^{-\frac{3}{2}}f^{\frac{3}{2}}_{R_{\star}}f^{-\frac{1}{2}}_{R}q\dot{q},$
(36) $\displaystyle P_{x}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\dot{x}}=12\mathcal{V}_{(3)}R_{\star}^{-\frac{3}{2}}f^{\frac{3}{2}}_{R_{\star}}f^{-\frac{1}{2}}_{R}q^{3}\dot{x}.$
(37)

Then, the corresponding Hamiltonian is

$\displaystyle\mathcal{H}=$
$\displaystyle-\mathcal{V}_{(3)}q^{3}\left(\frac{R_{\star}f_{R}}{f_{R_{\star}}}\right)^{-3/2}\bigg{\\{}f+6k\frac{R_{\star}}{f_{R_{\star}}}\frac{f_{R}^{2}}{q^{2}}$
$\displaystyle-
Rf_{R}+\frac{6R_{\star}^{3}}{(12)^{2}\mathcal{V}_{(3)}^{2}f_{R_{\star}}^{3}}\frac{f_{R}^{2}}{q^{4}}\left[P_{q}^{2}-\frac{P_{x}^{2}}{q^{2}}\right]\bigg{\\}}.$
(38)

For the canonical quantization procedure, we assume
$P_{q}\to-i\hbar\partial_{q}$ and $P_{x}\to-i\hbar\partial_{x}$. As a result,
the classical Hamiltonian constraint $\mathcal{H}=0$ becomes the modified
Wheeler-DeWitt (mWDW) equation for the wave function $\Psi$ of the Universe
[56, 57, 71]. Therefore, the mWDW equation reads

$\mathcal{\hat{H}}\Psi=0.$ (39)

After suitable rearrangements, the preceding expression can be cast in the
form of the hyperbolic differential equation [57]

$\Big{[}\hbar^{2}q^{2}\partial^{2}_{q}-\hbar^{2}\partial^{2}_{x}-V(q,x)\Big{]}\Psi(q,x)=0,$
(40)

where the effective potential is given by

$V(q,x)=\frac{q^{4}}{\lambda^{2}}\left(k+\frac{q^{2}}{6R_{\star}f_{R_{\star}}}(f-Rf_{R})e^{-4x}\right),$
(41)

with $\lambda\coloneqq R_{\star}/(12\mathcal{V}_{(3)}f_{R_{\star}})$. Note
that when the expression of the $f(R)$ is given, the variables $x$ and $q$ are
univocally fixed. Then, $f$ and $Rf_{R}$ must be rewritten in terms of $x$. In
the next section, we will again focus on the case $k=0$.

###  IV.2 Asymptotic wave function close to the LR

From here on, we consider the expression for $f(R)$ gravity found in Eq.
(III.2) by means of the reconstruction methods applied in Sec. III.2. We
recall that the group of metric $f(R)$ theories of gravity presented in Eq.
(III.2) lead to the same cosmological evolution as the relativistic model
filled with phantom DE described by the EoS (4). Hence, both models share a
common final fate: the LR abrupt event.

It is important to note that, for a given $f(R)$ expression, the analytic
inversion of the definition for the variable $x$ in Eq. (34b) is critical for
the computation of the mWDW equation. This is because the relation $R=R(x)$ is
crucial in order to express the effective potential (41) in terms of $x$.
Unfortunately, the term multiplying $c_{2}$ in Eq. (III.2) prevent us from
inverting relation (34b) due to the presence of the exponential integral
function Ei. Therefore, for the sake of simplicity, we consider $c_{2}=0$.
Thus, we study the subgroup of metric $f(R)$ theories of gravity with a LR
abrupt event given by

$\displaystyle
f(R)=c_{1}\left[27\beta^{4}+150\beta^{2}R-\beta\left(9\beta^{2}+12R\right)^{\frac{3}{2}}+2R^{2}\right].$
(42)

For this particular $f(R)$, the change of variables (34) reads

$\displaystyle q$
$\displaystyle=a\sqrt{\frac{2c_{1}R_{\star}}{f_{R_{\star}}}}\left(75\beta^{2}-9\beta\sqrt{9\beta^{2}+12R}+2R\right)^{\frac{1}{2}},$
(43) $\displaystyle x$
$\displaystyle=\frac{1}{2}\ln\left[\frac{2c_{1}}{f_{R_{\star}}}\left(75\beta^{2}-9\beta\sqrt{9\beta^{2}+12R}+2R\right)\right],$
(44)

with
$f_{R_{\star}}=2c_{1}\left(75\beta^{2}-9\beta\sqrt{9\beta^{2}+12R_{\star}}+2R_{\star}\right)$.
Concerning the value of $R_{\star}$, note that Eqs. (4), (5) and (12) imply

$\displaystyle R=$ $\displaystyle
4\rho_{0}\left[1+\beta\sqrt{\frac{3}{\rho_{0}}}\ln\left(\frac{a}{a_{0}}\right)\right]^{2}$
$\displaystyle+6\beta\sqrt{\frac{\rho_{0}}{3}}\left[1+\beta\sqrt{\frac{3}{\rho_{0}}}\ln\left(\frac{a}{a_{0}}\right)\right].$
(45)

Thus, following the spirit for a physical meaningful $R_{\star}$ discussed in
Ref. [57] (and also in Refs. [51, 52]), we define this constant as the value
of the scalar curvature evaluated at some future scale factor $a=a_{\star}$ on
which the description of the universe by means of DE only becomes appropriate.
For the sake of concreteness, we set this moment to occur at
$a_{\star}=100a_{0}$ 111At that point in the expansion, the matter content
will be diluted with respect to the present concentration by a factor of
$100^{3}$. Roughly speaking this is $\Omega_{M}=10^{-6}\Omega_{M,0}\approx
3.06\times 10^{-7}$ and, therefore, $\Omega_{DE}\approx 1$.. Thence,

$\displaystyle R_{\star}\coloneqq$ $\displaystyle
4\rho_{0}\left(1+\beta\sqrt{\frac{3}{\rho_{0}}}\ln 100\right)^{2}$
$\displaystyle+6\beta\sqrt{\frac{\rho_{0}}{3}}\left(1+\beta\sqrt{\frac{3}{\rho_{0}}}\ln
100\right).$ (46)

As already discussed, $A$ and $\rho_{0}$ are constrained by observations, see
Eq. (10) (and Ref. [36]), and, therefore, it can be seen that this definition
for $R_{\star}$ makes the change of variables in Eq. (34) well-defined.
Furthermore, for $c_{1}>0$ , the metric $f(R)$ gravity in Eq. (42) fulfils the
viability conditions for late-time cosmology in $f(R)$ alternative theories of
gravity [81]. Namely, have a stable high curvature regime and a positive
effective Newton’s constant, $f_{RR}>0$ and $f_{R}>0$ for all $R>R_{\star}$,
respectively.

Next in the quantization procedure, we compute the inverse of Eq. (44). This
is

$\displaystyle R(x)=$ $\displaystyle
84\beta^{2}+\frac{1}{4c_{1}}f_{R_{\star}}e^{2x}$ $\displaystyle\pm
54\beta\sqrt{\frac{1}{48c_{1}}f_{R_{\star}}e^{2x}+2\beta^{2}}.$ (47)

However, only the positive branch is compatible with $R$ being an increasing
function of $x$ and $R>R_{\star}$. Therefore, we choose the positive sign in
the preceding expression. Then, the effective potential in the mWDW equation
reduces to

$\displaystyle V(q,x)=-U(x)q^{6},$ (48)

where

$\displaystyle U(x)=$
$\displaystyle\frac{f_{R_{\star}}}{48c_{1}\lambda^{2}R_{\star}}\Bigg{\\{}1+1644\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}+205992\frac{c_{1}^{2}\beta^{4}}{f_{R_{\star}}^{2}}e^{-4x}$
$\displaystyle+36\beta
e^{-x}\left(1+336\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}\right)\sqrt{\frac{3c_{1}}{f_{R_{\star}}}\left(1+96\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}\right)}$
$\displaystyle-12\beta\sqrt{\frac{3c_{1}}{f_{R_{\star}}}}e^{-x}\left[1+330\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}+18\beta
e^{-x}\sqrt{\frac{3c_{1}}{f_{R_{\star}}}\left(1+96\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}\right)}\right]$
$\displaystyle\times\left[1+339\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}+18\beta
e^{-x}\sqrt{\frac{3c_{1}}{f_{R_{\star}}}\left(1+96\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}\right)}\right]^{\frac{1}{2}}\Bigg{\\}}.$
(49)

Since the main motivation of the present work is to evaluate the wave function
$\Psi$ at the LR abrupt event, it is not strictly necessary to solve the mWDW
equation for the whole potential but focus only in the configuration space
near the abrupt event. For that aim, note that the most important condition
for the occurrence of the LR is the divergence of the scalar curvature. This
corresponds to $q\to\infty$ and $x\to\infty$. Then, without loss of
generality, we can suppose the variable $x$ to be large (but finite) when
close to the LR. This assumption allows us to expand the above expression as

$\displaystyle U(x)\approx$
$\displaystyle\frac{f_{R_{\star}}}{48c_{1}\lambda^{2}R_{\star}}\left[1+24\beta\sqrt{\frac{3c_{1}}{f_{R_{\star}}}}e^{-x}+672\beta^{2}\frac{c_{1}}{f_{R_{\star}}}e^{-2x}\right.$
$\displaystyle\left.+1152\beta^{3}\left(\frac{3c_{1}}{f_{R_{\star}}}\right)^{\frac{3}{2}}e^{-3x}+\mathcal{O}\left(e^{-4x}\right)\right].$
(50)

Additionally, further approximations can be implemented. As the LR is
approached, the majority of the terms appearing in the latter expansion are
exponentially suppressed (since $x\to\infty$). Then, in order to analyse the
asymptotic behaviour of $\Psi$, we keep only the dominant part of the
effective potential222This simplification is also supported by the
observational constraints on the model [36]. As the subdominant terms are
proportional to some power of the parameter $A$ (remember
$\beta=\sqrt{3}A/2$), then they are naturally suppressed since $A$ takes small
values when observationally constrained, see expression (10).. Therefore, at
first order approximation, we consider

$V(q,x)\approx-Bq^{6},$ (51)

being $B\coloneqq\frac{f_{R_{\star}}}{48c_{1}\lambda^{2}R_{\star}}$.
Consequently, the mWDW equation simplifies to

$\Big{[}\hbar^{2}q^{2}\partial^{2}_{q}-\hbar^{2}\partial^{2}_{x}+Bq^{6}\Big{]}\tilde{\Psi}(x,q)=0.$
(52)

We have used the notation $\tilde{\Psi}$ to emphasize that this is the
asymptotic form of Eq. (40), where only the leading-order part of the
effective potential is taken into account. This differential equation can be
solved via a separation ansatz for the wave function of the Universe,

$\displaystyle\tilde{\Psi}(q,x)=\sum_{\tilde{k}}b_{\tilde{k}}\zeta_{\tilde{k}}(q)\theta_{\tilde{k}}(x),$
(53)

being $b_{\tilde{k}}$ the amplitude of each solution and $\tilde{k}$ standing
for an integration constant related to the associated energy. Please, not
confuse $\tilde{k}$ with the spatial curvature $k$, which has been neglected.
As a result of the separation ansatz, the mWDW equation (52) implies

$\displaystyle\hbar^{2}\diff[2]{\theta_{\tilde{k}}}{x}-{\tilde{k}}^{2}\theta_{\tilde{k}}$
$\displaystyle=0,$ (54)
$\displaystyle\hbar^{2}q^{2}\diff[2]{\zeta_{\tilde{k}}}{q}+\left(Bq^{6}-\tilde{k}^{2}\right)\zeta_{\tilde{k}}$
$\displaystyle=0.$ (55)

The former equation can be straightforwardly solved and leads to

$\displaystyle\theta_{\tilde{k}}(x)=d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right),$
(56)

being $d_{1}$ and $d_{2}$ arbitrary constants. The above solutions correspond
to trigonometric or exponential functions of $x$, depending on the sign of
$\tilde{k}^{2}$. On the other hand, the equation for $\zeta_{\tilde{k}}$ can
be solved in an exact way by means of Bessel functions, c.f. 9.1.53 of Ref.
[69]. The solution can be written as

$\displaystyle\zeta_{\tilde{k}}(q)=$
$\displaystyle\sqrt{q}\left[u_{1}J_{\frac{1}{6}\sqrt{1+4\frac{\tilde{k}^{2}}{\hbar^{2}}}}\left(\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+u_{2}Y_{\frac{1}{6}\sqrt{1+4\frac{\tilde{k}^{2}}{\hbar^{2}}}}\left(\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right],$
(57)

being $J$ and $Y$ the Bessel functions of first and second order,
respectively, and $u_{1}$ and $u_{2}$ integration constants.

Therefore, as $q$ and $x$ explode when approaching the LR, the solutions for
the $\theta_{\tilde{k}}$ part remain finite if the constant $d_{1}$ is set to
zero for the case of $\tilde{k}^{2}$ being positive. Whereas the solutions for
$\zeta_{\tilde{k}}$ vanish as $q$ diverges since, for large values of $q$, Eq.
(IV.2) reduces to

$\displaystyle\zeta_{\tilde{k}}(q)\approx$
$\displaystyle\sqrt{\frac{6\hbar}{\pi}}\frac{1}{B^{\frac{1}{4}}q}\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right],$
(58)

where $\tilde{u}_{1}$ and $\tilde{u}_{2}$ depends on $\tilde{k}$, c.f. 9.2.1-2
of Ref. [69]. Thus, the total wave function $\tilde{\Psi}$,

$\displaystyle\tilde{\Psi}(q,x)\approx$
$\displaystyle\sqrt{\frac{6\hbar}{\pi}}\frac{1}{B^{\frac{1}{4}}q}\sum_{\tilde{k}}b_{\tilde{k}}\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{B}}{3\hbar}q^{3}\right)\right]\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right.$
$\displaystyle\left.+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right],$
(59)

shrinks to zero as the LR abrupt event is approached. Hence we conclude that
the DW condition is satisfied if one of the integration constants is fixed to
zero. Thus, as it happens in GR, this result hints towards the avoidance of
the LR doomsday in $f(R)$ cosmology due to quantum gravity effects.

Qualitatively, the wave function (IV.2) has the same asymptotic form as the
exhibited in Ref. [52] for the wave function of the LSBR abrupt event in
$f(R)$ gravity. However, owing to the fact that both events are different at
the classical level, differences were expected to manifest in the shape of the
wave functions. In fact, since the effective potentials entering the mWDW
equation in each case are different, see Eq. (48) and the analogous expression
for the LSBR scenario given in Eq. (38) of Ref. [34], then, the similarity
between both solutions seems to be an artefact of the approximations performed
in order to solve the mWDW equations in the asymptotic limit. Hence, this
resemblance is no longer expected to hold when further terms in the expansion
(IV.2) are taken into account. Following this reasoning, in the next section
we propose a less restrictive approach that allow us to solve the mWDW
equation for a wider region in the configuration space.

###  IV.3 Born-Oppenheimer approximation for $\Psi$

In this section, we address the solution of the complete mWDW equation,

$\displaystyle\Big{[}\hbar^{2}q^{2}\partial^{2}_{q}-\hbar^{2}\partial^{2}_{x}+U(x)q^{6}\Big{]}\Psi(q,x)=0,$
(60)

where $U(x)$ is given by Eq. (IV.2). Owing to the fact that the potential term
in Eq. (60), that is $U(x)q^{6}$, now contains both variables, a separation
ansatz like (53) will no longer apply. Instead, we propose an adiabatic semi-
separability type ansatz for the wave function of the Universe. This is based
on the so-called Born-Oppenheimer (BO) ansatz, originally formulated in the
context of molecular physics [82]. In cosmological scenarios, this
approximation is often implemented such that the geometrical part of the total
wave function (usually depending on the scale factor $a$) is factored out from
the section encompassing the physical fields contained in the universe, which
in turn are considered to depend adiabatically on the background geometry, see
Refs. [83, 84, 85] (for a recent work see, e.g., Ref. [48]). In this way, it
may feel tempting to apply the BO quasi-separability by factorizing the wave
function $\Psi$ into a part depending only on $q$ and another related to both
$q$ and $x$, since the scale factor $a$ enters only in $q$ and, therefore,
disregarding the variable $x$ as carrying exclusively “matter” degrees of
freedom. However, this way of reasoning is misleading. In fact, this naive
separation results ultimately in both parts of the wave function depending on
both original variables, $a$ and $R$ (since $a$ and $R$ enters the variable
$q$). Consequently, following that procedure the interaction between $a$ and
$R$ will be present in both parts of the wave function. Thus, this would make
the quasi-separability ansatz pointless from the very beginning. The
misconception leading to such failure is the innocent consideration that the
mWDW is caring two degrees of freedom of different nature, those are: one
geometrical and another rather related to matter fields. Although this is
often true in cosmological scenarios, it is not longer the case for the mWDW
equation of $f(R)$ cosmology, that is Eq. (40). Here we have two genuinely
geometrical variables. These are the scale factor $a$ and the scalar curvature
$R$, both contained in the definition of $q$ and $x$. Hence, a different
formulation of the BO ansatz for solving Eq. (60) is needed.

For that purpose, it should be stressed that $R$ can be considered to be more
fundamental from a geometrical point of view than the scale factor. Hence, if
we are to conserve the spirit of the original BO ansatz, this is to quantize
the geometry at first place and, after that, the remnant physical fields
taking into account the back-reaction effects, then we should factorize $\Psi$
into a part depending only on $R$ (as the main geometrical variable) and
another depending on $R$ and $a$. Therefore we propose the following ansatz à
la Born-Oppenheimer:

$\displaystyle\Psi(q,x)=\sum_{\tilde{k}}b_{\tilde{k}}\chi_{\tilde{k}}(q,x)\varphi_{\tilde{k}}(x).$
(61)

We emphasise that $x$ depends only on $R$, whereas $q$ contains both $a$ and
$R$, see definitions in (34). In addition, $b_{\tilde{k}}$ stands for the
amplitude of each solution and $\tilde{k}$ is related with the associated
energy. As a result, the mWDW equation (60) reads

$\displaystyle\hbar^{2}q^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{q}-\hbar^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{x}-2\hbar^{2}\diffp{\chi_{\tilde{k}}}{x}\diff{\varphi_{\tilde{k}}}{x}$
$\displaystyle-\hbar^{2}\chi_{\tilde{k}}\diff[2]{\varphi_{\tilde{k}}}{x}+U(x)q^{6}\chi_{\tilde{k}}\varphi_{\tilde{k}}=0.$
(62)

Then, the contribution of the second and third terms can be neglected due to
the adiabatic assumption. (The validity of this approximations is justified in
the Appendix A.) Thus, Eq. (IV.3) implies the following equations

$\displaystyle\hbar^{2}\diff[2]{\varphi_{\tilde{k}}}{x}-\tilde{k}^{2}\varphi_{\tilde{k}}$
$\displaystyle=0,$ (63)
$\displaystyle\hbar^{2}q^{2}\diffp[2]{\chi_{\tilde{k}}}{q}+\left[U(x)q^{6}-\tilde{k}^{2}\right]\chi_{\tilde{k}}$
$\displaystyle=0.$ (64)

The former equation can be solved in the same fashion as Eq. (54). The
solutions are exponential and trigonometric functions, depending on the sign
of $\tilde{k}^{2}$,

$\displaystyle\varphi_{\tilde{k}}(x)=d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right),$
(65)

being $d_{1}$ and $d_{2}$ integration constants. On the other hand, due to the
adiabatic approximation, the potential term $U(x)$ appearing in Eq. (64) is
treated like a (quasi)constant parameter when solving for $\chi_{\tilde{k}}$.
Note that this assumption is supported on the fact that $U(x)$, given in Eq.
(IV.2), converges very quickly to a constant value when $\beta$ is
observationally constrained [see the expansion in Eq. (IV.2)]. Thence, the
most general solution for $\chi_{\tilde{k}}$ is, c.f. 9.1.53 of Ref. [69],

$\displaystyle\chi_{\tilde{k}}(q,x)=$
$\displaystyle\sqrt{q}\left[u_{1}J_{\frac{1}{6}\sqrt{1+\frac{4\tilde{k}^{2}}{\hbar^{2}}}}\left(\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+u_{2}Y_{\frac{1}{6}\sqrt{1+\frac{4\tilde{k}^{2}}{\hbar^{2}}}}\left(\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right],$
(66)

with $u_{1}$ and $u_{2}$ integration constants.

Near the LR abrupt event, expression (65) remains bounded for $\tilde{k}^{2}$
negative. However, when $\tilde{k}^{2}$ is positive, $\varphi_{\tilde{k}}$ is
finite if and only if the constant $d_{1}$ is set to vanish. In contrast, the
solutions for the $\chi_{\tilde{k}}$ function have all the same asymptotic
form. This is

$\displaystyle\chi_{\tilde{k}}(q,x)\approx$
$\displaystyle\sqrt{\frac{6\hbar}{\pi}}\frac{1}{U(x)^{\frac{1}{4}}q}\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right],$
(67)

for large values of $q$, c.f. 9.2.1-2 of Ref. [69], where the integration
constants $\tilde{u}_{1}$ and $\tilde{u}_{2}$ now depend on $\tilde{k}^{2}$.
Therefore, the asymptotic form of the total wave function $\Psi$ reads

$\displaystyle\Psi(q,x)\approx$
$\displaystyle\sqrt{\frac{6\hbar}{\pi}}\frac{1}{U(x)^{\frac{1}{4}}q}\sum_{\tilde{k}}b_{\tilde{k}}\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right.$
$\displaystyle\left.+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right]$
$\displaystyle\times\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right].$
(68)

As $U(x)$ tends to a constant value when $x$ explodes, the wave function
cancels at the LR abrupt event when one of the integrations constants is set
to zero, $d_{1}=0$ for $\tilde{k}^{2}$ positive, in accordance with the
results of the previous section. Thus, pointing towards the avoidance of this
fatal fate. Nevertheless, since $U(x)>B$ and $U(x)\to B$ asymptotically, then,
the rate at which the wave function shrinks is increased with respect to the
asymptotic approach performed in the previous section. Ergo, subdominant order
contributions to the effective potential speed up the vanishing rate of the
wave function $\Psi$. Furthermore, since the approximation presented in this
section is less restrictive than the asymptotic approach previously performed,
the wave function here obtained is valid in a broader region in the
configuration space.

##  V Conclusions

The LR abrupt event is a cosmic doomsday arising in some cosmological models
where the accelerated expansion of the universe is driven by a DE of phantom
nature. Since some of these models have been shown to be able to describe the
current cosmological observations [36], our own Universe may evolve towards
this singular fate. However, quantum effects can ultimately become significant
and prevent the occurrence of such a doomsday. In fact, for the case of the
background evolution being that provided by GR, it has already been
established that the DW criterion for singularity avoidance can be satisfied
for the particular LR model considered here. Subsequently, in this work we
have addressed the question whether this is still true when the classical
evolution of the universe is due to an $f(R)$ metric theory of gravity.

Hence, we have applied the so-called reconstruction methods to find a group of
$f(R)$ theories of gravity that produce the same expansion history as that of
a relativistic model filled with a DE fluid described by the EoS given in Eq.
(4). Thus, obtaining the group of metric $f(R)$ theories of gravity that
predict a classical fate à la little rip.

Thereafter, we have studied the quantum fate of the cosmos governed by one of
the metric $f(R)$ gravity theories obtained before. The quantum analysis was
performed within the framework of $f(R)$ quantum geometrodynamics, with the
mWDW equation playing a central role. We have solved the mWDW equation and
showed that the solutions satisfy the DW criterion when one of the integration
constants is set to zero. Thus, as it also happens in general relativity, the
fulfilment of the aforementioned condition hints towards the avoidance of this
doomsday in $f(R)$ cosmology.

Furthermore, we have performed different approaches when solving the mWDW
equation. On a first approximation, we have solved the asymptotic form of the
mWDW when only leading order terms in the potential are considered.
Afterwards, in a second approach, we have used a BO type approximation (61) in
order to analyse the behaviour of the wave function $\Psi$ for the complete
mWDW equation. As a consequence of this less restrictive approach, the
resulting wave function exhibits a richer behaviour. Indeed, as the compliance
region of the BO approach was shown to be wider than just the asymptotic
regime, this method for solving the mWDW equation can be useful for comparing
the wave functions of different events. Especially when those events share a
similar asymptotic regime, since the asymptotic approach to $\Psi$ will not
find any differences there.

It should be noted, however, that we have fixed to zero an integration
constant in order to find vanishing solutions at the abrupt event. Therefore,
we have disregarded a subgroup of solutions to the mWDW equation as
unphysical. If future investigations show the importance of the dismissed
solutions, then, it would be concluded that the DW criterion may not always be
fulfilled for solutions of physical interest.

## Acknowledgments

The research of T. B. V. and P. M. M. is supported by MINECO (Spain) project
PID2019-107394GB-I00 (AEI/FEDER, UE). T. B. V. also acknowledge financial
support from project FIS2016-78859-P (AEI/FEDER, UE) through grant
PAII46/20-08/2020-03, and from Universidad Complutense de Madrid and Banco de
Santander through grant CT63/19-CT64/19. The research of M. B. L. is supported
by the Basque Foundation of Science Ikerbasque. She also would like to
acknowledge the partial support from the Basque government Grant No. IT956-16
(Spain) and the project FIS2017-85076-P (MINECO/AEI/FEDER, UE).

##  Appendix A Validity of the BO approximation

During the application of the BO-type ansatz (61) performed in Sec. IV.3, we
have considered that $\chi_{\tilde{k}}(q,x)$ depends adiabatically on $x$.
Therefore, we have neglected the contribution of some parts in Eq. (IV.3).
This approach is valid as long as the corresponding solutions satisfy

$\displaystyle\hbar^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{x},$
$\displaystyle\
2\hbar^{2}\diffp{\chi_{\tilde{k}}}{x}\diff{\varphi_{\tilde{k}}}{x}\ll$
$\displaystyle\hbar^{2}q^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{q},\
\hbar^{2}\chi_{\tilde{k}}\diff[2]{\varphi_{\tilde{k}}}{x},\
U(x)q^{6}\chi_{\tilde{k}}\varphi_{\tilde{k}}.$ (69)

As a result of this approximation, the solutions for $\varphi_{\tilde{k}}$ and
$\chi_{\tilde{k}}$ are presented in Eq. (65) and Eq. (IV.3), respectively.
Then, the terms we keep in (IV.3) read

$\displaystyle\hbar^{2}q^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{q}\approx-U(x)q^{6}\chi_{\tilde{k}}\varphi_{\tilde{k}}\approx-\sqrt{\frac{6\hbar}{\pi}}U(x)^{\frac{3}{4}}q^{5}$
$\displaystyle\times\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right]$
$\displaystyle\times\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right.\left.+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right],$
(70)
$\displaystyle\hbar^{2}\chi_{\tilde{k}}\diff[2]{\varphi_{\tilde{k}}}{x}\approx\sqrt{\frac{6\hbar}{\pi}}\frac{\tilde{k}^{2}}{U(x)^{\frac{1}{4}}q}$
$\displaystyle\times\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right]$
$\displaystyle\times\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right].$
(71)

While, the neglected terms behave asymptotically as

$\displaystyle\hbar^{2}\varphi_{\tilde{k}}\diffp[2]{\chi_{\tilde{k}}}{x}\approx-\frac{1}{36}\sqrt{\frac{6\hbar}{\pi}}\frac{U^{\prime}(x)^{2}}{U(x)^{\frac{5}{4}}}q^{5}$
$\displaystyle\times\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)+\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right]$
$\displaystyle\times\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)+d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right],$
(72) $\displaystyle
2\hbar^{2}\diffp{\chi_{\tilde{k}}}{x}\diff{\varphi_{\tilde{k}}}{x}\approx\frac{i}{3}\sqrt{\frac{6\hbar\tilde{k}^{2}}{\pi}}\frac{U^{\prime}(x)}{U(x)^{\frac{3}{4}}}q^{2}$
$\displaystyle\times\left[\tilde{u}_{1}\exp\left(i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)-\tilde{u}_{2}\exp\left(-i\frac{\sqrt{U(x)}}{3\hbar}q^{3}\right)\right]$
$\displaystyle\times\left[d_{1}\exp\left(\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)-d_{2}\exp\left(-\frac{\sqrt{\tilde{k}^{2}}}{\hbar}x\right)\right].$
(73)

Note that, for $\tilde{k}^{2}$ positive, the constants $d_{1}$ must be zero in
order to have a vanishing wave function at the LR. Thus, to analyse the
validity of the performed approximation we compare the largest of the
neglected terms with the smallest of the saved ones. This is the ratio
$\varepsilon$,

$\displaystyle\varepsilon=\left|\frac{\hbar^{2}\varphi_{\tilde{k}}\partial^{2}_{x}\chi_{\tilde{k}}}{\hbar^{2}\chi_{\tilde{k}}\partial^{2}_{x}\varphi_{\tilde{k}}}\right|\approx\frac{U^{\prime}(x)^{2}}{U(x)}\frac{q^{6}}{36\tilde{k}^{2}}.$
(74)

Consequently, the approximation is valid as long as $\varepsilon\ll 1$. To
obtain the compliance region of this condition, note

$\displaystyle\frac{U^{\prime}(x)^{2}}{U(x)}\approx$ $\displaystyle
36\frac{\beta^{2}}{\lambda^{2}R_{\star}}e^{-2x}\left[1+\frac{40}{3}\sqrt{\frac{3c_{1}}{f_{R_{\star}}}}\beta
e^{-x}\right.$
$\displaystyle\left.+\frac{832}{3}\frac{c_{1}\beta^{2}}{f_{R_{\star}}}e^{-2x}+1152\sqrt{\frac{3c_{1}^{3}}{f_{R_{\star}}^{3}}}\beta^{3}e^{-3x}\right.$
$\displaystyle+\mathcal{O}\left(e^{-4x}\right)\Bigg{]},$ (75)

when $x$ is large. Then, in the configuration space near the cosmic event,

$\displaystyle\varepsilon\approx\frac{\beta^{2}}{\lambda^{2}\tilde{k}^{2}R_{\star}}e^{-2x}q^{6}.$
(76)

Finally, $\varepsilon\ll 1$ near the LR if $\beta$ is sufficiently small, i.e
for small value of $A$. Note that this corresponds, in fact, to the
observationally preferred situation [36]. [We recall that $A$ is of order
$10^{-28}$ when observationally constrained, see Eq. (10)]. Therefore, when
the parameters of the theory are observationally constrained, the
approximation is valid throughout the semi-classical regime towards the abrupt
event; where the variables $q$ and $x$ increase but not sufficiently rapidly
to compensate the small value of $\beta^{2}$. Hence, for the purpose of this
work, that is to analyse the fulfilment of the DW criterion in the
configuration space close to the LR, this approximation is valid.

## References

  * [1] A. G. Riess et al. [Supernova Search Team], “Observational evidence from supernovae for an accelerating universe and a cosmological constant”, Astron. J. 116 (1998) 1009, [arXiv:astro-ph/9805201].
  * [2] S. Perlmutter et al. [Supernova Cosmology Project Collaboration], “Measurements of Omega and Lambda from 42 high redshift supernovae”, Astrophys. J. 517 (1999) 565, [arXiv:astro-ph/9812133].
  * [3] S. Weinberg, “The cosmological constant problem”, Rev. Mod. Phys. 61 (1989) 1.
  * [4] R. R. Caldwell, R. Dave and P. J. Steinhardt, “Cosmological imprint of an energy component with general equation of state”, Phys. Rev. Lett. 80 (1998), 1582-1585, [arXiv:astro-ph/9708069 [astro-ph]].
  * [5] S. Tsujikawa, “Quintessence: A Review”, Class. Quant. Grav. 30 (2013), 214003, [arXiv:1304.1961 [gr-qc]].
  * [6] T. Chiba, T. Okabe and M. Yamaguchi, “Kinetically driven quintessence”, Phys. Rev. D 62 (2000), 023511, [arXiv:astro-ph/9912463 [astro-ph]].
  * [7] R. R. Caldwell, “A Phantom menace?”, Phys. Lett. B 545 (2002), 23-29, [arXiv:astro-ph/9908168 [astro-ph]].
  * [8] G. W. Gibbons, “Cosmological evolution of the rolling tachyon”, Phys. Lett. B 537 (2002), 1-4, [arXiv:hep-th/0204008 [hep-th]].
  * [9] T. Padmanabhan, “Accelerated expansion of the universe driven by tachyonic matter”, Phys. Rev. D 66 (2002), 021301, [arXiv:hep-th/0204150 [hep-th]].
  * [10] A. Y. Kamenshchik, U. Moschella and V. Pasquier, “An Alternative to quintessence”, Phys. Lett. B 511 (2001), 265-268, [arXiv:gr-qc/0103004 [gr-qc]].
  * [11] M. C. Bento, O. Bertolami and A. A. Sen, “Generalized Chaplygin gas, accelerated expansion and dark energy matter unification”, Phys. Rev. D 66 (2002), 043507, [arXiv:gr-qc/0202064 [gr-qc]].
  * [12] M. Li, “A Model of holographic dark energy”, Phys. Lett. B 603 (2004), 1, [arXiv:hep-th/0403127 [hep-th]].
  * [13] S. Nojiri and S. D. Odintsov, “Modified f(R) gravity consistent with realistic cosmology: From matter dominated epoch to dark energy universe”, Phys. Rev. D 74 (2006), 086005, [arXiv:hep-th/0608008 [hep-th]].
  * [14] S. Nojiri and S. D. Odintsov, “Gravity assisted dark energy dominance and cosmic acceleration”, Phys. Lett. B 599 (2004), 137, [astro-ph/0403622].
  * [15] G. Allemandi, A. Borowiec, M. Francaviglia and S. D. Odintsov, “Dark energy dominance and cosmic acceleration in first order formalism”, Phys. Rev. D 72 (2005), 063505, [gr-qc/0504057].
  * [16] M. H. Dehghani, “Accelerated expansion of the Universe in Gauss-Bonnet gravity”, Phys. Rev. D 70 (2004), 064009, [arXiv:hep-th/0404118 [hep-th]].
  * [17] S. Nojiri and S. D. Odintsov, “Modified Gauss-Bonnet theory as gravitational alternative for dark energy”, Phys. Lett. B 631 (2005), 1-6, [arXiv:hep-th/0508049 [hep-th]].
  * [18] T. Harko, F. S. N. Lobo, S. Nojiri and S. D. Odintsov, “$f(R,T)$ gravity”, Phys. Rev. D 84 (2011), 024020, [arXiv:1104.2669 [gr-qc]].
  * [19] G. R. Bengochea and R. Ferraro, “Dark torsion as the cosmic speed-up”, Phys. Rev. D 79 (2009), 124019, [arXiv:0812.1205 [astro-ph]].
  * [20] G. W. Horndeski, “Second-order scalar-tensor field equations in a four-dimensional space”, Int. J. Mod. Phys. 10 (1974), 363-384.
  * [21] R. Kase and S. Tsujikawa, “Dark energy in Horndeski theories after GW170817: A review”, Int. J. Mod. Phys. D 28 (2019) no.05, 1942005, [arXiv:1809.08735 [gr-qc]].
  * [22] D. Huterer and D. L. Shafer, “Dark energy two decades after: Observables, probes, consistency tests”, Rept. Prog. Phys. 81 (2018) no.1, 016901, [arXiv:1709.01091 [astro-ph.CO]].
  * [23] K. Bamba, S. Capozziello, S. Nojiri and S. D. Odintsov, “Dark energy cosmology: the equivalent description via different theoretical models and cosmography tests”, Astrophys. Space Sci. 342 (2012), 155-228, [arXiv:1205.3421 [gr-qc]].
  * [24] N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters”, Astron. Astrophys. 641 (2020), A6 [arXiv:1807.06209 [astro-ph.CO]].
  * [25] P. A. R. Ade et al. [Planck Collaboration], “Planck 2015 results. XIII. Cosmological parameters”, Astron. Astrophys. 594 (2016) A13, [arXiv:1502.01589 [astro-ph.CO]].
  * [26] G. Risaliti and E. Lusso, “Cosmological constraints from the Hubble diagram of quasars at high redshifts”, Nature Astron. 3 (2019) no.3, 272-277, [arXiv:1811.02590 [astro-ph.CO]].
  * [27] E. Di Valentino, A. Melchiorri, E. V. Linder and J. Silk, “Constraining Dark Energy Dynamics in Extended Parameter Space”, Phys. Rev. D 96 (2017) no.2, 023523, [arXiv:1704.00762 [astro-ph.CO]].
  * [28] E. Di Valentino, E. V. Linder and A. Melchiorri, “Vacuum phase transition solves the $H_{0}$ tension”, Phys. Rev. D 97 (2018) no.4, 043528, [arXiv:1710.02153 [astro-ph.CO]].
  * [29] R. R. Caldwell, M. Kamionkowski and N. N. Weinberg, “Phantom Energy: Dark Energy with $w<-1$ Causes a Cosmic Doomsday”, Phys. Rev. Lett. 91 (2003), 071301 [arXiv:astro-ph/0302506 [astro-ph]].
  * [30] M. Bouhmadi-López, P. F. González-Díaz and P. Martín-Moruno, “Worse than a big rip?”, Phys. Lett. B 659 (2008), 1-5, [arXiv:gr-qc/0612135 [gr-qc]].
  * [31] M. Bouhmadi-López, P. F. González-Díaz and P. Martín-Moruno, “On the generalised Chaplygin gas: Worse than a big rip or quieter than a sudden singularity?”, Int. J. Mod. Phys. D 17 (2008), 2269-2290, [arXiv:0707.2390 [gr-qc]].
  * [32] P. H. Frampton, K. J. Ludwick and R. J. Scherrer, “The Little Rip”, Phys. Rev. D  84 (2011) 063003, [arXiv:1106.4996 [astro-ph.CO]].
  * [33] P. H. Frampton, K. J. Ludwick, S. Nojiri, S. D. Odintsov and R. J. Scherrer, “Models for Little Rip Dark Energy”, Phys. Lett. B 708 (2012) 204, [arXiv:1108.0067 [hep-th]].
  * [34] M. Bouhmadi-López, A. Errahmani, P. Martín-Moruno, T. Ouali and Y. Tavakoli, “The little sibling of the big rip singularity”, Int. J. Mod. Phys. D 24 (2015) no.10, 1550078, [arXiv:1407.2446 [gr-qc]].
  * [35] I. Albarran, M. Bouhmadi-López and J. Morais, “Cosmological perturbations in an effective and genuinely phantom dark energy Universe”, Phys. Dark Univ. 16 (2017), 94-108, [arXiv:1611.00392 [astro-ph.CO]].
  * [36] A. Bouali, I. Albarran, M. Bouhmadi-López and T. Ouali, “Cosmological constraints of phantom dark energy models”, Phys. Dark Univ. 26 (2019), 100391, [arXiv:1905.07304 [astro-ph.CO]].
  * [37] J. D. Barrow, “Sudden future singularities”, Class. Quant. Grav. 21 (2004), L79-L82, [arXiv:gr-qc/0403084 [gr-qc]].
  * [38] K. Lake, “Sudden future singularities in FLRW cosmologies”, Class. Quant. Grav. 21 (2004), L129, [arXiv:gr-qc/0407107 [gr-qc]].
  * [39] J. D. Barrow, “More general sudden singularities”, Class. Quant. Grav. 21 (2004), 5619-5622, [arXiv:gr-qc/0409062 [gr-qc]].
  * [40] S. Nojiri, S. D. Odintsov and S. Tsujikawa, “Properties of singularities in (phantom) dark energy universe”, Phys. Rev. D 71 (2005) 063004, [arXiv:hep-th/0501025 [hep-th]].
  * [41] M. P. Dabrowski, “Are singularities the limits of cosmology?”, [arXiv:1407.4851 [gr-qc]].
  * [42] M. Bouhmadi-López, C. Kiefer and P. Martín-Moruno, “Phantom singularities and their quantum fate: general relativity and beyond—a CANTATA COST action topic”, Gen. Rel. Grav. 51 (2019) no.10, 135, [arXiv:1904.01836 [gr-qc]].
  * [43] M. P. Dabrowski, C. Kiefer and B. Sandhöfer, “Quantum phantom cosmology”, Phys. Rev. D 74 (2006), 044022, [arXiv:hep-th/0605229 [hep-th]].
  * [44] A. Kamenshchik, C. Kiefer and B. Sandhöfer, “Quantum cosmology with big-brake singularity”, Phys. Rev. D 76 (2007), 064032, [arXiv:0705.1688 [gr-qc]].
  * [45] S. Nojiri and S. D. Odintsov, “The Final state and thermodynamics of dark energy universe”, Phys. Rev. D 70 (2004) 103522, [arXiv:hep-th/0408170 [hep-th]].
  * [46] E. Elizalde, S. Nojiri and S. D. Odintsov, “Late-time cosmology in (phantom) scalar-tensor theory: Dark energy and the cosmic speed-up”, Phys. Rev. D 70 (2004), 043539, [arXiv:hep-th/0405034 [hep-th]].
  * [47] S. Nojiri and S. D. Odintsov, “Quantum escape of sudden future singularity”, Phys. Lett. B 595 (2004), 1-8, [arXiv:hep-th/0405078 [hep-th]].
  * [48] M. Bouhmadi-López, C. Kiefer, B. Sandhöfer and P. Vargas Moniz, “On the quantum fate of singularities in a dark-energy dominated universe”, Phys. Rev. D 79 (2009) 124035,[arXiv:0905.2421 [gr-qc]].
  * [49] I. Albarran, M. Bouhmadi-López, C. Kiefer, J. Marto and P. Vargas Moniz, “Classical and quantum cosmology of the little rip abrupt event”, Phys. Rev. D 94 (2016) no.6, 063536, [arXiv:1604.08365 [gr-qc]].
  * [50] I. Albarran, M. Bouhmadi-López, F. Cabral and P. Martín-Moruno, “The quantum realm of the ”Little Sibling” of the Big Rip singularity”, JCAP 11 (2015), 044, [arXiv:1509.07398 [gr-qc]].
  * [51] A. Alonso-Serrano, M. Bouhmadi-López and P. Martín-Moruno, “$f(R)$ quantum cosmology: avoiding the Big Rip”, Phys. Rev. D 98 (2018) no. 10, 104004, [arXiv:1802.03290 [gr-qc]].
  * [52] T. Borislavov Vasilev, M. Bouhmadi-López and P. Martín-Moruno, “Classical and quantum fate of the little sibling of the big rip in $f(R)$ cosmology”, Phys. Rev. D 100 (2019) no.8, 084016, [arXiv:1907.13081 [gr-qc]].
  * [53] S. Nojiri, S. D. Odintsov and D. Sáez-Gómez, “Cyclic, ekpyrotic and little rip universe in modified gravity”, AIP Conf. Proc. 1458 (2012) no.1, 207-221, [arXiv:1108.0767 [hep-th]].
  * [54] A. N. Makarenko, V. V. Obukhov and I. V. Kirnos, “From Big to Little Rip in modified F(R,G) gravity”, Astrophys. Space Sci. 343 (2013) 481, [arXiv:1201.4742 [gr-qc]].
  * [55] A. N. Makarenko, S. D. Odintsov and G. J. Olmo, “Little Rip, $\Lambda$CDM and singular dark energy cosmology from Born-Infeld-$f(R)$ gravity”, Phys. Lett. B 734 (2014), 36-40, [arXiv:1404.2850 [gr-qc]].
  * [56] B. S. DeWitt, “Quantum Theory of Gravity. I. The Canonical Theory”, Phys. Rev. 160 (1967) 1113.
  * [57] A. Vilenkin, “Classical and Quantum Cosmology of the Starobinsky Inflationary Model”, Phys. Rev. D 32 (1985) 2511.
  * [58] H. S̆tefanc̆ić, “Expansion around the vacuum equation of state - Sudden future singularities and asymptotic behavior”, Phys. Rev. D 71 (2005) 084024, [arXiv:astro-ph/0411630 [astro-ph]].
  * [59] S. Nojiri and S. D. Odintsov, “Inhomogeneous equation of state of the universe: Phantom era, future singularity and crossing the phantom barrier”, Phys. Rev. D 72 (2005), 023003, [arXiv:hep-th/0505215 [hep-th]].
  * [60] M. Bouhmadi-López, “Phantom-like behaviour in dilatonic brane-world scenario with induced gravity”, Nucl. Phys. B 797 (2008), 78-92, [arXiv:astro-ph/0512124 [astro-ph]].
  * [61] S. Nojiri and S. D. Odintsov, “Unified cosmic history in modified gravity: from F(R) theory to Lorentz non-invariant models”, Phys. Rept. 505 (2011), 59-144, [arXiv:1011.0544 [gr-qc]].
  * [62] S. Capozziello, V. F. Cardone and A. Troisi, “Reconciling dark energy models with $f(R)$ theories”, Phys. Rev. D 71 (2005) 043503, [arXiv: astro-ph/0501426].
  * [63] Á. de la Cruz-Dombriz and A. Dobado, “A $f(R)$ gravity without cosmological constant”, Phys. Rev. D 74 (2006) 5, [arXiv:grqc/0607118].
  * [64] S. Nojiri and S. D. Odintsov, “Modified gravity and its reconstruction from the universe expansion history”, J. Phys. Conf. Ser. 66 (2007) 012005, [arXiv:hep-th/0611071].
  * [65] S. Nojiri, S. D. Odintsov, A. Toporensky and P. Tretyakov, “Reconstruction and deceleration-acceleration transitions in modified gravity”, Gen. Rel. Grav. 42 (2010), 1997-2008, [arXiv:0912.2488 [hep-th]].
  * [66] P. K. S. Dunsby, E. Elizalde, R. Goswami, S. Odintsov and D. S. Gomez, “On the LCDM Universe in f(R) gravity”, Phys. Rev. D 82 (2010), 023519 [arXiv:1005.2205 [gr-qc]].
  * [67] S. Carloni, R. Goswami and P. K. S. Dunsby, “A new approach to reconstruction methods in $f(R)$ gravity”, Clas. Quant. Grav. 29 (2012) 135012, [arXiv:1005.1840 [gr-qc]].
  * [68] J. Morais, M. Bouhmadi-López and S. Capozziello, “Can $f(R)$ gravity contribute to (dark) radiation?”, JCAP 1509 (2015) no.09, 041, [arXiv:1507.02623 [gr-qc]].
  * [69] M. Abramowitz and I. A. Stegun, Handbook of Mathematical Functions, Dover Publications (1980), ISBN 978-0486612720.
  * [70] C. Kiefer and B. Sandhöefer, “Quantum Cosmology”, [arXiv:0804.0672 [gr-qc]].
  * [71] C. Kiefer, Quantum gravity, 3rd edition. Oxford University Press, Oxford (2012), ISBN 978-0199212521.
  * [72] K. V. Kuchar and M. P. Ryan, “Is minisuperspace quantization valid?: Taub in mixmaster”, Phys. Rev. D 40 (1989), 3982-3996.
  * [73] J. A. Wheeler, “On the Nature of quantum geometrodynamics”, Annals Phys. 2 (1957) 604.
  * [74] M. Bouhmadi-López and C. Y. Chen, “Towards the Quantization of Eddington-inspired-Born-Infeld Theory”, JCAP 11 (2016), 023, [arXiv:1609.00700 [gr-qc]].
  * [75] M. Bouhmadi-López, C. Y. Chen and P. Chen, “On the Consistency of the Wheeler-DeWitt Equation in the Quantized Eddington-inspired Born-Infeld Gravity”, JCAP 12 (2018), 032, [arXiv:1810.10918 [gr-qc]].
  * [76] I. Albarran, M. Bouhmadi-López, C. Y. Chen and P. Chen, “Quantum cosmology of Eddington-Born–Infeld gravity fed by a scalar field: The big rip case”, Phys. Dark Univ. 23 (2019), 100255, [arXiv:1811.05041 [gr-qc]].
  * [77] P. W. Higgs, “Quadratic lagrangians and general relativity”, Nuovo Cim 11 (1959) 816 .
  * [78] G. V. Bicknell, “Non-viability of gravitational theory based on a quadratic lagrangian”, J. Phys. A: Math. Nucl. Gen. 7 (1974) 1061.
  * [79] B. Whitt, “Fourth-order gravity as general relativity plus matter”, Phys. Lett. B 145 (1984) 176.
  * [80] J. D. Barrow and S. Cotsakis, “Inflation and the Conformal Structure of Higher Order Gravity Theories”, Phys. Lett. B 214 (1988) 515.
  * [81] L. Pogosian and A. Silvestri, “The pattern of growth in viable f(R) cosmologies”, Phys. Rev. D 77 (2008), 023503; Phys. Rev. D 81 (2010), 049901, [arXiv:0709.0296 [astro-ph]].
  * [82] M. Born and J. R. Oppenheimer, “On quantum theory of molecules”, Ann. Physik 84 (1927) 457.
  * [83] C. Kiefer, “Continuous measurement of mini-superspace variables by higher multipoles”, Class. Quantum Grav. 4 (1987) 1369.
  * [84] R. Brout, G. Horwitz and D. Weil, “On the onset of time and temperature in cosmology”, Phys. Lett. B 192 (1987) 318.
  * [85] C. Kiefer, “Wave packets in minisuperspace”, Phys. Rev. D 38 (1988) 1761.